Moës

M Le Cren

A Martin

P Massin

N Moës

A robust 3D crack growth method based on the eXtended Finite Element Method and the Fast Marching Method

Keywords: Level sets, Signed distance functions, Fast Marching Method, Crack propagation

  L'archive ouverte pluridisciplinaire

Introduction

The level set method coupled with X-FEM [START_REF] Belytschko | Elastic crack growth in the finite elements with minimal remeshing[END_REF][START_REF] Moes | A finite element method for crack growth without remeshing[END_REF]) is very effective to simulate 2D [START_REF] Stolarska | Modelling crack growth by level sets in the extended finite element method[END_REF]) and 3D [START_REF] Moes | 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] crack growth. The initial crack is represented by two level set functions: the faces of the crack belong to the zero level set of a first level set function while the second level set function is defined in such a way the intersection of the zero level set of the two level set functions describes the crack front. The two level set functions are signed distance functions. The two level set functions are orthogonal in the sense that their gradients are orthogonal. The growth of the crack is discretized in increments. Linear elastic fracture mechanics is used to describe the displacement of a point of the crack front from its position at the beginning of the growth to its position at the end of the growth increment. At the end of each growth increment, the level set functions are updated to describe the new crack position.

The udpate of both level set level set functions is based on the level set method introduced by [START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. The procedure proposed by [START_REF] Gravouil | Non-planar 3D crack growth by the extended finite element and level sets -Part II: Level set update[END_REF] is as follows. In a first step, called the extension step, the displacement field is first extended to a neighborhood of the zero level set. In a second step, called the update step, the Hamilton-Jacobi equation describing the evolution of each level set function is integrated. The intersection of the zero level set of the two functions obtained after the update step describes the new position of the crack front. These functions are not expected to be true signed functions. The last step consists of two interlocked steps, called the reinitialization step and the orthogonalization step. The reinitialization step builds two new level set functions, from the functions obtained at the end of the update step. The orthogonalization step ensures the two updated level set functions are orthogonal.

Diffrent approaches to perform each of the four steps of the update of the two level set functions have been explored. [START_REF] Gravouil | Non-planar 3D crack growth by the extended finite element and level sets -Part II: Level set update[END_REF] proposed to integrate Hamilton-Jacobi equations to perform these four steps. [START_REF] Sukumar | Extended finite element method and fast marching method for threedimensional fatigue crack propagation[END_REF][START_REF] Sukumar | Threedimensional non-planar crack growth by a coupled extended finite element and fast marching method[END_REF] and [START_REF] Shi | Abaqus implementation of extended finite element method using a level set representation of three-dimensional fatigue crack growth and life predictions[END_REF] formulated the extension and the reinitialization steps in terms of eikonal equations and used the Fast Marching Method, introduced by [START_REF] Sethian | A marching level set method for monotonically advancing fronts[END_REF]. [START_REF] Prabel | Level set X-FEM non matching meshes: application to dynamic crack propagation in elastic-plastic media[END_REF] performed the four steps of the update of the two level set functions on a dedicated regular grid, in order to use an upwind scheme to integrate Hamilton-Jacobi equations. [START_REF] Colombo | Fast and robust level set update for 3D non-planar X-FEM crack propagation modelling[END_REF] introduced a geometric approach to perform the extension and update steps. [START_REF] Colombo | An implicit geometrical approach to level sets update for 3D non planar X-FEM crack propagation[END_REF] extended this geometric approach to perform also the reinitialization and orthogonalization steps, in order to obtain a more robust method.

The developments are implemented in code aster finite element package, developped by EDF (1989EDF ( -2021)). Three methods are available to perform the reinitialization and orthogonalization steps. Nevertheless, none of these methods fit our desire of a robust tool, able to perform the mixed mode propagation of a crack in an industrial structure and, for example, the Brokenshire test (see [START_REF] Barr | Torsion fracture tests[END_REF]) is still a challenge.

In this paper, we propose a new approach in which the geometric approach proposed by [START_REF] Colombo | An implicit geometrical approach to level sets update for 3D non planar X-FEM crack propagation[END_REF] is used for the orthogonalizaton step and the Fast Marching Method for the reinitialization step. Our implementation of this approach is designed for triangulated meshes. We then extend this method to all types of (linear) volume elements available in a standard finite element library.

Level set update to model the growth of a crack

The growth of the crack is discretized in growth increments. At the begining of a growth increment k, the crack Γ k is described by two level set functions : the normal level set function Φ k n and the tangential level set function Φ k t . The crack faces are described by the set of points:

Γ k = {M, Φ k n (M ) = 0} ∩ {M, Φ k t (M ) < 0}. ( 1
)
The crack front Γ k 0 is discretized by linear segments (cf. Fig. 1). Let P i be a point on the crack front. The local frame (T i , N i , B i ) is computed from the level set functions Φ k n and Φ k t :

B i = ∇Φ k n (P i ) ∇Φ k n (P i )
,

T i = ∇Φ k t (P i ) × ∇Φ k n (P i ) ∇Φ k t (P i ) × ∇Φ k n (P i ) , N i = B i × T i .
(2)

The standard finite element interpolation of the level set functions is used to compute ∇Φ k n (P i ) and ∇Φ k t (P i ). The mechanical fields are computed by means of X-FEM. We typically use an energetic approach to compute the energy release rate G i and stress intensity factors K i I , K i II and K i III . G i is then computed using the domain integral method [START_REF] Destuynder | Sur une Interprétation Mathématique de l'Intégrale de Rice en Théorie de la Rupture Fragile[END_REF][START_REF] Li | A comparison of methods for calculating energy release rates[END_REF] while the stress intensity factors are computed using interaction integrals [START_REF] Gosz | An interaction energy integral method for computation of mixed-mode stress intensity factors along non-planar crack fronts in three dimensions[END_REF]. We can also use the displacement jump extrapolation technique [START_REF] Chan | On the finite element method in linear fracture mechanics[END_REF] to compute K i I , K i II and K i III and apply Irwin's formula [START_REF] Irwin | Analysis of stresses and strains near the end of a crack traversing a plate. Journal of Applied Mechan-(a) (b) Fig. 22 Fracture surface observed experimentally (Barr and Brokenshire 1996) (a), and obtained by means of the Fast Marching Method[END_REF]) to obtain G i . These quantities are used to determine the crack growth size and crack growth direction.

Following [START_REF] Gravouil | Non-planar 3D crack growth by the extended finite element and level sets -Part II: Level set update[END_REF] and [START_REF] Sukumar | Threedimensional non-planar crack growth by a coupled extended finite element and fast marching method[END_REF], we assume a fatigue crack growth law and we use a modified Paris' law [START_REF] Paris | A Critical Analysis of Crack Propagation Laws[END_REF], in which G plays the role of the stress intensity range ∆K. The maximal crack growth size during a growth increment ∆a max is an input parameter of the simulation. The value of ∆a max corresponds to the typical element size in the neighborhood of the crack. The crack growth size ∆a i associated to the point P i on the crack front is:

∆a i = G i G max m ∆a max , (3)
where G max is the maximum of the energy release rate G, for all the points of the crack front, and C and m are the parameters of Paris' law.

The maximum hoop stress criterion is used to find the crack growth direction. The angle β of the crack growth direction with respect to the plane tangent to the crack is obtained by [START_REF] Erdogan | On the Crack Extension in Plates Under Plane Loading and Transverse Shear[END_REF]:

β = 2 tan -1   1 4   K I K II - K II |K II | K I K II 2 + 8     . (4)
A planar crack growth will correspond to assume β = 0. Finally, the displacement of the point P i on the crack front is given by the vector:

∆a i = ∆a i (cos β i N i + sin β i B i ), (5)
where β i is the angle β evaluated at point P i . At the end of the growth increment, the point Q i , defined by:

OQ i = OP i + ∆a i , (6)
lies on the new crack front Γ k+1

0

. The problem is now to build new level set functions Φ k+1 n and Φ k+1 t so that the new crack front Γ k+1 0 can be described as the intersection of the zero level set of Φ k+1 n and the zero level set of Φ k+1 t .

PDE-based approaches: Simplex and Upwind Methods

We now expose how the extension, update, reinitialization and orthogonalization steps are implemented for Simplex and Upwind Methods. The general scheme follows the approach proposed by [START_REF] Colombo | Fast and robust level set update for 3D non-planar X-FEM crack propagation modelling[END_REF].

Extension step

The extension step extends the displacement defined on the crack front to a neighbourhood of the crack front. This neighbourhood is the domain in which the level set functions will be updated during the update step. The simplest approach consists in updating the level set functions in the whole domain. One can also restrict the updated domain, in order to save computational time.

The first step of the extension consists in mapping each point of the domain to a point on the crack front, by projecting the point of the domain onto the crackfront. This first step is always performed for all the nodes in the domain, because it enables the computation of the distance from each node of the mesh to the crack front. This distance can then be used to restrict the domain in which the level set functions will be udpdated to a torus surrounding the crack front, following the procedure described in [START_REF] Colombo | Fast and robust level set update for 3D non-planar X-FEM crack propagation modelling[END_REF].

Let M be a node in the domain. The projection algorithm finds point P on the crack front Γ k 0 such that point M belongs to the plane (P, N, B). The algorithm computes point P and segment [P i P j ] to which P belongs. The displacement of point P reads:

∆a P = ∆a(cos βN + sin βB). (7)
The computation of crack growth size ∆a, angle β and local frame (T, N, B) at point P , from the data known at points P i and P j , is discussed in detail in [START_REF] Colombo | Fast and robust level set update for 3D non-planar X-FEM crack propagation modelling[END_REF].

Displacement at point M is computed from displacement at point P . Displacement ∆a P n is associated to the displacement of the zero level set of Φ n and displacement ∆a P t is associated to the displacement of the zero level set of Φ t . Displacement at point M is decomposed as:

∆a M = ∆a M n + ∆a M t , (8)
where ∆a M n is associated to the displacement of the zero level set of Φ n and ∆a M t is associated to the displacement of the zero level set of Φ t . Displacement ∆a M n is zero where Φ k t is negative to ensure Γ k ⊂ Γ k+1 at the end of the growth increment. Displacement ∆a M n where Φ k t is positive is computed assuming ∆a n varies linearly from ∆a n = 0 where Φ k t = 0 (i.e. on the crack front Γ k 0 ) to ∆a n = ∆a sin β where Φ k t = ∆a cos β.

Update step

We use the geometrical approach proposed in [START_REF] Colombo | Fast and robust level set update for 3D non-planar X-FEM crack propagation modelling[END_REF]. The updated level set functions at point M read:

Φ n (M ) = Φ k n (M ) -∆a M n • ∇Φ k n (M ), Φ t (M ) = Φ k n (M ) -∆a M t • ∇Φ k t (M ). (9)
Updated level set functions Φ n and Φ t are neither expected to be signed distance functions nor expected to be orthogonal.

Reinitialization and orthogonalization steps

The reinitialization and orthogonalization steps are interlocked. Actually, in the case of a non-planar crack, it is not possible to fulfill the orthogonality and unit gradient conditions at the same time, except in the neighborhood of the crack front. We use the following procedure:

∂Φ n ∂τ = - Φ n |Φ n | ( ∇Φ n -1), (10)
with initial data:

Φ n (τ = 0) = Φ n . ( 11
)
Let I n be the set of nodes belonging to cells which contain the zero level set of Φ n . Level set function Φ k+1 n is built so that the zero level set of Φ k+1 n coincides with the zero level set of Φ n . Thus, level set function Φ k+1 n at a given node of set I n is the distance from this node to the zero level set of Φ n , where Φ n is positive, and the opposite of the distance from this node to the zero level set of Φ n , where Φ n is negative. The computation of the distance from each node of I n to the zero level set of Φ n is based on a projection algorithm proposed in [START_REF] Colombo | Fast and robust level set update for 3D non-planar X-FEM crack propagation modelling[END_REF].

Orthogonalization. A level set function Φ t orthognal

to Φ k+1 n is computed as the steady-state solution of the following evolution problem:

∂Φ t ∂τ = - Φ k+1 n |Φ k+1 n | ∇Φ k+1 n ∇Φ k+1 n • ∇Φ t , (12)
with initial data:

Φ t (τ = 0) = Φ t . (13)
We use the result of the The difference between the Simplex Method and the Upwind Method is in the way the Hamilton-Jacobi equations ( 10), ( 12) and ( 14) are integrated. Simplex Method is dedicated to unstructured meshes (triangles or tetrahedra) and implements the method proposed in [START_REF] Barth | Numerical schemes for the Hamilton-Jacobi and Level Set equations on triangulated domains[END_REF]. Upwind Method is dedicated to regular grids (quadrangles or hexahedra) and implements the upwind scheme proposed in [START_REF] Prabel | Level set X-FEM non matching meshes: application to dynamic crack propagation in elastic-plastic media[END_REF]). Both methods use a one step explicit time integration scheme (explicit Euler) to integrate the Hamilton-Jacobi equations. Therefore the time step must be less than a critical value in order to satisfy the CFL condition.

Geometric Method

We now expose the fully geometric method proposed by [START_REF] Colombo | An implicit geometrical approach to level sets update for 3D non planar X-FEM crack propagation[END_REF]. The extension step discussed in section 2.1.1 is first performed. Thus, point P on the crack front Γ k 0 is such that point M belongs to the plane (P, N, B) which is known. Displacement ∆a P and angle β at point P are also known.

Let Q be the point on the crack front Γ k+1 0 corresponding to point P on the crack front Γ k 0 . Since ∆a P is the displacement at point P , we have:

PQ = ∆a P = ∆a(cos βN + sin βB). ( 16
)
Since β is the angle of the crack growth direction to the plane tangent to the crack, local frame (N Q , B Q ) at point Q is obtained by rotating local frame (N, B) by an angle β with respect to T. The level set values Φ k+1 n (M ) and Φ k+1 t (M ) form a cartesian coordinate system in plane (Q, N Q , B Q ) and we have:

QM = Φ k+1 n (M )B Q + Φ k+1 t (M )N Q . ( 17
)
The level set functions at the end of the growth increment are then computed as:

Φ k+1 n (M ) = QM • B Q , Φ k+1 t (M ) = QM • N Q . ( 18
)
This method condenses the update, reinitialization and orthogonalization steps in one step. No property of the underlying mesh is exploited so that this method can handle all types of elements.

Table 1 gives a concise comparison of the methods to update the level set functions. How each method performs extension, update, reinitialization and orthogonalization steps is shown.

Numerical tests

A mode I-II crack growth

We propose to compare the results obtained using the Simplex Method, the Upwind Method and the Geometric Method on a simple case. We consider the cracked plate depicted on Fig. 2. The length of the plate is L = 8 m, its height is H = 18 m and its thickness is e = 1 m. The initial crack is inclined with respect to the loading symmetry plane (xOy) by an angle of 45°. The initial crack front is located at a quarter of the

G (J/m 2 ) K I (MPa/ √ m) K II (MPa/ √ m) 45.9905 2.72750 1.38547
Table 2 Tabulated valued of the energy release rate and the stress intensity factors.

length and at half the height. The initial crack length is a 0 = 2 √ 2 m. The plate is submitted to a symmetric uniform traction t applied on the top face and the bottom face (see Fig. 2). Since the pre-crack is inclined with respect to the loading symmetry plane, a mode I-II crack growth is expected.

We consider two meshes. The first one is a rectangular grid with twenty eight elements along the length, sixty elements along the height and four elements in the thickness. The second one is obtained by splitting each hexaedron of the first mesh in six tetrahedra. Both meshes share the same nodes, so that the Geometric Method gives the same results with both meshes. The Upwind Method is applied to the first mesh, made of hexaedra. The Simplex Method is applied to the second mesh, made of tetrahedra. The energy release rate and the stress intensity factors are smoothed to make them uniform along the crack front. This way, all points on the crack front travel the same distance. Each simulation consists in one growth increment and the distance travelled by each point of the crack front is ∆a max = 0.4 m. Paris' law exponent is m = 1.

In order to test only the level set update, we use tabulated values of the energy release rate and of the stress intensity factors (cf. Table 2). These values were obtained using the mesh made of tetrahedra, a Young's modulus equal to 205000 MPa, a Poisson's ratio equal to 0.3 and a load amplitude t = 1 MPa. The energy release rate and the stress intensity factors are computed using the energetic approach in a domain restricted to a torus surrounding the crack front of radius R S = 1.16 m, which corresponds approximately to four times the size of an element.

The domain in which the normal and tangential level set functions are updated is restricted to a torus surrounding the crack front. The radius of this torus is automatically computed by the software from a minimal radius given by the user. We chose a minimal radius R = R S +∆a max = 1.56 m. The radius of the torus computed by the software is 2.31 m.

In this test the problem does not depend on the x coordinate. Consequently, we expect the normal and tangential level set functions to be uniform along axis Ox.

Figure 3 shows the results we obtained using the Simplex Method, the Upwind Method and the Geometric Method to update the level set functions. We computed level sets of the normal and tangential level set functions and we represented the projection of these level sets onto the plane (yOz).

The projection onto the plane (yOz) of the level sets of the level set functions obtained using the Simplex Method (cf. Fig. 3(a)) and the Upwind Method (cf. Fig. 3(b)) do not result in curved lines but surfaces. The level set functions obtained using the Simplex Method and the Upwind Method are thus not uniform along axis Ox. Conversely, the projection onto the plane (yOz) of the level sets of the level functions obtained using the Geometric Method (cf. Fig. 3(c)) result in curved lines so that the level functions obtained using the Geometric Method are uniform along axis Ox.

The level sets computed from the normal level set function obtained using the Geometric Method are not parallel to each other. Actually, due to the presence of the bifurcation of the crack, level sets of the normal level set function computed using the Geometric Method are discontinuous, since the normal level set function is not updated where the initial tangential level set function is negative. The post-processing is not able to represent these discontinuities in the general case. However, it is able to represent the kink of the zero level set of the normal level set function. Taking into account this artifact, we conclude the level sets computed from the normal level set function obtained using the Geometric Method are parallel to one another. The level sets computed from the tangential level set function obtained using the Geometric Method are parallel to one another. Furthermore, the level sets of the normal level set function are orthogonal to the level sets of the tangential level set function, where the normal level set function has been updated. This analysis emphasizes the difficulty to obtain accurate results by integrating the Hamilton-Jacobi equations to perform the reinitialization and orthogonalization steps. The level set functions obtained using the Geometric Method are quite good. Nevertheless, this method faces other issues, as we shall see in the analysis of the Brokenshire test.

Brokenshire test

Brokenshire test is a torsion test first proposed by [START_REF] Barr | Torsion fracture tests[END_REF]. The initial work of Barr and Brokenshire consisted in an experimental analysis. The setup of Brokenshire test is depicted on Fig. 4.

Several authors investigated this test with various approaches. [START_REF] Jefferson | Three dimensional finite element simulations of fracture tests using the Craft concrete model[END_REF] used a plastic-damagecontact model for concrete. [START_REF] Su | Finite Element Modelling of Complex 3D Static and Dynamic Crack Propagation by Embedding Cohesive Elements in Abaqus[END_REF] used cohesive elements embedded between solid elements. Baydoun and Fries (2012) used X-FEM and studied dif- We consider an unstructured tetrahedral mesh. In order to accurately calculate the energy release rate and the stress intensity factors the mesh is refined at each propagation step in a region surrounding the crack front. To do so we use HOMARD, an automatic mesh adaptation tool distributed along with code aster (Nico-las et al. 2016). The element size in the initial mesh is h 0 = 20 mm. HOMARD splits the edges such that the element size is divided by two in the refined region. We call HOMARD three times in a row, so that the element size is h = h 0 /2 3 = 2.5 mm in a torus surrounding the crack front of radius R S = 5h = 12.5 mm. The stress intensity factors and the energy release rate are computed using the displacement jump extrapolation technique. We extract the displacement of points lying on the crack faces, within a distance 4h = 10 mm from the crack front. The maximal crack growth size during a growth increment is ∆a max = 4h = 10 mm and the Paris' law exponent we arbitrarily chose is m = 1. The normal and tangential level set functions are updated in the whole domain.

Simplex and Upwind Methods fail to converge during the first increment step (the maximum number of iterations is fixed to five hundred). Geometric Method is able to perform at least the first two growth increments. The zero level sets of the resulting normal and tangential level set functions are shown in Fig. 5.

After the second growth increment, oscillations of the zero level set of the tangential level set function are observed (cf. Fig. 5(b)). As a result, the crack front becomes discontiunous and the computation stops.

We conclude none of the three approaches we presented so far is robust enough. We thus propose a new implementation of the Fast Marching Method and we expect it to be robust enough to be able to perform the simulation until the total failure of the specimen is reached. 3 Fast Marching Method for crack growth [START_REF] Sukumar | Extended finite element method and fast marching method for threedimensional fatigue crack propagation[END_REF][START_REF] Sukumar | Threedimensional non-planar crack growth by a coupled extended finite element and fast marching method[END_REF] proposed to couple the Fast Marching Method and X-FEM to simulate planar and non-planar crack growth. They used the Fast Marching Method to perform the extension and reinitialization steps. In this section, we propose a new approach in which only the reinitialization step is performed by means of the Fast Marching Method.

A level set function Φ satisfies the following eikonal equation:

∇Φ = 1. ( 19
)
The Fast Marching Method proposes to solve this eikonal equation, assuming the values of Φ in the neighborhood of the zero level set of Φ are known. Let d = |Φ| be the distance to the zero level set of Φ. The idea of the method is to propagate the information outward from smaller values of d to larger values. To ensure a monotonic propagation of the information, we split the process in two steps:

1. The function d = Φ is computed in the domain {M : Φ(M ) > 0} from known values of d in the neighborhood of the zero level set of Φ. We then take Φ = d. Each step propagates the (positive) distance function d from smaller to larger values. The propagation of the information during the Fast Marching Method process is illustrated in Fig. 6.

Γ Φ = 0 Φ = -1 Φ = -2 Φ = 1 Φ = 2
In practice, the neighborhood of the zero level set of Φ refers to a set of nodes I computed as follows. We loop over the cells and identify the ones which contain the zero level set of Φ. If the sign of the level set function Φ is not the same at all the nodes of a given cell, this cell contains the zero level set of Φ and the nodes belonging to this cell are added to I.

We now expose our new approach:

1. Extension step. The extension step discussed in section 2.1.1 is performed. 2. Update step. The update step discussed in section 2.1.2 is performed. We recall the updated level set functions Φ n and Φ t obtained at the end of this step are neither expected to be signed distance functions nor expected to be orthogonal. 3. Reinitialization of the normal level set function. The normal level set function at the end of growth increment Φ k+1 n is computed as the solution of the eikonal equation:

∇Φ n = 1, ( 20
)
by means of the Fast Marching Method. We recall I n is the set of nodes belonging to cells which contain the zero level set of Φ n . Initial data consists in the known values of level set function Φ k+1 n at nodes of set I n . Level set function Φ k+1 n is built so that the zero level set of Φ k+1 n coincides with the zero level set of Φ n . Thus, level set function Φ k+1 n at a given node of set I n is the distance from this node to the zero level set of Φ n , where Φ n is positive, and the opposite of the distance from this node to the zero level set of Φ n , where Φ n is negative. Computation of the distance from each node of I n to the zero level set of Φ n is based on the projection algorithm proposed in [START_REF] Colombo | Fast and robust level set update for 3D non-planar X-FEM crack propagation modelling[END_REF]. 4. Reinitialization of the tangential level set function.

The tangential level set function at the end of growth increment Φ k+1 t is computed as the solution of the eikonal equation: During each growth increment, the normal level set function and the tangential level set function are reinitialized. Each time a level set function is reinitialized, the zero level set of this level set function is altered (an illustration of this phenomeom is given in [START_REF] Chopp | Some improvements of the fast marching method[END_REF]). This behavior does not depend on the method used to solve the reinitialization problem and may occur using the Fast Marching Method. The numerical tests reported in section 3.3 show that despite the alteration of the zero level set of the level set functions induced by this phenomenon, our new approach gives acceptable results for our applications.

∇Φ t = 1, ( 21
)
The Fast Marching Method travels across the nodes of the mesh and updates the level set value of each visited node. The way the Fast Marching Method selects the next node to be updated and the way this node is updated depend on the mesh structure. In the following sections, we expose an approach dedicated to unstructured grids, and we then extend this approach from tetrahedra to all types of elements.

V 4 1 V V 2 V 3 Fig. 7 Representation of the tetrahedron V 1 V 2 V 3 V 4 .

Fast Marching Method applied to unstructured grids

We discuss in this section our implementation of the Fast Marching Method dedicated to unstructured grids, which consists in an extension to the tridimensional case of the method proposed by [START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF]. The distance function d represents the distance to the zero level set of a given function Φ, which can be either the normal level set Φ k+1 n or the tangential level set Φ k+1 t . The computation is performed for a given set of nodes N , which can correspond to the domain where Φ is positive or the domain where Φ is negative.

Estimation of the distance function at a given node

We consider a non degenerate tetrahedron in an unstructured grid. The vertices of this tetrahedron are denoted V 1 , V 2 , V 3 and V 4 (see Fig. 7). The value of the distance function d at vertex V i is denoted d i . Distances

d 1 = d(V 1 ), d 2 = d(V 2 ) and d 3 = d(V 3 ) are known and distance d 4 = d(V 4
) is unknown. The first order Taylor series expansion of distance function d around the vertex V 4 evaluated at point P reads:

d(P ) ≈ d(V 4 ) + V 4 P • ∇d(V 4 ). ( 22
)
We assume that this approximation is exact for P = V 1 , for P = V 2 and P = V 3 . We obtain the following system of three equations:

d 1 = d 4 + V 4 V 1 • ∇d(V 4 ), d 2 = d 4 + V 4 V 2 • ∇d(V 4 ), d 3 = d 4 + V 4 V 3 • ∇d(V 4 ). ( 23
)
Since the tetrahedron is not degenerate, (

V 4 V 1 , V 4 V 2 , V 4 V 3 ) is a basis of R 3 . Thus we can write: ∇d(V 4 ) = n 1 V 4 V 1 + n 2 V 4 V 2 + n 3 V 4 V 3 , (24)
where n 1 , n 2 and n 3 are the coordinates of ∇d(V 4 )

over (V 4 V 1 , V 4 V 2 , V 4 V 3 )
. System (23) can then be rewritten as follows:

   d 1 d 2 d 3    =    d 4 d 4 d 4    + G n, (25)
where n = (n 1 , n 2 , n 3 ) T and matrix G reads:

G =   V 4 V 1 • V 4 V 1 V 4 V 1 • V 4 V 2 V 4 V 1 • V 4 V 3 V 4 V 2 • V 4 V 1 V 4 V 2 • V 4 V 2 V 4 V 2 • V 4 V 3 V 4 V 3 • V 4 V 1 V 4 V 3 • V 4 V 2 V 4 V 3 • V 4 V 3   .
Since d is a distance function, the gradient of the distance at vertex V 4 is a unit vector. We thus have :

∇d(V 4 )•∇d(V 4 ) = 1. Matrix G represents the metric tensor in the basis (V 4 V 1 , V 4 V 2 , V 4 V 3 ).
Since the metric tensor represents the scalar product, we have:

n T G n = ∇d(V 4 ) • ∇d(V 4 ) = 1. ( 26
)
Equations ( 25) and (26) represent a system of four equations involving four unknowns: distance d 4 at vertex V 4 and coordinates n 1 , n 2 and n 3 of the gradient of distance at vertex V 4 . We define two criteria to accept a solution of this system of equations:

1. the consistency condition d 4 ≥ max(d 1 , d 2 , d 3 ); 2. the monotonicity condition n 1 ≤ 0, n 2 ≤ 0 and n 3 ≤ 0.

The consistency condition states that, since we want the update to propagate the information outward from smaller values to larger values, we want d 4 to be larger than d 1 , d 2 and d 3 . The monotonicity condition states that -∇d(V 4 ) lies in the tetrahedron, i.e. -∇d(V 4 ) ∈ C, where C is the convex cone defined as follows :

C = {α 1 V 4 V 1 + α 2 V 4 V 2 + α 3 V 4 V 3 , α 1 ≥ 0, α 2 ≥ 0, α 3 ≥ 0}.
We suppose the solid angle at vertex V 4 is acute, so that we have:

V 4 V 1 • V 4 V 2 ≥ 0, V 4 V 1 • V 4 V 3 ≥ 0 and V 4 V 2 • V 4 V 3 ≥ 0.
As a consequence, if the monotonicity condition is fulfilled, each component of the vector G n is negative so that system (25) gives:

d 1 -d 4 ≤ 0, d 2 -d 4 ≤ 0, d 3 -d 4 ≤ 0. ( 27
)
The monotocity condition then appears as a sufficient condition to fulfill the consistency condition. Nevertheless, we accept a solution of the system of equations if, and only if, these two conditions are fulfilled simultaneously.

The way to solve the system of equations, depends on a a priori hypothesis on ∇d(V 4 ):

1. The tridimensional update corresponds to the case where -∇d(V 4 ) is supposed to lie in the interior of convex cone C.

The bidimensional update corresponds to the case

where -∇d(V 4 ) is supposed to lie in a face of convex cone C. 3. The unidimensional update corresponds to the case where -∇d(V 4 ) is supposed to lie on an edge of convex cone C.

Tridimensional update In this case, -∇d(V 4 ) is supposed to lie in the interior of convex cone C. System (25) is rewritten as follows:

d -d 4 1 = G n, (28)
where d = (d 1 , d 2 , d 3 ) T and 1 = (1, 1, 1) T . Since matrix G represents a meric tensor, G is symmetric and positive-definite. We deduce from (28):

G -1 (d -d 4 1) = n, (29)
(d -d 4 1) T = n T G. (30)
We thus have:

(d -d 4 1) T G -1 (d -d 4 1) = n T G n = 1. (31)
Finally, d 4 is solution of the following quadratic equation:

ad 2 4 -2bd 4 + c = 0, (32)
where:

a = 1 T G -1 1, b = 1 T G -1 d, c = d T G -1 d -1. ( 33
)
We assume the discriminant ∆ of ( 32) is non-negative, so that this equation admits two solutions δ 1 and δ 2 , such that δ 1 ≥ δ 2 . We consider only the solution d 4 = δ 1 to enforce the consistency of the update. Since d 4 is known, we can solve (28) to obtain coordinates n 1 , n 2 and n 3 . If distance d 4 fulfills the consistency condition and if coordinates n 1 , n 2 and n 3 of the gradient of distance ∇d(V 4 ) fulfill the monotonicity condition, we accept distance d 4 . Otherwise, we try the bidimensionnal update.

Bidimensional update In this case, -∇d(V 4 ) is supposed to lie in a face of convex cone C, i.e. the intersection of convex cone C and one of the faces

V 1 V 2 V 4 , V 2 V 3 V 4 or V 3 V 1 V 4 .
One of the coordinates n 1 , n 2 and n 3 is set to zero :

n 3 is set to zero if we consider face V 1 V 2 V 4 , n 1 is set to zero if we consider face V 2 V 3 V 4 and n 2 is set to zero if we consider face V 3 V 1 V 4 .
As a consequence, one of the equations of system (25) is removed. Without loss of generality, we consider face V 1 V 2 V 4 , so that n 3 is set to zero and the third equation of system (25) has to be removed. System ( 25) is rewritten as follows:

d -d 4 1 = G n , ( 34
)
where

G = V 4 V 1 • V 4 V 1 V 4 V 1 • V 4 V 2 V 4 V 2 • V 4 V 1 V 4 V 2 • V 4 V 2 , d = (d 1 , d 2 ) T , 1 = (1, 1)
T and n = n 1 , n 2 T . We deduce from ( 34):

(d -d 4 1 ) T G -1 (d -d 4 1 ) = n T G n . ( 35
)
Since n 3 = 0, equation ( 26) gives:

n T G n = n T G n = 1. ( 36
)
Finally, d 4 is solution of the following quadratic equation:

a d 2 4 -2b d 4 + c = 0, ( 37
)
where a , b and c are obtained by replacing 1, d, G with 1 , d , G in (33). We proceed similarly to the tridimensionnal update to determine distance d 4 for face V 1 V 2 V 4 but we try also to compute distance d 4 from faces V 2 V 3 V 4 and V 3 V 1 V 4 and we keep the smallest obtained value. If a value of distance d 4 cannot be computed from any face, we try the unidimensionnal update.

Unidimensional update In this case, -∇d(V 4 ) is supposed to lie on an edge of convex cone C, i.e. the intersection of convex cone C and one of the edges V 1 V 4 , V 2 V 4 or V 3 V 4 . Two of the coordinates n 1 , n 2 and n 3 are set to zero : n 2 and n 3 are set to zero if we consider edge V 1 V 4 , n 1 and n 3 are set to zero if we consider edge V 2 V 4 and n 1 and n 2 are set to zero if we consider edge V 3 V 4 . As a consequence, two of the equations of system (25) are removed.

System ( 25) is rewritten as follows:

d 1 -d 4 = V 4 V 1 • V 4 V 1 n 1 . ( 38
)
The gradient of distance ∇d(V 4 ) reads:

∇d(V 4 ) = n 1 V 4 V 1 . ( 39
)
Since ∇d(V 4 ) is a unit vector, the coordinate n 1 fulfilling the monoticity condition is:

n 1 = - 1 V 4 V 1 . ( 40
)
Distance d 4 fulfilling ( 38) is then:

d 4 = d 1 + V 4 V 1 . (41)
A value of d 4 can be computed from edges

V 4 V 1 , V 4 V 2 and V 4 V 3 .
We select the smallest distance from the three possible values:

d 4 = min {d 1 + V 4 V 1 , d 2 + V 4 V 2 , d 3 + V 4 V 3 } . ( 42
)
A pathological case for which the tridimensional update fails Fig. 8 illustrates a case for which the consistency condition is fullfilled but not the monotonicity condition. The coordinates of the vertices over the reference frame (O, e x , e y , e z ) are V 1 (0, 1/2, 1), V 2 (0, -1/2, 1), V 3 (1, 0, 1) and V 4 (-1/2, 0, 2). The distances associated to vertices V 1 , V 2 and V 3 are respectively d 1 = 1, d 2 = 1 and d 3 = 1. The largest root of (32) in this case is

d 4 = 2 > 1 = max(d 1 , d 2 , d 3 )
, so that the consistency condition is fullfilled. The coordinates of ∇d(V 4 ) over the basis

(V 4 V 1 , V 4 V 2 , V 4 V 3 ) are n 1 = -3/4
, n 2 = -3/4 and n 3 = 1/2. Since n 3 > 0, the monotonicity condition is not fullfilled. The tridimensional update thus fails. We now focus on face V 1 V 2 V 4 . Fig. 9 illustrates this problem. The largest root of (37) in this case is d 4 = 1+ √ 5/2 > 1, so that the consistency condition is fullfilled. The coordinates of ∇d(V 4 ) over basis (

V 4 V 1 , V 4 V 2 , V 4 V 3 ) are n 1 = -1/ √ 5, n 2 = -1/ √
5 and n 3 = 0, so that n 1 ≤ 0, n 2 ≤ 0, n 3 ≤ 0 and the monotonicity condition is fullfilled. The bidimensional update is thus successful.

Sweeping process

The initial data consists in the known values of the distance function d at nodes of a given set I. The set of nodes for which a candidate value of distance d is known is denoted K and the set of nodes with a minimal distance is denoted P. Algorithm 1 outlines the way the nodes are visited by the Fast Marching Method applied to an unstructured grid.

In Algorithm 1 D refers to the set of candidate values for distance d j . There are three cases: Traditional exposition of the Fast Marching Method involves three sets of nodes [START_REF] Sethian | A marching level set method for monotonically advancing fronts[END_REF][START_REF] Sethian | Fast Marching Methods[END_REF]Kim-mel and Sethian 1998;[START_REF] Sethian | Fast methods for the Eikonal and related Hamilton-Jacobi equations on unstructured meshes[END_REF][START_REF] Bronstein | Weighted distance maps computation on parametric threedimensional manifolds[END_REF]):

-The distance function at a node in the set of "accepted" nodes is known and frozen. -The distance function at a node in the set of "tentative" nodes (the so called "Narrow band") is known but can be recomputed. -The distance function at a node in the set of "distant" nodes is unknown.

V 4 e x e z V 2 e y V 1 V 3 O -n
Fig. 8 Illustration of a configuration for which distance d 4 can not be computed using the tridimensional update. The monotonicity condition is not fullfilled, since

-n = -∇d(V 4 ) does not lie in convex cone C. V 1 V 2 V 3 V 4 e z e y e x

O

-n' Fig. 9 Illustration of the result of the bidimensional update for face V 1 V 2 V 4 in the configuration for which the tridimensional update fails. Distance d 4 can be computed because the monotonicity condition is fullfilled, since -n = -∇d(V 4 ) lies in convex cone C.

With our notations, the set of "accepted" nodes is I ∪P, the set of "tentative" nodes is K\(I ∪ P) and the set of "distant" nodes is N \K. Our algorithm does not distinguish the update of the distance function at a "tentative" node and the computation of an initial distance function at a node going from the set of "distant" nodes to the set of "tentative" nodes. In the same way, our algorithm does not perform the initialization of the nodes in the set of "tentative" nodes in a distinct block. Actually, while the node with the smallest value for the distance in K\P belongs to I, the algorithm populates the initial set of "tentative" nodes. When the initial set of "tentative" nodes is fully populated, the node with the smallest value for the distance in K\P belongs to the set of "accepted" nodes.

3.2 Extension of the Fast Marching Method from tetrahedra to all types of elements.

The method we just exposed is designed for a tetrahedral volume mesh. The simplest way to handle any type of element is to build a tetrahedral mesh by splitting the element of the given mesh in a preprocessing step. This solution introduces a new step, a new tool and a new source of problems. We thus propose in this subsection a way to extend the method to arbitrary volume meshes.

In the algorithm, the mesh structure is used to iterate over the tetrahedra adjacent to the node to update. The idea is to build a set of tetrahedra from an arbitrary cell. The algorithm then consists in two nested loops: the outer loop iterates over the cells adjacent to the node to update, while the inner loop iterates over the tetrahedra.

We recall the method computes an updated value of the distance function at the given node, from a given tetrahedron adjacent to the node, by computing an approximation of the gradient of the distance function. We thus propose to iterate over all the tetrahedra we can build from the vertices of the cell, so as many approximated gradients of the distance function as possible are explored. Obviously, some of the tetrahedra built this way can not be used to compute an updated value of the distance function at the given node, since the value of the distance function at some vertices of the cell can be not known yet. Such tetrahedra are simply omitted, as will be explained in the following.

We consider three types of volumic cell, except the tetrahedron: the pyramid, the pentahedron and the hexaedron. We now expose how we compute a set of tetrahedra from each of these cell types.

Pyramid Let {1, 2, 3, 4, 5} be the connectivity of the reference pyramid (cf. Fig. 10). It is possible to build the four following tetrahedra from the pyramid:

{1, 3, 4, 5}, {1, 2, 3, 5}, {1, 2, 4, 5}, {2, 3, 4, 5}. (43)
where {a, b, c, d} is the connectivity of the tetrahedron based on vertices a, b, c and d.

There are three different cases in which we can use the pyramid to update the value of the distance function at some vertex of the pyramid:

1. If the distance function is known at the four vertices of the base (i.e. vertices 1, 2, 3 and 4) and we want to update the value of the distance function at the apex (i.e. vertex 5), we iterate over the four tetrahedra of (43). Fig. 10(a) illustrates this case. 2. If the distance function is known at four nodes of the pyramid, including the apex, we iterate over only three tetrahedra. Pentahedron Let {1, 2, 3, 4, 5, 6} be the connectivity of the reference pentahedron (cf. Fig. 11). We first remark we can build the six following pyramids from the pentahedron:

{1, 2, 3, 4, 5}, {1, 2, 3, 4, 6}, {1, 4, 5, 6, 2}, {1, 4, 5, 6, 3}, {2, 3, 5, 6, 4}, {2, 3, 5, 6, 1}.

(44)

From each of these pyramids, we can build four tetrahedra. We then obtain the twenty four tetrahedra listed in Table 3. Each column in Table 3 refers to one of the six pyramids. In each column, the first row gives the connectivity of the pyramid, while the other rows give the connectivities of four tetrahedra corresponding to that pyramid. Some of these twenty four tetrahedra can be built from two pyramids. For example, we can build terahedron {1, 2, 4, 5} from both pyramids {1, 2, 3, 4, 5} and {1, 4, 5, 6, 2}. By removing the duplicates, we end up with the following list of twelve tetrahedra: {1, 2, 4, 5}, {1, 2, 4, 6}, {1, 2, 5, 6}, {1, 2, 3, 5}, {1, 2, 3, 6}, {2, 4, 5, 6}, {2, 3, 4, 5}, {2, 3, 4, 6}, {1, 3, 5, 6}, {1, 3, 4, 5}, {1, 3, 4, 6}, {3, 4, 5, 6}.

(45) Hexahedron Let {1, 2, 3, 4, 5, 6, 7, 8} be the connectivity of the reference hexahedron. We remark we can build twelve pentahedra from the hexahedron, by splitting the faces two by two. Figure 12(a) illustrates a first way to build four pentahedra from the hexaedron, Fig. 12(b) illustrates a second one and Fig. 12(c) illustrates a third Table 3 The tetrahedra obtained from the six pyramids one can build from the reference pentahedron. The duplicates are in red. one. Since we showed a pentahedron generates twelve tetrahedra, one hundred forty-four ( 144) tetrahedra (including duplicates) can be built.

In practice, we obviously do not store the connectivities of all these tetrahedra. We implemented a function which returns the connectivities of the twelve tetrahedra one can build from a pentahedron. We call this function for each of the twelve pentahedra one can build from a hexahedron to dynamically build the connectivities of the one hundred forty-four tetrahedra.

Numerical tests

A mode I-II crack growth

We propose to apply the Fast Marching Method to the mode I-II crack growth problem we already investigated using the Simplex Method, the Upwind Method and the Geometric Method in section 2.3.1. The results are depicted on Fig. 13.

The results obtained using the mesh made of tetrahedra or the mesh made of hexaedra are almost identical, which validates the extension of the method from tetrahedra to other types of elements. The level sets computed from the normal level set function obtained using the Fast Marching Method are parallel to each another. We remark the level sets of the normal level set function computed using the Fast Marching Method are smooth where the level sets of the normal level set function computed using the Geometric Method are discontinuous. The level sets computed from the tangential level set function obtained using the Fast Marching Method are parallel to each another. Furthermore, the level sets of the normal level set function are orthogonal to the level sets of the tangential level set function, where the normal level set function has been updated.

We conclude that the Fast Marching Method is as robust as the Geometric Method, since both methods lead to similar results where the level sets of the normal level set function are continuous. However the Fast Marching Method has the advantage of building smoother level set functions.

Convergence study

Our new approach based on the Fast Marching Method involves approximations at each step. In order to evaluate a posteriori the convergence of our new approach (involving the four steps of the update of the two level set functions), we propose a convergence study. We consider the uniform planar (β = 0) growth of a circular planar crack. The structure occupies domain Ω = [-1, +1]×[-1, +1]×[-1/8, +1/8]. The radius of the initial crack is R = 1/2. The initial crack occupies domain Γ 0 = {M (x, y, z) ∈ Ω : x 2 + y 2 ≤ R 2 , z = 0}. Each simulation consists in three growth increments and the distance travelled by each point of the crack front during one growth increment is ∆a max = ∆R = 1/20. Paris' law exponent is m = 1. At the end of growth increment k, normal level set function Φ k n and tangential level set function Φ k t read:

Φ k n (M ) = z, (46)
Φ k t (M ) = x 2 + y 2 -(R + k∆R). (47)
A tetrahedron for which the solid angle at each vertex is acute is called an acute tetrahedron. A tetrahedron which is not an acute tetrahedron is called an obtuse tetrahedron. A procedure to compute the solid angle at each vertex of a tetrahedron, based on the formula given by [START_REF] Writh | Relations between edge lengths, dihedral and solid angles in tetrahedra[END_REF], has been implemented. The caracteristics of the unstuctured meshes used in this convergence study are reported in table 4. Each mesh used in the convergence study includes a small part of obtuse tetrahedra.

The convergence of the Fast Marching Method is ensured only if the triangulation is acute, i.e. each tetrahedron of the triangulation is acute [START_REF] Sethian | Ordered upwind methods for static hamilton-jacobi equations: Theory and algorithms[END_REF]. To handle general triangulations, one can either use a method to split the obtuse tetrahedra [START_REF] Bronstein | Weighted distance maps computation on parametric threedimensional manifolds[END_REF][START_REF] Fu | A fast iterative method for solving the eikonal equation on tetrahedral domains[END_REF] or the Order Upwind Methods proposed by [START_REF] Sethian | Ordered upwind methods for static hamilton-jacobi equations: Theory and algorithms[END_REF].

We want to determine if a specific treatment of the obtuse tetrahedra is needed for our applications. We thus perform two convergence studies. In the first one, each time we need to compute the absolute value of a level set function at a vertex, from a tetrahedron such that the solid angle at the vertex is obtuse, we apply the Geometric Method instead of the procedure discussed in 3.1.1 (cf. Fig. 14). In the second one, no specific treatment is applied to obtuse tetrahedra (cf. Fig. 15). With or without the treatment of obtuse tetrahedra, we obtain the same order of convergence for both level set functions. We conclude no specific treatment of obtuse tetrahedra is needed for our applications.

In our approach the Geometric Method is used to initialize the Fast Marching Method. As a consequence, we have to consider two convergence mechanisms: 1-the convergence of the distance function at nodes used to intialize the Fast Marching Method and 2-the convergence of the distance function at the other nodes of the domain. The first convergence mechanism caracterizes the convergence of the Geometric Method while the second one caracterizes the convergence of the Fast Marching Method. Since the approximation of the distance at a given node in the Fast Marching Method is based on a first order Taylor series expansion, we expect the error in terms of L 1 norm to converge at order 1, as reported by [START_REF] Fu | A fast iterative method for solving the eikonal equation on tetrahedral domains[END_REF]. Using the Geometric Method in the whole domain, the error in terms of L 1 norm converges -1 e φ t -step 1 (slope = 2.07) φ t -step 2 (slope = 2.04) φ t -step 3 (slope = 2.00) Fig. 16 Error in terms of L 1 norm for the tangential level set function as a function of the normalized mesh size, obtained using the Geometric Method in the whole domain. at order 2 (cf. Fig. 16). For each growth step, the convergence order of the whole process is in between these two values (1 and 2), which seems consistent.

A three-point bending test

A cracked plate subjected to three-point bending is depicted Fig. 17. The length of this beam, associated to direction x, is L = 260 mm. Its thickness, associated to direction z, is t = 10 mm and its width, associated to direction y, is W = 60 mm. The initial crack length, associated to direction y, is a = 20 mm. The initial crack is inclined with respect to plane (yOz) by an angle of 45°. We have AB = 2L e = 240 mm. Line segments AA and BB correspond to the supports. The loading consists in a compressive force distribution applied on line segment CC . This configuration has been investigated experimentally by [START_REF] Lazarus | Comparison of predictions by mode II or mode III criteria on crack front twisting in three or four point bending experiments[END_REF] and numerically by [START_REF] Citarella | Comparison of crack growth simulation by DBEM and FEM for SEN-specimen undergoing torsion or bending loading[END_REF]. It has already been used by [START_REF] Colombo | Fast and robust level set update for 3D non-planar X-FEM crack propagation modelling[END_REF] and [START_REF] Geniaut | A simple method for crack growth in mixed mode with X-FEM[END_REF] to challenge crack growth methods available in code aster. Since the pre-crack is inclined with respect to the loading symmetry plane, the crack is subjected to a mixed mode loading condition and a twisted propagation is observed. We expose in this section the behaviour of the Fast Marching Method applied to this three-point bending test.

The displacement along direction y is prescribed to zero on both line segments AA and BB . The displacement along direction z is prescribed to zero at both points A and B. The displacement along direction x is prescribed to zero at point C. The resultant of the load applied on line segment CC is F = 2.4 kN. The material is PMMA with Young's modulus E = 2800 MPa and Poisson's ratio ν = 0.38.

We consider an unstructured tetrahedral mesh. In order to accurately calculate the energy release rate and the stress intensity factors we use HOMARD to refine the mesh at each propagation step in a region surrounding the crack front. The element size in the initial mesh is h 0 = 2.5 mm. We call HOMARD five times in a row, so that the element size is h = h 0 /2 5 = 0.078125 mm in a torus surrounding the crack front of radius R S = 5h = 0.390625 mm. The energy release rate and the stress intensity factors are computed using the energetic approach in a domain restricted to this torus.

The maximal crack growth size during a growth increment is ∆a max = 4h = 0.3125 mm and the Paris' law exponent we arbitrarily chose is m = 1. The domain in which the normal and tangential level set functions are updated is restricted to a torus surrounding the crack front. We asked the radius of this torus to be greater than R = R S + ∆a max = 0.703125 mm.

At the end of the simulation, the total failure of the specimen has almost been reached (cf. Fig. 18). We compare the evolution of the crack front position during the simulation we obtained using the Fast Marching Method with the same evolution of the crack front position reported in [START_REF] Citarella | Comparison of crack growth simulation by DBEM and FEM for SEN-specimen undergoing torsion or bending loading[END_REF]. Fig. 19 shows the projection of the crack fronts on plane (yOz) while Fig. 20 shows the projection of the crack fronts on plane (xOz). Citarella and Buchholz extracted the crack fronts at the beginning of some growth increments. We extracted from our data the crack fronts corresponding to the ones reported by Citarella and Buchholz. Each curve is indexed by the growth increment to which it corresponds. The growth increments we extracted do not correspond to the ones reported by Citarella and Buchholz because they used a value of the maximal crack growth size during a growth increment ∆a max = 1 mm, which corresponds to about three times the value we used. One also remarks Citarella and Buchholz used criteria to compute the direction of the crack growth direction and the growth size different from the ones we used. Nevertheless, the solution we obtained with the Fast Marching Method is, qualitatively, in good agreement with the one reported in [START_REF] Citarella | Comparison of crack growth simulation by DBEM and FEM for SEN-specimen undergoing torsion or bending loading[END_REF].

Brokenshire test

We expose in this section the behaviour of the Fast Marching Method applied to the Brokenshire test we introduced in 2.3.2. The zero level sets of the normal and tangential level set functions computed after the first two growth increments are depicted Fig. 21.

Unlike the result we obtained using the Geometric Method, no oscillations are observed on the zero level set of the tangential level set function after the first two growth increments. The simulation went on after the second growth increment and the total failure of the specimen was reached. We conclude our implementation of the Fast Marching Method is more robust than the one of the Geometric Method.

The fracture surfaces observed experimentally and obtained by means of the Fast Marching Method are depicted on Fig. 22. They are similar. This qualitative comparaison allows us to conclude our implementation of the Fast Marching Method is able to reproduce a complex crack path.

Conclusion

We proposed in this paper a new crack growth method based on X-FEM and Fast Marching Method. We used a Fast Marching Method designed for tetrahedral volume meshes and extended to all types of linear elements (pyramids, pentahedra and hexahedra). Thus, our approach allows us to use the same mesh to solve the mechanical problem and to update the level set functions. The price to pay for this ease of use is an order one discretization scheme of the Eikonal equation used to update the level set functions. For their part, [START_REF] Sukumar | Extended finite element method and fast marching method for threedimensional fatigue crack propagation[END_REF][START_REF] Sukumar | Threedimensional non-planar crack growth by a coupled extended finite element and fast marching method[END_REF] chose to update the level set functions on an auxiliary rectangular grid to benefit from a second order scheme. We proposed a new implementation of the Fast Marching Method applied to tetrahedral volume meshes, which consists in an extension of the method proposed by [START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF]. Our implementation is suitable for a finite element software. We also proposed a new presentation of the algorithm which is straightforward to implement.

We proposed several numerical examples to show the capabilities of our approach. We used a simple mode I-II crack growth to show we are able to compute smooth level set functions in the presence of kinks. We performed a convergence study of the proposed approach on a crack growth problem for which an analytical solution, in terms of the level set functions, is known. We investigated the three point bending test introduced by [START_REF] Lazarus | Comparison of predictions by mode II or mode III criteria on crack front twisting in three or four point bending experiments[END_REF] to demonstrate the approcah is able to handle non planar cracks. Finally, this approach was able to simulate the Brokenshire test [START_REF] Barr | Torsion fracture tests[END_REF] until the total failure of the specimen while all other methods available in code aster stopped after two growth increments.

We showed the approach we propose is able to compute smooth level set functions. We now want to get rid of the robustness issues coming from accumulated errors in the position of the crack front during the crack growth. The position of the crack front is computed, at the end of each growth increment, from the energy release rate and the stress intensity factors. Consequently, our next objective is to robustify the computation of the energy release rate and the stress intensity factors.

Fig. 1

 1 Fig. 1 Discretized crack front and local frame associated with each point on the crack front.

Fig. 2

 2 Fig. 2 Geometry and loading for the mode I-II crack growth.

Fig. 3

 3 Fig. 3 Representation of the level sets of the normal and tangential level set functions computed at the end of the crack growth using different implementations for the update step: (a) Simplex Method, (b) Upwind Method and (c) Geometric Method. The zero level sets of the two level set functions are shown in red. The update of the level set functions is restricted to a torus surrounding the crack front.

Fig. 4

 4 Fig. 4 Geometry and loading of Brokenshire test (Su et al. 2010).

  (a) After one growth increment.(b) After two growth increments.

Fig. 5

 5 Fig. 5 Representation of the zero level set of the normal level set function (in grey) and of the zero level set of the tangential level set function (in purple) during the growth of the crack, obtained using the Geometric Method.

Fig. 6

 6 Fig. 6 Propagation of the information during the Fast Marching Method process applied to the computation of the signed distance function Φ to the interface Γ.

  2. The function d = -Φ is computed in the domain {M : Φ(M ) < 0} from known values of d in the neighborhood of the zero level set of Φ. We then take Φ = -d.

1.

  If distance d j can be computed from the tetrahedron, D is a singleton. 2. If distance d j can be computed from at least one of the faces, D contains from one to three candidate values. 3. If distance d j must be computed from the edges, D contains three candidate values.

Fig. 10

 10 Fig. 10 Spliting a pyramid into tetrahedra. Orange nodes represent known level sets values and red nodes represent level sets values to be determined.

  Fig. 10(b) illustrates the case of the update of the distance function at vertex 3 from known values of the distance function at vertices 1, 2, 4 and 5. In this case, we iterate over tetrahedra {2, 3, 4, 5}, {1, 3, 4, 5} and {1, 2, 3, 5}. 3. If the distance function is known at three vertices of the pyramid, we consider only one tetrahedron. Fig. 10(c) illustrates the case of the update of the distance function at vertex 1 or 4 from known values of the distance function at vertices 2, 3 and 5. If we want to update the distance function at vertex 1 we consider only tetrahedron {1, 2, 3, 5} while if we want to update the distance function at vertex 4 we consider only tetrahedron {2, 3, 4, 5}.

Fig. 11

 11 Fig. 11 Spliting of a pentahedron into pyramids.

Fig. 12

 12 Fig. 12 Spliting of a hexahedron into pentahedra.

Fig. 13

 13 Fig. 13 Representation of the level sets of the normal and tangential level set functions computed at the end of the crack growth using the Fast Marching Method for the update step and applied to: (a) the mesh made of tetrahedra and (b) the mesh made of heaxaedra. The zero level sets of the two level set functions are shown in red.

  Fig. 14 Error in terms of L 1 norm for the level set functions as a function of the normalized mesh size, obtained with the treatment of obtuse tetrahedra: normal level set function (top) and tangential level set function (bottom).

Fig. 15

 15 Fig. 15 Error terms of L 1 norm for the level set functions as a function of the normalized mesh size, obtained without the treatment of obtuse tetrahedra: normal level set function (top) and tangential level set function (bottom).

  Fig. 173D cracked plate subjected to a three-point bending test (the initial crack is in red).

  Fig. 173D cracked plate subjected to a three-point bending test (the initial crack is in red).

Fig. 18

 18 Fig. 18 Crack path obtained at the end of the simulation of the three-point bending test.

Fig. 19 Fig. 20

 1920 Fig. 19 Crack front plots, coordinate y as a function of the coordinate z: (a) results obtained using the Fast Marching Method, (b) results obtained by[START_REF] Citarella | Comparison of crack growth simulation by DBEM and FEM for SEN-specimen undergoing torsion or bending loading[END_REF]. The legend gives the growth increment corresponding to each position of the crack front.

  (a) After one growth increment.(b) After two growth increments.

Fig. 21

 21 Fig. 21 Representation of the zero level set of the normal level set function (in grey) and the zero level set of the tangential level set function (in blue) during the growth of the crack, obtained using the Fast Marching Method.

Table 1

 1 Comparison of the methods to update the level set functions studied in this article: Simplex Method, Upwind Method, Geometric Method and the new approach based on the Fast Marching Method.

			Simplex Method	Upwind Method	Geometric Method	Fast Marching Method
	Extension	step	Projection of the nodes onto the crack front.	Projection of the nodes onto the crack front.	Projection of the nodes onto the crack front.	Projection of the nodes onto the crack front.
	Update	step	Geometrical approach.	Geometrical approach.	Implicit in the fully geometric method.	Geometrical approach.
	Reinitialization	step (Φ n )	Integration of a time-dependent PDE; solver dedicated to triangulated meshes; boundary condition: Φ n = Φ k+1 n at nodes	Integration of a time-dependent PDE; solver dedicated to regular meshes; boundary condition: Φ n = Φ k+1 n at nodes	Implicit in the fully geometric method.	Integration of a stationnary PDE; solver dedicated to triangulated meshes; boundary condition: Φ n = Φ k+1 n at nodes
			of I n .		of I n .		of I n .
	Orthogonalization	step	Integration of a time-dependent PDE; solver dedicated to triangulated meshes; boundary condition: Φ t = Φ t at nodes of I n .	Integration of a time-dependent PDE; solver dedicated to regular meshes; boundary condition: Φ t = Φ t at nodes of I n .	Implicit in the fully geometric method.	Replaced by the computation of Φ k+1 t at nodes of I t , using the fully geometric method.
	Reinitialization	step (Φ t )	Integration of a time-dependent PDE; solver dedicated to triangulated meshes; boundary condition: Φ t = Φ k+1 t at nodes	Integration of a time-dependent PDE; solver dedicated to regular meshes; boundary condition: Φ t = Φ k+1 t at nodes	Implicit in the fully geometric method.	Integration of a stationnary PDE; solver dedicated to triangulated meshes; boundary condition: Φ t = Φ k+1 t at nodes
			of I t .		of I t .		of I t .
	1. Reinitialization of the normal level set function. The
		normal level set function at the end of growth incre-
		ment Φ k+1 n	is computed as the steady-state solution
		of the following evolution problem, with respect to
		a virtual time τ :

  is the distance from this node to the zero level set of Φ t , where Φ t is positive, and the opposite of the distance from this node to the zero level set of Φ t , where Φ t is negative. The computation of the distance from each node of I t to the zero level set of Φ t is based on the projection algorithm presented in the paragraph dedicated to the reinitialization of the normal level set function.

	growth increment Φ k+1 t	is computed as the steady-
	state solution of the following evolution problem:
	∂Φ t ∂τ	= -	Φ t |Φ t |	( ∇Φ t -1),	(14)
	with initial data:
	Φ t (τ = 0) = Φ t .	(15)
	Let I t be the set of nodes belonging to cells which
	contain the zero level set of Φ t . The level set func-
	tion Φ k+1 t	is built so that the zero level set of Φ k+1 t
	coincides with the zero level set of Φ t . Thus, the
	level set function Φ k+1 t	at a given node of the set
	I t This procedure ensures the level set functions Φ k+1 n Φ k+1 t are signed distance functions.	and
						projection algorithm pre-
						sented above to compute the value of the level set
						function Φ t at nodes of I n . The tangential level set
						function Φ t is first interpolated at the vertices of the
						triangulation of the zero level set of Φ n . Let M be a
						node in I n . The projection algorithm associates to
						M a triangle T and a point P T in T . The value of
						the level set function Φ t at M is computed as the
						linear interpolation of the tangential level set func-
						tion values computed at the vertices of the triangle
						T evaluated at point P T .
						3. Reinitialization of the tangential level set function.
						The tangential level set function at the end of the

  by means of the Fast Marching Method. We recall I t is the set of nodes belonging to cells which contain the zero level set of Φ t . Initial data consists in the known values of level set function Φ k+1 at a given node of set I t cannot be deduced from the distance from this node to the zero level set of Φ t . We thus apply the Geometric Method discussed in section 2.2 to compute level set function Φ k+1 t at nodes of set I t . This way, the reinitialization of the tangential level set function and the orthogonalization step are condensed in one step.Table1gives a concise comparison of our new approach based on the Fast Marching Method, with respect to the three other methods to update the level set functions (Simplex Method, Upwind Method and Geometric Method).

		t	at nodes
	of set I t . Since level set function Φ t and level set
	function Φ k+1 n Φ k+1 t	are not orthogonal, level set function

Table 4

 4 Statistics of the unstructured meshes used in the convergence study.

	mesh size normalized mesh size (h) # nodes # tetrahedra # obtuse tetra. obtuse tetra. (%) max. angle (°)
	1.47E-1	1.00E+0	892	2893	20	0.69	104.19
	7.92E-2	5.40E-1	4771	19520	30	0.15	115.89
	4.10E-2	2.79E-1	30564	152574	119	0.08	120.95
	2.09E-2	1.43E-1	212917	1183921	622	0.05	128.34

Acknowledgements We sincerely thank the two anonymous reviewers whose suggestions helped us to improve this manuscript.

Algorithm 1: Fast Marching Method on an unstructured grid INITIALIZATION: K = I. P = ∅.

PROCESSING:

While K\P = ∅ Find the node i 0 with the smallest value for the distance in K\P. Add node i 0 to P.

For each tetrahedron T adjacent to the node i 0 For each node j of the tetrahedron T If j / ∈ I and j / ∈ P Then If all the other nodes adjacent to T are in K Then

The vertices of 32) is non-negative Then Compute d 4 as the largest root of (32).

Compute the coordinates of ∇d(V 4 ) by solving (28).

If the consistency condition and the monotonicity condition are fulfilled Then