A complex 3D fracture simulation of quasi-brittle material in quasi-static is still hard to tackle nowadays. Many methods and models propose partial solutions to this problem. The Thick Level Set (TLS) model, which uses an approach mixing damage mechanics and explicit crack representation, provides an easy fracture initiation, a complex crack growing capability (coalescing or branching) and an accurate tortuous fracture path. In this thesis we will demonstrate that the implementation of this model in a parallel 3D context provides an accurate and versatile tool that potentially scales. Regarding the accuracy, a novel tool called the "double cut algorithm" has enhanced the existing TLS implementation by letting pass a straight fully damaged zone in a mesh element without conditions on its size. This tool also brings a way to optimize the discretization by coarsening the mesh in a crack front wake. This adaptation reduces the size of the discrete mechanical problem and therefore the effort for the linear algebra resolution. As far as the scaling is concerned, the bottleneck is the linear algebra resolution time and its associated memory consumption. The parallel solving strategy developed in this thesis to tackle this problem starts first with a basic approach. Then by switching to a method close to domain decomposition and later to a two-scale method, it permits increasing scalability. The other TLS tasks are also partially parallelized. The principal concern is to obtain a tool that either runs faster or can treat a more significant problem if we provide much more computational units. Finally, some test cases illustrate the obtained results with a parallel 3D implementation.

Chapter I

Context and objectives

This thesis was written at "Ecole Centrale de Nantes" at the Gem Laboratory and is related to the complex 3D fracture simulation of quasi-brittle material. The term complex implies here that the simulations can initiate one or more cracks which develop by describing a set of tortuous paths which can coalesce or branch out Three categories composed of strong, weak or discrete discontinuity approaches cover such kind of simulations where the cracks opening must be represented somehow. The first one, associated with the fracture mechanics is able to reproduce the crack discontinuous displacement behavior by using an explicit discretization of the crack. It, therefore, gives a correct crack path. However, it needs a priori crack initiation and coalescing or branching are hard to achieve. The second, associated with the damage mechanics, gives an easy crack initiation, branching or coalescing. Nevertheless, it needs a regularization and does not provide an explicit crack path (and the associated discontinuous displacement). The third one, related to the discrete element method, uses beam or trust failure to represent the crack which leads to some mesh dependency. The recent Thick Level Set model used in this thesis sets a bridge between these first two discontinuity approaches by mixing a damage model with an explicit crack representation. Used in the 2D simulation it already offers a valuable tool in terms of results quality and computational consumption. In this thesis the objective is to maintain such quality for a 3D simulation and, in particular, to get a reasonable simulation computation time. This chapter develops the fracture simulation context and details the meaning of what "maintains such quality" and "reasonable simulation computation time" mean.

I.1 General context

In many situations the robustness of the objects or constructions that one manufactures is a security problem. Sometimes, a crack initiation in a structure is not affordable. Satellite, for example, is not expected during the launching phase to be deteriorated as it cannot be repaired. Sometimes having a way to predict a crack creation and its development in a structure helps to design it. For example, to evaluate the security risk after a blast in an Atomic power plant, one would try to evaluate the leaks of the contaminated fluid that may pass through cracks and then define the structure of the wall in accordance. In other cases, the crack path knowledge may help to set the appropriate reinforcements to avoid these cracks. On the contrary, in some situations one wants to crack objects. For example, in hydraulic cracking, people are eager to optimize water injection to obtain the best-fractured ground. In any case, the cracks and their propagation are phenomena that one wishes to control.

To control such phenomena, one needs to reproduce the cracks to be able to act on them. A perfect tool to achieve this goal is numerical simulation. It reproduces the complex behavior of objects in a low-cost and accurate way. More precisely we will focus on the quasi-static failure simulation of quasi-brittle material in 3D structures.

The general vocabulary and the characteristic lengths usually associated with a fracture problem are presented in figure I.1. The characteristic lengths, given in this figure, help compare the problem dimension together. Characteristic lengths do not represent the exact dimension but only its order of magnitude. L is associated to object sizes. In this object, a crack is described by a and B, respectively associated with its length and its tip width. This crack may develop at a distance associated with W , from the edge of the object. It may open (i.e., its lips may separate from each other) under some loading conditions. It may grow, with crack tips advancing in the material (a and/or B may increase). ρ is associated with the crack path curvature. Finally, h is associated with the size of the spacial problem discretization (illustrated here with a finite element mesh). In this document h will specifically be named h t , h w or h o for respectively mesh size in the crack tip vicinity, mesh size in the crack wake region or mesh size outside both regions.

With these conventions, the following chapter introduces the numerical model retained in this work.

I.2 Fracture and damage mechanics

Fracture modeling implies considering the discontinuity in the displacement field, which is introduced by the cracks (lips separate from each other). The numerical methods to simulate such problems can be classified into three approaches: strong, weak and discrete discontinuity.

In the strong discontinuity approach, the crack introduces a discontinuous displacement in the solution field explicitly. In the weak one, the displacements stay continuous but the material strength drops in a region where the crack is expected to be located. It provides a crack opening behavior without a clear identification of lips. In the discrete approach, inherited from the Discrete Element method devised by [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], the failure of the discrete elements, such as lattices of trusses and beams with Voronoi tessellation [START_REF] Bolander | Fracture analyses using spring networks with random geometry[END_REF]), is interpreted as a discontinuous phenomenon (a crack). This last approach is not studied here. However, the two others are the cement for the model used in this work.

The strong discontinuity approach, based on [START_REF] Griffith | The Phenomena of Rupture and Flow in Solids[END_REF] and Ir-win [1957] fracture mechanics theory uses the finite element method with an explicit cracks discretization. In other words, it means that both crack lips are meshed independently; they may have their own life because of the allowed discontinuous displacement across the crack.

A fine discretization to capture singularity is also mandatory near the crack tip. It is possible then to use the assumption of fracture mechanics to compute the crack growth. This assumption is based on the stress intensity factor computed at the crack tip. The two problems related to this approach are:

• The size of the crack compared to the dimension of the object studied is so small (a L and b L ) that the mesh size needs to be adapted around it. When the crack curvature (ρ) is small, the path also needs a fine discretization. It introduces a considerable meshing effort (with robustness issues) and gives a rather large system to solve.

• When a crack grows, the mesh must be almost completely recreated.

Many strategies were used to minimize those issues. Among them there is the crack box. It corresponds to a strategy where a region (a box) is assumed to encapsulate the crack during all the simulation. This zone is finely discretized to capture the crack behavior and is connected somehow to a coarser mesh. It is supposed to reduce the global problem size and to ease re-meshing by limiting it to the box region. Samcef (from Siemens [2019]) formerly uses "glue" to connect the box to the coarser mesh. Franc3D (from Fracture Analysis Consultants [2019]) creates a box region by extracting a portion of a 3D finite element model. A re-meshing of this portion with mesh boundary constraints is made so that the connection to the model remains identical.

Alternatively, the boundary element method (used in the first Franc3D version) proposes a more straightforward discretization approach. Only the crack and the studied object boundary are meshed.

The 3D crack progression is simply associated with a 2D meshing tool. However, the system matrix is dense with this method. It leads to unacceptable linear algebra computational effort.

Still chasing for a solution without meshing issues, researchers propose meshfree methods. In those techniques [START_REF] Belytschko | Crack propagation by element-free Galerkin methods[END_REF], [START_REF] Bordas | Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment[END_REF]), the use of element-free Galerkin methods gives among other things, an easy way to propagate a crack by setting moving particles to follow the crack path.

Nevertheless, the real jump in performance was introduced by the use of the eXtended Finite Element Method (X-FEM) in Mo ës et al. [1999]). The definition of the crack no longer relies on a fine mesh discretization. The mesh is now of a uniform fixed size, and the level set tools (J.A. [START_REF] Sethian | Level Set Methods and Fast Marching Methods Evolving Interfaces in Computational Geometry[END_REF]) are used to represent crack and place it where it has to be in the object. Enrichment technique is then used to describe the discontinuity (crack opening) and the particular strain field at the crack tip. The crack growth, still based on fracture mechanics, is achieved only by moving the level set with the mesh remaining fixed. The main drawback is that the crack initiation needs to be done manually (there must be an existing crack somewhere to start the computation). Moreover, the crack coalescence or branching is hard or even impossible to achieve.

The weak discontinuity approach associates, on its side, irreversible microscopic phenomena such as micro-cracks or micro-voids to the material rigidity loss. It is done through the use of an internal variable, which represents a micro-default density. In its scalar form1 , called further d, it ranges from 0 (sound material) to 1 (completely damaged material without any rigidity). For instance the free energy density of linear elastic damage law can be written using Hooke's law tensor C, strain tensor and d as follows:

ϕ( , d) = 1 2
(1d) : C :

(I.1)
The damage energy release rate Y associated with this potential (Y = -∂ϕ ∂d ) is then used to settle the evolution laws that are used to make damage evolving in the studied part. Compared to the strong discontinuity approach, the damage mechanics provide an easy way to initiate a crack from a sound situation. It also gives access to a complex path with coalescence or branching. However, this local treatment of the damage induces some spurious localizations with a pathological mesh dependency. To overcome this problem, different models were proposed to regularize the damage evolution. They share the introduction of a characteristic length that tunes the averaging process used for regularization.

Integral-type or gradient-type model families are commonly proposed in the scientific literature to smooth the quantity Y . The former, the integral-type family, proposes a spatial weighted averaging process at the studied material point. It is based on the integration of some weighted quantity, which is performed over a domain of fixed sizes. The possible smoothed quantities can be directly Y or the associated strain as proposed in Pijaudier-Cabot and [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF], [START_REF] Bažant | Nonlocal Integral Formulations of Plasticity and Damage : Survey of Progress[END_REF]. The latter, the gradient-type family, devised by [START_REF] Peerlings | Gradient-enhanced damage modelling of concrete fracture[END_REF], is based on the use of an averaged Y implicit gradient formulation. This regularized Y is obtained from a solving of a partial differential equation where Y is the source term and where arbitrary null boundary condition is imposed. In terms of the computational cost, the former adds an extra non-local averaging task for each integration at Gauss point and the latter adds a resolution of an auxiliary Laplacian problem.

Other approaches add to the free energy, a term based on gradient of damage [START_REF] Fr | Damage, gradient of damage and principle of virtual power[END_REF], [START_REF] Lorentz | Gradient damage models: Toward full-scale computations[END_REF]) that penalizes the sudden variations of damage.

To conclude this non-exhaustive review, a more recent model assimilates d to a crack phase-field such as in [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF], [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF]. The regularization of d is transformed into a minimization problem of a diffusive crack topology.

The TLS model has taken the best of both discontinuity approach.

Like in damage mechanics, it uses a damage variable d but it considers its treatment in a configurational mechanics manner. A surface (front) separates sound from damaged material. The damage evolution is then linked to the front movement. It leads to compute Y on this front. Furthermore, to again avoid spurious localization on the front, a regularization such as provided in the integral-type model, performs a Y weighted averages on segments 2 originated from front. The length of the segments is not fixed during time and can evolve from zero to a maximum length l c . The orientation of these segments is al-ways aligned with the damage gradient. This weak discontinuity treatment provides to the TLS model, its capacity to initiate a damaged zone from a sound situation. The damage progression, linked to the front progression, is made very versatile (tortuous path, merging, forking, ...). Moreover, the model also adds automatically crack lips in the damaged zone where d = 1. Using the enrichment techniques from the strong discontinuity approach, those lips become independent. It adds to the model a precise discontinuous displacement location that lacks in damage mechanics. This model is detailed in the dedicated chapter II.

I.3 Thesis objectives and outlines

The global goal of this thesis is to illustrate the capabilities of the TLS model in the quasi-static 3D simulation context of quasi-fragile materials. From an implementation point of view, switching from 2D to 3D has not been a significant effort.The TLS library created during the work of [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF] and incorporated in eXlibris later on, was from scratch designed to deal with 2D or 3D problems. However, due to 3D discretization, computation consumption becomes an issue. Behind this computation consumption, the causes need to be distinguished from the consequences. The causes are the discretization choices induced by the model, which can be optimized (which mesh size, where to adapt, which polynomial order to choose for interpolation and where to apply it, ... ). The consequences are the resolution method choice to handle very large system of equations. In this work, both are addressed (causes in chapters III and IV, consequences in chapter V).

In both cases, the parallel paradigm needs to be considered. Parallelism is nowadays provided by almost all computational hardware.

From a laptop to an enormous cluster computer, it offers a way to reduce simulation elapsed time and to overpass memory consumption issues. This thesis is focused on the scalability of the TLS problem resolution. The strong scaling will indicate if significant problems are achievable in a reasonable time with the parallelism. And the weak scaling will indicate if more significant problems can be solved using the parallelism.

The outlines of the thesis are the following:

Chapter II starts by describing the TLS model used in the context of this thesis. Chapter III introduces an improvement of the description of this model based on the "double cut" algorithm. This tool meets the thesis goal by opening the spatial discretization, treated in chapter IV, to coarser meshes. The main TLS resolution aspects are treated in chapter V. Then a set of test cases in chapter VI provides an insight on results that can be obtained in a quasi-static quasi-fragile 3D simulation with TLS model using some of the numerical improvements proposed in this thesis. Finally, chapter VII is the conclusion of what has been demonstrated. It also describes as a prospect, what remains to be done to value this work.

Chapter II

The TLS model a bridge between fracture and damage

The Thick Level Set model sets a bridge between fracture and damage mechanics by mixing a damage model with an explicit crack representation. The former offers complete freedom on how crack may initiate, where it progresses or forks and how it coalesces. The latter, using enrichment techniques, adds a strong discontinuous behavior for the crack opening . This chapter presents the TLS ingredients:

• φ the signed distance function used, with the help of level set tools, to split the domain in sound, damaged or fully damaged region,

• d(φ) the function provided by the user to associate d damage variable to φ,

• l c the characteristic length that is used, such as in other damage models, to regularize the damage energy release rate but here, only along the φ gradient segment.

It is also used to automatically place the crack lips when they do exist, giving support for the enrichment and the discontinuous behavior.

In this chapter, the governing equations are also provided, showing that d is used in a rather conventional way, with the asymmetric constitutive law from [START_REF] Mazars | A description of micro-and macroscale damage of concrete structures[END_REF]. Moreover, as mentioned above, the damage evolution law is regularized to avoid some spurious damage localization. In the end, the resolution scheme is summarized in the presented staggered algorithm.

18

This chapter slightly modified in this thesis have been published in [START_REF] Salzman | On use of the thick level set method in 3d quasi-static crack simulation of quasi-brittle material[END_REF] II.1 The Thick Level Set model components

In the TLS framework, as detailed in Mo ës et al. [2014], the domain Ω of interest is decomposed into three zones (possibly empty): a zone Ω -in which damage, d, evolves in a local manner, a localizing zone Ω + in which the evolution is non-local and, finally, a zone Ω c in which the material is completely damaged (d = 1). Nevertheless, in this thesis, we consider constitutive damage models without hardening that leads to a damage localization as soon as the damage starts. So, the damage is always zero in Ω -. This choice is done in [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF] and has been used in most of our TLS model implementation.

The interface between Ω + and Ω -is denoted Γ 0 and it is located by the iso-zero level set of a signed distance function φ (counted negative in Ω -). The "distance function" characteristic is achieved by :

∇φ(x) = 1 (II.1)
The damage is then expressed in terms of φ by a scalar function.

d = d(φ(x)) (II.2)
The function d(φ) is a material data that must satisfy a set of requirements, among which the fact that the damage reaches one, only beyond a distance l c . The interface between Ω + and Ω c is thus the iso-l c1 and is denoted

Γ c .
The use of a level set is a powerful aspect of the method which offers support for complicated Γ 0 , with potential branching and coalescence, as it will be seen in chapter VI.1. Another powerful aspect is the introduction in the model of this characteristic length l c which achieves two goals:

• The automated introduction of a crack, considering the iso-l c values of the level set, as the crack lips. In this thesis a more accurate discretization of Γ c is proposed in chapter III.

• The introduction of a maximum length over which an averaging of damage is performed.

l c Ω -: d = 0 Ω + : 0 d < 1 Ω c : d = 1 l c l c
Cut plan • A is a damage initializing zone where all points are at a distance less than l c from Γ 0 (in blue) and for which d < 1.

Ω + : 0 d < 1 l c l c A B C (a)
• B is a simple crack represented by Γ c (in red). In this case, one crack tip is located at a distance greater or equal than l c from Γ 0 tip. It illustrates a complex damage growth pattern in front of the crack tip.

• C is illustrating a branching situation where the crack first develops "horizontally"2 and then moves "vertically" 2 after a change of the loading conditions.

The damage shape function d(φ) must be

• increasing,

• continuous,

• bounded in [0, 1],

• null outside the damage zone, in Ω -(when φ(x) 0),

• equal to 1 in the fully damaged zone Ω c (when φ(x) l c ) Specific choices of d(φ) have permitted some comparisons with the phasefield method in

II.2 Governing equations

As in [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF] work the asymmetric constitutive law from [START_REF] Mazars | A description of micro-and macroscale damage of concrete structures[END_REF] is used here. This asymmetric behavior is mandatory to avoid damage growth in the compression direction. In this thesis, however, we do not deal with contact on the crack faces.

In Ω, the equilibrium and kinematics equations are:

           ∇.σ = 0 on Ω σ. n = f on ∂Ω N = 1 2 ∇ u + (∇ u) t on Ω u = ū on ∂Ω D (II.3)
where σ is the Cauchy's stress tensor, the strain tensor, f the external loading on part ∂Ω N of Ω boundary , n the outgoing normal vector of the domain, ū the prescribed displacements on part ∂Ω D of Ω boundary and u the displacement field. Here, the small strains and displacements are assumed.

The stress and the damage energy release rate Y are derived from the free energy density as follows:

σ = ∂ϕ ∂ Y = -∂ϕ ∂d (II.4) with: ϕ( , d) = λ 2 (1 -αd) tr( ) 2 + µ 3 i=1 (1 -α i d)Λ i 2 (II.5)
where:

• λ and µ are the Lam é elasticity coefficients

• Λ i are the eigenvalues of the strain tensor

• α i and α are coefficients introduced to take into account an asymmetric behavior in traction/compression:

-α i = β if Λ i < 0 -α i = 1 if Λ i 0 -α = β if tr( ) < 0 -α = 1 if tr( ) 0
• β is a parameter, bounded in [0, 1], to drive the asymmetry.

When β is equal to 1, α i and α terms are always equal to 1 and the free energy density is then a linear elastic damage potential (see (I.1)) When β < 1, the potential becomes non-quadratic. The damage participation to free energy density becomes negligible in compression when β tends to 0. When β is null, no damage growth is obtained in compression because, in this case Y = 0. The intermediate β values are allowed to fit the material properties where irreversible degradation slightly modifies the stiffness in compression.

For the damage evolution law, such as detailed in [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF], the local damage growth model:

Y Y c , ḋ 0, (Y -Y c ) ḋ = 0 (II.6) is regularized into its non-local counterpart: Ȳ Ȳc (II.7a) ḋ 0 (II.7b) Ȳ -Ȳc ḋ = 0 (II.7c)
where Ȳ and Ȳc are the regularized damage energy release rate and the regularized critical damage energy release rate, respectively.

An averaging operator is introduced to fulfill a non-local estimation of the energy (critical energy) release rate on the damage front. By construction, (II.1) with (II.2) impose damage along the φ gradient segment in Ω + (figure II.2). Thus, any point of Ω + may not have their damage value modified without changing the damage values of all gradient segment points belonging to it. It

II.3 Staggered algorithm

The TLS solver used for this work is given by algorithm 1 where µ i is the load factor corresponding to load step i.

Algorithm 1 Schematic staggered algorithm scheme.

Starting from u i , d i , and

µ i repeat Find d i+1 such that d i+1 = g µ i , u i (I) Find ( u i+1 , µ i+1 ) such that K d i+1 , u i+1 u i+1 = F µ i+1 (II) max k f k µ i+1 , u i+1 = 0 (III)
until Complete failure or user given load level

The g operator is related to the damage growth model (II.7). The K and F operators are related to the structural equilibrium. Finally, the f k operators are related to the damage criteria (II.7a).

In algorithm 1, (I) (corresponding to (29) in [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF]) uses an explicit resolution (as it is based on regularized quantities, the averaging operator resolution 3 has to be added to this step). In (II) the solution u i+1 is obtained through a Newton-Raphson algorithm, because of the non-linear nature of K. The condition (III) is solved with an explicit resolution. Solving (III) is eased by the fact that u i+1 linearly depends on µ i+1 (as long as it does not change sign).

The displacement field u i is discretized in space using a classical finite element approximation as well as a ramped Heaviside enrichment [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF]) to represent the existing cracks. The determination of Γ c is greatly improved compared to [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF] by the use of the so-called "double cut algorithm" detailed in chapter III.

The solution ( u i , d i , µ i ) at any load step is thus in equilibrium (in the finite element sense) and fully satisfy (II.7a) and (II.7b) even if (II.7c) is slightly violated (see [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF] for discussion).

The user does not give the load step. It is part of the solution process. Only the maximum damage (or φ to be more precise) increase over each load step is provided. This increase falls within the definition of g. The g operator is needed both to grow the existing damage zone and to initiate the 3 Averaging operator computation was not examined in this thesis despite that the slightly modified version of [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF], used in chapter VI simulations, adds two linear system to solve. In [START_REF] Moreau | Explicit dynamics with a nonlocal damage model using the thick level set approach[END_REF] the introduction of an explicit modal space replace one linear resolution by a projection. And in [START_REF] Moreau | Concurrent development of local and non-local damage with the Thick Level Set approach: Implementation aspects and application to quasi-brittle failure[END_REF], the use of Fast Marching techniques remove the other linear resolution. So in this work, even if this averaging operator resolution does have a computational cost, the use of those last techniques will hopefully alleviate it. new ones 4 . Regarding the initialization of the staggered scheme, we find the elastic loading for which the damage occurs first. It gives ( u 0 , d 0 = 0, µ 0 ). If needed, a non-zero initial damage field may be prescribed.

In algorithm 1, we have revealed an explicit dependence of K on u i+1 . Indeed, the asymmetric damage model, introduced in (II.5), yields a nonlinear elastic problem. This dependence of K with respect to u i+1 is, in fact, very limited in space since it affects only Ω + . Thus, the tangent matrix used to solve structural equilibrium is implemented, to spare CPU time, as two parts. A fixed part, which is assembled once per load step. And a varying part, which is reassembled at each Newton-Raphson iteration. The same strategy is also used for the Newton-Raphson residual vector, where the linear part of the internal forces is simply obtained by multiplying the fixed tangent matrix part by u. This early optimization does not change, as presented in chapter V.1, the fact that the non-linear problem resolution is the major computational bottleneck in 3D. The chapter V will address this problem.

Finally, the finite element mesh requirement for this algorithm in 3D is discussed in chapter IV.

Chapter III

The double cut algorithm or the art of describing a complex crack pattern On the one hand, the level set method allows representing a curve (2D) or a surface (3D) using an iso-value of a function. Generally this function, which is discretized as a field on the studied mesh, is used to identify a possible iso-value position over all edges. Those points that are also called edge cut gives, the iso-curve or iso-surface mesh according to an element pattern meshing strategy. On the other hand the Thick Level Set model uses such iso-curve/surface (iso-l c level set obtained from signed distance function φ) to represents the crack lips. Those lips can be very close compared to the element size and can cut twice an edge. In this case standard level set technology is faulty since it provides only one cut per edge. The new approach presented in this chapter allows cutting an edge twice. Firstly, a versatile tool, the vector distance function, that comes from the work of [START_REF] Gomes | The vector distance functions[END_REF] is reused. In this work a signed form of this tool is proposed and leads to the signed vector distance function concept (abbreviated to SVDF). It introduces extra information compared to the standard level set tool. Secondly, adapted to SVDF, a new algorithm, the "double cut" algorithm or DCA, gives a way to cut an edge twice. Associated with a database of element cut pattern (which is introduced for computation performance), this new algorithm allows passing a fully damaged zone in a mesh element without any condition on its size. It leads to a crack lips representation quite insensitive to the mesh size, allowing one to discretize over bigger elements.

This chapter slightly modified in this thesis have been published in [START_REF] Salzman | On use of the thick level set method in 3d quasi-static crack simulation of quasi-brittle material[END_REF] III.1 Double vs single cut algorithm

Usual treatment of the level set in the context of XFEM starts with the determination of the iso-zero location. This definition is then used to generate the integration cells (called sub-element hereafter) embedded in the elements crossed by the iso-zero. Depending on the type of problem, some or all those sub-elements are used to integrate the problem using the approximation of elements crossed by iso-zero. Some additional enrichment aspects may also be added and depends on this initial treatment. In the TLS framework this scheme is followed for both Γ 0 and Γ c . In this chapter we are going to focus on the iso-l c improvement which represents the crack discontinuity.

From a numerical point of view, if no enrichment is introduced around the iso-l c , the crack lips will move apart until an element is entirely in the damaged zone (figure III.1a). It allows obtaining the expected discontinuous displacement by having a fully damaged element inside the lips of the crack. However, depending on the size of the elements, more energy can be dissipated than wanted since the fully damaged zone depends on the mesh definition.

To avoid this issue, the ramped Heaviside enrichment described in [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF], offers a numerical mechanism to introduce the displacement discontinuity by adding extra degrees of freedom. These latter are added along the iso-l c of entities, having their support divided into two parts at least by the fully damaged zone. This condition is directly related to the Γ c definition that depends on the level set representation.

Since the level set is defined by algebraic values at the mesh nodes, all what can be done with them on an edge is their linear interpolation to find cuts. It implies that edges can only be cut once by the iso-l c . This restriction implies that the elements can not be cut in two parts by the fully damaged zone. In consequence, the enrichment will only occur if at least two elements are cut by the fully damaged zone (one element per crack lip). It imposes that the crack zone grows (or shifts) enough to start having a discontinuous displacement. It can again dissipate more energy than in reality (or can lead to a wrong crack path) if the mesh size is too big (figure III.1b).

The new approach allows passing over the restriction imposed by the standard level set method. The vector distance function, coming from the work of [START_REF] Gomes | The vector distance functions[END_REF], is reused and presented in a signed form, the signed vector distance function, in chapter III.2. It introduces at mesh nodes, a vector vision instead of a scalar vision of the mathematical function used to define the level set. It allows with the new algorithm called the "double cut" algorithm (DCA hereafter) to cut an edge twice. This allows a fully damaged zone to cut a single element (figure III.1c). This double cut algorithm first cuts all edges and then splits all elements. It is explained in chapter III.3.1 and III.3.2.

III.2 The Signed Vector Distance Function concept

Given a manifold, a vector distance function (VDF) gives at each surrounding point a vector pointing out to the closest point on the manifold [START_REF] Gomes | The vector distance functions[END_REF]). The manifold is also denoted as the VDF iso-zero in what follows. The manifold of interest in this chapter is Γ c .

To keep the domain partitioning available, a sign has to be associated with vectors (so we use the concept of signed vector distance function (SVDF) and not only VDF): a node is either in the negative or in the positive domain delimited by SVDF iso-zero.

In the TLS framework, the SVDF replaces the offset level set for the definition of the iso-l c . Indeed, the offset level set technology can not find multiple cuts within an element as depicted in figure III.1c.

In this context no intrinsic evolution of the SVDF is considered as it is the case in [START_REF] Gomes | The vector distance functions[END_REF]. It is only used as an enhanced tool to represent the Γ c manifold, reset at every load step using a new φ definition. It follows the same sign convention as the level set it replaces (i.e. positive for fully damaged zone and negative for all other parts of the domain). The SVDF can be constructed for each mesh node, from the level set that defines the front, by the following steps, illustrated in Figure III.2:

• Find the closest point (CP) on Γ 0 from the node (N).

• Compute the vector -→ V cp from N to CP .

• From this vector, construct the vector

-→ V lc = sign lc -→ Vcp -→ Vcp starting from CP
where sign is 1 if level set sign for N is "-" and -1 if level set sign for N is "+".

• Signed vector distance function vector

---→ V svdf for N is then -→ V cp + -→
V lc starting from N and ending at the point called L.

• Signed vector distance function sign for N is "+" if level set sign for N is "+" and

-→ V cp > lc. Otherwise, it is "-".
This procedure is the one used in this thesis. However, we can notice that a fast marching tool, introduced in eXlibris recently by Nicolas Chevaugeon, Γ c and SVDF iso-zero φ iso-zero (i.e. Γ 0 ) gives a way to obtain the SVDF more accurately using φ and its gradients computed during the marching.

+ + - ---→ V svdf -→ V cp -→ V lc + - - N N N CP CP CP L L L Figure III.

III.3

The Double Cut Algorithm (DCA) III.3.1 First step of DCA: cut mesh edges

From the SVDF information, one first has to compute the Γ c cutting points for all mesh edges. The conditions for having an edge AB cut, come from a logical consideration of the position of the edge with respect to the parts of the domain (table III.1).

The algorithm uses then the following general geometrical predicates:

• Geometrical SVDF iso-zero passes through the ending point L.

• Geometrical SVDF iso-zero tangent plane at point L is orthogonal to ---→ V svdf .

• Geometrical SVDF iso-zero tangent plane at point L corresponding either to node A or B of an edge may cut it.

- The algorithm to obtain the discrete SVDF iso-zero points on edges (i.e. Γ c points) is based on these geometrical predicates. Tangent plane intersection(s) is(are) considered as SVDF iso-zero point(s) under the following conditions illustrated in Figure III.3 (taking also into account the logical consideration of table III.1) :

• one cut cases (+-/-+) -The cut is first identified with the negative SVDF information (III.3b).

It agrees with the assertion that the quality of the negative SVDF information is potentially better than the one of the positive SVDF information (that are closer to skeleton location which is hard to track).

-If the above fails, the cut is based on the positive SVDF information (III.3c).

-If the above fails once again, a linear interpolation of the magnitude of the SVDF vector taken as level set values is done(III.3d).

• zero or two cut cases (--)

-If tangent planes do not cut the segment then there is no cut (III.3e)

-If tangent planes cut the segment twice but the order of cut point is not "A, cut from point A, cut from point B, B" then there is no cut (III.3f). This condition avoids some topological inconsistency (unwanted intersections) when following the SVDF iso-zero.

-If tangent planes cut the segment twice in the right order (see above) there are 2 SVDF iso-zero points(III.3g).

An important aspect to consider in this process is what should be done when the identified cut point lies close to A or B. First, a definition of "close" must be given. A cut point is close to an edge node when its abscissa1 s in the edge coordinate system is such that s ∈ [0, ε metric [ or s ∈]1ε metric , 1] where ε metric is an arbitrary small number (we choose 10 -5 in this work).

With a single cut algorithm, a simple treatment that modifies the level set is usually enough so that its iso-zero surface passes through the node. Nevertheless, it is no longer possible with the double cut algorithm because there is no way to identify the crack lips if we merge information at nodes. Thus, the cut node is placed on an edge node if it is close to it, but it still belongs to the edge from a topological point of view.

Regarding the double cut position on an edge, a similar proximity consideration must be considered. Do two cutting nodes on edge are close enough to be considered coincident? Once again the same ε metric is used to compare the abscissa (in the edge coordinate system ) s 1 and s 2 : if |s 1s 2 | ε metric then both cut points are metrically considered at the same location ( s1+s22 ) and topologically distinct.

Another case is also taken into account during the coincident node treatment. Independently of what might give the above plan cut algorithm, if the SVDF magnitude at a node divided by the length of the emanating edge is small (less then ε metric ) then a cut is considered to exist and is located at the node. To avoid different treatment around a node due to edge length fluctuation, we first compute the mean length of edges connected to a node. It gives a consistent metric reference across edges to evaluate the magnitude of the SVDF.

All these overlapping position identification are then used during the construction of the sub-elements for the integration and during the identification of the enrichment (see annex A and chapter III.4.1 for further details).

III.3.2 second step of DCA: cut mesh elements

Once the edge cut positions and topologies are determined, the construction of the geometric domain using the following assertions (also available in 2D) can be achieved:

• The convex hull generated by the cut points and positive node-set corresponds to a + domain (Ω c in the TLS framework). It is a polytope2 (polyhedron or polygon depending on the case).

• The subtraction (geometric operation sense) of this + domain to the original element gives the -domain polytope (Ω + in the TLS framework).

• The common boundary of the + and -domains gives the SVDF iso-zero (Γ c in the TLS framework).

The element cut pattern cases can be reduced to a small set of unique cases using the elementary topological rotations and symmetries. Depending on the way one chooses to implement the cutting algorithm, generation of the null volume parts (coming from the cut point merged with other nodes) can be eliminated from the integration or not. It leads to a more significant number of element cut pattern cases or not.

All those cases are depicted for the 3D case in annex A. This annex also presents the cutting database introduced to accelerate the DCA computations. Moreover, it presents a pure geometric illustration of the performance of the double cut algorithm with SVDF compared to the level set single cut algorithm. Note that the library that provides DCA has also been implemented in a distributed version, not presented here, but used in chapter V.4.

III.4 Discussion about the new Γ c representation

The use of the "double cut" algorithm for the discretization of Γ c in the TLS framework has two impacts:

1. On the way ramped Heaviside enrichment is computed. It will be discussed in chapter III.4.1 2. On the way φ values are evaluated in the Γ c vicinity. Chapter III.4.2 briefly comments on this aspect.

III.4.1 Enrichment impact

The fact that an element can be cut twice by the same iso-l c is of small influence on the process that finds the dof to be enriched in an element as described in [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF]. The general guidance is still to consider for a given dof its support and to count the unconnected parts from Ω + created while splitting it by Ω c (figure III.4).

As soon as there is more than one Ω + part, the enrichment must be used to allow the crack opening.

The "double cut algorithm" presented in this work handles with care the cases where the crack runs close to the nodes (chapter III.3.1 and annex A). It has an impact on the way Ω + parts are identified and counted. In figure III.4, the zoom shows a basic example of Γ c passing on a mesh node. In this case, a topological Ω + part exists but it is of zero measure. Those parts do not count for the enrichment decision. Subtle impact on Ω + parts identification is illustrated in figure III.5. Let's consider the mesh given in figure III.5a and A one of its node . In figures III.5c and III.5f, where Γ c does not run through node A, either we have 2 sane parts (III.5c) or 3 sane parts (III.5f) depending on whether the node A lies outside Ω c or not. In figure III.5d, Γ c runs through node A (from the above), there are two sane parts (such as in figure III.5c). In figure III.5e, Γ c runs through node A (from below) and there are 3 sane parts (like in figure III.5f). We remind that even though Γ c runs through a node, the information about which edges are cut around the node by Γ c (which enables the tracking of Ω + part definition) is kept.

III.4.2 φ in Γ c neighborhood

The "double cut algorithm" uses a signed vector distance function instead of an offset level set for Γ c construction. As exposed in the previous chapter and annex A, it gives a suitable Γ c discretization. Nevertheless, it introduces some differences in Γ c vicinity, between the φ value and the crack front location. discretization at a distance of 0.499. This value is in contradiction with the 0.51 φ value. It is possible because Γ 0 is a discretization of the φ iso-zero that can induce a difference with the true iso-zero (in dot). The consequence in this case is that the SVDF sign at this vertex is negative! Consequently, the DCA uses the '--' rule of table III.1 and there are only cuts at the bottom edge. The Γ c appears then on it. Now, if considering φ value instead of -→ V cp the upper vertex would have been in Ω c . This problem is treated modifying the φ value on this vertex (for example, by using neighborhood values) so that it is not greater than l c as shown in figure III.6d. Note again that the SVDF fast marching alternative raised in chapter III.2 would suppress such difference (and would give a better vector orientation).

The other issue (not linked to the SVDF construction) is related to integration. Figure III.6c presents in green a sub-element introduced to take into account the Γ c presence. As the integration is made at the Gauss point G for this sub-element, one must know the φ value at this point. There are two ways to interpolate φ on G : Clearly in both cases, using the φ level set will not give a computed value corresponding to the presence of Γ c . For example in figure III.6c, at the Γ c edge nodes locations, φ should be 0.49 and 0.494 from geometric linear interpolation of element level set values.

The idea, to be always consistent with DCA construction, is to consider that the level set values on Γ c nodes are imposed to a l c value. Then, the computation of φ on G, using the geometric linear interpolation of the subelement modified level set values, will be in accordance with the presence of Γ c as shown in figure III.6d.

These two modifications lead to the fact that φ is no longer a distance function in the vicinity of Γ c . It implies that somehow d() function response is not anymore what we expect in this region (It is locally another d() function). But the damage is now 1 on Γ c and no spurious fully damaged location does appear anymore.

Chapter IV

Mesh adaptation strategy

Let us recall h t , h w and h o definitions as respectively the mesh size in the crack tip vicinity, the mesh size in the crack wake region and the mesh size outside damage band. With the Thick Level Set model, the discretization requirements are of two kinds: the accurate capture of the mechanical behavior at the crack tips (that imposes h t to be roughly equal to lc 5 for a right approximation of the strain field) and the right discretization of the crack path (that imposes h w to be between h t and h o depending on the crack wake curvature for a correct integration computation). h o must only fulfill the general mechanical constraint over the domain (such as boundary condition, fillet, hole , ...). In a 3D simulation, the lazy option that imposes h t everywhere or in a rather large defined region, leads to a huge system of equations that will be hard or impossible to solve on a computer (because of lack of memory or irrelevant computation time). Some specific computational strategies can help as exposed in chapter V. But adapting the mesh size is almost essential to reduce the effort of the linear algebra computation. A first easy adaptation consists of tracking the damaged band and imposes within it, h t (with h w = h t ) and h o elsewhere. It has been studied in this chapter with an octree, unstructured and two-scale embedded mesh. This last mesh is related to the twoscale strategy exposed in chapter V where adaptation is made at a fine-scale level.

An alternate adaptation strategy that sticks to the TLS requirement given above, consists of using h w such that h t < h w ≤ h o . In particular, the "double cut" (chapter III) algorithm capability will be advantageously used to fulfill this new requirement in the octree mesh context.

40

In the TLS model, Ȳ computed from the strain field, needs to be of excellent quality to provides the right crack progression. As the strain field is singular in the crack tip vicinity, a suitable discretization in this region is mandatory. That is the reason why h t must be roughly equal to lc 5 to capture the strain singularity accurately [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF]).

In the crack wake, only a correct material integration and a proper enrichment do mater. Both are related to ρ and crack complexity. If the crack curvature is small, the discretization must be rather precise to capture lips geometry correctly: h w is close or equal to h t . If ρ is large then the crack is almost straight and h w can be equal to h o . The same holds for branching or merging scenario where the lips are described with sharp angles.

In the region outside the damaged band, h o needs only to fulfill constraint imposed by the mechanical problem: local Neumann boundary condition, holes in the part, fillet ... In this work, h o is used with a constant order of magnitude. However, in reality, the mesh size is adapted around the holes, the fillet, ... This adaptation will be considered as a fixed input in this chapter.

The TLS mesh size requirement can then be summarized as h t ≈ lc 5 h w h o with h w being a parameter to adapt.

On the one hand, one may implement the TLS model without dynamic mesh size adaptation. In this case, he will be forced to impose h t everywhere or in a rather large defined region. Nevertheless, no re-meshing adaptation nor φ transport has to be done. In 3D, it will lead to a system size difficult or impossible to solve on a computer ( lack of memory or irrelevant computation time). However, all simulations of chapter VI have followed this principle. It has required a specific solving strategy, called active zone (see V.3), to partially bypass memory issues and time consumption.

On the other hand, if the mesh adaptation tools are available, a good choice of h o and h w can reduce the size of the system. The next chapter, proposes a mesh adaptation strategy that tunes h w taking h o (fixed adaptation) has an optimal input.

IV.1 Adapting mesh in damaged band

It corresponds to the simplest but not the most efficient strategy. h w is arbitrarily taken equal to h t . It is simple because it doesn't require algorithmic efforts to separate the crack tip vicinity from the crack wake. Only damaged band localization is requested and φ gives an easy way to do so.

IV.1.1 Octree context

It is the TLS historical adaptation strategy . It uses eXlibris Octree library that provides an effortless coarsening/refinement of cells in a cube or arrangement of cubes. The octree implementation uses the Linear octree concept [START_REF] Gargantini | An effective way to represent quadtrees[END_REF]). It stores the cell status in a linear array instead of a tree data structure. The array size is fixed by the maximum refinement level imposed by the user. At level 0 only one cell is defined: the original hexahedron (or set of hexahedron). At Level 1, 8 cells are corresponding to the first cube split in 8 equal hexahedron. At Level 2 there are 8 2 cells, and so on.

The library adaptation functionality takes as input a set of quadruplet: a level set, a level to reach in the octree, an offset from iso-zero of the level set and the level set side to use for refinement. The octree cells at the appropriate level are activated following those parameters with an imposed 2:1 constraint: a cell at a level level i can not have a neighboring cell with a level higher than level i+1 or lower than level i-1 . It ensures a priori smooth mesh size transition and only one hanging node on edges or faces (annex C).

The TLS implementation uses this library via an octree/unstructured mesh conversion tool since all cutting and enrichment features are acting only on tetrahedrons. The unstructured mesh here, as inherited from an octree mesh, is mostly a regular mesh turned into simplex. The conversion tool transports back and forth the level set function from an unstructured mesh to/from an octree mesh. When passing from an octree to an unstructured mesh, the hanging nodes are treated either with supplementary linear kinematic equations or additional connecting elements (annex C). Of course, the level set function used for the octree adaptation is φ. In Ω + Ω c , we impose h t and enlarge the refined zone using an l c percentage offset. This offset is essential in this strategy since the adaptation does cost (mesh creation, octree definition, level set function transfer, ... ). By having this margin around Γ 0 , the front can progress during several load steps without the need to reshape the unstructured mesh. The only constraint is that Γ 0 remains in the refined area. As soon as this condition is not fulfilled anymore, a new adaptation is required. In the future, the offset has to be tuned. For now, it is taken all around Γ 0 which is unnecessary most of the time. An "Only where Γ 0 will progress" criteria can tune the offset appropriately and can reduce the problem size slightly. However, this "Only where Γ 0 will progress" criteria is complex to tackle because it implies to know where the crack will progress before computing its progression.

Many simulations, not presented here, did successfully use this adaptation but the octree strategy is only available for parallelepiped shaped objects. Mentioned by members of the team, to bypass this limitation, a level set function can be used to define objects embedded in a cube. However, it only gives a discrete definition of its boundary (the geometric representation depends on the maximum octree refined level). It also requires a complex implementation to deal with an element cut by several level sets. Alternatively, using a mapping from/to a cube to/from a specific shape can work for a simple object. It has to be tested. In the end, coping with this shape limitation is investigated in the next chapter, transposing the octree refinement on an unstructured mesh.

IV.1.2 Unstructured context

The same principle as the octree refinement with a 2:1 constraint (see previous chapter) is used in this unstructured distributed mesh adaptation tool. Annex D provides the full parallel algorithm 5 and the details of the used refining procedure. This tool takes as input a mesh and a criterion indicating whether an element should be refined or not. In the context of damaged band mesh adaptation, this criterion is a function that asks to refine the element if it is in or partially in Ω + Ω c and if its size is greater than h t . As in the octree context, adding a smart offset in this criterion would be a good mechanism to reduce the number of adaptation. However, it as not been investigated yet.

In terms of scalability load balancing is a crucial point (as illustrated in annex D with a dedicated test). If the initial mesh is always distributed so that each process gets elements to split, the load may be balanced and the scalability will be present. Besides, in this rather crude tool, the splitting history is not kept. It limits the operations only to the refining from the original mesh: neither successive adaptation nor coarsening are possible. A future development can provide them and will also potentially ameliorate the scalability. However, for now this implementation provides an experimental adaptation tool to quickly test the assumptions formulated in this work using the appropriate mesh. Furthermore, the next chapter shows, in particular, how this unstructured adaptation tool has been used in the two-scale context.

IV.1.3 Two-scale context

The two-scale method is presented in chapter V.4 but it is interesting to recall, in this part, its discretization mechanism. At this point, all we need to consider for the two-scale method is that it needs to have a fine mesh embedded into a coarse mesh. A fine mesh element must know in which coarse mesh element it is embedded. Reciprocally, a coarse mesh element must know all its embedded fine mesh elements. This two-scale method applied to the TLS can use the adaptation criterion suggested in this chapter: h w = h t . In fact, the unstructured adaptation tool of chapter IV.1.2 may be used with a coarse element set to create all the attendees described above. The procedure may be summarized as follows:

1. Selection of the coarse elements that cover the damaged band. It will be, hereafter, named the C element set. An SVDF description of φ on coarse mesh is needed for that. Eventually paired with DCA, it will allow determining if the damaged band crosses a coarse element or not. Note that a small isolated damaged band can be untraceable by the DCA. It has to be investigated in the future.

2. Distribution of C using ParMetis weighted partitioning (an element that covers a more significant part of the damaged band gets higher weight). It is mandatory for a proper load balancing of the adaptation(as mentioned in annex D.3) as well as of the two-scale resolution (chapter V.4.3)

3. Use of C to create a new mesh, the SuperPatch, that will correspond to the fine scale. A connection between the C elements and their equivalent in SuperPatch is also created.

4. Adaption of the SuperPatch using the tool and criteria of chapter IV.1.2.

A transport function has to be added to maintain a connection between the fine mesh elements newly created and their ancestors in C.

5. Creation of the reverse connection from C to SuperPatch using the connection between the SuperPatch elements and their ancestor in C.

Compare to the adaptation in chapter IV.1.2, the main difference, in this case, is the presence of two data structures to store the initial and final (refined) meshes instead of only one transformed mesh . It has two consequences. Firstly, it allows launching a new adaptation from a coarse level easily: only the fine data structure has to be cleaned. Secondly, it provides different results since the 2:1 constraint does not propagate in the same way for the two approaches. In this two-scale context, the SuperPatch is limited to the C set. Hence, the neighboring elements of C are not impacted by 2:1 constrain as they are in the adaptation of chapter IV.1.2.

IV.2 Coarsening mesh in the crack wake

This strategy proposes h t in the crack tip vicinity, h o outside the damaged band and h w in the crack wake with: h t h w h o . Three questions have to be associated with this new approach. The first one, fairly obvious, is what value should be specified for h w and should it be uniform in the crack wake? This question is not addressed in this work. Nevertheless, considering

h w ≈ max h t , min h o , ρ L
should be a good starting point. The second one is related to the crack tip vicinity definition itself. Where should vicinity stop and where does wake start? This question is treated in chapter IV.2.1. The third one is what tools are mandatory to represent the TLS framework correctly in the coarsened crack wake? This last question is related to the Γ 0 and Γ c description. If we want to keep both of them in the coarsened crack wake, a standard level set tool is not appropriate. The DCA (chapter III) will let Ω c pass inside an element. It lets h w grows up to a limit where an edge is cut more than twice (2 times by Γ c and one or two times by Γ 0 ). So with current tools, we must have h w < l c without further treatment.

IV.2.1 Identifying the crack tip vicinity

A priori, the crack tip vicinity can only be described if the tip is located. The Γ c discretization provides the tip and lips. However, where is the tip on this curve (in 2D) or surface (in 3D)? For now, no topological information is given by the tool that creates Γ c . Investigation, not presented here, based on graph manipulation and zoning, gives a bad approximate of the tip localization on Γ c . Besides, the computational effort is important.

Alternatively, since only the crack tip vicinity does matter, why trying to isolate the tip accurately? Knowing some points, close to the crack tip, can be enough. In this work, the enrichment information around Γ c are used to find those points as presented in annex E. A quick analysis based on the current mesh and the TLS framework location, gives the elements which are close (at least in contact) to the crack tip. The points sought are then only the centers of gravity of these selected elements.

Figures IV.1a and IV.2a show, in a 2D context, two situations where point locations are given by enrichment information around Γ c (see annex E for the steps that select those points). It is then possible to define a tip neighborhood by: constructing spheres (3D) or circles (2D) using those points as centers (figures IV.1c and IV.2c); taking their convex hull (see the areas with the white stripes in figures IV.1d and IV.2d). The sphere (or circle) radius is intended to be large enough so that the damaged band front is covered. It imposes a radius at least greater than l c . But reducing the adaptation number, as in chapter IV.1, is possible by taking a larger radius value. It allows staying with the same mesh for several times steps as long as Γ 0 remains in the refined zone. Once having those refined regions, it is not difficult to remove them from the damaged band. Figures IV.1e and IV.2e show the resulting regions where a coarsening will be done with a h w element size target. The adaptation tool will ensure that 2:1 constraint is valid. Furthermore, the DCA will permit the construction of Γ c in those mesh size transitions. Finally, figures IV.1f and IV.2f show the mesh that will not be adapted (excepted marginally to ensure the 2:1 constraint). Note that this zoning is better than the offset strategy presented in chapter IV.1, even without coarsening. To limit the phases of adaptation, the additional refined elements are applied only in the vicinity of the crack tip.

Note that some care is required in the absence of cracks (during the initiation phase).In this case, no point close to the tip can be located because there is no tip. The workaround solution is to switch to the damaged band adaptation presented in chapter IV.1.

Finally, as discussed in annex E, the forking case presented in figure IV.2 shows that the proposed algorithm also identifies a region where the crack lips are not close to each other, as a tip vicinity. In this specific example, it is The first one corresponds to a plate (2D) under a shear loading. The boundary conditions correspond to tearing the plate in a specific direction (30 degrees from the vertical figure axis). It leads to a crack path not aligned with any edges of the mesh. The simulation is started with a seed (little damaged circle) at the plate center since the problem is initially symmetric. Then, the damaged band grows without an explicit crack. The octree adaptation imposes h t in the band. When the damage front grows enough, Γ c is introduced and the crack tip identification starts. On the first adaptation following this event, the tips at both ends are too close and the refined zone covers all the band. Then, later on in the simulation, both tips are far enough so that a coarsening region appears as depicted in figure IV.3.

The second simulation is a three-points bending test with two initial notches at the bottom of the beam. The load is applied at the top center of the beam. One notch is centered (just below loading ) and the other is shifted from the center. In the beginning,both notches, represented by the TLS framework, are of the same length. The first adaptation covers both. The two notches start growing and the shifted one stops obviously while the centered one continues toward the top of the beam. In figure IV.4, the tip of the centered crack is far enough from the bottom of the beam so that coarsening occurs.

During those two simulations, slightly modified Γ 0 and Γ c positions have been observed in the coarsened regions. It is related to level set transport during the adaptation and this aspect must be investigated in the future.

The work presented in this chapter is more a proof of concept than a clean quantitative evaluation. This one has been postponed to focus on the a priori more profitable solutions studied in chapter V.

IV.2.3 Other experiments

Note that an anisotropic mesh adaptation [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF]) has been tested with the TLS model by other team members. The idea was to replace the isotropic coarsening in the crack wake by an anisotropic meshing. It was done with the ICI-lib from ICI-Tech suites (ICI laboratory at ECN). The control of ρ and more precisely the curvature of Γ 0 , give a parameter to adapt the mesh in all regions. At the tip, Γ 0 has a high curvature and the elements are kept isotropic. In the crack wake, Γ 0 is more straight and the elements can be extended in the crack path direction reducing their numbers in this direction.

The tests were encouraging but some level set transport from one library to the other added numerical issue. It would make sense to investigate this path again.

Chapter V

Parallel and multi-scale solving strategy: distributed TLS approach

We show, in this chapter that in the implementation of the Thick Level Set model, the computational throughput is mainly linked to the performance of linear system resolution, as in many finite element applications. The system creation and the damage front advance are, respectively, second and third CPU time consumer. Two general approaches are commonly used to solve a linear system: The "frontal" strategy that takes the problem as it is and uses direct or iterative solver to obtain a solution. The second approach corresponds to the "indirect". In this case, this strategy tries upstream manipulations to reduce the problem size before using a direct or iterative solver. In this chapter, both strategies have been investigated and parallelism has been chosen to chase High-Performance-Computing with a high scalability objective. Regarding the "indirect" strategy, among all available possible solutions, two numerical methods have been studied: the domain decomposition method with a specific application called the "active zone" and the two-scale or GFEM gl computational framework. The first one reduces the system size using a linear algebra condensation technique. The second one keeps low the system sizes by working on a global discretization (relatively small) which is enriched by some local fine-scale information. Both have been chosen for their potential scalability. Even though every thing is not fully stated with solving strategies proposed in this chapter, positive preliminary results comfort the numerical and implementation choices. V.1 Initial profiling and strategies for improvement

Sequential TLS profiling shows that the greatest computational effort is related to the nonlinear resolution ((II) in algorithm 1). This iterative resolution (Newton-Raphson) can be split into two main tasks: the system of equation creation (definition of sparse matrix structure, integration and assembly of tangential rigidity matrix and right-end side residual vector) and the system of equation resolution (direct solver factorization and solving at each Newton-Raphson iteration). The computational performances of these tasks and of the damage updates (solving of equation (II.7) and φ updating) are presented in table V.1. They correspond to a chalk test case simulation (chapter VI.2 ) in a sequential computation context. It is representative of 3D elapsed time consumption. This table gives results for load steps 2 and 84. In both cases, the damage front advance computation is negligible. At load step 2, the simulation begins and a small damaged zone is only initiated at the chalk edge. The Newton-Raphson resolution is done in 2 iterations and the linear algebra resolution is proportionally the most important consumer. The assembly does not cost much because the non-linear part (damaged band) is small and is only computed twice. At load step 84, the crack has reached the chalk center and the Newton-Raphson resolution is done in 12 iterations. Consequently, the assembly task is proportionally more important than in step 2. There are more elements to assemble and this task is done more time. However, the linear algebra resolution remains the hot spot. Note that this last assertion was made on the basis of the resolution of linear systems with a direct solver.

To decrease the time consumption, one first solution consists in reducing the number of Newton-Raphson iterations by "simplifying" non-linearity. This point, not investigated here, is related to the asymmetric non-linear potential(II.5). Smoothing this law at tension/compression switch may ease the nonlinear resolution (team proposal). Now keeping the same non-linearity,the Newton-Raphson time consump-tion can be reduced by modifying the update frequency of the tangent matrix. With the modified Newton-Raphson or other variants, the tangent matrix factorization is done only a few times (depending on the modification applied to the original Newton-Raphson algorithm). Then, between these factorizations, only cheap forward and backward substitutions are used to solve the linear systems. It has been tested in a three-points bending simulation with an update of the tangent matrix every arbitrary N iterations. It is not present here but the results were not positive. The nonlinear resolution takes more iteration to converge decreasing the benefits given by reducing the number of factorizations. Nevertheless, it might be interesting to continue digging this track.

If the nonlinear resolution is not modified, it remains to decrease the computational time of the linear systems resolution. Two modus operandi have to be considered for that. The former is call "frontal" since it takes the linear system as it is and tries to find the fastest algorithm to solve it. The latter is called "indirect" because it first tries to reduce the problem size and then uses an adapted solver to solve the reduced system. The "frontal" approach is described in this chapter as the frontal strategy presented in chapter V.2.

For "indirect" approach, many strategies are available. The domain decomposition methods are good candidates. As stated in [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF], these techniques "revolve around the principle of divide and conquer". Stepping in the resolution process very early, they divide the initial problem in a set of smaller ones. Those smaller problems, called domains, contribute to the global resolution by providing computed information individually. This domain treatment independency gives excellent scalability properties to those methods. In this thesis, non-overlapping domains are considered (eliminating Schwarz alternating method which seems complex with enrichment) and a element-base partitioning is chosen (the vertex or edge-base partitioning are also eliminated). This last choice is strongly related to the eXlibris tools capabilities in a distributed context. A first basic direct Schur complement resolution with two domains is proposed for the TLS model in Chapter V.3. Then in chapter V.5, an iterative Schur complement approach is benchmarked (however not with the TLS model) against the other strategies proposed in this chapter.

Alternatively, an "indirect" approach can reduce the problem size by the use of multi-scaling strategies. The most basic form has been investigated in previous chapter IV where the mesh adaptation copes with the different problem scale phenomena. Besides, mentioned in the introductory chapter I.2, a crack box can be seen as a specific scale where the discontinuous crack behavior is well captured. The rest of the part, embedding the box, corresponds to another scale only able to reproduce a overall mechanical behavior. Another technique, the multi-model and multi-scale Arlequin frame-work [START_REF] Dhia | Problemes mecaniques multi-echelles: La methode Arlequin[END_REF]), introduces a generic way to superimpose a fine-scale model to a large-scale model combined at the energy level by coupling operator (often Lagrange Multiplier). The first attempt with this technique proposes to follow the crack evolution with a fine mesh covering it. Then, using XFEM, in Ben [START_REF] Dhia | On the use of XFEM within the Arlequin framework for the simulation of crack propagation[END_REF], only the crack tips vicinity is tracked by a fine-scale problem. At the global level, the enrichment with a Heaviside function ensures the discontinuous behavior in the crack wake1 . A similar approach can be applied to the TLS. Despite the complexity linked to the Arlequin/TLS mixture, it is not clear whether the above proposals would have had the right scalability property. By nature, the Arlequin system structure is well suited for the domain decomposition-like solvers. But that implies that the crack tip vicinity (in 3D it's a pretty big problem in itself) has to be treated efficiently (in many domains in parallel) which is not obvious from bibliographic study. Around the same period [START_REF] Pereira | A two-scale approach for the analysis of propagating three-dimensional fractures[END_REF] propose, for the same kind of simulation, a multi-scaling with the use of enrichment to "connect" scales. The associated two-scale (TS) or GFEM gl computational framework (first introduced in Duarte and [START_REF] Duarte | Analysis and applications of a generalized finite element method with global-local enrichment functions[END_REF]) links a global scale to a set of overlapping local (or fine) scale problems ( called patches). This last point offers a priori good scalability properties [START_REF] Kim | Parallel simulations of threedimensional cracks using the generalized finite element method[END_REF]) despite suffering from redundancy. This justifies studying how the TLS model can be introduced in this two-scale framework. Chapter V.4 presents this work.

V.2 Frontal strategy

This strategy may be summed up as choosing the fastest appropriate linear solver for the TLS non-linear resolution. To make the appropriate choice, the linear system matrix nature must be analyzed. The tangent matrix is a sparse symmetric positive definite matrix. This is interesting because it involves the use of a simpler and faster resolution algorithm compared to an indefinite matrix. However, the matrix conditioning is not very good due to the ramp Heaviside enrichment. In this thesis, this last point has contributed to eliminate the iterative solvers based on Krylov subspace from the possible solutions. Those types of sparse solvers, although they provide good scalability properties, are very sensitive to the matrix condition number. A vast research field is dedicated to preconditioning techniques that improve the solved matrix properties and then the iterative solver performance. A specific work on that aspect would be necessary to find the best preconditioner for the TLS model. Once done, other investigation about Newton-Krylov methods would be fruitful to really gain performance at the non-linear resolution level. Note that the multi-grid method may also play a role in this kind of solver. Finally, mixing the direct and iterative solvers in Newton-Raphson resolution might give some results easily (the direct solver computes, at some non-linear iterations, the costly tangent matrix factorization (eventually to be fast with low precision : MUMPS BLR) and this factorization is used as a pres-conditioner for all linear resolution done with iterative solver. It was successfully tested a few years ago at EDF R&D by Olivier Boiteau (no citation here sorry, we just talked together)).

After eliminating the iterative solver, the direct solver, less sensitive to poor conditioning is then the key solution. Many libraries are providing good algorithm for such sparse symmetric defined matrix: Su-perLu (https://portal.nersc.gov/project/sparse/superlu/), Taucs (https: //www.tau.ac.il/ ~stoledo/taucs/), Pardiso (https://www.pardiso-project. org/), PaStiX (https://gitlab.inria.fr/solverstack/pastix), MUMPS (http: //mumps-solver.org/), ... The results of table V.1 show that a sequential approach is not adapted to the kind of problem addressed in this work. So parallel versions of those codes are a must. However, different parallel paradigms are proposed. For example, Pardiso proposes a multi-threaded solution that implies a computer with lots of unified memory connected to many cores (even though Pardiso developers switch now to hybrid direct/iterative distributed resolution). The message passing protocol with or without multithreading is, on his side, more hardware flexible ( it can be used indifferently on a set of personal computers networked together or a cluster interconnected with high bandwidth network or ...). For example our target platforms (annex J.2) are more adapted (especially Liger) to the message passing protocol. In this regard, MUMPS, PaStiX or SuperLU are better candidates. They all propose a version using MPI with multi-threaded or GPUs hybrid computation. Mumps has been retained for its well-known performance, its large set of features and its maintainers' group size (guarantee of longevity). Now, how is it possible to introduce the MPI parallel solver into the TLS framework sequential implementation? A first try was to isolate the sequential code in process 0 and span, with a master/slave paradigm, solver steps in all processes. It gives quickly a way to scale the linear resolution but it has a limited action: the system creation becomes the bottleneck. The MUMPS solver offers a distributed input capability. The next move was therefore to fully distribute the TLS resolution using an element-based partitioning. It allows the system creation tasks parallelization which well scales (quite ideally if load is well-balanced). This move has not been achieved yet. Because many algorithms used in the TLS framework, take time to be adapted to MPI communication (such as the double cut algorithm of chapter III, for example). Nevertheless, this goal is slowly getting closer. This strategy though, for the non-linear resolution, is only relying on the MUMPS scalability to perform. A recent version of this solver proposes multithreading capabilities as well as block low rank (BLR) approximation [START_REF] Amestoy | Improving multifrontal methods by means of block low-rank representations[END_REF]): they push both the performance at a high level. A point in figure H.8c gives insight of those promising feature, with the uniform benchmark test (with 64 MPI process, elapsed time for MUMPS resolution has been divided by 3.5 through the multi-threading). However, for sure, BLR must be studied at some point. Chasing access to higher performance, the condensation techniques to reduce the problem size have been tested in the next chapter V.3 with the active zone strategy.

V.3 Active zone strategy V.3.1 General concept

In the crack wake, the damage evolution is not expected unless the boundary condition does impose drastic modifications or unless another crack joins from elsewhere. If the damage does not generally change in those regions, it means that the mechanical properties are locally constant (the lips are at the same location and d is then unchanged). This said, one can try to isolate the crack tips to focus the computation effort on these variable zones (called Active Zones or AZ hereafter). It meets the concern of chapter IV.2: finding the crack tip vicinity. So naturally, the same enrichment base algorithm 6 (annex E) is used to locate the points in the crack tip vicinity. The simple example of an artificial TLS framework applied to a 2D mesh is presented in figure V.1. As the adaptation, a generic envelope is constructed base on a circle (2D) or a sphere (3D) (figure V.1c) centered on those points (figure V.1a). In figure V.2, this envelop is presented in a 3D simulation context. Using this envelope, the elements are separated into two groups :

• the AZ group ( orange elements in figure V.1d which are cut or inside the AZ envelope)

• the fixed group ( green elements in figure V.1d)

The idea is to construct an AZ group for a certain amount of load steps so that the damage front may progress into it. As soon as the Γ 0 moving part reaches the AZ boundary, a new active zone has to be considered. However, during all this progression, information of the fixed zone (φ,d,...) are considered constant and their computational cost can be optimized.

The larger is the AZ group, the higher is the number of load steps that benefit from the optimization of the fixed group (but the smaller is this group and the greater is the computation for the AZ).

In this work, the fixed and AZ groups are considered as two domains where decomposition techniques are applied. The Schur complement contribution and the condensed right-hand sides for both domains are computed. Then, Schur complement reduced linear system is solved. It represents a part of algorithm 1 resolution II. However, during all the AZ lifetime, the fixed group is considered constant. So this domain can be condensed only once at AZ construction. Then fixed domain linearized form, the Schur complement contribution, is simply added to reduced problem at resolution. In this way, the factorization effort for the fixed group is amortized over many load steps and non-linear iterations. The AZ domain itself remains non-linear and is subject to the φ variation. Hence, its condensation is systematic for each Newton-Raphson iteration during its lifetime. It gives the new stagger algorithm 2 where the ˜corresponds to the condensed space and G is the expansion operator (from the condensed to the full domain).

Algorithm 2 Modified staggered algorithm scheme with AZ.

Starting from u i , d i , and

µ i repeat Find d i+1 such that d i+1 = g µ i , u i Eliminate fixed group from problem for any new AZ Find ( u i+1 , µ i+1 ) such that        K d i+1 , ũi+1 ũi+1 = F µ i+1 u i+1 = G ũi+1 max k f k µ i+1 , u i+1 = 0
until Complete failure or user given load level

In algorithm 2,it is in the application of the g operator that the Γ 0 evolution is checked against the AZ boundary. Any progression, close to this boundary, launches a new AZ construction. Note also that when a damage initiation outside the damaged band occurs a small damaged sphere is added to φ and a new AZ must be created to encapsulate this extra zone.

Like in frontal strategy (chapter V.2) message-passing paradigm is chosen to gain computational performance with parallelism. In the current AZ strategy, only two domains are used. So parallelism will not be, like in many domain decomposition methods, related to domains parallel treatment. Instead, the parallelism will be employed to reduce the algebraic resolution cost of both domains and the reduced system. MUMPS is again chosen since it provides excellent performance for Schur complement creation by simply adapting its parallel multi-frontal factorization to an incomplete block factorization. It also nicely gives Schur complement dense matrix in a distributed cyclic block format compatible with the Scalapack library. It simplifies the manipulation of reduced problem and drives to use the Scalapack distributed dense solver to resolve it. The load balancing induced by parallelism is treated in this thesis, with a specific partitioning that tries to stick as much as possible to the AZ concept. It is depicted in annex F.1.

Note that an extra feature not presented in this thesis has been tested. It corresponds to the use of a mixed shape function approximation order. One may want to use a higher-order to solve the problem. With the concept of focusing the computational effort on crack tip, it is rather natural to restrict this higher-order to the AZ group. An order 1 for the fixed group is enough as any way it is only used in a linearized way. Compared to a full higher-order problem, gains in computational time are expected. Though this requires that the size of the common frontier between the AZ and the fixed groups is kept at order 1 (i.e. order transition must be done in the AZ boundary element) to get the smallest possible size of dense problem.

V.3.2 Feed back

This strategy has been used in many test cases of chapter VI. Note that none of these simulations uses a mesh adaptation strategy from chapter IV: a fixed region using h t mesh size covers a priori the crack path. This choice was guided by the wish to eliminate any mesh adaptation interference in 3D TLS model validation ( independent damage initiation, no level set transport interference during the adaptation, ...). The AZ strategy itself has also justified it. The condensation and the mesh adaptation strategy are similar: they concentrate the work load in the crack tip vicinity and amortize it over few load steps. Hence, all optimizations that the mesh adaptation can provide are only used in the fixed group. The cost of this domain condensation being amortized over many load steps/Newton-Rapshon iterations, the mesh adaptation optimization is less worthwhile.

Not using the mesh adaptation also simplifies the implementation. This one, like the frontal strategy in chapter V.2, starts from the TLS framework sequential implementation. It must be parallelized to integrate the AZ strategy. A full distributed implementation would be the best. However, for the same reasons seen in the frontal approach, a mid-term solution has been adopted. It loads the full mesh in each process. Then, the distributed algorithm is only introduced for key tasks used in the AZ strategy. All other steps remain sequential (but done on each process).

A full description of obtained performance on these test cases is given in annex I. What has to be retained is first that the enrichment base algorithm 6 (as said in annex E.2) provides excellent to moderate locations for the AZ envelope creation. The algorithm is faulted when a not fully detached crack lips zone exists or when the crack region includes a full element. In both cases, the AZ group becomes larger than expected and the computation slows down. The crack tip vicinity localization is to be improved. Nevertheless, when the AZ envelop is well identified, this strategy gives, as shown in table I.5, a good time consumption reduction. However, the scalability performances are hard to investigate with this mixed sequential/parallel implementation. What can be said is that, like the frontal strategy, the scaling of nonlinear resolution depends mostly on MUMPS parallel performance (Scalapack also comes into play). And if a 3D crack tip becomes very long, the associated AZ can have a size that is no more appropriate to a direct solver treatment of the incomplete block factorization. Extending domain decomposition approach, with the AZ divided into many domains, can allow regaining scalability. However, more domains imply a larger Schur complement (even larger if the fixed domain is also split into many domains). The reduced matrix direct resolution with Scalapack can then become an issue. Switching to an iterative solver for reduced system resolution would then be a natural solution. Nevertheless, as evoked in V.2 (even if Schur complement does improve the conditioning ), the matrix conditioning will impose to find an ad hoc TLS model preconditioner. All this work has been postponed to try a completely different approach. It is based on a multi-scaling strategy which should give better scalability. This new strategy is presented in the next chapter V.4. However, the idea of this AZ enhanced version (with many domains) as been evaluated indirectly in the benchmark of chapter V.5. until s ≥ smx or

V.4 Two-scale strategy

Us-U s-1 U s-1 < ε V.4.

General concept

The two-scale (TS) or GFEM gl computational framework (first introduced in Duarte and [START_REF] Duarte | Analysis and applications of a generalized finite element method with global-local enrichment functions[END_REF]) uses interdependent mechanical problems for the simulation. The global level problem treats the entire domain with a finite ele-ment discretization which is only able to capture the global phenomena. The related fine2 level problems capture localized specific behaviors with a fine discretization. A fine level problem definition corresponds to the region covered by a coarse entity support that includes the localized specific behaviors to capture. Consequently, fine level problems overlap.

The local behavior can be of different natures such as welds spot in Li and Duarte [2018], crack in [START_REF] Pereira | A two-scale approach for the analysis of propagating three-dimensional fractures[END_REF], thermo-structural effects in [START_REF] Plews | Generalized finite element approaches for analysis of localized thermo-structural effects[END_REF], ... In all cases, the same approach is used for the interaction between the global problem and the set of fine problems:

1. The global level solution is used to impose the boundary conditions of the fine level problems.

2. The fine level solutions are used directly or not, as enrichment functions to enrich the global space using the partition of the unity framework.

The enrichment function construction from the fine level solution is explained more in detail in chapter V.4.4.

Different algorithms are conceivable using those two interactions. For simple fixed phenomena, a single loop between the scales provides the global solution as shown in algorithm 3. For evolving behavior such as crack growth, a first approach would be to use the algorithm 3 to solve the mechanical problems at each dissipating load step that makes the crack advances.

But with evolutionary behavior problems, the scale loop can also be mixed with a general solving process. This is what has been done in [START_REF] Pereira | A two-scale approach for the analysis of propagating three-dimensional fractures[END_REF] where the authors consider that the growth of the cracks is sufficiently small so that the global solution evolves smoothly. It leads to the general algorithm 4 where the TS enrichment functions are updated using the global solution of the previous load step. Intermixing the load steps and the scale loops is the point that can make the difference compared to other resolution strategies. Because, in this case, a part of the resolution process cost is amortized between load steps.

V.4.2 Integrating the TLS model into the two-scale method

Following the general concept presented in chapter V.4.1, the discontinuous problems represented by the TLS model will be treated at the fine level in the two-scale framework. Here, one possible approach would be to focus on the crack tips only for the TLS computation at the fine scale. Moreover, in the crack wake, a Heaviside enrichment function at the global level would [START_REF] Pereira | A two-scale approach for the analysis of propagating three-dimensional fractures[END_REF] with the GFEM gl method and in Ben [START_REF] Dhia | On the use of XFEM within the Arlequin framework for the simulation of crack propagation[END_REF] and [START_REF] Giovanardi | A hybrid XFEM -Phase field (Xfield) method for crack propagation in brittle elastic materials[END_REF] with other methods. However, to maintain the TLS model versatility (late fork in the crack wake, merging, ...) and to respect the crack path during a simulation ( small ρ ), the whole TLS framework has to be discretized. It does not imply that in the crack wake there is a need to pass by the TS enrichment function to obtain a discontinuous displacement at the global level. Using at the global level, a Ramp Heaviside enrichment function, built on top of the fine-scale discretization, is cheaper.

The alternative is then to maintain at the fine level an appropriate discretization for the TLS model in the whole damaged band. And at the global level, one must use enrichment with the TS function only3 if the patch may not be treated by a Ramp Heaviside function. Note also that the weak discontinuity of the TLS model (damaged material) is taken into account at the global level: macro elements are integrated using the fine-scale elements so that Ω -, Ω + and Ω c are distinguished in these global element stiffnesses. This solution is illustrated in figure V.3, in a fictitious 2D simulation context. known somehow (it could be a mathematical definition at the beginning of the simulation , or a level set defined in another mesh obtained from a previous simulation step, ...). It is represented here just for visualization purposes since it is not defined at this scale.

A first task is to identify the nodes4 that do have their support covering the damaged band (it is illustrated in figure V.3b with 4 supports of different color). All nodes and their support covering the damaged band are visible in figure V.3c. Then, as mentioned above, the fine-scale problems have to be created for each of them. By recalling that in the TLS model, φ is discretized on nodes of the mesh, this field should logically be defined at the fine scale within the TS approach (to provide a correct model description). Now does the φ definition have to be replicated on each patch? And, does each patch must have its own fine mesh discretization? It looks more handy to have the level set defined on a unique mesh representing the fine scale. This unique mesh called the "SuperPatch" (similar to the master-local domain in [START_REF] Li | A two-scale generalized finite element method for parallel simulations of spot welds in large structures[END_REF]) would correspond to the union of fine-scale elements of all patches if they were all discretized in the same way. This last point (uniform discretization) provided by the use of a unique fine-scale mesh is essential for the integration: a fine-scale element is unique for its covering supports . Hence, at its Gauss points, all associated enrichment functions are defined (it would not be the case with a different discretization per patch) when integrating the matrix at the global level. Using a SuperPatch also simplifies (I) and (III) of the algorithm 1. Moreover, it allows running only once the cutting algorithm (in particular, the DCA presented in chapter III). The SuperPatch construction, by itself, follows the unstructured adaptation strategy of chapter IV.1.3. It is illustrated in figure V.3e for the simple 2D example where C corresponds to the colored elements of figure V.3c (i.e. the support union).

Having the SuperPatch, the level set φ is transported on this mesh and Γ 0 and Γ c are computed ( shown on the same figure V.3e even if they are drawn far too precisely compared to the non-realistic refined mesh that we see).

Once done, one can analyze the node supports of global mesh, to count, using the fine-scale definition of Ω c , the number of parts of these supports divided by the crack. If there is only one part, the TS enrichment is used (red node in figure V.3f). Otherwise, a ramp Heaviside enrichment function is chosen (green node in figure V.3f). Figure V.3f visualizes these parts for 3 supports. For the ramp Heaviside, no extra computation is required. However, for the TS enrichment the enriched function has to be computed for each patch.

The patches are created by extracting the SuperPatch elements related to their support as illustrated in figure V.3g. Regarding the patch boundary conditions, they come from the global level displacement and strength fields solutions, which are imposed as Robin boundary conditions5 . Here, some preliminary tests lead to the sole use of loosely imposed displacement by a penalization. Other global boundary conditions are inherited. The TLS framework is naturally present in the patches since the elements are the SuperPatch one. In particular, the damage and the ramp Heaviside enrichment (at fine level, not to be confused with global level enrichment) are already created in the SuperPatch. The resulting patch displacement solution is then used as a TS enrichment function at the global level as detailed in chapter V.4.4 and scale loop can be used.

The scale loop integration with the TLS non-linear resolution (II of algorithm 1) is still under analysis. In a very preliminary work [START_REF] Salzman | On performances of the thick level set method in 3d quasi-static fracture simulation of quasibrittle material[END_REF]), the solution of embedding Newton-Raphson loop into a scale loop has been tested with a sequential TLS/TS framework version. Fine-scale problems were interpreted as linear problems where the strains used in (II.5) that give traction/compression states are taken from the global scale. When entering the scale loop, the previous load step solution is used to compute fine-scale problems (boundary conditions and strains). Then, the problem at the global scale is solved (Newton-Raphson iterations) with a fixed TS enrichment function. At the convergence, the fine-scale problems are again computed with the new global solution. And the non-linear resolution is launched again. The scale loop is stopped when the fine-scale updates do not change the non-linear solution. The test shows that the convergence is slow with this scheme. It could be improved by solving at fine-scale the non-linearity. Alternatively, using a completely different approach, the scale loop can be reduced to only update the TS enrichment functions from the previous load step. It corresponds to disregard the error introduced by not looping again after the Newton-Raphson resolution in the current load step (as proposed in algorithm 4). This is left as a prospect.

Finally, as φ progress, two adaptation criteria should be taken into account to update the SuperPatch: the conventional element size criterion (chapter IV.1.2) and the fact that Γ 0 passes in a support not activated yet. In this thesis, the SuperPatch adaptation has not been tested.

V.4.3 Two-scale method in a distributed context

In the TLS approach, the crack tip in 3D may generate a significant amount of patches that do overlap. This redundancy may be masked only if the patches represent tiny problems or if their treatment is done in parallel. In [START_REF] Kim | Parallel simulations of threedimensional cracks using the generalized finite element method[END_REF], the authors use a multi-threaded parallelism to compute each patch. However, from the TLS simulations done with an AZ strategy (chapter V.3), millions of dofs can easily be encountered in the case of long crack tips. Using only multi-threaded in this context is expected to be limited by memory issues. So, as in the frontal approach (chapter V.2), the idea is to distribute the global mesh on processes relying on MPI, to exchange useful information amongst themselves. In this strategy, an element is only present in one process. This choice impacts the TS method treatment:

As fine-scale problems are concerned, how do we treat a patch that has its global elements in more than one process6 ? One possible solution is to consider that this kind of patches are duplicate somehow on each process. Then, each process computes the same patch but this is too costly. An intermediate solution is to consider that the patch is dispatched as its macro element on several processes. However, only one process holds the full definition of this patch and compute it. The obtained solution must then be scattered on the other processes related to this patch. It is complicated, because part of the mesh (or the contribution to the matrix system from this mesh part) must be duplicated on the process performing the resolution of the patch. It does not respect the rule "one element is held only by one process" and implies many developments in eXlibris. It also introduces load balance asymmetries by leaving only one process of the patch to do most of the work.

The other solution is to dispatch the patch as their macro elements and to assume that the patch problem itself is treated in a distributed manner. It introduces some complexity in the patch treatment as explained below, but it offers a double-level of parallelisms. It can be useful in the TLS context. With the mesh adaptation from chapter IV.1.3, the patches are not of the same size: only the damaged band is locally refined so a patch entirely inside the band will be bigger (more dofs) than the one which is partially in the band. As soon as there are differences, one can arrange the coarse mesh partitioning (dispatching of coarse elements on all available processes) so that largest patches are treated in parallel and the smallest one in sequential. The load balancing is then expected to be good. A further argumentation about this double-level of parallelism, related to the size jump between scales, is given in the benchmark test analysis (of chapter V.5). However, the greatest complexity related to this double-level parallelism comes from the fact that the patch treatment can no longer be done in random order. Two patches defined on the process 1 and 2 for the first one and 1 and 3 for the second one, may not be computed at the same time on process 1 (even though it is possible on process 2 and 3). The patch treatment has to be ordered across all processes. This sequencing algorithm, quickly presented here, retake the order proposed in [START_REF] Kim | Parallel simulations of threedimensional cracks using the generalized finite element method[END_REF] ("bigger" patch treated first). Besides that, it imposes sequences where the possible distributed patches are chosen so that the greatest number of them are computed in the same sequence. If in a process, no distributed patch can be computed in the current sequence, a local patch is taken. And if no local patch is available, nothing is done (which is the worst case).

When all distributed patches are treated, the remaining local patches are computed independently, following the local ordering rule "biggest firstly" (processes are not synchronized anymore). For that, in one process, the natural sequential computation can be used but a multi-threaded computation can also be a solution. This last choice would turn this application into a hybrid MPI/multi-threaded simulation code. For now, many eXlibris components (assembly, sparse solver ) are not thread-safe. That is the reason why a hybrid TS was not implemented but this is not impossible from a conceptual point of view. In terms of efficiency, it has to be compared with MPI: does the usage of many threads with few MPI processes (and fewer distributed patches) better/same/worse than the usage of only many MPI processes? In any case, it would give great flexibility to the user. This question is unresolved and left as a future work prospect.

Apart from the patch processing, one must not forget that the assembly must also be well distributed. As integrating global elements leads in some cases to review the associated fine-scale elements, the assembly cost is not the same across the global-scale elements. Hence, to sum up, the load balancing and the associated global mesh partitioning is mainly driven by :

• Number of fine-scale elements per global element, used for the integration.

• Number of enriched vertices per global element.

• Cost of the patch resolution (as said above the "biggest" patches should themselves be distributed) which, in general is related to its number of fine-scale elements.

• Use of the TS or the ramp Heaviside enrichment at the global level.

In this work, the load balancing has not been investigated in detail and this question is left as a perspective. A simple crude partitioning is used in the current testing implementation ( based on the number of patches associated with a global element)

The TS method implementation in the distributed context proposed in this chapter has been carefully evaluated in chapter V.5 (but in a non TLS con-text). It has also been used successfully (giving correct results in parallel) in various tests not presented here.

V.4.4 Enrichment strategy

Let us consider the standard first order finite element shape functions N i ( x), i ∈ I g = {1, 2, ..., N }, associated with a global node x i , with a support given by the union of all finite elements sharing the node x i , in a global mesh discretization with N nodes. In this chapter, only the TS enrichment function is retained at the global level. This function that can be in scalar or vector form is considered here in its vector form to simplify the presentation. Let us call it

F k ( x), k ∈ I g e = {1,
U ( x) = i∈I g N i ( x) U i + k∈I g e A k N k ( x) F k ( x) (V.1)
At the fine-scale level for a given patch k, we consider the classical firstorder finite element shape functions n j ( x), j ∈ I l k = 1, 2, ..., n k , associated with the fine-scale node x j , with a support given by the union of all finite elements sharing the node x j , in the fine mesh discretization with n k nodes. The TLS scalar ramp Heaviside function associated withe the k patch is called RH k r ,r ∈ I l k e = 1, 2, ..., t k . It is associated with the global node x r and has the same support as n r ( x). I l k e represents the t k fine-scale nodes enriched by the TLS model. Kinematic equation used to describe the fine-scale displacement field is then:

u k ( x) = j∈I l k n j ( x) u k j + r∈I l k e n r ( x)RH k r ( x) a k r (V.2)
In scientific publications, many solutions are proposed to create F k ( x) with the displacement u k ( x), solution of the fine-scale problem k. The simplest solution is to use the fine-scale displacement without further manipulations:

F k ( x) = u k ( x), k ∈ I g e (V.
3) It is the solution that has been used in chapter V.5 (where all nodes are enriched). However, when the blending elements 7 are present, a more accurate solution is possible. Proposed in [START_REF] Gupta | Stable GFEM (SGFEM): Improved conditioning and accuracy of GFEM/XFEM for three-dimensional fracture mechanics[END_REF] for 3D fracture simulation with GFEM (similar to XFEM), the Stable GFEM (or SGFEM) method removes from the enrichment function its projection on linear standard finite element space. In the TS context, it corresponds to a fine-scale displacement field projection over the global-scale space removal. The enrichment function becomes:

F k ( x) = u k ( x) - p∈I g s k N p ( x) u k ( x p ), k ∈ I g e (V.4)
where I g s k is the set of the global nodes used to discretize the global node k support. A test, not presented here, has been done with the problem described in annex H. Similarly, as in annex D.3, the plate has been refined in its center to have some blending elements. It confirms that results with a non-TLS framework are improved with (V.4). Nevertheless, in this test as well as in the TS literature, u k ( x p ) always exists at x p nodes. With the TLS model, it is not anymore the case. The node x p may be in Ω c and as such eliminated from the fine-scale displacement field. This issue, that can be addressed by taking u k s ( x p ) = U s-1 ( x p ) (with s being the currently considered scale loop and s -1 the previous one, see algorithm 3), remains an open question.

Otherwise, note that F k ( x) as it uses u k ( x) is not anymore a dimensionless quantity such as the Heaviside function values. Depending on the chosen unities, the tangent matrix can have some terms with very high or low values compared to other standard non-enriched dofs. The conditioning is then impacted. A scaling of F k ( x) can be a solution. For (V.4), the scaled enriched function becomes:

f or k ∈ I g e                  G k ( x) = u k ( x) - p∈I g s k N p ( x) u k ( x p ) F k ( x) = G k ( x) max j∈I l k max i∈{1,2,3} G k ( x j ). e i (V.5)
with ( e 1 , e 2 , e 3 ) an orthonormal basis of 3D Cartesian space.

Using the same test that check (V.4), the matrix conditioning has been improved with this scaling. Discovered later in the literature, a scaling SGFEM has already been proposed in [START_REF] Sillem | The Orthonormalized Generalized Finite Element Method-OGFEM: Efficient and stable reduction of approximation errors through multiple orthonormalized enriched basis functions[END_REF] where the authors call it the scaled SGFEM or sSGFEM. An illustration of this enrichment function is given in annex G. Note that this scaling factor is a good indicator concerning the enrichment pertinence. If it is close to zero it means that u k ( x) can already be represented in the global-scale space. Moreover, in this situation, the global node k must not be enriched to avoid a singular problem (independently of the fact that dividing by zero for scaling is not possible). It has allowed switching off enrichment automatically, when using (V.5).

Introducing the TLS model in the TS method is still a work in progress at time of writing this thesis. The interaction of the ramp Heaviside enrichment at a coarse level (not parallelized yet) with the TS enrichment is still subject to questions. Those two enrichments overlap. It can even be chosen to superimpose them (by forcing the full overlapping). The idea would then be to remove, not the projection (such as in SGFEM), but a piece-wise linear projection, to u k ( x) so that the contribution of the global-scale ramp Heaviside enrichment is removed from F k ( x). This idea has been proposed with a different approach (not a piece-wise linear projection but a double Heaviside function application) in Sanchez-Rivadeneira and [START_REF] Sanchez-Rivadeneira | A stable generalized/eXtended FEM with discontinuous interpolants for fracture mechanics[END_REF]. Note that such a method can alleviate the issue pointed above with the problematic global nodes in Ω c .

V.4.5 Feed back

The first thing to mention is that finding a reference to compare TS enrichment functions and test ramp Heaviside/TS overlapping is not easy. Using the firstorder frontal approach as a reference is difficult: first, since the SuperPatch is built on a subset of global elements, the refined mesh can differ in the frontal and the TS problem; second, with the SGFEM enrichment a second-order approximation is introduced by this function that first-order reference will not be able to represent. Moreover, no simple mathematical reference has been found yet. Nevertheless, many tests (one of them is used as an illustration in annex G) have been conducted using a symmetric version of II.5 (i.e. β = 1.) and using a fixed φ (no evolution). These tests permit to check:

• That the TLS ingredients are well introduced in the TS framework using the SuperPatch.

• That the scale loop does converge to a similar solution to the one of an equivalent problem treated by a frontal strategy.

• That the choice of the TS enrichment functions does influence this convergence. Finding the right function is an essential point to be able to progress on this subject.

• That the parallel features proposed in chapter V.4.3 does work.

This last point is illustrated in the next chapter.

V.5 Uniform test case benchmark

This test is an artificial benchmark dedicated to comparing the strategies performance of chapter V. All detailed information for this test case are given in annex H. It corresponds to a simple plate under a known complex loading that gives a way to compare scalability of:

• The frontal strategy (chapter V.2) called "frontal" in this benchmark.

• Some form of the active zone strategy ( in fact it corresponds to the enhanced version with many domains proposed in chapter V.3.2) called "domain" in this benchmark.

• The TS strategy (chapter V.4) called "two-scale" in this benchmark.

This plate, more or less refined, has to be seen as a 3D crack tip with a far simpler strain field without discontinuous displacements. The simulation is intended to evaluate strategies performances when the crack tips become spatially very wide. The overall part and the crack wake discretization are de-facto neglected in this study. Their real impact on the performance depends on the strategy. Let us hereinbelow use the term "reality" to refer to a complete TLS simulation.

With the TS strategy, neglecting the part and the crack wake is rather true because the global level is in charge of them. And normally it is expected to remain a fast step even if the crack wake global-scale integration is impacted by the fine-scale discretization. So, in this benchmark, the "two-scale" results are rather close to the reality.

Regarding the frontal strategy, adding overall part and refined/coarsened crack wake discretization implies a larger system size to be solved. So, in this benchmark, the "frontal" results are optimistic compared to the reality.

For the active zone strategy, what has been tested is the decomposition of the active zone itself. Each process holds a domain which is condensed. The resulting distributed condensed dense system is solved with a home-made parallel Conjugate Gradient iterative solver, with a Block Jacobi preconditioner. Note that using such iterative resolution is possible in this benchmark due to the right problem conditioning. In reality, this conditioning aspect, as already mentioned, can be an issue for an iterative resolution. Apart from this point, the Schur complement size (i.e. condensed problem size) should be increased because of the interface between the AZ group and the fixed group (corresponding precisely to what has been neglected in this benchmark). Besides that, the time for the fixed domain condensation should be added. So, in this benchmark, the "domain" results are optimistic compared to the reality.

Note that for the "domain" approach, fixing only one domain per process looks rather unfair: reducing domain sizes is only done by increasing the number of processes, which can slow down the iterative resolution at some point. Maybe, a more efficient implementation would have used an independent number of domains and processes (i.e a process may handle one or more domains). It would have permitted, for example, to have many domains treated in sequential, and obtain an efficiency for a constant domain number. Due to the complexity of the implementation as well as lack of time, only the one domain per process solution has been tested. It already gives some fruitful information.

What we can retain from results given in annex H.3 is the following:

• The TS strategy both for strong scaling and weak scaling has better efficiency than the two other strategies.

• Modifying the active zone, using many domains looks promising because of an excellent strong scaling efficiency compared to the original strategy that can be assimilated to the "frontal" results in this benchmark (i.e. MUMPS is used in parallel to factorize matrices, in both cases).

A multi-domains per process strategy may change the observed results positively.

• The TS strategy even with small problem, still has a good strong scaling capability.

• The most critical case shows that the TS strategy has a good strong scaling efficiency up to the point where:

-The parallel global-scale resolution slows down the computation.

The resolution scalability corresponds to the performance of the direct sparse solver. Due to the global level moderate problem size, using the same distribution for both scales is then not efficient.

Gathering a global level problem matrix on a smaller set of processes can provide a way to use more efficiently the direct solver with hopefully a small computational effort spent in communicating.

-All the patches are themselves distributed which looks like playing a role int the performance.

• As expected, the sequential TS strategy is slower than other strategies even if we can notice that when the fine-scale discretization grows, this assertion becomes less valid. Up to the point where lack of memory becomes an issue. And in this case only the TS strategy can obtain a sequential resolution. The same ascertainment may be done with a moderate parallelism and small problems: with few processes, a frontal strategy defeats, in absolute terms, the TS strategy.

• For a given strategy, the energy error which is computed in this benchmark is independent of the number of processes used. This evidence only proves that there is no bug in the implementation.

• As expected, the size jump between scales in the TS strategy drives the computational performance. If the jump is high then the patch problem size will be large and can reach the computational effort of a global frontal approach. That is the reason why in [START_REF] Kim | Parallel simulations of threedimensional cracks using the generalized finite element method[END_REF], the authors use a multi-scaled strategy so that the size jump between levels remains small: A global level has "small" patches, themselves considered as the global level of new TS strategies, themselves having "small" patch ... In this thesis, we only use one scale level because:

-The distributed context with a distributed patches resolution is expected to be equivalent to the use of many scale levels if the biggest patches (i.e. those with a large problem size) are distributed over more than one process. This equivalence need still to be proved. However, we can consider that this proof is almost already given through this benchmark. Using extra scales or using distributed patches resolution is similar, respectively to the use of the TS strategy or the use of the frontal strategy. Regarding small to moderate size problems with few processes it turns in favor of a frontal solution.

-It simplifies this first testing implementation in a distributed context. Maintaining local and distributed patch at all scales adds a lot of extra work to afford the right load balance.

• The amortization of the TS strategy computation over several loading steps has been studied in this benchmark by recovering the time elapsed by iteration of the scale loop. This quantity is not counting the factorization step of fine-scale problems resolution that can be the same between non-linear iterations (chapter V.4.2 ) or load steps (φ constant in the patch). Considering this measure reveals the potential of the TS strategy compared to other strategies even for small size problems.

In conclusion, this benchmark has comforted the idea that using a TS strategy appears as a valuable solution for the TLS 3D simulations targeted by this thesis.

Chapter VI

Simulation gallery: confrontation to complex 3D problems

This chapter illustrates the capabilities of the Thick Level Set model in 3D. However, a full simulation campaign (with convergence rate analysis, l c testing, ... ), already well studied in 2D [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF],...), was not carried here in 3D. Only a few qualitative and quantitative performance of this method compared to some experiments and other numerical models are presented. A first simulation corresponding to a cube with spherical holes, under traction, starts presenting the TLS model features, such as the automatic damage initiation or the merging capability. Then, a twisting chalk simulation gives a first clue to the correct quality of this method when compared qualitatively to the experiment. After, a more quantitative comparison to experimentations is given by a set of notched beam simulations. With a different setting of the TLS parameter, those beams are loaded under three points bending conditions. Next, an L-shape simulation with modes I+III loading gives a way to compare with another numerical method qualitatively. Finally, a spiral bevel pinion gear of a helicopter transmission system, is studied under a critical loading. It allows some good comparison with tests and simulations carried out with a more complicated loading. Note that, most of those simulations have been used to illustrate the computational performances of the AZ strategy exposed in chapter V.3. This chapter slightly modified in this thesis have been published in Salzman et al. [2016] All material characteristics and simulation settings are given in annex J. Annex I gives the computational cost of those tests. Note that most of those tests use the AZ strategy (chapter V.3) to solve the problem. Moreover, none of those simulations use the mesh adaptation strategy of chapter IV as explain in section V.3.2. The DCA (chapter III) is used in all simulations to obtain Γ c .

VI.1 A cube with spherical holes

This test case shows the capabilities of the TLS framework in terms of initiation and merging of crack. A cube with 4 embedded spherical holes (Figure VI.1 ) is in tension with a uniform loading condition applied on 2 faces (perpendicular to y-axis).

The simulation starts with an undamaged model. In the first initial step, a search is performed to find the location of the maximum damage criterion. A spherical iso-zero level set is put at this location introducing a small damaged zone. It is added around the point A (figure VI.1). Then at each step, a search for an additional location that locally violates the damage criterion (Y Y c ), is done outside the already initialized damaged zone. As the load factor is still slightly increasing during this part of the simulation, extra local zones The material section between the spherical hole and the cube face is so thin that its local degradation does not affect the overall stiffness too much. After 5 steps, the load starts decreasing and at step 10, a crack appears (i.e. Γ c appears). From this point, there is almost no more additional location violating the damage criterion found during the simulation. During the propagation, one can observe the merging capabilities offered by the TLS. Figure VI.3 exhibits the behavior when, after turning around the first spherical hole, the fronts are joining together. First iso-zeros of both sides of the front merge (top-right view). As the distance from the hole is less than l c , the propagation continues with a rather steady location of Γ c ( bottom-left view). The propagation continues and the front gets far enough from the hole (more than l c ) which generates a quick coalescence of the Γ c front (bottom-right view). Here, the fronts were colliding in a rather straight- The final crack location is given in figure VI.5 after 386 load steps. The cube is split into two components.

θ - → F = - → F max . sin θ 2 L = 200mm

VI.2 Chalk under torsion

This test case was studied in [START_REF] Bordas | Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment[END_REF] with a meshfree method. A cylindrical chalk bar is twisted with two opposite torques at its ends. Geometry and loading are presented in figure VI.6. The loading is applied along the tangential direction of the chalk surface. To follow [START_REF] Bordas | Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment[END_REF], the loading magnitude depends on the angle θ defined in figure VI. As the object is axisymmetric, a defect must be introduced to start at a defined location. In this work, a simple initial damage is set by a Γ 0 small sphere with a center over the chalk surface at L/2 along the cylinder axis. It allows reducing the time consumption by using a mesh with a refined slice where a crack will start and is expected to develop. Figure VI.7 presents the result of the simulation after 103 load steps, when the chalk is completely split into two parts.

One can observe a good agreement of the crack shape between the experiment and the simulation: the development of a helicoid (bottom-left and middle points of view) is followed by a blunt finish (bottom-right point of view) where the crack lips shape is straighter.

VI.3 Quantitative comparison : Three-point bending test of a notched beam

In this chapter, we try to provide some insight into the quantitative quality of the TLS method by making a comparison with an experimental test. We choose one of [START_REF] Hoover | Comprehensive concrete fracture tests: Description and results[END_REF] campaign which provides data for a size effect study on a concrete beam under three points bending conditions. The general beam geometry and loading are depicted in figure VI.8. Two (out of 18) size combinations are chosen from [START_REF] Hoover | Comprehensive concrete fracture tests: Description and results[END_REF] campaign : the Bc and Cb cases (dimensions given in table VI.1).

The force displacement curves given in figure VI.9 show the experimental response. The CMOD (crack mouth opening displacement) is measured at the bottom of the beam between two points symmetrically located ac- cording to the notch (and separated by 137 mm and 59mm in the test case, respectively Bc and Cb.).

The Cohesive Zone Model (CZM) simulations corresponding to the test case are given in [START_REF] Hoover | Cohesive crack, size effect, crack band and work-of-fracture models compared to comprehensive concrete fracture tests[END_REF] and are reproduced in figure VI.9.

The idea is to compare the TLS simulations with these two results sets (experiment and CZM). In order to model properly the concrete behavior, the brutal softening regime of figure II.3 can not be used. In order to adopt a different softening regime, we consider the work of Parrilla G ómez et al. [2015] where Y c now depends on d :

Y c (d) = Y 0 c h(d) (VI.1)
where Y 0 c is a constant.

Using [START_REF] Parrilla G Ómez | Comparison between Thick Level Set (TLS) and cohesive zone models[END_REF], the CZM Hoover and Bažant [2014] parameters are transformed into TLS parameters ( h(d) and Y 0 c ).

The results given in figure VI.9 show that the TLS simulations are rather close to the test and CZM simulations. The test peak load is obtained with an error of 2.5% and 4.5%, for the Bc (top) case and for the Cb (bottom) case, respectively (Note that the CZM simulations do give errors of 1% (Bc) and 4.1%(Cb) with respect to the test peak load). The last part of the post-peak is in discrepancy with the test. In [START_REF] Parrilla G Ómez | Comparison between Thick Level Set (TLS) and cohesive zone models[END_REF], the load-CMOD curves on another test case show a more accurate post-peak region. In [START_REF] Parrilla G Ómez | On the capability of the Thick Level Set (TLS) damage model to fit experimental data of size and shape effects[END_REF], a full comparison with [START_REF] Hoover | Comprehensive concrete fracture tests: Description and results[END_REF] campaign in 2D, have investigated this aspect and the authors observed better coherence.

In figure VI.10, Γ 0 and Γ c are presented for points A, B and C. Those points, shown in figure VI.9, correspond to the curve peak (A), the first appearance of Γ c (B) and the end of simulation (C). One can observe that the point B, for both cases, happens quite late in the process. The damage first develops in a very long process zone. In VI.10e, we see that Γ c appears in a non-uniform manner. It is related to the lengthy process zone where, on the skeleton of φ, the damage is almost 1. The transition to a damage equal to 1 with an automatic introduction of Γ c discontinuity is then hard to softly achieve because a large part of the skeleton can vary abruptly. In those simulations, we chose to slow down the front advance (which gives this non-uniform intermediate state). However, the new TLS model version introduced in L é et al. [2018] is certainly a better solution to this problem and could also have some effects in the last part of the post-peak curves too.

VI.4 L shape -mode I+III

In [START_REF] Lorentz | Gradient damage models: Toward full-scale computations[END_REF], the authors modify the standard L Shapemode I test case by introducing an additional lateral imposed displacement. TLS results are iso surface at the end of computation, Γ 0 in blue, Γ c in red.

Consequently, the mode III is activated and the crack path is slightly modified.

The geometry and loading are presented in figure VI.11. The mesh is refined around the corner A up to the edge B.

Two types of computations have been conducted: one with an automatic damage initiation and the other one with a forced damage initiation (along the corner, setting an initial cylindrical φ of radius 1.05 × l c ). This second case is needed to improve the comparison with [START_REF] Lorentz | Gradient damage models: Toward full-scale computations[END_REF] results where the authors obtain, at the corner, a rather straight damage distribution. The force-displacement curves are given in figure VI.12 and the final situations are presented in figure VI.13. The comparison in figure VI. 13 shows that there is little to no difference between the final crack shapes between the two damage initiation strategies, except in the corner vicinity. Compared to [START_REF] Lorentz | Gradient damage models: Toward full-scale computations[END_REF], the crack is twisting the same way but the simulation was conducted further in this work as compression is considered in the free energy density (II.5) in Ω + . This test illustrates the interest of an automatic damage initiation which here allows obtaining a more complex crack path in the corner vicinity.

The displacement field at load step 120 in figure VI.14 illustrates a correct enrichment and representation of Γ c . The opening of the crack faces is evident.

In conclusion, this simulation shows that the TLS model gives equivalent results when compared to another numerical tool which uses a damage gradient method.

VI.5 A spiral bevel gear

This test case mimics an industrial study on a spiral bevel pinion gear of a helicopter transmission system. It first appeared in a Nasa report [START_REF] Spievak | Simulating Fatique Crack Growth in Spiral Bevel Gears[END_REF]) and was published in [START_REF] Spievak | Simulating fatigue crack growth in spiral bevel gears[END_REF]. Later [START_REF] Ural | Three-dimensional, parallel, finite element simulation of fatigue crack growth in a spiral bevel pinion gear[END_REF] again use this problem with a new computational technique. Those studies were conducted to help the pinion gear designers by giving them a crack path coming from the simulation. The purpose of this test case is to show that the TLS method, with a complex geometry, gives right information about the crack shape without specific mesh refinement and any initial crack placement.

From [START_REF] Ural | Three-dimensional, parallel, finite element simulation of fatigue crack growth in a spiral bevel pinion gear[END_REF], [START_REF] Spievak | Simulating fatigue crack growth in spiral bevel gears[END_REF], an approximate geometry has been rebuilt from scratch with only three teeth. In this work no fatigue study with any varying complex loading is done since the testing implementation only deals with a quasi-static loading. The load location corresponds to the highest point of a single tooth contact (HPSTC see [START_REF] Spievak | Simulating Fatique Crack Growth in Spiral Bevel Gears[END_REF]). An ellipse E on an arbitrary CAD plane P1 is projected on the tooth to follow the Hertz contact shape. In figure VI.15c, it appears as a yellow zone where the mesh over the surface is following its boundary. The loading magnitude F in this zone is computed with the following expression:

F = H c .F 0 where: • F 0 is the maximum load in this zone • H c coefficient is H c = 1 - x a 2 - y b 2
where x and y are the projected coordinates of a point in the loading zone, on plane P with Cartesian coordinates corresponding to the center and axes of ellipse E. a and b correspond to the major and minor radius respectively of the ellipse E. This simulation produces the following results:

• The (automatic) damage initiation is found in the concave part of the fillet at almost mid-distance (8% shift) from the toe and heel sides (see VI.15d for sides location). It is where [START_REF] Ural | Three-dimensional, parallel, finite element simulation of fatigue crack growth in a spiral bevel pinion gear[END_REF] put the crack. Only a small damaged zone is put at this location which makes no assumption about the future crack path.

• The first loading step extends this zone, and at the 9 th load step a crack (Γ c ) is automatically put. It emanates from the damage initiation location (figure VI.16).

• The crack shape is consistent with the experimental results given in [START_REF] Ural | Three-dimensional, parallel, finite element simulation of fatigue crack growth in a spiral bevel pinion gear[END_REF] and reported in figure VI.18.

The general situations after 246 and 380 load steps are depicted respectively with iso-surfaces in figure VI.17 and with the displacement field in figure VI.19.

Note that in the performed simulation, the chosen length l c is rather large with respect to the process zone size of the actual material. However, the chosen material strength and length l c do combine (Y c × l c ) to produce the right order of magnitude for the toughness of the material.

Chapter VII

Conclusions and perspectives

90

VII.1 Thesis assessment

The simulation of cracks initiation and propagation in 3D is nowadays still hard to tackle. On the one hand, the mechanical models have their pros and cons regarding the accuracy and the features they provide in a such context. So making the right choice is not simple. On the other hand, computational consumption remains, most of the time, an issue for 3D simulation. In this thesis, we have addressed both concerns.

This thesis first demonstrates that the TLS model is well-designed for the quasi-static simulations of quasi-brittle material in 3D. The presented simulations confirm that the versatility and quality of the model are well-kept in 3D. The new "double cut" algorithm added in this work participates in this quality, giving the right representation of the crack lips without any extra discretization cost.

Regarding the computational performance, this thesis again demonstrates that TLS model is well-designed for 3D simulations since it imposes a very localized constraint for the problem spatial discretization. It has been exploited in the mesh adaptation strategies, proposed in this work, to reduce the problem size and consequently its computational cost. Despite those efforts, the linear system resolution that arises from the non-linear problem treated in the TLS algorithm remains a computational bottleneck. Different strategies have been suggested in this work to overcome this problem. All of them lead to consider parallelism as a natural remedy. In particular, the message passing protocol has been adopted for its performance and hardware flexibility. It induces in all those strategies an element-based partitioning with a distribution of the problem on all available processes. Among the three proposed solutions, this thesis shows that the two-scale strategy offers the best potential in terms of scalability and algorithmic integration within the TLS solver.

Note also that regarding computational throughput, the "double cut" algorithm implementation has exhibited excellent performance thanks to its cutting database. Moreover, a distributed mesh, which implies distributed integration and assembly, has also greatly contributed to reduce the elapsed computation time.

VII.2 Perspectives

An obvious perceptive is to finish the TLS integration in the distributed twoscale method. For that, one will need to determine the right enrichment to use at the TS global level to obtains the best results convergence. Then, one will have to optimize the scale loop integration in the staggered algorithm 1. All this will be possible only if all components of this algorithm are well parallelized (in particular if the new modal concept proposed in Moreau et al.

[2017] is adopted). Regarding the parallelism, hybrid MPI/multi-threaded paradigm should be benchmarked. Finally, the scale size jump, load balancing and SuperPatch update (in particular, the tracking of small isolated damaged zone will deserve some attention) will have to be addressed. On top of all previous steps, one may wish to adapt the SuperPatch mesh using the proposed coarsening strategy of chapter IV.2. This will meet the mesh adaptation perspectives. On this topic, the crack tip vicinity identification proposed in this work has to be consolidated. The tool that provides Γ c might give this additional topological information in the future. The offsetting concept should also be improved. As far as the coarsening is concerned, this type of adaptation has been possible, thanks to the "double cut" algorithm. However, some transport issue has appeared in this adaptation, and needs to be addressed in the future. Finally the unstructured mesh adaptation tool created in this work will have to be enhanced to provide better throughput and more functionalities.

In the end, as far as valorization is concerned, continuing the transfer in eXlibris of the code developed in this thesis (distributed "double cut" algorithm version, distributed TS implementation, ...) is expected. If the TLS model integration in the distributed TS method comes to an end, publishing an article on the subject would definitively value the preliminary work proposed in this thesis.

A.1 Element Cutting

Following the scheme given in chapter III.3.2, using [START_REF] Hert | 3D convex hulls[END_REF] for the convex hull computation and [START_REF] Hachenberger | 3D Boolean operations on Nef polyhedra[END_REF] for the subtractions, one can obtain, in 3D, with tetrahedron element, 17 reference patterns of element cut as shown in figure A.2. A cut pattern corresponds to a set of edge cut points and a tetrahedron node sign: it leads to a unique topological element cut. The 17 patterns do not take into account the specific metric (single cut are put at the edge middle and double cuts at one-third and two-thirds of the edge).

We can group the 17 patterns into four different categories:

• Patterns 1 to 3 correspond to a simple cut that one can obtain with a standard scalar level set.

• Patterns 3, 5, 6 and 9 have potential Γ c warped surface (four points surface) where the choice of diagonal is arbitrary.

• Patterns 11 to 14 give a non-convex negative domain polytope (Ω + ) which has to be subdivided in convex polytopes to generate subelement tetrahedrons correctly. It corresponds to the extra blue edges on some element faces. It describes the sub-cut used to split the polytopes into convex ones. See figure A.1a for illustration of pattern 13 where the negative domain polytope is split in 3 convex ones.

• Patterns 15 to 17 are termination patterns. Γ c is only present on the tetrahedron boundary and no positive domain is present. For example pattern 17 may terminate pattern 6 bottom face of the positive zone. Those choices are arbitrary.

In table A.1, the number of possible permutations from each pattern is given. Without close node treatment, a total of 111 permutations has to be handled. With a close node treatment, some polytopes change and are not giving anymore null sub-element for the integration. See figure A.1b for an illustration where one of the negative domain parts reduces to a single tetrahedron producing 1 instead of 3 tetrahedrons for the integration. This is interesting when some specific extra computation is made with those subelements (such as some fluid flow in the crack) and when the null volume is an issue. For this work we could have used those null sub-element. However, creating a tool to handle double cut algorithm, we chose to be a little more general from the start, and eliminate null volume sub-elements. It leads to handle 8399 permutations with the close node treatment. From an implementation point of view, a natural choice is to use a database to hold all those patterns. The cutting procedure is then reduced to cut the edges and query the database with the edge cut pattern as a search key. Database gives then the integration cell and the topological relation for the enrichment identification (number of independent parts of supports).

A.2 Geometrical comparison between single and double cut algorithms

To illustrate how single and double cut algorithms perform, we consider a hammer head shark surface ( from Lutz Kettner home page at Max Planck institute: https://people.mpi-inf.mpg.de/ kettner/proj/obj3d/) which is discretized with triangles and plunged into unstructured meshes of increasing number of elements (figure A.3). For each mesh, a level set and a signed vector distance function are computed from the shark surface. An associated cutting algorithm is then used. Fins appear sooner with the double cut algorithm, as expected, since this algorithm detects layers even if they are embedded inside an element. Less expected, is that, even in a rather simple zone such as the shark's body, the use of signed vector distance function gives better accuracy. Two reasons can be pointed out: a more accurate cut locations with the SVDF and a possibility of warped iso-contour inside the elements (patterns 3, 5, 6 and 9).

C.2 Treatment

To get a continuous field description, the hanging nodes have to be eliminated.

Kinematic linear relation is a solution that can be applied in all contexts (octree, unstructured, ...). For example, when computing the mechanical displacement field of figure C.1a in 2D, the use of additional equation, given in figure C.1b, eliminates the displacements of B and D. This solution has the following disadvantages: It is costly when elementary matrices are assembled, the final system to solve is densified which increases its resolution cost and it is less obvious to use with higher-order field.

When using a simplex mesh, two solutions can be used, at least. Both are using the simple principle of meshing elements that do have hanging edge or face to connect the hanging nodes with extra elements. The former illustrated in figure C.1c adds a node at the center of the element for that. Then, a simple loop on edges ( or faces in 3D) creates one triangle (or 2 tetrahedrons in 3D) or two triangles (or 8 tetrahedrons in 3D) if a hanging node is present. The pros are that the meshing algorithm is simple and works perfectly well in 3D. The cons are the addition of many elements (compared to the second solution). Moreover, in an unstructured context, it gives elements with less interesting aspect ratio. The latter illustrated in figure C.1d uses an ad hoc meshing adapted to the hanging situation. The pros are that it gives fewer elements with a better aspect ratio than the first solution. The cons are that treating all possible cases in 3D is difficult. In [START_REF] Schroeder | Compatible triangulations of spatial decompositions[END_REF], they were forced to use a specific constrained Delaunay procedure to achieve this ad hoc meshing.

D

Distributed unstructured mesh adaptation 108

D.1 General concept

The inputs given by the user to this adaptation tool are a distributed tetrahedron mesh, an optional transport function set and a function that returns true or false if the element given as an argument must be split or not. Split an element in 3D is the action of:

• Creating nodes on the middle of its edges (if they are not already created)

• Replace the tetrahedron by eight smaller tetrahedrons joining original vertex and mid-edge nodes. The four tetrahedrons at the middle (not connected to original vertex) can be created in three different ways depending on the new internal edge choice. An aspect ratio criterion (see B) is used to select the best solution among the three.

• A set of actions to handle the hanging nodes treatment and their associated edges or faces.

The same kind of actions are possible in 2D with triangles. The user provides the transport function set if he needs to transport information to the newly created element. For example, an initial mesh coloring can be propagated to the new refined embedded mesh by such transport function.

In figure D.1, the global idea of the algorithm in a 2D sequential context, is described with a small example (figure D.1a presents the initial mesh and the fictitious refinement location). The adaptation is made by successive passes that split the elements fulfilling the criterion. The first pass selects all initial mesh elements fulfilling the criterion requirement (in orange in figure D.1b). Then those elements are split and pass 1 gives the mesh in figure D.1c. The algorithm enters a second pass in figure D.1d with the mesh created from the pass 1. It selects the elements based on the criterion (in orange) and the elements that need to be split to fulfill the 2:1 constraint (in yellow). These elements are split and give the mesh in figure D.1e. Same thing in pass 3, selection in figure D.1f and splitting in figure D.1g. Then, in this example, the pass 4 does not detect any new element fulfilling the criterion which stops the pass loop. The adaptation ends, selecting all elements having, at least, one hanging node (figure D.1h) and uses a particular splitting scheme to obtain the final unstructured refined mesh (figure D.1i).

D.2 Parallel algorithm

The algorithm 5 as mentioned in D.1 works by successive passes. Compared to the sequential vision given in figure D.1, several communications have to be added to propagate the 2:1 constraint and the hanging situation. For each pass, or level, a set of element ( S ) is created based on the criterion (user function) that can have been itself updated by the new mesh definition from the previous pass. This set is filled with the elements that need to be split

Algorithm 5 Adaptation algorithm.

Starting from a distributed mesh with a set of elements M for each process N bpass = 0 repeat update criterion against new distributed mesh Let S ⊂ M be the initial ordered set of element verifying criterion repeat repeat Let O ⊂ M \ S be the set of elements that must be split to respect one level rule if elements of S are split. The first split elements are then not concerned by criterion but their treatment allows splitting those that are concerned, without breaking the 2:1 constrain. This first loop ends when no more elements are in O. The second loop which embedded the first one, works at inter-process level. On process boundary, the elements that are put in S, must inform their adjacent remote element that they will be split. After the communication, a set Q of new elements to be split in order to respect 2:1 constraint on process boundary, is created. It is again drained and added first in S to respect the dependency. This second loop will end if Q is empty for all processes.

S = S ∪ O, O being placed first in S until O = ∅ Let Q ⊂ M \ S
After setting S, the effective splitting is done following S order to re- Then, all hanging nodes are removed using a the central node strategy. The idea is to split the connecting elements (elements that have, at least, one hanging edge or face) by adding in their center of gravity, a node used to construct the internal tetrahedrons appropriately connected to each of their faces. Unfortunately, this is not preserving the mesh quality and a case by case algorithm (where all cases are optimized so that tetrahedrons can fill them without the addition of an extra central node ) may give better elements. In this thesis, we did not push further investigations on this topic (see also annex C).

The sentence "update criterion against new distributed mesh" in pass loop of algorithm 5 means that the user must provide an extra function to cope with the elements newly created in the mesh. Depending on what is used by the criterion function of the user a criterion transport can be mandatory for newly created elements from a previous pass. Otherwise, the criterion function can not give the right answer for those elements. It can be a source of problem if the criterion is not easily transportable.

D.3 Performance

To test the performance of algorithm 5, a parallelepiped is meshed with tetrahedrons. It has the same geometry as the case presented in figure H.1. A uniform refinement would have had no interest. No unbalancing would oc- In this test, two different initial distributions are used. Figure D.2 presents those two distribution in the context of 8 processes. In figure D.2b the initial mesh is distributed using the conventional ParMetis partitioning. In figure D.2a, the initial mesh is distributed by a horizontal roughly equal strip (basic algorithm far from perfect). This last case is expected to give a better load balancing. The refined zone, perpendicular to those strips, is expected to be dispatched in all processes.

The obtained results are presented in figure D.3. Figure D.4 presents the strong scaling results for both initial distributions. The elapsed time curves presented in figure D.4c show that adding process does reduce the time consumption up to a point where the consumption increase. The striped distribution gives a more regular decrease than the Parmetis distribution. In figure D.4a and D.4b, the speedup and the efficiency confirm that the Parmetis distribution gives bad results which is quite natural. The load is not balanced in this case. We can check this aspect in figure D.4d and D.4e that presents the initial (before refinement) and final (after refinement) number of elements per process. The curves in those figures correspond to the averaged value across all processes. It turns to be the same for both distributions, because it corresponds to the global numbers of elements divided by the number of processes. More interesting is this quantity range (minimum and maximum). In figure D.4d, the initial min/max zone is close to the curve with the ParMetis distribution. ParMetis is correctly doing its job by giving a quite uniform dispatching of elements. For the strip distribution, up to 32 processes, the dispatching is correct. However, for 64 process, the stripes are incorrectly created (some process have no element). In figure D.4e, the ParMetis distribution exhibits its weakness with a large distribution unbalancing of refined element.This imbalance begins immediately upon the use of 2 processes. With 8 processes the unbalancing grows explaining a complete loss of efficiency at this point. For the strip distribution, the unbalancing appears later with 32 processes. However, in fact with 16 processes the stripes are already not perfect. It explains why the efficiency is maintained up to 40% until 16 processes are used. A better strip dispatching based on a ParMetis weighted partitioning may give a better efficiency.

In conclusion, as shown by those curves, the right strong scaling of algorithm 5 is difficult to obtain. Nevertheless, communication is not an issue. With 512 processes, the elapsed time remains in order of magnitude of the best obtained values. The intrinsic algorithmic aspects are certainly to be revisited. However, an appropriate initial mesh distribution choice is essential from a performance point of view for this kind of tool.

E.1 General concept

The damaged band zoning is used in chapter IV.2 and chapter V.3. In both cases the crack tip vicinity has to be located. ). This property is related to DCA which almost always stops a crack on an edge (chapter III). When it is not the case, the 3D mesh generally provides a close location where it is the case. In a simple crack wake, the elements usually have, at least, one node enriched. So, this classification brings the appropriate topological information to identify the crack tip. The final step, presented in E.1d, creates the points at the gravity center of the identified element. All those steps are summarized in algorithm 6.

E.2 Limitations

In figure E.2 a more complicated case is presented. It is a forking situation where Ω c has grown up to the point of fully embedding elements. In this case, as it can be seen in E.2d, the elements not at the crack tip are marked. They all do have their nodes not enriched but do not correspond to the tip location. It can be a limitation or not, depending on which application is using this algorithm. For the chapter IV.2 coarsening strategy, it is positive since the refinement must be maintained in such a forking zone. For the chapter V.3 active zone strategy, there is no need to integrate such a region in the AZ group. However, that is not a stumbling block. It will just slow down the resolution.

Figure E.3 illustrates a case with a problem that as been observed during the simulations. It corresponds to the presence of a zone of not completely detached crack lips. In this case, the algorithm locates few elements in this zone (E.3d). Generally, those locks break as the crack progress but it can take some load steps. During this period, things are not optimal. For chapter IV.2 strategy, it can alleviate the coarsening itself. It is not perfect and those connections may (numerically) disappear with larger elements. For chapter V.3 active zone strategy, it is far from perfect as it slows down the resolution during all the periods where they are present. In current AZ implementation, MUMPS is used with its distributed entries capability and is also tuned to use all processes for Schur matrix block-cyclic distribution. As MUMPS does its own matrix distribution for the factorization, all that has to be considered for the load balancing is the distributed matrix creation. It basically comes down to balancing the assemblies. So, with algorithm 2, four assemblies appear.The AZ and the fixed group are separated in linear (in Ω -) and nonlinear (in Ω + ) elements group. Those four groups of elements are then partitioned to obtain a appropriate load balancing for the assembly computations.

A group partitioning is obtained with ParMetis. Figure F.1 depicts the process with a 2D example. On the right, the fixed zone is only treated when a new AZ is computed. The linear and nonlinear assemblies are done separately, so they need an independent partitioning. On the left, the AZ, here made by 2 circles around the tips of the oblique crack , is also split into linear and non-linear independent assemblies. So, an independent partitioning is also conducted for those two groups. During the load steps, those 2 partitioning sets are kept unchanged considering that the re-partitioning costs more than the possible degradation of the load balancing. When φ evolves some elements can pass from the linear group to the non-linear group and then, the partitioning may become less optimal. However, this is expected to be of sufficient quality during all load steps until the next AZ creation.

A more ad hoc partitioning would be to start from a partitioned Γ 0 and by some coloring algorithm, expand this initial partitioning to the AZ and the fixed zone. It would give a priori a more equilibrated partitioning during a load step. It also opens the AZ zone to a natural domain decomposition.

H.1 Description

The test case is a 3D plate with the geometry depicted in figure H.1. We can consider that it represents the tip crack zone in terms of computational effort (number of dofs). The coarser mesh called level zero or L0 (figure H.2a) is unstructured, with a maximal aspect ratio of 3.47 (figure H.3). Different refined meshes are used. At each refinement step, all edges are split into two edges and the tetrahedrons are split accordingly. Each step is called level and named Lx, where x denotes the level number ranged from 0 (no refinement) to 5 (initial element size is divided by 2 5 , i.e, 32).

The mechanical properties are those of a simple uniform elastic isotropic material. The approximation order is 1 for all fields.

H.2 Analytical formula

We want to test a known cubic displacement field with this plate in order to obtain an accurate error computation (i.e. comparison with an analytic formula). The following -→ u (x, y, z) displacements are used as analytic reference:

- → u (x, y, z) = (ν + 1)(1 -2ν)F E    (x 2 ( L 2 -x 3 ) + 2ν(y 2 -z 2 )) -4νxy 4νxz    (H.1)
where :

• E is the Young's modulus

• ν is the Poisson's ratio Considering small displacements, the strain is given by :

ij = 1 2 ( u i,j + u j,i ) (H.2)
and with the displacement field coming from (H.1), we have:

[ (x, y, z)] = (ν + 1)(1 -2ν)F x E    (L -x) 0 0 0 -4ν 0 0 0 4ν    (H.3)
Hooke's law gives the stress tensor:

σ(x, y, z) = C : (x, y, z)

(H.4)
where C is the elasticity tensor corresponding to an isotropic material defined by E and ν.

It results in the following expressions using C and (H.3):

[σ(x, y, z)] = xF    (ν -1)(x -L) 0 0 0 ν(L -x -4(1 -2ν)) 0 0 0 ν(L -x + 4(1 -2ν))   
(H.5) The equilibrium equations are:

-→ div σ(x, y, z)

+ - → f v = 0 on Ω σ(x, y, z). - → n = - → f s on ∂Ω (H.6)
From (H.6) and (H.5) we obtain the load to impose on the test case that gives the (H.1) displacement field :

                                                 - → f v =    (1 -ν)F (2x -L) 0 0    - → f s (x=0) = - → f s (x=L) = - → 0 - → f s (y=H) =    0 xF ν(L -x -4(1 -2ν)) 0    - → f s (y=0) = - - → f s (y=H) - → f s (z=T ) =    0 0 xF ν(L -x + 4(1 -2ν))    - → f s (z=0) = - - → f s (z=T ) (H.7)
The analytic equation (H.3) is used as a reference for the numerical computations of the energy error associated with the different meshes and/or methods:

error( - → u computed ) = Ω ( ( - → u computed ) -(x, y, z)) : C : ( ( - → u computed ) -(x, y, z)) dΩ Ω (x, y, z) : C : (x, y, z) dΩ (H.8)

H.3 Results

In this benchmark, the plate (described in annex H.1 using boundary conditiond from equation (H.7)) has to be considered as a fictitious 3D crack tip with a far simpler strain field. In this context, all the strategies of chapter V have been compared, with for each of them:

• Only the creation, assembly and resolution of the mechanical problem system have been monitored.

• The crack wake, and overall part definition have been neglected (it corresponds to different approximations for each strategy described in chapter V.5).

• The hybrid parallelism is almost never used in the simulations for a fair comparison. For example, using Mumps solver in some approaches with OMP capability and multi-threaded Blas, is unfair compared to none multi-threaded TS strategy. Moreover, multi-threading would have also complexify the creation of the curves (how many threads are needed per MPI process for each strategy, how to consider MKL Blas that adapts the number of threads when computing speedup,...).

• For a given number of processes, the coarse mesh is partitioned on them and splits uniformly up to the target refinement level. This refined mesh is used in the TS strategy (V.4) for the fine scale level. The same partitioned refined mesh is used for the comparison of the frontal (V.2) and active zone (V.3) strategies.

• The simulations have been conducted on the Liger computer (annex J.2) without interference from any other applications at node level (it is not completely the case for the network usage).

• For TS strategy, the benchmark follows the simple algorithm 3. Nevertheless, instead of a test based on the relative variation of U s to stop the loop, the energy error given by (H.8) is compared to the one obtained with the frontal strategy or the modified AZ strategy. The loop stops as soon as TS energy error is less than the one obtain with other strategies.

This test provides a rough estimate of what can be expected from all strategies if they were finalized. It is made without extra complexity (the mesh, the mesh partitioning, ... are the same between the strategies). Recalling the naming convention of chapter V.5, hereafter, "frontal" will designate the frontal strategy, "domain" will designate the modified AZ strategy and "twoscale" the TS strategy. • In (c), the curves of elapsed time for mechanical problem resolution (in second). It corresponds to an average elapse time across processes. For the frontal results, it counts field declaration, system assembly and resolution. For domain results, it counts local domain field declaration and associated system assembly and condensation, preconditioner factorization, iterative resolution of the dense system,and decondensation phase. For TS results, it counts patch field creation, patch boundary condition integration, patch resolution, global enriched field declaration, global enriched system assembly and resolution, and for most of the same items their accumulated contribution while looping on scales.

H.3.1 Strong scaling

• In (a), the speedup curves based on the elapsed time (c). For each approach, the sequential reference time is the one corresponding to its strategy, except for the domain results based on the frontal computation (with one process, there is only one domain without any condensation; it is therefore not possible to obtain a result with our implemen-tation). The ideal speedup is plotted in black (if we use x processes we have to be x times faster than with a single process).

• In (b) the curves of strong scaling efficiency. They correspond to the speedup curves (a) divided by the number of used process and transformed in %. 100% corresponds to an ideal speedup (if we use x processes the gain obtained, compared to sequential time, is the maximum expected).

• In (d), the dofs curves. They show the solved system dimension for all strategies and their different involved scale. For the frontal, coarse and TS global enriched results, the dimensions remain constant when the number of processes changes. With the TS strategy "per patch" gives the average size of the fine-scale problems. Due to a uniform refinement, all patches got quite the same size on average and the parallelism does not change the situation. It roughly leads to a constant size. With the domain results, two curves are given: the first one ("domain(Schur)") indicates the size of the condensed system. The second one ("per domain") gives the average size of the domain problems.

• In (f) the information curves about the TS patch number. "total" gives the global number of patches for the problem. "per proc" gives the averaged number of patches treated by a process. "dist per proc" gives the averaged number of distributed patches treated by a process.

• In (e), the curves give the ratio between the TS or domain and the frontal elapsed times for a given number of processes. If this ratio is 1, the considered strategy has the same performance as the frontal approach. If lower than 1, it is slower. If greater than 1, it is faster. With the TS strategy, the "iter" curves represent the ratio if we only consider one scale iteration without counting the patch creation. This situation is expected when TS is used with algorithm 4.

• All the TS results got a "from x" indication that indicates from which global level, the method starts. "x" correspond to Lx. "from 1", in L3 figure H.6, means that the global level corresponds to the mesh L1 (coarser L0 refined once) and that the fine level is obtained by refining this mesh twice which provides an equivalent L3 mesh (3=1+2).

• It is hard to ensure that the presented curves are representing a continuous processing task increase. Take the example when going from 16 to 32 cores. On the one hand (16 tasks), the inter-process messages are exchanged inside a node (shared memory mechanism) and do not pass through the network. On the other hand (32 tasks), some messages are remaining in the node and some are passing through the network layers since a Liger node is limited to 24 cores. It forces, at least, to use 2 nodes for the computation. A fair comparison, but non-realistic from a practical point of view, is to use a distribution with one process per node. In this context, increase the number of tasks would systematically raise the network load. No additional simulation has been conducted to clarify this point.

Of the L5 results, only those that have been successfully executed are shown in figure H.8. With this kind of problem size (9965229 dofs), memory becomes an issue. Liger computer permits only 128GB per node. Some observations show that the frontal strategy needs memory between 1Tb to more than 3TB (mainly due to the system factorization). For the domain approach, as only one domain is considered per process, the memory consumption remains high with few processes. For the TS method, only around 300GB is needed. It explains why it was impossible to obtain a sequential solution for any approach1 . With the TS method, the smallest possible number of processes to be used is 300 127 ∼ 2.37 which corresponds to 4 processes in power of 2 progression. To obtain a speedup and an efficiency, a perfect scaling between 1 and 4 processes has been considered. It allows extrapolating the elapsed time of the sequential computation. And it explains why we have a perfect strong scaling efficiency with 4 processes. For the other strategies, the first possible computation starts with 32 processes (1 core per node). Extrapolation for the elapsed time of the sequential computation is too hazardous in this case. The speedup and efficiency are then not provided. For the TS method, extra speedup and efficiency curves are provided with a name ending by "ideal global". They correspond to data where the elapsed time of the global level problem resolution has been kept constant to its minimal value (with 128 processes) when minimum has been reached (i.e. for 256, 512 and 984 processes).

As mentioned above, increasing the number of MPI process with the same network loading is not easy to achieve, especially for these L5 results. The first curve points, due to memory consumption, imply a rather natural increase of nodes per computation. However, when reaching 900 cores, the amount of work per core starts to be dominated by communication. Then the results, when passing from 984 cores on 41 nodes (41x24) to 1024 cores on 43 nodes, are too much perturbed by the network load increase to be relevant to compare the methods. This assertion comes from a test, not presented here, where the same 984 cores are used for the computation on 41, 42 and 43 nodes. A drastic drop in performance is observed. It leads to using a maximum of 984 cores on 41 nodes. Finally, as the L5 results are concerned, note that a multi-threaded Blas (with MKL on-demand free strategy) has been tested with 64 MPI processes (named with "th"). It explains why in H.8e, the abscissa is given in number of MPI processes and not in number of cores. With the frontal and domain results, the Mumps multi-threaded capability has also been used with 4 threads per MPI process.

H.3.2 Weak scaling

Figure H.9 presents the "frontal", "domain" and "two-scale" weak scaling results. Three initial refinement Levels are considered (L1,L2,L3). For each initial level, two sets of curves are given: one corresponding to the computational elapsed time in seconds, and one corresponding to the weak scaling. The abscissa is the problem size in number of dofs. The points on those curves correspond to different problem sizes maintaining a constant number of dofs per process. The starting point corresponds to the initial level size treated in sequential. That is why all the weak scaling curves start at 100% since this scaling is the ratio of elapsed time in sequential divided by the elapsed time in parallel. For the frontal and domain results, nothing can be changed except the starting point. For the TS, the global level may change."two-scale from X above" curves (with X=1, 2, or 3) represents the problem size obtained by always starting from a global level X level above it. I.e. for X=1 with L1, the computed points are:

• L1(3324 dofs) from L0 with 1 process The other curves "two-scale from X" (with X=1 or 2) represent the obtained problem size by always starting from the same global level X. For example, for X=1 with L2, the computed points are:

H.3.3 Task dispatching

Figure H.10 provides some information about the origin of the computational costs associated with the "frontal", "domain" and "two-scale" results for L5 simulation. The elapsed time for each task is accumulated. For each strategy, the identified task covers different items:

frontal: "field declaration" represents the dofs numbering (involving communication), the creation of the local matrix graph and the elimination of the local Dirichlet boundary conditions. "integration and assembly" corresponds to the local elementary integration and the elementary matrix assembly to the local part of matrix A. "Resolution" is the parallel direct sparse solver factorization and the resolution phases.

domain: "field declaration" represents the domain dofs numbering (no communication involved since the domains are restricted to only one process), the Schur dofs numbering (involving communication), the creation of the domain matrix graph and the elimination of the local Dirichlet boundary conditions. "integration and assembly" corresponds to the local elementary integration and the elementary matrix assembly to domain matrices. "Resolution (dense)" is the parallel iterative resolution of the distributed Schur problem. "Resolution " is the incomplete factorization of the direct sparse solver for each domain (creating their Schur contribution), the block Jacobi preconditioner factorization, the right-hand side b condensation and the de-condensation of the solution.

two-scale: "field declaration" represents the patch dofs numbering, the treatment of the patch boundary conditions, the integration and assembly (in sequential or parallel) of the patches, the sequential or parallel resolution of the patches using a direct sparse solver, the global level dofs numbering (involving communication) and the elimination of the Dirichlet boundary conditions. "field declaration update" only covers the patch boundary conditions integration and the sequential or parallel resolution of the patches with an already factorized matrix. It is the update phase of the enrichment functions that happen during the scale loop. "integration and assembly" corresponds to a local elementary integration of SuperPatch on coarse elements and an elementary matrix assembly in distributed local part of matrix A. "Resolution" is the parallel direct sparse solver factorization and the resolution phases at global scale.

H.3.4 Computational precision

In figure H.11, the energy error (from (H.8)) is plotted for the frontal strategy and the theoretical slop is given. The TS curves correspond to the initial error obtained using a global non-enriched field as boundary conditions for the patches. These curves give, for the TS resolution, the error of the "starting point". The scale iteration reduces those errors to be less or equal to the frontal one. In terms of performance in computation time, the state-of-the-art about the elapsed times for all the test cases of chapter VI is given in Table I.1, except for notched beams tests1 (chapter VI.3). Note that the older Spherical Holes test case was not re-computed with AZ. The elapsed times given in the table correspond to a pre-AZ version without condensation.

The AZ was very useful for the spiral bevel gear test case where the model is large (453 824 nodes) from the start. Notably, in the case of the first load steps, the gain is high thanks to the small AZ envelope. As mentioned in annex E.2, the AZ determination can be problematic in some cases, such as in this simulation. Instead of obtaining a reduction in the size of the envelope at the end of the computation as damage front is reaching domain boundaries (like in figure V.2 with the chalk example), the AZ is very large due to several spurious "close crack tip points" detection. The times given in table I.1 correspond to the computation up to load step 246 (figure VI.17) where the AZ envelope is almost the one expected. Additional 4.4 days were used to reach loading step 385. Once more, here, the AZ is not ideal during those last load steps and future work on crack tip detection will hopefully reduce those time consumption.

To have an insight on which part of the computation is consuming CPU The first table presents an overall distribution of the computation time. First, regarding the TLS update task, which is detailed in table I.4, it participates for 34.59 % of the total simulation elapsed time, which is not negligible. The actual main consumers are the computation of Ȳ and φ update. Both of those tasks are not well parallelized and correspond to a rather old algorithm. More recent fast marching techniques used in eXlibris will replace those computations and in particular, the algebraic system resolution attached to them.

One can note that the new double cut algorithm proposed in this thesis does not consume much computation time. It is embedded in the update domain item, which represents only 7.17 % of the total simulation elapsed time without being parallelized in these tests. Secondly, regarding the overall time distribution, the efforts of reducing the mechanical problem resolution via the parallelism and AZ techniques, are well justified because the corresponding task remains to be the most important consumer. As detailed in table I.3, it is mostly penalized by the algebraic resolution (78.91 %: factorization, condensation, solving, ...). The matrix creations (elementary creation and assembly), are not greedy in terms of time consumption despite the use of sub-elements for integration (XFEM technique). The parallelism associated with a proper load balancing (annex F.1) kept the consumption low anyway. Remaining sequential computation such as dof updates will have to be also treated in parallel in the future. .5). Comparison between the sequential and parallel versions with or without AZ. The ratio corresponds to a comparison per column, line and diagonal.

Without AZ With

To conclude, the AZ performances are again illustrated by the spiral bevel pinion gear test case which is an intensive simulation. In table I.5, the elapsed times after 80 load steps are given for different configurations. Upper left results correspond to a sequential computation with no AZ nor condensation. It is more or less the implementation of [START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasi-static loadings[END_REF] in terms of mechanical problem resolution. It is somehow the starting point of this work. On the one hand, the introduction of AZ and condensation (upper right) reduces the time consumption by 6.35. It validates by itself the interest of this technique. On the other hand, applying the parallel computation on an intensive task (lower left) reduces the time consumption by 8.56. Now applying the parallel computation on the AZ technique (lower right) divides the time consumption by 3.25. If the CPU time of mechanical problem computation is extracted from the global time, it is a ratio of 8.08, which is in fact obtained in this case. It illustrates the right choice of the association of parallel technique with the dense condensed resolution. We can add that the memory consumption becomes an issue with the condensation if no parallel distribution of the dense system is made. In parallel, the use of the AZ technique reduces the time consumption by 2.41, which is less than sequential ratio but still interesting. Compared to a non-AZ sequential version, adding both strategies as proposed in this work, reduces time consumption by 20.64.

J

J.1 Material characteristics and parameters

The material characteristics and parameters of the model used in chapter VI are presented in table J.1. For the test case in VI.3, the critical damage energy release rate is a function of the damage. This function is built using an equivalence with a bi-linear cohesive zone model, as stated in [START_REF] Parrilla G Ómez | Comparison between Thick Level Set (TLS) and cohesive zone models[END_REF]. The four extra information mandatory to describe Y 0 c h(d) are coming from the description of the bi-linear softening law given in terms of stress versus opening displacement (figure J.1). [START_REF] Hoover | Cohesive crack, size effect, crack band and work-of-fracture models compared to comprehensive concrete fracture tests[END_REF].

A | C | D | E | G | H | M | P | R | S | T | W | X A aspect ratio
Criterion that gives the finite element quality (annex B).

AZ

Active zone: Solving strategy using condensation techniques presented in chapter V.3.

C core

Basic computational unit of the CPU: A core can only run one program or thread of execution at a time (except with a hardware multithreaded core like Intel chipset).

CPU

Central processing unit also called the main processor.

GFEM

Generalized Finite Element Method.

H hanging node

In Quadtree or Octree context, the hanging nodes exist in a cell boundary when they have no counterpart in neighborhood cells due to cell scale change. See annex C.

h o
Mesh size outside the TLS damage band.

h t
Mesh size in the crack tip vicinity.

h w Mesh size in the crack wake region.

M MPI

Message Passing Interface: normalized parallel protocol that uses independent process on eventually independent computers which communicate by exchanging messages.

P patch

In the TS method, we will define a patch of an entity as the independent fine-scale problem which is embedded into this entity support.

R ρ

Characteristic length associated with the curvature of a crack path . See figure I.1b.

S scalability

The scalability (also referred as the scaling efficiency) indicates how efficient an application is, when using an increasing number of parallel processing elements.

SGFEM

Stable Generalized Finite Element Method.

speedup

The speedup of a parallel program can be defined as the ratio of this program execution time ran with one core by the program execution time ran with N cores.

sSGFEM Scaled Stable Generalized Finite Element Method.

strong scaling

The strong scaling of a parallel program can be defined as its scalability when starting from a given problem size, we increase the numbers of process used. It represents the application capacity to face a huge CPU-time consumption.

SuperPatch

In the TS strategy, it corresponds to a mesh that groups all patches finescale elements.

support

The support of an entity is the set of elements over which the approximation function associated with this entity is non-zero..

SVDF

Signed Vector Distance Function.

T thread-safe

A program or a function is thread-safe if, when it is launched in a multithreaded environment, it does not do concurrent uncontrolled writing access to a unique memory location.

TLS

Thick Level Set.

TS

Two-scale or GFEM gl method.

W weak scaling

The weak scaling of a parallel program can be defined as its scalability when starting from a given problem size, we increase both the numbers of process and the problem size so that the amount of work (dofs in FEM) per computing unit remains constant. It represents the application capacity to face huge memory consumption.

X XFEM

Extended Finite Element Method.

  Figure I.1: General fracture problem

  Figure II.1: Model parameters.

Figure

  FigureII.1a shows the above-defined ingredients. It represents different situations, using a planar cross-section to see inside Ω + , where:

  Figure II.2: graphical representation of Ȳ (in 2D for better understanding but all principles also hold in 3D) : In the small dashed line, φ skeleton. In the large dashed line, some intermediate iso-values of φ. In blue Γ 0 . In red Γ c . Cyan segments are aligned with the gradient of φ. Over these segments, Ȳ is uniform and is the average of the underlying Y field.

  Figure III.1: Location of Ω c (orange) depending on the technique used to construct Γ c and enrichment scheme.

Figure

  Figure III.4: Selection of nodes to be enriched by the ramped Heaviside. Star nodes are enriched once and the 'sun' node is enriched twice. Note that K nodes are not enriched because both cracks lips run on the boundary of their support.

Figure III. 6

 6 Figure III.6 shows the issues with an element extracted from a mesh in a 2D context (with l c = 0.5).Figure III.6a depicts the following. The vertex element

  Figure III.6 shows the issues with an element extracted from a mesh in a 2D context (with l c = 0.5).Figure III.6a depicts the following. The vertex element φ values are : 0.48, 0.497, 0.51. The SVDF construction of triangle upper vertex (based on figure III.2 steps) has discovered the closest point CP on the Γ 0 discretization at a distance of 0.499. This value is in contradiction with the 0.51

  Figure III.5: Decomposition of the support of a node A (a) into sane Ω + parts. Two sane Ω + parts appear in (c) and (d) whereas three sane Ω + parts appear in (e) and (f). In (d), Γ c runs through node A while cutting edges above. In (e), Γ c runs through node A while cutting edges below.

Figure III. 6 :

 6 Figure III.6: Reshaping φ around Γ c . III.6a SVDF construction and edges cut processing, III.6c error on Gauss points G computation corresponding to original φ, III.6d corrections applied to φ . lc = 0.5, 2D

  Figure IV.1: Simple zoning example in 2D

  Figure IV.2: Complex zoning example in 2D

  Figure V.1: Simple AZ envelope example in 2D

Algorithm 3

 3 Scale algorithm scheme: smx maximum number of scale iterations, ε arbitrary numerical convergence criteria, U s solution at global level and scale iteration s, u k∈I g e s set of fine-scale problems solution at iteration s, I g e set of enriched global mesh entities identification number. Starting from the global problem solution U s=0 repeat s ← s + 1 for k ∈ I g e do Impose fine-scale problem k boundary conditions with U s-1 Solve fine-scale problems k to obtain u k s solution end for Compute global enriched problem with u k∈I g e s as enrichment functions to give U s solution.

  Figure V.3: TLS framework imbrication with the TS method presented in a fictitious 2D problem context

  Figure VI.1: Spherical holes test case (dimensions in mm).

Figure VI. 2 :

 2 Figure VI.2: Spherical holes test case: initiation locations at step 4.

Figure VI. 3 :

 3 Figure VI.3: Spherical holes test case: merging (inside view with only skin mesh).

Figure VI. 4 :

 4 Figure VI.4: Undocumented test case (one crack diving under the other): presenting of the correct partial merging of Γ 0 fronts but not Γ c fronts.

Figure

  Figure VI.5: Spherical Holes test case: iso surface after 336 load steps, Γ 0 in blue, Γ c in red.

Figure VI. 6 :

 6 Figure VI.6: Chalk test case.

Figure VI. 7 :

 7 Figure VI.7: Chalk comparison of the final state between simulation and experiment presented in Bordas et al. [2008] (Reprinted from Bordas et al. [2008] Copyright (2007) with permission from Elsevier).

  Figure VI.8: Three-point bending test of a notched beam: generic geometry and boundary conditions.

  Figure VI.9: Load displacement curves for notched beam case Bc and Cb.

Figure

  Figure VI.10: Notched beam case Bc (a,c,e) and Cb (b,d,f) : iso surface at points A B and C of figure VI.9, Γ 0 in blue, Γ c in red.

  Figure VI.12: L shape: force-displacement curves (along Y and Z) with automatic (a) or forced (b) damage initiation.

  Figure VI.13: L shape: TLS results (a and b), gradient method results (c). TLS results are iso surface at the end of computation, Γ 0 in blue, Γ c in red.

Figure

  Figure VI.14: L shape: displacement field (scaled) at load step 120 (with automatic damage initiation).

  Figure VI.15: Spiral bevel pinion gear: mesh, boundary conditions and loading.

Figure

  Figure VI.16: Spiral bevel pinion gear at load step 9: in blue Γ 0 and the little red spot is the initial crack location.

Figure

  Figure VI.17: Spiral bevel pinion gear at load step 246: Iso surfaces Γ 0 (blue) and Γ c (red).

Figure

  Figure VI.15d shows, in red, the clamped face (ring at the end of the long shaft) and the connection sliding pivot (cylinder over the surface of the smaller shaft). Figures VI.15a and VI.15b show general and focus views of the mesh. The mesh is overall much finer on the central tooth. The element size in the refined zone is at least lc 4 . The mesh contains 453 824 nodes.

  Figure VI.18: Spiral bevel pinion gear at load step 385: comparison with Ural et al. [2005].

Figure VI. 19 :

 19 Figure VI.19: Spiral bevel pinion gear at load step 380: displacement field.

  (a) Tetrahedron negative domain polytope splitting: out of initial domain, 3 convex polytopes are created for pattern 13 of figure A.2 (b) Pattern 5 of figure A.2 with 2 cut nodes from the same edge collapsing and one cut node for two edges collapsing on one node of the tetrahedron.

Figure A. 1 :

 1 Figure A.1: Convexity and close node treatment illustrations.

Figure A. 2 :

 2 Figure A.2: Tetrahedron topological cut patterns: the red zone corresponds to a positive domain. On each edge, the red segments are in the positive domain, the blue segments are in the negative domain and the grey segments are in the iso-l c .

Figure C. 2 :

 2 Figure C.2: Hanging situation example in 3D with hanging nodes on edges and faces

  Figure D.1: Unstructured mesh adaptation illustration in 2D sequential context. Elements are split because they fulfill user criterion , 2:1 constraint or hanging node treatment . The fictitious refined location is represented in pink.

  be the set of elements that must be split to respect one level rule if elements of remote S are split. communication S = S ∪ Q, Q being placed first in S until Q = ∅ Split elements of S following S order. Create N bc the maximum number of elements refined across all process communication Clean and update distributed mesh and M communication N bpass = N bpass + 1 until N bc = 0 and N bpass < M axN bP ass Remove hanging node by specific splitting communication to respect the 2:1 constraint if element of S are split. Two loops do it. The first one, local to each process, adds a new element to S via an intermediate set O filled by inspecting the neighbor of S.The drained elements from O are added first in S so that elements are treated in reverse order of dependence.

Figure D. 2 :

 2 Figure D.2: Test mesh dispatching on 8 process. Each color corresponds to a process owning this part of the mesh

  Figure D.3: Refined mesh from both dispatching on 8 process. Each color corresponds to a process owning this part of the mesh

Figure D. 4 :

 4 Figure D.4: Plate refinement strong scaling curves

  Figure E.1 introduces, in a 2D context, the steps to obtain the coordinates of the points close to a crack tip. The principle is to use a previously adapted mesh where the TLS computation has been done (figure E.1a). The nodes, with a support that touches or intersects Ω c , are classified into two categories. Nodes that have support separate in two or more parts by Ω c are designated as enriched nodes. The others, designated as not enriched, are those whose support is not divided by Ω c (only touched or reduced). This classification is presented in figure E.1c. The main idea is to consider that at the tips of the crack, it always exists, at least, one element made only of non-enriched nodes (shown in yellow in figure E.1d

Algorithm 6 Figure E. 2 :

 62 Figure E.2: Complex example of points localization in 2D

  Figure E.3: Problematic example of points localization in 2D

  Figure F.1: Partitioning: a 4 processes example in 2D, each color represents a set of elements treated by a process.

Figure

  Figure G.1: SuperPatch mesh (colored per processes id) plunged in a global-scale mesh.

Figure G. 2 :

 2 Figure G.2: Iso-zero and iso-l c surface created in SuperPatch and colored per processes Id.

Figure

  Figure G.1 shows the mesh of each scale: the global-scale mesh (i.e. the whole beam) and the union of all patches mesh (i.e the SuperPatch). The dis-

Figure G. 4 :

 4 Figure G.4: Damage field for patch 106.

  Figure G.5: Boundary of patch 106: G.5a with SuperPatch colored per process id, G.5b from part boundary colored per process id with interprocess boundary in blue.

  Figure G.7: Global-scale displacement field solutions.

  Figure H.1: Benchmark test: simple plate under volume and surface loading (see (H.7) ). A, B, and C are the corner points used to fix the rigid body modes.

  Figure H.2: Plate meshes

Figure H. 4 ,

 4 Figure H.4, H.5, H.6, H.7 and H.8 present strong scaling results respectively for the final refined level L1,L2,L3,L4 and L5. For each figure we have:

  Figure H.4: Cubic field, Level 1, strong scaling

Figure H. 8 :

 8 Figure H.8: Cubic field, Level 5, strong scaling.

•

  Figure H.9: Cubic field, weak scaling for different starting problem sizes (L1,L2 and L3).

Figure

  Figure H.10: Cubic field, Level 5, Elapse time dispatching

  Figure J.1: Bilinear cohesive zone softening law in terms of stress versus opening displacement.

0

 0 crack tipCrack location where both lips are connected. It is also called the crack front since the crack grows from it. It corresponds to points (in 2D) or lines (in 3D). See figure I.1b . crack wake Crack location where both lips are not connected. It is assimilated to the crack path. See figure I.1b . scalability: for the strong scaling, it corresponds to the ratio of the speedup by the numbers of used cores. For the weak scaling it corresponds to the ratio of the program execution time ran with one core by the program execution time ran with N cores (equivalent to the speedup but in the context of a constant amount of work per process). eXlibris Set of C++ libraries providing FEM, XFEM, TLS, linear algebra, cutting algorithm, tensors, .... tools. This thesis relies on those libraries. G Γ Damage front. iso-0 level set of the signed distance function φ. Γ c Crack lips. iso-l c level set of the signed distance function φ.

  Algorithm 4 Scale algorithm integrated into a general growth loop: tmx maximum number of dissipation load step, CS t crack status at load step t , U t solution at global level and load step t, u k∈I g t e t set of finescale problems solution at load step t, I g t e set of enriched global mesh entities identification number at load step t. To draw a parallel with TLS , ( †) correspond to step (I) and ( ‡) more or less to step (II) and (III) of algorithm 1, and CS t is related to φ level set Starting with a crack status CS t=0 compute global enriched problem solution U t=0 with algorithm 3 repeat t ← t + 1 ( †) Compute new crack shape CS t from U t-1 and CS t-1 for k ∈ I g t e do Set or reset fine-scale problem k with CS t Impose fine-scale problem k boundary condition with U t-1 Solve fine-scale problems k to obtain u k

	t
	end for ( ‡) Compute global enriched problem with u k∈I g t e t	to give U t solu-
	tion.
	until t ≥ tmx or crack stop growing
	do the job. This is what has been proposed in

  2, ..., T }. It is associated with the global node x k , it has the same support as N k ( x) and I g e represents the T global nodes enriched by the TS method. Kinematic equation used to describe the global displacement field is then:

  6.

	Label	L	D	a	S	n	T
	Bc	516	215. 16.125 467.84 1.5 40
	Cb	223.2 93.	13.95 202.37 1.5 40
	Table VI.1: Three-point bending test of a notched beam (dimension in
	mm).

Table I .

 I 2: Computation times dispatching, in hours, of the spiral bevel pinion gear (the 246 first load steps).

	Test case	Elapsed time	Number of processes
	L shape (VI.4)	6.17h	16
	Chalk (VI.2)	5.43h	24
	Spherical Holes (VI.1)	9.25h	24
	Spiral bevel gear (VI.5) 1.6d (up to step 246)	48
	Table I.1: Test elapsed time in hours (h) or days (d) (run on Titan cluster,
	see annex J.2)
	Item	Elapsed time % of total
	Total	37.93	100
	Mechanical problem res-	23.68	62.42
	olution Algorithm
	1 II
	TLS update Al-	13.12	34.59
	gorithm 1 I+III
	AZ creation	0.32	0.84
	Other	0.82	2.16

Table I .

 I 3: Detailed computation times in hours for mechanical problem resolution of the spiral bevel pinion gear(the 246 first load steps).

	Item			elapsed time % of total % of task
	Mechanical	problem	23.68	62.42	100
	resolution Algorithm 1 II
	System creation (Inte-	0.63	1.66	2.65
	gration + assembly)
	Algebraic system resolu-	18.68	49.25	78.91
	tion
	Dof update and results	2.21	5.82	9.33
	output
	Matrix structure cre-	2.16	5.69	9.11
	ation,	parallel task,
	nonlinear resolution, ...
	Item		elapsed time % of total % of task
	TLS update Algo-	13.12	34.59	100
	rithm 1 I+III
	Ȳ computation	5.24	13.8	39.9
	chapter III.4
	Update domain	2.72	7.17	20.74
	(Γ0 and Γc creation)
	Update φ and	5.16	13.61	39.35
	other

Table I .

 I 4: Detailed computation times in hours for the TLS update of the spiral bevel pinion gear ( the 246 first load steps).time in a rather intensive simulation, the test case of the spiral bevel pinion gear was profiled. The results for 246 first load steps are given in tableI.2, I.3 and I.4.

Table I .

 I 5: Elapsed times computation (in hours) up to 80 th load step for the spiral bevel pinion gear (chapter VI

			AZ	ratio
	1 process	113.18	17.82	6.35
	48 processes	13.22	5.48	2.41
	ratio	8.56	3.25	20.64

Table J .

 J Table J.2 gives the values fitted by Hoover and Bažant [2014] and used in this paper. 2: Bilinear cohesive zone softening law coefficients from Hoover

		Values
	σ f	3.92MPa
	σ k 0.588MPa
	w 1	25.3µm
	w f	94.8µm

only this scalar form will be considered in this work

The iso-lc is determined by a -lc offset of φ iso-zero.

in figure perspective

In the current implementation, the initiation corresponds to the introduction of a new damaged location. It is made by merging φ with a level set which is describing a little sphere encapsulating the new default.

s is dimensionless and varies from 0 (edge node origin location) to 1 (other edge node location).

see[Coxeter, 1973, p. 118] for a general definition

Note that in weak discontinuity approach community, the same performance concerns are also addressed. Moreover, the same mixing approach is proposed in[START_REF] Giovanardi | A hybrid XFEM -Phase field (Xfield) method for crack propagation in brittle elastic materials[END_REF] with XFEM in the wake and phase-field patch at the crack tip.

In this thesis most of the time "fine" is used in place of "local" original GFEM gl designation

Future work may prove that some enrichment overlapping can be mandatory to get better results. In this case some nodes can be enriched with both TS and ramp Heaviside function.

Note that in this thesis, the TLS model integration into the TS method has been mainly evaluated in the context of first-order approximation for fields. This leads to talk only of enriched nodes.

note that no boundary condition is imposed in Ωc

Illustrated in Figure V.3d with the global mesh of the 2D simple fictitious example distributed over 4 processes. In this figure, green and blue nodes correspond to patches respectively entirely embedded in one process or split across two or more processes

not shown here but roughly located in front of the studied tooth, locally parallel to its face

solver out of core capability has not been considered. The speedup, in this case, would mainly have measured hardware capability instead of the algorithmic performance

For this notched beams test the AZ envelop is not well computed. This problem combined with too small elements size (ht ≈ lc/20), makes irrelevant the usage of the AZ implementation. Therefore, the elapsed times are not representative.

Acknowledgments

The author is grateful for the research support of the European Research Council through ERC starting Grants ERC-XTLS N°291102. "Ecole Centrale de Nantes" school through the GeM laboratory is also acknowledged for its financial help. The author also wants to warmly thanks Nicolas Mo ës who has allowed this work to exist, to come to fruition and to succeed. The whole GeM MS team is also to be associated with this work for their constant help, advice and support.

A Double cut algorithm details

where length edgei is the length of the i th edge of the tetrahedron and height j is the tetrahedron height from the j th corner. For a regular tetrahedron of edge length a, the height is a 2 3 and then the criterion is 1. The higher this criterion, the worse the quality of the elements (element is rather flat, having edges of very different sizes).

C.1 Origin

If a mesh must be adapted, the element splitting process implies adding nodes in the middle of its edges. It is done independently for each mesh element. However, mid-nodes are merged for edges common to split elements. A hanging node appears when in an edge adjacency, there exists, at least, one element split and one not split. In this case, the mid-node of this edge is said to be hanging. This node is only related to some mesh elements, but not to all the adjacency, which leads for example, to a discontinuous displacement field in mechanics. In 2D, only node related to edge are involved as depicted in figure C.1a. The term hanging edge is then commonly used to talk about an edge that is facing two edges connected to a hanging node. Title : Thick Level Set model implementation in 3D parallel context Keywords : TLS, fracture, parallel, two-scale, mesh adaptation

Glossary List of table and figures

Abstract :

A complex 3D fracture simulation of quasi-brittle material in quasi-static is still hard to tackle nowadays. Many methods and models propose partial solutions to this problem. The Thick Level Set (TLS) model, which uses an approach mixing damage mechanics and explicit crack representation, provides an easy fracture initiation, a complex crack growing capability (coalescing or branching) and an accurate tortuous fracture path. In this thesis we will demonstrate that the implementation of this model in a parallel 3D context provides an accurate and versatile tool that potentially scales. Regarding the accuracy, a novel tool called the "double cut algorithm" has enhanced the existing TLS implementation by letting pass a straight fully damaged zone in a mesh element without conditions on its size. This tool also brings a way to optimize the discretization by coarsening the mesh in a crack front wake.

This adaptation reduces the size of the discrete mechanical problem and therefore the effort for the linear algebra resolution. As far as the scaling is concerned, the bottleneck is the linear algebra resolution time and its associated memory consumption. The parallel solving strategy developed in this thesis to tackle this problem starts first with a basic approach. Then by switching to a method close to domain decomposition and later to a two-scale method, it permits increasing scalability. The other TLS tasks are also partially parallelized. The principal concern is to obtain a tool that either runs faster or can treat a more significant problem if we provide much more computational units. Finally, some test cases illustrate the obtained results with a parallel 3D implementation.