Raphaëlle Chaine

Christophe Geuzaine

Lévy, who saw in me the researcher I could be but tried to ignore and made me discover the wonderful, full of surprises world of computational geometry and mesh generation. I would like to thank Jean-François Remacle in whose brilliant team I had the pleasure to work and where I learned so much about meshing for real engineering problems while tackling some bewildering theoretical questions.

I would like to warmly thank the students I was lucky to work with, Pierre Anquez, Célestin Marot, Kilian Verhetsel, their contributions are the core of this habilitation.

Finally I would like to

Résumé

Les maillages sont omniprésents en informatique car ils permettent de définir la géométrie des objets virtuels, de les visualiser, ainsi que de prédire les propriétés physiques de leurs contreparties réelles. En ingénierie, les maillages ont un rôle central dans la plus grande partie des méthodes de calcul scientifique et sont notamment un prérequis pour obtenir des simulations numériques fiables et efficaces. Les maillages non-structurés, particulièrement flexibles, permettent des variations de taille, d'orientation et de forme de leurs éléments qui sont nécessaires pour certaines simulations. L'objectif de mes travaux de recherche est de développer des algorithmes de génération de maillages non-structurés en trois dimensions, sujet qui combine des défis de géométrie algorithmique, de physique computationnelle et d'informatique. Dans une première partie, je présente mes contributions au pré-traitement des modèles surfaciques dans l'espace 3D ainsi qu'une méthode de remaillage basée sur le diagramme de Voronoï. La deuxième partie se concentre sur la génération de maillages volumiques. Des algorithmes pour générer deux types de maillages non-structurés sont décrits: les maillages tétraédriques avec une approche type Delaunay et les maillages éléments finis hybrides avec une approche dite indirecte. Ces recherches m'ont conduite à l'étude de problèmes théoriques sous-jacents à la génération de maillage, travail présenté dans la troisième partie de ce manuscrit avec mes contributions à l'énumération des subdivisions combinatoires des polyèdres en tétraèdres ou hexaèdres. Le dernier chapitre conclut ces travaux de génération de maillage et propose un ensemble de perspectives combinant une approche théorique et pratique.

Abstract

Meshes are ubiquitous in computer science, they define the geometry of the virtual representations of objects and permit to study the physical properties of their real world counterparts. In engineering, meshes have a key role in most scientific computation methods and are a prerequisite to reliable and efficient numerical simulations. Unstructured meshes allow size, orientation and shape variations, that are necessary for some numerical methods. The objective of my research work is to develop algorithms to generate unstructured meshes in three dimensions, a topic that combine computational geometry, computational physics and computer science challenges. In a first part, I describe my contribution to the pre-processing of surface meshes in the 3D space as well as a remeshing method based on the Voronoi diagram. The second part focuses on 3D mesh generation. Algorithms to generate two types of meshes are described: tetrahedral meshes using a Delaunay based approach and hybrid mesh generation using an indirect approach. This work in 3D led me to study on the underlying theoretical questions of computational geometry which are presented in the third part where my contributions to the enumeration of the combinatorial subdivisions of polyhedra into tetrahedra or hexahedra are described. The final chapter draws some conclusions and presents perspectives of this research work in combining a practical and theoretical approaches in mesh generation.

Introduction 1.1 An Introduction to Mesh Generation

What is a mesh? When watching the latest 3D animation movie or playing your favorite video game, you are probably not thinking about the meshes hidden behind the beautiful landscapes or your beloved characters. But without meshes, you would not be able to see these landscapes or characters whose actual 3D shapes they define. Meshes are virtual geometrical representations of objects, and are ubiquitous in computer science. In engineering, meshes are used for the prediction of the physical behavior of manufactured, or natural objects by numerical simulations. When running physical simulations, meshes are at the forefront, and you should look at them, because they have an important impact on the results. If a mesh is not good enough, computations based on it will fail or the results will be incorrect. Meshes have indeed a central role in most scientific computation methods. The great variability of the objects studied in engineering (e.g. sand grain, mechanical piece, wind turbine, subsurface model) as well as the great variability of physical processes (e.g. flows, wave propagation, mechanics, electromagnetics, heat transfer) and numerical methods (e.g. finite element, finite volumes, spectral elements) are a major challenge for mesh generation methods which should in all these cases generate the best mesh possible, or at least, one good enough.

The mesh generation challenge Generating a mesh may be an easy task: draw a circle, place ten points on this circle, and finally connect the neighboring points by segments. The result is a mesh of the circle. More points may be added to improve the circle shape approximation by the segments. In the general case, the mesh generation problem is significantly more difficult, but the main steps remain the same: (1) process the input shape, (2) place the mesh vertices, (3) connect the vertices, and (4) optimize the mesh. Depending on the constraints on the size, shape, and number of elements, the problem can be easily solved, rather difficult, extremely difficult, or have no solution. Trivial in one dimension, mesh generation is fairly easy in two dimensions. The high difficulty of generating three dimensional meshes is however surprising and widely underestimated. Mesh generation is considered by some practitioners as the most expensive step of engineering analysis in terms of human interaction time. NASA considered in 2014 that mesh generation is one of the six technologies whose development is critical for future numerical simulations stating that "Current meshing methods are neither reliable enough, nor sufficiently robust, and cannot automatically produce high quality meshes at the desired resolution for complex configurations" [START_REF] Slotnick | CFD vision 2030 study: a path to revolutionary computational aerosciences[END_REF].

Defining our terms

In mesh generation, it is particularly important to define precisely the terms we use. Indeed, the difficulty, and sometimes the feasibility, of one mesh is dependent on the use, definition, and understanding of a few words. Given an object defined by its boundaries, generating a mesh of this object consists in generating a set of elements (points, segments, quadrilaterals, triangles, tetrahedra, pyramids, prisms, hexahedra, etc.) whose union is an adequate approximation of the object and whose interiors do not intersect. Depending on the application, a mesh should also enforce a set of constraints on the cell types, total number of cells, their shape, size, and orientation.

The mesh type is defined from the types of its elements and their connectivity. Structured meshes have a regular connectivity, i.e. all their vertices have the same number of neighbors. They are generally constituted of quadrilateral or hexahedral elements. To the contrary, unstructured meshes do not have a regular connectivity. When unstructured 2D mesh are extruded along the third dimension, the mesh is semi-structured (also called 2.5D). The diversity of unstructured meshes is large: the mesh cells can be general polyhedra, or specified polyhedra (tetrahedra, hexahedra, prisms, pyramids). When all the elements have the same type, the mesh type is defined from that element: quadrangular mesh, triangular mesh, hexahedral mesh, tetrahedral mesh, prismatic mesh, etc. When the mesh is constituted of elements of different types, it is said to be hybrid or mixed. A mesh is conformal when the non-empty intersection any of two of its elements is an element common the boundary of these two elements. The mesh is constrained if some of its elements (edges and facets) are imposed to be part of the final mesh before its construction (very often elements the domain boundaries).

Objectives The objective of all mesh generation methods is to generate a valid mesh that fulfills user expectations. The challenge is that this adequate mesh, one that ensures reliable and efficient numerical simulation results for one application, will probably not be adequate for another application. Indeed, meshing is always an exercise in tradeoffs between the characteristics of the mesh that allows the numerical simulations to be faithful to the physics at the desired precision and the efficiency of the computation (in other words the mesh quality). Mesh requirements are typically incompatible and generating a mesh is always a compromise. Geometrical characteristics of the objects to mesh (low angles, thin features and thin sections, high curvature) are in direct conflict with the geometrical characteristics desired for the mesh. The increasing complexity and size of the models, as well as the increasing complexity of the physics simulated are major challenges for mesh generation. This is without even mentioning the question of the existence (and feasibility) of an optimal mesh whose characteristics cannot be known in advance. It is important to keep in mind that (1) mesh generation is intrinsically an ill-posed geometrical problem, and that (2) the difficulty of the problem is infinitely superior in three dimensions than in two dimensions. In that context, it is essential to develop tools with solid theoretical basis that are quite general, but which remain motivated by given applicative needs.

The meshing workflow Despite major differences between applications and desired properties for optimal meshes, practical mesh generation generally follows four steps:

1. Input geometry pre-processing 2. Generate the mesh 3. Optimize the mesh quality (remove a priori problematic elements, fit size map) 4. Adapt the mesh to the solution (refine or coarsen the mesh using a posteriori error evaluation)

Depending on the type of mesh and on the application, the relative importance and difficulty of each of these four steps vary a lot. For visualization purposes, e.g. in computer graphics applications, the primary interest is to automatize the procedure, fidelity to the actual physics not being the crucial point. The trend is therefore to automate the correction of the input models and to develop fast provable meshing algorithms that generate automatically meshes. Such approaches are not (and may never be) reliable options for the engineering oriented meshing community. Practical meshing methods should be robust, applicable to a large range of real models, and allow control to the user. When evaluating complex physical phenomena about complex geometries, strong approximations are not acceptable, e.g. when simulating air flows around a plane wing. In complex cases, generating the mesh is an iterative process based on its adaptation to capture the small crucial features of the numerical solution. The actual objectives of the computational engineering and computer graphics communities are different, and the developed approaches, though complementary, deeply differ in terms of methods and applicability.

Main Contributions

Mesh generation has been the focus of my research work since my PhD thesis. I am more particularly interested in the development of algorithms to generate unstructured meshes in three dimensions. During my PhD thesis I developed mesh generation methods that allowed modifications of the input model, an approach motivated by geological models with very thin layers that vanish laterally. More recently, I have worked on robust meshing volumetric methods which leave the input mesh unchanged: tetrahedral meshing, hybrid meshing and hexahedral meshing. These practical topics lead me to study on the underlying theoretical questions that may be the key to new meshing algorithms.

An abstract point of view on meshing

The problem of mesh generation may be subdivided into two interdependant steps: (i) determining the coordinates of the mesh vertices (ii) connecting the vertices to build valid mesh cells. Given a domain Ω and a set of points V, I study the problem of the construction of a set of cells C that connect all the points of V and enforce the two fundamental properties of a mesh:

1. The union of the cells is equal to the domain: C = Ω (Union property) 2. The intersection of two cells is either empty or is a face/edge/vertex common to the two cells (Intersection property).

These two properties have two components: one combinatorial and one geometrical whose interdependance is one of the main challenge in 3D mesh generation. In theory, the set of meshes that can be built from V, denoted H V , is finite, but is in general impossible to evaluate in 3D. Each mesh with vertices V is a combination of cells in C that verify the union and intersection properties. Each mesh is then a maximal clique of the compatibility graph G of C in which there is one node for each cell of C and one edge if two cells are compatible, i.e. the intersection property is valid for this pair. This approach is however far from being practical, since even for small polyhedra the set of cells C has a prohibitive size.

This abstract view of the meshing problem may however be used to analyse and classify my contributions to unstructured mesh generation methods. Practical methods heavily rely on the Delaunay triangulation and its geometrical dual the Voronoi diagram. The Delaunay triangulation of a point-set in general position, V, is the unique triangulation such that the interior of the circumsphere of any tetrahedra is empty. The Voronoi diagram is defined from its Voronoi cells: the subsets of the space closer to a point p ∈ V than to any other point of V.

Mesh generation from the Voronoi diagram

We propose meshing methods for geological models, Ω that are defined by their triangulated boundary surfaces δΩ, using the intersection of the Voronoi diagram of points V with these surfaces. Our contribution is the definition of strategies to build a mesh whatever are the topology of the intersections between the surfaces and the Voronoi cells. Indeed, depending on the point sampling resolution, a Voronoi cell can cut several model entities and define several connected components. To control the number of cells in the output mesh, our algorithms allow to locally modify the model to reduce the difference between the model Ω level of detail and the point set V resolution. The generated meshes are constituted of different types of cells: quadrilaterals and triangles for surfaces, tetrahedra, prisms, and pyramids for volumes (Section 3.2 and Chapter 5).

Mesh generation from the Delaunay triangulation

We propose an efficient implementation of the sequential 3D Delaunay triangulation algorithm that is about three times faster than existing open-source implementations and is able to triangulate a million points in about 2 seconds. Our main contribution is a scalable parallel version of the Delaunay triangulation algorithm devoid of heavy synchronization overheads. The domain is partitioned using the Hilbert curve and conflicts are detected with a simple coloring scheme. On a AMD EPYC 64-core machine, we have been able to generate three billion tetrahedra in less than a minute. The Delaunay triangulation is the base of a tetrahedral mesh generation process where the input is the boundary surfaces of the domain to mesh (Chapter 4).

Indirect hybrid mesh generation

We propose a method to build hexahedra, prisms and pyramids by merging subsets of cells of a tetrahedral mesh T of vertices V. Our contribution is an algorithm that build the set H T of all possible combinations of tetrahedra of T in hexahedra, i.e. the subsets t ⊂ T that define a topological ball and whose boundary is a triangulation of the cube. The sets of prisms P T and of pyramids S T are determined similarly to obtain H = H T ∪ P T ∪ S T ∪ T . The final mesh is built by determining a clique of the compatibility graph G where there is a node for each cell of H and one edge if the cells are compatible. (Chapter 6).

Enumerating tetrahedral meshes Using the indirect meshing algorithm, we discovered subdivisions of the hexahedron into 8 and more tetrahedra that contradicted previous beliefs on that problem. Given V = {1, 2, 3, 4, 5, 6, 7, 8}, we further investigated the enumeration of the combinatorial tetrahedral meshes of the hexahedron. The set of tetrahedral cells C has 8 4 = 70 elements. The set of all parts of C, denoted P(C) has 2 70 elements. A naive brute force approach testing the validity of each combination succeeds but is extremely time consuming.

To show that there are 174 combinatorial triangulations of the hexahedron up to isomorphism with 5 to 15 tetrahedra, we consider the triangulations of the 3-sphere with nine vertices determined by discrete mathematicians in the 1970's. The second contribution is to demonstrate that exactly 171 of these 174 triangulations have a geometrical realization in R 3 . Our proof consists in (1) solving the underlying combinatorial problem of the point configurations that can realize a given triangulation T , before (2) finding point coordinates in R 3 corresponding to one configuration (Chapter 7).

Enumerating hexahedral meshes We now consider the extremely constrained hexahedral meshing problem, the generation of hexahedral meshes whose boundary facets match a input quadrilateral mesh. For a given number n of vertices V, the set of all possible hexahedral cells C has n 8 elements. We first propose a backtracking search to enumerate the subsets of C that constitute a valid combinatorial hexahedral meshes of a given polyhedron.

We then propose to enumerate directly a subset of all possible hexahedral meshes. The key idea is to exploit the equivalence between quadrangle flips and hexahedron insertion in a mesh. The tree of all sequences of flipping operations is explored, searching for a path that transforms the input quadrangulation Q to the boundary of a cube. For larger meshes the search stops when a sequence transforming Q to a set of pre-computed hexahedral meshes is found. This is the first practical algorithm to build constrained combinatorial hexahedral meshes. We compute small hexahedral meshes of quadrangulations for which the previously known best solutions could only be built by hand or contained thousands of hexahedra.

Our major theoretical result is to prove that an arbitrary ball bounded by n quadrangles can be meshed using only 78 n hexahedra. This significantly lowers the previous upper bound of 5396 n (Chapter 8). Surface Mesh Models

Among the various possibilities to construct and represent geometrical objects in 2D and 3D, boundary representations are one of the most flexible. The idea of the boundary representation (or BRep) is (not a surprise) the representation of a solid by its boundary. In 2D these boundary are segments and curves and in 3D planes or more complex surfaces. Boundary representations are an ubiquitous representation scheme for solids in most geometric modelers, since they are able to represent in a portable light way solids exhibiting a high level of geometric complexity. They are used in many industrial applications as well as to produce computer animations in movies. Dedicated methods and software have been developped for dedicated applications. Software for automotive, shipbuilding, oil and gas, or aeronautics industries may share algorithms but developments are mostly driven by the specific applications. My research work was first motivated by geological applications where the surfaces of the solid models are meshed with triangles (Figure 2.1b). Surface models are particularly adapted to geological models that are generally built from the surfaces delimiting rock layers (horizons, faults, unconformities, etc).

In this chapter, we define the 3D models and related notions used throughout this manuscript, their components, validity, implementation. I describe my contributions to the management of geological models outside geological modeling software.

Surface Models

In a boundary representation, 3D solid models are constituted of a set of volumes (or regions) bounded by a set of surfaces, each surface is bounded by a set of lines, and each line is bounded by two end points. The choice of a datastructure is key to manipulate efficiently these models and the simplest representation is often the best choice. A surface model is built as sets of entities. The base entities define completely the topology and geometry of the model, while additional physical entities define the features that have a meaning for the user (Figure 2.3). In geological modeling, those are rock layers. When implementing the library RingMesh we adopted a simple representation that revealed itself to be similar to the one of Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in preand post-processing facilities[END_REF]. Several alternative computer representations for these models have been proposed in the literature, see for instance [START_REF] Caumon | Building and editing a sealed geological model[END_REF]) and the references therein.

Entities

Base entities Four base entities constitute the surface model: Corner (dimension 0), Line (dimension 1), Surface (dimension 2), and Region (dimension 3). Each base topological entity is restricted to be a simply connected manifold of arbitrary genus. Simply connected means that the entity has a unique connected component (any two points of the entity can be connected by a path contained in the entity), manifold means that the entity does not contain non-manifold point (for example points at a T intersection). Arbitrary genus means that an entity may have holes and internal boundaries. A specific structure stores the model extension. Adjacencies in the model are represented with a bi-directional data structure: each entity stores a list of upward adjacencies and a list of downward adjacencies.

Corner (1 or 2) (1..n) Line (0..n) (1..n) Surface (1 or 2) (1..

n) Region

• A Region is bounded by a set of oriented Surfaces.

• A Surface is either closed (no boundary) or bounded by a set of Lines and is incident to one or two Regions.

• A Line is either closed or bounded by one or two Corners (closed Lines are cut at one of their vertices), and is incident to at least one Surface.

• A Corner is incident to at least one Line.

Physical entities Physical features, that are typically boundaries between different materials (in geology: faults, horizons, fault-horizon contact, stratigraphic layers, etc.), are represented by physical entities. Each physical entity is a group of base entities to which are associated a name and a feature recording its role in the model (in geology: normal fault, reverse fault, horizon, unconformity, boundary of the model). Each Contact corresponds to a group of Lines, each Interface to a group of Surfaces, and each Layer to a group of Regions. Adjacency relationships between physical entities are not stored.

Meshes

In the models considered in this work, the geometry of all base entities is represented by a mesh of the same dimension as the entity. Each Corner corresponds to a 3D point. Each Line is represented by a set of adjacent segments, and each Surface by a polygonal surface mesh. Regions can either be defined by their oriented boundary Surfaces, or by a volumetric mesh. To ensure a correct definition of the model, the geometry of the entities should enforce two conditions: (i) the boundary of each entity is a union of entities of lower dimension and (ii) the intersection of two entities is a union of entities of equal or lower dimension. These conditions are the same than those required to have a Piecewise Linear Complex (PLC), the BRep representation used in Tetgen [START_REF] Si | TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator[END_REF].

Validity

Before any manipulation or processing of a surface model, one must be sure that this model is consistent. Without this verification, problem identification and debugging of algorithms applied on the models is difficult and cumbersome if not impossible. Indeed the visual inspection of a 3D model is generally not sufficient to check its validity. The validity of a model is tightly linked to the chosen representation, to the application and to the applied algorithms. On the meshing side, the most efficient methods will stop or crash because of bad input.

Entity validity Before all, a valid model should be constituted of valid individual entities. When the input models are already meshed, these meshes must be valid. The entity meshes should be valid and conformal meshes, i.e. they are defined by a set of elements (vertices, edges, facets, and cells) such that (i) the interior of each element is not empty, and (ii) the intersection of two elements is either empty or is an element common to their boundaries (Figure 2.4c.&d.). Empty elements are typically edges, facets, or cells incident twice to the same point (Figure 2.4b.). Second, a base entity must not contain any non-manifold point (Figure 2.4e.), and must have a unique connected component. Finally, Regions defined by their boundary surfaces are to be watertight closed volumes. A model that looks good may be invalid, e.g. the model of the Annot sandstones has few local mesh issues break its validity (Figure 2.5).

Model validity

The solid model must have a finite extension, i.e. the region defining the exterior of the model must be tight and closed. Second, the boundary of a model entity can only be a set of entities of the model. All the vertices, edges, or facets on the boundary of a Line, Surface, or Region have to be part of a Corner, Line, or Surface. The third validity condition ensures the consistency of the stored entity adjacencies with the geometrical representation: two distinct entities should intersect exclusively along entities of their common boundary. This implies that there is no intersection between two entities except at points that are on their (geometric) boundaries and that are part of a (topological) entity in their boundary. Some additional validity conditions may be imposed by the application. for example in geological modeling, only fault or fracture surfaces might end in a volumetric region and two distinct stratigraphic interfaces cannot cross one another.

It is important to note that if it may be relatively easy to implement some corrective functions to fix small defects of the mesh and remove duplicated vertices/edges/facets, isolated vertices, degenerate edges/facets [START_REF] Levy | Geogram[END_REF], it is extremely difficult to fix multiple defects in the general case in a robust way.

Extended geometrical validity Recently, in the PhD thesis of Pierre Anquez we extended the definition of validity to capture holes, gaps and thin features in geological models by using entity exclusion zones. We define for each model entity (Corner and Line, Surface ) exclusion zones to materialize volumes that cannot overlap and cannot include any other model entity. Intersections between exclusion zones define invalid features of the input model Section 3.1.1.

Data Structure

A significant part of my mesh generation research was performed for geological modeling purposes.

In geomodeling, the BRep geological models (also called structural models) are most of the time discrete models that are built in dedicated commercial software following specific modeling techniques [START_REF] Caumon | Surfacebased 3D modeling of geological structures[END_REF]). These commercial software do not ease (or try to prevent) interoperability with non-geological modeling tools. That is a major hassle for the geomodeling scientific community since the use, testing and improvement of more general methods (meshing algorithms, numerical solvers, visualization software) is most active in non-commercial software.

There is simply no standard file format for mesh storage. Difficulties are technical, and solutions are known, but when we started RingMesh no implementation was available for this specific type of models. In a joint effort to factorize in a common library a simple and efficient data model and to ease the tedious tasks of reading and writing files describing geological models in various formats we created RingMesh. I developed the very first version of that code to load the surface models of Skua-Gocad (Paradigm 2016) and compute complexity measures (Pellerin, Caumon, et al. 2015). In 2014, building on that very first version, we started the project with Arnaud Botella and I was one of the main developper till september 2016, writing about half of the code. In addition to file format conversion we added functionalities to check models validity, fix mesh errors, plug in meshing software, access to information on model topology and geometry, manage attributes, and visualize the models. A convenient functionality was the reconstruction of the volumetric regions of model for which only the surfaces are available in input. The documented code is open-source and distributed under the modified BSD license.

• Pellerin, J., Arnaud [START_REF] Pellerin | RINGMesh A programming library for geological model meshes[END_REF].

"RINGMesh A programming library for geological model meshes". In: Proc. 17th IAMG conference. Freiberg, Germany

• Pellerin, J., A. Botella, F. Bonneau, A. Mazuyer, B. Chauvin, B. Lévy, and G. Caumon (2017). "RINGMesh: A programming library for developing mesh-based geomodeling applications". en.

In: Computers & Geosciences 104, pp. 93-100

• GitHub repository: https://github.com/ringmesh/RINGMesh RingMesh has been used for the development of meshing methods dedicated to geological models by the team who developed it [START_REF] Pellerin | Automatic surface remeshing of 3D structural models at specified resolution: A method based on Voronoi diagrams[END_REF][START_REF] Pellerin | Automatic surface remeshing of 3D structural models at specified resolution: A method based on Voronoi diagrams[END_REF][START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF], and was integrated in several workflows to: find appropriate boundary conditions in 3D geomechanical restoration [START_REF] Chauvin | Defining Proper Boundary Conditions in 3-D Structural Restoration: A Case Study Restoring a 3-D Forward Model of Suprasalt Extensional Structures[END_REF], to perform homogenization for seismic wave propagation [START_REF] Cupillard | Homogenization of 3d geological models for seismic wave propagation[END_REF], and to determine far-field stress conditions from several borehole observations [START_REF] Mazuyer | Stress estimation in reservoirs by a stochastic inverse approach[END_REF].

The file format nightmare Meshes are compact flexible representation of potentially extremely complex object in 2D and 3D. Software choices to generate and use meshes are tightly linked to the application domain, to the available functionalities and to pragmatic reasons such as user proficiency, code availability, time, or budget. Implementing readers, writers, converters of mesh file formats conversion has been and is still a tedious part of my work. When developing RingMesh (see Section 2.1) we gave in to the temptation of writing our own file format. This was an error because there were equivalent, more general, more largely used file formats.

Synthetic 3D Geological Models

During my PhD thesis, no simple geological model was available to to develop, test, explain research methods. I built a suite of synthetic models to evaluate the impact of typical geological features on meshing algorithms (Section 2.3). These models are simple, most certainly not exactly consistent from a geologist point of view, but they have been extremely useful for testing new methods operating on geological surface models (Figure 2.7). When building them, I never suspected that they would still be used today (to the contrary of some of my meshing algorithms). They are available freely on the website of the Ring consortium. I proposed a suite of six simple models derived from a simple cylindrical anticline composed of three gently folded horizons where each layer has a constant thickness called model A1. In model A2, two regional normal faults affect the anticline. These faults are planar, parallel one to another, cut the whole volume of interest, and have dips close to 60 degrees toward the west. Moreover, they have a constant total slip, corresponding to parallel horizon cutoff lines. In model A3, the regional faults are restricted to an ellipsoid shape and terminate in the model. These faults do not compartmentalize the domain and fault slips vary from a maximum near fault centers to zero at fault tips. In model A4, one fault is regional while the eastern fault dies out to the north. Fault displacements increase to the south. As a result, horizon cutoff lines intersect with a small angle. The faults of model A5 intersect along a branch line, resulting in a Y configuration. Fault slips are regular and the model can be restored to model A1 by rigid block motion. However, the slip on the west fault is close to the thickness of the top layer, which generates thin features in the Allan diagram. Model A6 is obtained by cutting model A4 with a topography surface. This results in numerous isolated layer parts, some of them very small when compared to the model dimensions.

The three other models illustrate challenges arising in other contexts. Model B corresponds to a compressive fault-propagation fold. In the lower part, the fault has a low dip and branches onto a horizontal décollement level. The thrust dip changes to a medium angle (ramp) in the upper part and stops in the upper layer in which the shortening is accommodated by internal layer deformation. Model C is built from the folded basal horizon of model A1 overlaid by onlapping horizontal layers deposited at a low angle. The diapiric dome of model D intrudes and cuts three subhorizontal layers. Except in the intrusion influence area, the horizons are only slightly deformed. Dimensions for model B and C are similar to the ones for model A.

Comparing Geological Models: A Meshing Standpoint

In this section, I present a contribution that is a collaboration with Guillaume Caumon, Charline Julio, Pablo Mejia-Herrera and Arnaud Botella. Our primary objective was to evaluate the relative complexity of several models to allow a fair comparison of algorithms. We analyzed the complexity of 3D geological models by describing and quantifying entities which contribute to the geometrical complexity of the geological structures.

All 3D models, in geological modeling or in other fields, are, by definition, an idealized and simplified vision of reality. How far should each model represent the complexity of the actual object? At what scale? What is exactly meant when a theory, method or algorithm is claimed to address complex domains? The term complex is indeed generally used in a qualitative and subjective way depending on the person making that statement, on his/her education and experience, on the means at his/her disposal, and on the purpose of his/her work.

The measures proposed could be useful for five purposes: (i) to help understand by quantifying our perception of structural model complexity. (ii) to compare geological models (models of the same zone at different resolutions, or models of two different zones); (iii) to improve the objective comparison of geomodeling algorithms and computer codes (an objective comparison of algorithms requires an objective comparison of the models on which they succeed or fail); (iv) as an indicator of the required effort to process these models; (v) parameterize cleaning and meshing algorithms (thin feature size, mesh size).

• Pellerin, J., G. [START_REF] Pellerin | Elements for measuring the complexity of 3D structural models: Connectivity and geometry[END_REF]. "Elements for measuring the complexity of 3D structural models: Connectivity and geometry". In: Computers & Geosciences 76. Publisher: Elsevier, pp. 130-140

Previous work

In geosciences, several papers propose model complexity evaluations using a subdivision of the models into different cells. [START_REF] Andrle | Complexity and scale in geomorphology: Statistical self-similarity vs. characteristic scales[END_REF] uses radius-varying circles to compute a resolution-adapted angle measure and evaluate the complexity of geomorphic lines. [START_REF] Lindsay | Geodiversity: Exploration of 3D geological model space[END_REF] propose a local connectivity complexity measure for models represented as regular Cartesian grids. For each cell, the number of stratigraphic layers sampled by the cell and its adjacent cells is determined, and the average value of this local measure is computed over the model. This local counting principle is also used by box-counting methods to compute fractal dimension [START_REF] Kruhl | Fractal-geometry techniques in the quantification of complex rock structures: A special view on scaling regimes, inhomogeneity and anisotropy[END_REF].

In other modeling fields using a similar 3D representations, several papers propose global evaluations of model complexity (J. [START_REF] Rossignac | Shape complexity[END_REF][START_REF] Sukumar | Towards understanding what makes 3D objects appear simple or complex[END_REF][START_REF] White | Meshing complexity: predicting meshing difficulty for single part CAD models[END_REF][START_REF] Quadros | Finite Element Mesh Sizing for Surfaces Using Skeleton[END_REF]). In geometric modeling, (J. [START_REF] Rossignac | Shape complexity[END_REF] discusses several aspects which contribute to the complexity of 3D models and influence design, implementation, stability, and performance of 3D modeling systems. In visualization, [START_REF] Sukumar | Towards understanding what makes 3D objects appear simple or complex[END_REF]) identify surface variations, symmetry, number of parts, details and topology as the six aspects that influence visual shape complexity. To compute a mesh sizing function, [START_REF] Quadros | Finite Element Mesh Sizing for Surfaces Using Skeleton[END_REF]) identify the small geometrical characteristics in a model from the proximity between two surfaces (distance to medial surface), the proximity between two lines (distance to medial axis), surface curvature, length and curvature of the surface boundary lines.

My Model Is More Complex Than Yours

Global Measures

(1) Number of features The complexity of a model depends on the number physical entities it is representing. In geological modeling they are continuous layers (or rock units), the discontinuities affecting these layers. The more entities, the more difficult and time consuming it is to check and guarantee their validity.

(2) Interactions between features Another obvious statement, is that the spatial distribution of features inside the model has a direct impact on the modeling, meshing, and simulation steps. The more interactions (closeness, intersections) there are between the model entities, the more complex that model is. (3) Geological specificities Conformable layers may be challenging if one of them is locally very thin (Figure 2.8a). Abrupt thickness variations may also raise difficulties when building the model and will also be a challenge to mesh. Vertical relationships between layers controlled by unconformities are very common in stratigraphic domains which host most of the world's natural resources. They often involve thin layers and low angles between horizons (Figure 2.8c) that are extremely challenging when meshing the model, and when simulating physical processes.

Faults are discontinuities inducing a displacement of the rock units localized along a surface (Figure 2.8d). They are often difficult to characterize from subsurface data and introduce significant complexity due to their connectivity, shape, and specific properties [START_REF] Jolley | Structurally complex reservoirs: an introduction[END_REF]). Intersections at small angles of fault-horizons contact lines are extremely challenging for meshing software. Blind faults raise similar challenges because they stop in the model and the nil displacement at their tips is linked to very small angles between fault-horizon contacts in the Allan diagram (Figure 2.8e). When considering a fault network, the total complexity depends on the variability of the fault orientations and on the number, angle, length, and orientation of branch lines. These features have a major impact on the approximations required for meshing, e.g. Y-shaped configurations (Figure 2.8f).

Local Measures

We can simply defined a model complexity as the sum of the complexity of its components. To evaluate the complexity of the components we propose local and global measures of their connectivity (topology) and geometry.

(4) Number of components This simple measure counts the number of entities of the model (regions, surfaces, lines, and corners) excluding the ones defining the volume of interest, and giving the same weight to each entity.

(5) Variability of the component sizes This measure evaluates the complexity of the entities of a given type as a statistic on their sizes. We choose the coefficient of variation (mean over standard deviation), that characterizes the relative distribution of entity sizes and evaluates scale changes for one type of entity.

(6) Geometry We saw in the previous section that sizes, thicknesses, curvatures and angles are four important points of the geometrical complexity of structural models. We then propose to compute the geometrical complexity of each entity as the sum of four measures characterizing (1) its size C s ; (2) its shape C f ; (3) its thickness C t ; and (4) its angles C a . Corners are not explicitly taken into account in this measure, but have an indirect impact on the values obtained for the lines and surfaces. Each of the elementary measures is chosen so that it has values between 0 and 1. The measure are made relatively to a characteristic size h given by the usage measures and a characteristic angle α (Figure 2.9) and are defined as:

• C s compares the size of a component (line length, surface area, or region volume) to a characteristic size h given by the user.

• C t relates to the thin features of the entity. It is defined as the percentage of the boundary length for which the thickness of the entity is smaller than h.

• C f globally evaluates the line and surface deviation from a linear object. It is taken equal to 1 minus the size of the projection of tjhe entity on its mid-segment or mid-surface divided by the size.

Figure 2.9: Geometrical measure computations in the plane. Thickness measure C t and angle measure C a for a region R and shape measure C f for one of its boundaries B 2 .

• C a evaluates the percentage of the line length (respectively corners) where the angle between two surfaces (respectively two lines) is inferior to α. Model complexity can then be evaluated from statistics (mean, coefficient of variation, interquartile range, percentiles, entropy) on these values.

Comparing the Synthetic Models

The comparison of the models constructed in Section 2.2 using the different measures gives different classifications of the models (no surprise there) by increasing complexity:

• Number of entities: A1, C, D, B, A2, A3, A4, A5, A6

• Geometry complexity with h = 100m, α = 20 degrees: A1, C, B, A2, D, A4, A3, A5 and A6

• Size variation measure: A1, D, B, C, A2, A4, A3, A5, A6

• Volume sampling with 10,000 cells: B, D, A1, C, A3, A2, A4, A6, A5

One of the main advantages of local sampling measures is the possibility to understand the spatial organization of the complexity and estimate the extension of the zones where a given geomodeling method will fail from the cells in which the total number of entities is above a given number. These cells are the ones closest to the thin features of the models, but also to their corners and contact lines. Combining the number of cells that contains more than 6 elements and the maximum number of entities in all cells permits to establish a visual classification of the models (Figure 2.10).

These measures are far form being perfect and the question is which alternative do we have that would be able to measure volumetric features, keeping local signal like minimal distance that are crucial for meshing algorithms, and be globally consistent. Chapter 3

Automatic Correction

In general, the input surface models designed during long hours with a CAD software (or another modelling software) are not valid inputs for meshing algorithms. And when they are, you can be sure that the level of detail of the surface model is not the one needed to generate the mesh. This is the eternal challenge faced by engineers. In practice, most models are manually modified by skilled practitioners who are able to ensure that the generated mesh will fulfill the project requirements. Most automatic strategies simplify the model before generating its mesh, but it is extremely challenging to ensure the validity of the output model or its consistency for the application.

Correction before Remeshing of 2D Models

This section presents the work of the PhD of Pierre Anquez that I co-advised with Guillaume Caumon and Bruno Lévy. We proposed an automatic method to repair and simplify two-dimensional (2D) geological models: geological maps and cross-sections. The output model enforces specific practical quality criteria on the model topology and geometry to allow the generation of a mesh enforcing minimum angles and minimum edge size constraints. To detect the features of the input model to correct, we define exclusion zones around input model corners and lines.

• Anquez, P., J. Pellerin, M. Irakarama, P. [START_REF] Anquez | Automatic correction and simplification of geological maps and cross-sections for numerical simulations[END_REF]. "Automatic correction and simplification of geological maps and cross-sections for numerical simulations". en. In: Comptes Rendus Geoscience 351. 1, pp. 48-58 • Anquez, P. (2019). "Correction et simplification de modèles géologiques par frontières : impact sur le maillage et la simulation numérique en sismologie et hydrodynamique". PhD Thesis

Automatic Correction in 2D

Overview The method is based on a graph to encode the input model entities, their connectivity and the invalidities (Figure 3.1). This global topological information ensures the model consistency.

The method operates on that graph to keep track of all the modifications to build a new model in which invalid features are removed. Our strategy disconnects the graph edition from the actual geometrical model editing. The advantage is that all corrections are performed once all invalid features are identified, and all decisions are taken when operating on the geometry. The input model is taken as reference when building the final geometry. Avoiding mixing topological and geometrical modifications is an advantage, because it leads to instabilities or inconsistencies (e.g. [START_REF] Euler | A new tool to seal a 3D earth model: a cut with constraints[END_REF][START_REF] Pellerin | Automatic surface remeshing of 3D structural models at specified resolution: A method based on Voronoi diagrams[END_REF].

The input models are edited in three steps (Figure 3. Build the model graph The colored graph we use is an incidence graph to which we added edges representing invalid configurations. To each entity of the input model ( Corner and Surface) correspond a node of the model graph. There is a blue edge connecting two nodes if corresponding entities are incident to one another. There is a red edge connecting two nodes if the exclusion zones of the corresponding entities intersect (Figure 3.1). No red edges, means that the model is valid.

The construction of the graph is straightforward. To consider topological issues in the input model (cf. Section 2.1.3), such as horizon free borders, we use the graph to detect them and compute the nearest model entity so that a red edge may be built between the two. We implemented the exclusion zone of Lines as the union of the disks centered on all line points. We implemented the exclusion zones of corners as disks centered at the point. Intersections between exclusion zones can be detected by computing convex shape intersections, e.g. using the GJK algorithm [START_REF] Gilbert | Computing the distance between general convex objects in three-dimensional space[END_REF]. We identify all the entity points (not only its vertices) whose corresponding disks compose the intersections between exclusions zones. This geometrical information is stored alongside the graph. The detection of intersections can be implemented using exact predicates (e.g. Lévy 2015) for increased robustness.

Fix topology before geometry The next step is edit G to remove all red edges. This is the key idea of the method: operate on the connectivity of the model, encoded in the graph before modifying the defining the geometry of its entities. We perform successive graph elementary operations that correspond to operations on the geological model which will be actually performed at the next step, see Figure 3.3. The two geometrical operations that do not change the model topology aim at increasing the distance between entities. The three topological operations, merging, splitting and removal, modify the number of model entities and their connectivity; they aim at deleting thin model features by collapsing them.

(a) Edge deletion. Removing the edge connecting to entities corresponds to a modification of the model geometry that increase the distance between them and remove intersection between the excelusion zones.

(b) Edge contraction. Merging the two entities is only possible if both entities are of the same type. This corresponds to locally merge of two model entities.

(c) Node removal. Removing one entity (and subsequent red or blue edges) corresponds to the removal of a minor entity or a duplicated model feature. Prior input knowledge is required to determine which entity is removed.

(d) Node splitting in two nodes. This operation on Lines corresponds to its subdivision in two connected parts, e.g. to fix a hanging nearby corner.

The choice between the possible operation is constrained by (i) entity types (ii) a global strategy to preserve the topology or one that largley allows modifications, (iii) the physical meaning of the entities.

We first process red edges connecting a Line and a Corner, then the other edges, because this is the typical configuration for a horizon free-border invalid edge. Because the collision of Line exclusion zones may not concern the whole Lines, operations on Line-Line invalid edges like edge contraction or node removal may need, as a pre-process, a node splitting.

Reconstruction of the geometrical model

To determine the geometrical represention of the model entities encoded in the graph, we rely on the input model entity on which operations were performed. We are first set Corners, their output position is determined from the input model entities listed alongside the graph. It is set at the barycenter of input Corners, if Lines are involved that barycenter is projected on each of the Line part, and the barycenter of these projections is chosen. To set Lines barycenter points of the Lines tomerge are computed. The final model can be meshed using various tools, e.g., Triangle (J. R. [START_REF] Shewchuk | Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator[END_REF], Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in preand post-processing facilities[END_REF]) or mmg2d1 .

Correction of 2D Models for Numerical Simulations

Sealing of a 2D geological model By repairing non-watertight vertical cross-sections, the method allows to mesh the model and run numerical simulations on it. A cross section in the 3D geological surface-based model built by [START_REF] Collon | 3D geomodelling combining implicit surfaces and Voronoi-based remeshing: A case study in the Lorraine Coal Basin (France)[END_REF]) is given as an example Figure 3.4. The input model is not watertight; there are small gaps and overlaps at contacts between coal bed (pseudo-vertical lines) and fault (oblique lines) and boundary.

To seal the 2D model, minimal local entity size and minimal angle values are set to zero because the objective is to repair without simplification. We use additional geological knowledge: the coal veins are conformable one to another (they can not intersect). This constraints our method by preventing mutual branching of vein Lines. As a consequence, vein free borders can only be connected to the fault (red), to the Permo-Triassic base horizon (green) or to a boundary line. The generated output 2D model (Figure 3.4b) is watertight. The execution time for repairing this cross-section is 0.1 7 s on a laptop with a 2.50 GHz processor Intel R Core TM i7-6500U.

Accelerating numerical simulations An application presented in the manuscript of Pierre Anquez [START_REF] Anquez | Correction et simplification de modèles géologiques par frontières : impact sur le maillage et la simulation numérique en sismologie et hydrodynamique[END_REF] shows how automatic correction of 2D geological cross sections can accelerated the numerical simulations for seismic risk evaluation. This work is a collaboration with Nathalie Glinsky, Diego Mercerat, and Paul Cupillard. Site effects are amplification phenomena of seismic waves in geological layers close to the surface in basins. They can have major consequences in terms of human lives when a seism occurs, e.g. the Michoacán seism on September 19, 1985 in Mexico. In south-east France, Nice is localized in an area where seisms are a risk and site effects can appear. The geological cross-section considered here, is a 2D profile crossing the alluvial bassin of the lower Var valley. Site effects are computed with the simulation of planar waves with a Discontinuous Galerkin solver on a triangle mesh [START_REF] Peyrusse | A nodal discontinuous Galerkin method for site effects assessment in viscoelastic media-verification and validation in the Nice basin[END_REF]. This solver is time explicit, viscous-elastic, and the time step is constrained by smallest triangle height over the P-waves velocity to enforce the CFL condition. The reference model is a cross-section built from the 3D model of [START_REF] Rohmer | Combining borehole log-stratigraphies and ambient vibration data to build a 3D Model of the Lower Var Valley, Nice (France)[END_REF]. A triangular mesh of this reference cross-section is generated with Mmg2d. The simulation of planar shear wave propagation takes 18 days. The details on the model construction and on the parameterization of the simulation are available in [START_REF] Anquez | Correction et simplification de modèles géologiques par frontières : impact sur le maillage et la simulation numérique en sismologie et hydrodynamique[END_REF].

We compare this reference simulation result with the simulation on two models automatically simplified with our method. The first (M5) is obtained using a minimal size of de 0.18m and a minimal angle of 5 degrees. Complementary manual operations to remove small components were performed. The model was remeshed keeping the mesh size of the reference model. The second model (M6) is a further simplification of the first using the set of modifications to edit the model geometry at a fixed topology with a length of 0.18m, but changing the angle to 10 degrees. On Table 3.1, we can see that the simulation time is reduced from 18 days to less than 13 hours for (M5) and to 8 hours for M6. The simplifications targeting the small features increased the triangle height, impacting directly the time step required for the simulation. The comparisons of the results for the simplified models Figure 3.5 show that the results are slightly different but that the site effects are similar. The same frequencies are amplified at the same level leading to the same lengthening of the signals in the basin.

Models simplified automatically have strong advantages compared to the reference models or to models simplified by hand: (1) computational timings goes from 18 days to 6 -12 hours (2) simplifications are controlled by objectives geometrical criteria, they are parsimonious and local.

Simultaneous Surface Remeshing and 3D Model Correction

A more risky alternative to correct models is to authorize model topological modifications while remeshing. During my PhD thesis, I proposed such a method to remesh the surfaces of 3D model surfaces. From a set of conformal triangulated surfaces that are in contact along given lines and at given points, the method computes a set of surfaces meshed with triangles as equilateral as possible.

The method relies on a global Centroidal Voronoi optimization to place the vertices of the final surfaces combined with combinatorial considerations to either recover or simplify the surfaces, lines and points of the input. When the final resolution is sufficient, the input contact lines, and points are also contact lines and points of the final model. However, when dealing with models with complex contacts, resolution may be insufficient and instead of a refinement strategy that may lead to too many points, we propose to locally merge some features of the input model. This ability to simplify the input model is particularly interesting when the model is to be volumetrically meshed.

• Pellerin, J., B. [START_REF] Pellerin | Topological control for isotropic remeshing of nonmanifold surfaces with varying resolution: application to 3D structural models[END_REF]. "Topological control for isotropic remeshing of nonmanifold surfaces with varying resolution: application to 3D structural models". In: Proc. 13th IAMG conference. Salzburg, Austria Previous Work Whether it be in geological modeling or in other applications, the classical strategy to remesh a possibly high number of surfaces and their contact lines is to do so one after another. First, remesh each contact line, then parameterize each surface part to place it in a 2D space and remesh it with adequate strategy is a robust worklow implemented for example in Gmsh. The topology of the model is fixed and no modification is allowed, that is a key advantage when the input surface model is valid and at the correct level-of-detail, but may be a limitation if that is not the case. The second limitation of that strategy is that as all entities are meshed independently even when they are extremely close in space, the volumetric gridding is not the best possible.

An alternative approach was proposed in the computer graphics community to to generate high quality meshes by optimizing the coordinates of a fixed number of points to minimizing a Centroidal Voronoi Tessellation objective function [START_REF] Lloyd | Least squares quantization in PCM[END_REF][START_REF] Du | Centroidal Voronoi Tesselations: Applications and Algorithms[END_REF][START_REF] Alliez | Centroidal Voronoi diagrams for isotropic surface remeshing[END_REF][START_REF] Valette | Generic remeshing of 3D triangular meshes with metric-dependent discrete voronoi diagrams[END_REF][START_REF] Yan | Isotropic Remeshing with Fast and Exact Computation of Restricted Voronoi Diagram[END_REF]Levy and Y. Liu 2010b;[START_REF] Chen | Isotropic Surface Remeshing Using Constrained Centroidal Delaunay Mesh[END_REF]. Then, they locally recover theconnected components of the model (surface parts, contacts, and triple points) by adding one point per connected component. This simple strategy is very efficient to remesh thin layers and parallel contact lines, but fails when the resolution is not sufficient to capture other small-scale features in particular close triple points. We chose to modify these features, i.e. modify the model topology.

Remeshing methods that allow topological changes can be grouped in two main categories: those doing local modifications and those operating globally. The local modifications proposed by [START_REF] Garland | Surface simplification using quadric error metrics[END_REF]) generalize edge contractions and perform both close vertex contraction and edge contraction at the same time. Vertex clustering, introduced by (J. R. [START_REF] Rossignac | Multi-resolution 3D approximations for rendering complex scenes[END_REF], was used to remesh discrete fractures network by [START_REF] Mustapha | An efficient method for discretizing 3D fractured media for subsurface flow and transport simulations[END_REF]. A fixed size box scans the initial model in the three directions and the vertices inside it are replaced by a single vertex placed at the center of the box. Methods operating globally consider the entire model for the simplifications and subdivide into cells whose dimensions control the degree of simplification. For closed surfaces, [START_REF] Andujar | Topology-reducing surface simplification using a discrete solid representation[END_REF]) use an octree, flag each cell as inside or outside the surface, and reconstruct a simplified closed surface from the remaining cells. Considering globally the model and working with a volumetric subdivision of the space permit to analyze locally the relationships between its different parts.

Overview

Input Triangulated surface models following the definition given in Section 2.1. The method is resilient to input model errors (non-water tightness, holes, non-conformal triangulation of the input surfaces along contact lines) but does not guarantee to fix them.

Output Our method computes a global remeshing of all the surfaces. The triangles of the output mesh are as equilateral as possible. Their quality does not depend on the input mesh quality and the method is resistant to degenerated triangles (skewed elements). The conformal contacts between the surfaces remain conformal, the non-conformal ones are kept non-conformal. Some modifications are done to simplify the model where contacts lines are too close.

Principle We use a Centroidal Voronoi Tessellation to place adequately a given number of points near the surfaces and contact lines of the model. A topological control is then used to determine the vertices and triangles of the final surface from the intersection of the Voronoi diagram of these points with the model while recovering the model components.

Fundamental Geometrical Objects

There are two fundamental geometrical objects that are ubiquitous in computational geometry and mesh generation (and in this manuscript): the Voronoi diagram and the Delaunay triangulation

Voronoi Diagram and Restricted Voronoi Diagram

A Voronoi diagram V D(S ) is defined relatively to set of points S To each point correspond a Voronoi cell: the subset of the space R d closer to a point p ∈ S than to any other point of S (e.g., [START_REF] Aurenhammer | Voronoi diagrams-a survey of a fundamental geometric data structure[END_REF] In R d , if we consider the Euclidean distance . , the Voronoi cell of p is defined as V p = {x ∈ R n | px ≤ qx , q ∈ S }. Voronoi cells are convex closed polyhedra that may bounded or not, and cover the space without overlapping (Figure 3.6). The Voronoi diagram was initially defined by [START_REF] Voronoï | Nouvelles applications des paramètres continus à la théorie des formes quadratiques. {D}euxième mémoire. {R}echerches sur les parallélloèdres primitifs[END_REF].

A Voronoi diagram subdivides the space and all the objects Ω included in that space. For a set of sites S and an object Ω the Restricted Voronoi diagram is defined as the intersection of the Voronoi diagram of S with Ω. The intersection of a Voronoi cell, V p , with the object Ω is the restricted Voronoi cell of p to Ω and is defined by V p∩Ω = V p ∩ Ω. Note that to the contrary of Voronoi diagram elements, the elements of a restricted Voronoi diagram can have several connected components (Figure 3.9).

Remark Many useful Voronoi diagram generalizations have been developed by changing the distance measure between two points, changing the site natures etc, see the reviews of [START_REF] Okabe | Spatial tessellations: concepts and applications of Voronoi diagrams[END_REF]) and [START_REF] Aurenhammer | Voronoi diagrams-a survey of a fundamental geometric data structure[END_REF]. It is also possible to work in non Euclidean spaces, for example the one defined by a surface embedded in a 3D space.

Delaunay Triangulation and Restricted Delaunay Triangulation

A triangulation T (S ) of the n points S = {p 1 , . . . , p n } ∈ R d is a set of non overlapping simplices that covers exactly the convex hull Ω(S ) of the point set, and leaves no point p i isolated. The Delaunay triangulation DT (S ) of a point set S has the fundamental geometrical property that the interior of the circumsphere of any d-simplex is empty, i.e. it does not contain any point of S . If the point set S is in general position (in R d if it contains no group of d + 1 coplanar points and no group of d + 2 cospherical points) then its Delaunay triangulation exists and is unique. To manage robustly degenerated cases is mandatory in mesh generation. These degeneracies disappear with an extremely small perturbations of site positions, that may only be symbolic [START_REF] Edelsbrunner | Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms[END_REF][START_REF] Lévy | Robustness and Efficiency of Geometric Programs The Predicate Construction Kit (PCK)[END_REF].

The mathematical properties of Delaunay triangulation, initially proposed by [START_REF] Delaunay | Sur la sphere vide[END_REF]), made it a favored object of study in computational geometry and meshing, see e.g. the excellent review in the book [START_REF] Dey | Delaunay Meshing of Piecewise Smooth Complexes without Expensive Predicates[END_REF]. Its well known geometric dual relationship with the Voronoi diagram is the base of many meshing algorithms Figure 3.7. In 3D, a Voronoi point correspond to a Delaunay triangulation cell, an edge to a face, a face to an edge, and a cell to a point.

The Restricted Delaunay Triangulation (Edelsbrunner and Shah 1997), is the geometric dual of the restricted Voronoi diagram (Figure 3.8). It is the subset of the elements of the Delaunay triangulation of the points S such that the dual restricted Voronoi cell/face/edge/point is not empty.

Surface Remeshing Using the Restricted Delaunay Triangulation

The restricted Delaunay triangulation of given points S to an input surface remeshes this surface. However, it does not necessarily have the same topology as the input surface (not homeomorphic), see for example Figure 3.8. To ensure that the restricted Delaunay triangulation is homeomorphic to an input manifold surface it is sufficient to verify the topological ball property [START_REF] Edelsbrunner | Triangulating topological spaces[END_REF]): all the restricted Voronoi cells (respectively facets and edges) of S to Ω are topological disks (respectively segments and points) and all the restricted Voronoi cells (respectively facets and edges) of S to the boundary of Ω are topological segments (respectively points and the emptyset). A geometrical mean to enforce the topological ball property is to have an ε-sampling of the input surface [START_REF] Amenta | Surface reconstruction by Voronoi filtering[END_REF]: if for each point x of the surface Ω, there is a sampling point of S at a distance smaller than ε × l f s(x), where ε < 0.3 and l f s denotes the local feature size (distance to the medial axis of Ω).

Surface Remeshing à la Voronoi

This section details the two main steps of the surface remeshing method: the optimization of a given number of sampling points on the model (Section 3.2.3.1) and the building of the final mesh from the connected components of the restricted Voronoi diagram of the sampling points (Section 3.2.3.2).

Optimization of the Mesh Vertex Coordinates

First, a fixed number of points are placed so that they are a good sampling of the surface model. Each point samples the model in the sense that it represents the part of the model closest to it than to any other point. The input number of points then determines the resolution at which the model will be remeshed, it may be computed from the square root of the model area divided by the target edge length. To have a good sampling the idea I developed was to optimize the point positions such that their restricted Voronoi diagram to the model is centroidal.

Centroidal Voronoi Diagram

The space subdivision determined by a Voronoi diagram of sites randomly distributed is random, and its optimization appears in numerous statistics, image processing, or mesh generation problems. The goal is to optimize point placement to reach a specific objective, for example, minimize the distances between each one of the sites and points inside their Voronoi cell. This optimization tends toward a specific Voronoi diagram: a centroidal Voronoi diagram.

The Voronoi diagram of points S is centroidal if each site is at the centroid p * of its Voronoi cell V p , with ρ the density function we have:

p * = V p yρ(y)dy V p ρ(y)dy (3.1)
Let's consider the problem of the computation of the partition of a domain Ω in k regions Ω i and the positions of k points s i that minimizes the function:

F((s i , Ω i ) i=1...k ) = k i=1 y∈Ω i ρ(y)||y -s i || 2 dy (3.2)
This function evaluates the sum of the square distance between points s i and points of the region Ω i with the same index. [START_REF] Du | Centroidal Voronoi Tesselations: Applications and Algorithms[END_REF] show that, to minimize this function, it is necessary that the regions Ω i are the Voronoi cells of points s i and that each point must be at the centroid of its Voronoi cell. They also show that this function has the same minimums than the function in which Ω i = V i . Function F parameters are only point positions, the integration being done on the Voronoi cells V i . To compute a centroidal Voronoi diagram it is then sufficient to minimize this function of the site positions. In practice, it is very difficult to obtain a global minimum, and a local minimum is often considered satisfactory. Note also that, in the same space, for a given number of sites, there are multiples centroidal Voronoi diagrams, all of them minimizing the function F.

Restricted Centroidal Voronoi Diagram

The restricted voronoi diagram of a point set S to a domain Ω is centroidal if each site p is at the centroid of its restricted voronoi cell. By integrating function F only on the restricted cells we have:

F Ω = k i=1 y∈V i ∩Ω ρ(y)||y -s i || 2 dy (3.3)
When points s i belong to domain Ω, the restricted centroidal voronoi diagram is constrained [START_REF] Du | Constrained centroidal Voronoi tessellations for surfaces[END_REF]. A site optimization example on a sphere is given Figure 3.10. After the optimization restricted Delaunay triangles are almost equilateral. The centroidal Voronoi diagram is one of the key elements of the variational meshing methods that consist in optimizing the positions of all final mesh vertices before building their (restricted) Delaunay triangulation (Figure 3.10) [START_REF] Alliez | Variational tetrahedral meshing[END_REF][START_REF] Tournois | Interleaving Delaunay refinement and optimization for practical isotropic tetrahedron mesh generation[END_REF][START_REF] Tournois | Optimisation de maillages[END_REF][START_REF] Dardenne | Variational tetrahedral mesh generation from discrete volume data[END_REF]Levy and Y. Liu 2010b). The two features distinguishing this approach from more classical tetrahedron meshing methods is that (1) the number of vertices can be fixed and that (2) tetrahedra shape and quality are globally optimized by the objective function.

Implementation The three steps to perform the optimization of a given number of points over a model Ω are:

1. The initial random placement of the sampling points on the model surfaces is done using the algorithm given by [START_REF] Lévy | Variational Anisotropic Surface Meshing with Voronoi Parallel Linear Enumeration[END_REF].

2. The computation of the restricted Voronoi diagram is done using the fast parallelized method also described by [START_REF] Lévy | Variational Anisotropic Surface Meshing with Voronoi Parallel Linear Enumeration[END_REF].

3. The contributions of each cell of the restricted Voronoi diagram to the objective function and to its gradient is computed following [START_REF] Yan | Isotropic Remeshing with Fast and Exact Computation of Restricted Voronoi Diagram[END_REF] for the CVT energy and the above for the boundary term. The minimization of the objective function F is done with a L-BFGS algorithm [START_REF] Nocedal | Updating quasi-Newton matrices with limited storage[END_REF]. The optimization can be stopped when the norm of the gradient is inferior to a given value. From our experience, convergence is very fast, and in practice we stop the optimization process after 100 iterations. Specific convergence rates are discussed by (Y. [START_REF] Liu | On Centroidal Voronoi Tessellation -Energy Smoothness and Fast Computation[END_REF]. In all the cases we considered the input mesh resolution does not impact the convergence while increased feature density slightly decreases the convergence.

Once the points have been optimally distributed, we compute their restricted Voronoi diagram to the structural model and compute the connected components of the restricted Voronoi cells to determine the vertices and triangles of the output mesh.

Connecting the Dots

Surface remeshing Let's first consider a (non-geological) model in which the different surface connected components do not intersect and have no boundary, e.g. two nested spheres Figure 3.11. The two surfaces are sampled by 100 points whose optimized positions are between the spheres. To recover the input surface parts, we set one point for each connected component of the restricted Voronoi cell. There is then one triangle to build for each point shared by three restricted Voronoi cells (Figure 3.12). The obtained mesh is a dual of the connected components of the restricted centroidal Voronoi diagram, it is closer to the input mesh than the restricted Delaunay triangulation (one sphere in this case). The multi-nerve theorem gives that it is homotopy equivalent to the input model (Colin de Verdière et al.

2012).

Line remeshing Now, consider a surface with a boundary. Similarly to what happens for surfaces, the boundary is not correctly remeshed if the number of points is too small.. To remesh them correctly, each point is replaced by as many points as the number of connected components of the intersection between its Voronoi cell and the boundary lines (Figure 3.13).

These additional points challenge the triangle building step since one connected component of one restricted Voronoi cell may correspond to several points. As a consequence, the dual of points shared by three cells may not be a triangle and additional polygons, dual of the edges shared by two cells and intersecting twice the boundary, should be built between close boundary lines (see the gray segments on Figure 3 .13c & d).

The more intersections between the restricted Voronoi cell and the model boundary, the more vertices in this cell. This may lead to configurations where the polygons to build intersect (Figure 3.14a). To avoid this, if there are more than two points for a connected component of a restricted Voronoi cell, they are merged. This makes our method more robust, but at the cost of modifications of the surface connections that are questionable from a geological point of view and depend highly on the optimized point positions.

Corner remeshing The last entities to account for are Corners (triple points), i.e. points at the intersection of at least two contact lines. To recover all the corners, one should put one point for each of them. So, to fully reconstruct the input model, we need to have for each restricted Voronoi cell one point per triple point, one point per contact, and one point per surface part.

When there are more than one Corner in a connected component of a restricted Voronoi cell to the boundary, i.e. the final resolution is not sufficient, we merge them (Figure 3.14b). When this merging operation is done, the previously described merging is also performed. This way, each restricted Voronoi cell part has 1 or 2 points and the quads or triangles to build with these points do no intersect. The last modification is the merging of the vertices that correspond to close features, close meaning that the distance between them is inferior to a specified input value (Figure 3.14c & d). This is a way to make the model easier to mesh and simplify very small features by removing small fault throws and joining fault tips close to another fault.

Implementation The input of the method is a restricted Voronoi diagram, a polygonal surface obtained from the intersection of the input model with the Voronoi diagram of the optimized points. Each polygon of the restricted Voronoi diagram is the intersection of one triangle of the input model and one Voronoi cell whose ids are known.

1. Each restricted Voronoi cell and its connected components are built using that information.

2. The vertices to put for each restricted Voronoi cell are computed.

3. The last step is to build the cells linking these vertices.

Discussion

Some of the remeshing results of this Voronoi based method are given on (Table 3.2, Figure 3.17, Figure 3.22). For more models, see the original paper [START_REF] Pellerin | Automatic surface remeshing of 3D structural models at specified resolution: A method based on Voronoi diagrams[END_REF]. With input mesh containing between several thousand triangles to almost one million triangles, computation times to optimize the point positions (100 iterations) and build the final mesh range between 13s and 150s on a 8-core laptop (in 2013).

Advantages Looking back, these results are quite good. The method is completely automatic and performs in several minutes or less on real data geological surface models. The method produces a result whatever the input, the potentially extremely bad quality of the input triangles and the nasty features that might have stopped another method are not an issue. These small features like small holes, gaps, thin features, whose size is smaller than the resolution of the final mesh are altered.

The final global resolution in terms of number of vertices is controlled by the user. The output quality of the triangles does it depend on the quality of the input mesh. The output triangle quality

Q = 6 S /( √ 3 h max p)
where S is the area of the triangle, h max the length of its longest edge, and p its half perimeter [START_REF] Frey | Surface mesh quality evaluation[END_REF] is good.

Drawbacks The method has however major drawbacks that makes it as it impracticable for remeshing of industrial models. First of all, no mesh size can be respected, that is one of the key job of a meshing software in engineering and despite some attempts Voronoi based methods are bad at it. Another major problem for engineers (less so for computer graphics applications) is the fact that there is absolutely no guarantee whatsoever on which of the input entities (corners, lines, points) will be preserved, the ones that are modified are not tracked and no one knows what happened or how to measure it. No guarantee is provided for the output model validity, when using an valid input surface model the output could be invalid. there is no guarantee either on the distance between the input and output models (link with mesh size). Some incomplete post-processing was implemented to fix some defects in the final meshes (1) to improve the quality and ensure conformity along contact lines when triangles are the dual of small triangular restricted Voronoi cell (2) to fix triangle intersections that can appear when making some simplifications. The real issue with this type of method is that there is the consistency between local modifications and the global model consistency is not adequately managed (Is it possible?) Locally modifications can be made that have a impact on a global scale, or at least on the neighboring local areas. The worst is probably that the modifications performed depend on the intersection of a Voronoi cell with the surface model. Move slightly that cell and the modifications may change significantly.

One click meshing is a bad dream. In general, performing topological modifications while keeping a surface model valid is extremely difficult, and this is without mentioning the remeshing part. This is why, in practice, most of the model modifications are checked by the user that can either use semi-automatic methods or even do it by hand (component suppression, mesh repair, etc.) My current point of view (in 2020) is that mesh generation and model modification/simplification/defeaturing must be separated for robustness puposes. Clicking on a button to get a mesh may be relevant and may become true for computer graphics applications. But I do not believe in this approach for engineering applications where validity, control, and quality must be guaranteed to the user (and not in theory). This is why in the work of Pierre Anquez we focused on the simplification of the model, and left the subsequent meshing to an adequate software.

Part II

Volumetric Meshing

Chapter 4

Tetrahedral Meshing à la Delaunay Motivations The Delaunay triangulation is a fundamental geometrical object that associates a unique triangulation DT (S ) to a given point set S in general position. This triangulation and its geometrical dual, the Voronoi diagram, have locality properties that makes them extremely useful to analyse points sets and of course to generate meshes. Generating robustly and efficiently the Delaunay triangulations is one key ingredient to many efficient robust and efficient tetrahedral meshing. It is important to remember that, for meshing purposes, the Delaunay triangulation is not a mesh but a mean to compute a mesh because its boundary is the convex hull of the points and not the boundary of the object to mesh.

Delaunay triangulation algorithms face a major challenge: the availability of more and more massive point sets. LiDAR or other photogrammetry technologies are used to survey the surface of entire cities and even countries, like the Netherlands [START_REF] Martinez-Rubi | Taming the beast: Free and open-source massive point cloud web visualization[END_REF]. The size of finite element mesh is also growing with the increased availability of massively parallel numerical solvers. It is now common to deal with meshes of over several hundred millions of tetrahedra [START_REF] Rasquin | Scalable fully implicit finite element flow solver with application to high-fidelity flow control simulations on a realistic wing design[END_REF][START_REF] Ibanez | PUMI: Parallel Unstructured Mesh Infrastructure[END_REF][START_REF] Vázquez | Alya: Multiphysics engineering simulation toward exascale[END_REF]).

Can we parallelize tetrahedral meshing? Parallelizing 3D Delaunay triangulations and Delaunaybased meshing algorithms is however very challenging because the algorithm is essentially sequential. The most part of today's clusters have two levels of parallelism. Distributed memory systems contain thousands of nodes, each node being itself a shared memory system with multiple cores. In recent years, nodes have seen their number of cores and the size of their memory increase, with some clusters already featuring 256-core processors and up to 12tb of RAM [START_REF] Nystrom | Bridges: A Uniquely Flexible HPC Resource for New Communities and Data Analytics[END_REF][START_REF] Fu | The Sunway TaihuLight supercomputer: system and applications[END_REF]. As many-core shared memory machines are becoming standard, Delaunay triangulation algorithms designed for shared memory should not only scale well up to 8 cores but to several hundred cores.

Contributions This chapter presents one of the main result of the PhD of Célestin Marot whom I co-advised during my post-doc in UCLouvain in Belgium.

Before parallelizing, we ensure that our sequential implementation of the Delaunay triangulation algorithm in 3D is efficient. It is about three times faster than existing open-source implementations and is able to triangulate a million points in about 2 seconds Section 4.1. Our implementation is similar to those, being based on the incremental Delaunay insertion algorithm, but the gain in performance is to ascribe to the details of the specific data structures we have developed, to the optimization of geometric predicates evaluation, and to specialized adjacency computations. We then propose a scalable parallel version of the Delaunay triangulation algorithm devoid of heavy synchronization overheads (Section 4.2). The domain is partitioned using the Hilbert curve and conflicts are detected with a simple coloring scheme. On a AMD EPYC 64-core machine, we have been able to generate three billion tetrahedra in less than a minute (about 10 7 points per second). We finally demonstrate how this efficient Delaunay triangulation algorithm can be easily integrated in a tetrahedral mesh generation process where the input is the boundary surfaces of the domain to mesh (Section 4.3).

A first version of this work was presented in a conference research note in 2017, the final results were published in 2018:

• Marot, C., J. Pellerin, J. Lambrechts, and J.-F. [START_REF] Remacle | A two-level multithreaded Delaunay kernel[END_REF]. "Toward one billion tetrahedra per minute". In: 26th International Meshing Roundtable, Research Notes. Barcelona, Spain

• Marot, C., J. Pellerin, and J. Remacle (2018). "One machine, one minute, three billion tetrahedra". en. In: International Journal for Numerical Methods in Engineering 117.9, pp. 967-990

The corresponding reference implementation is open-source and available in Gmsh 4 at http:// gmsh.info. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement ERC-2015-AdG-694020).

Sequential Delaunay

The fastest serial algorithm to build the 3D Delaunay triangulation DT (S ) is the Bowyer-Watson algorithm, which works by incremental insertion of points in the triangulation. The Bowyer-Watson algorithm was devised independently by [START_REF] Bowyer | Computing Dirichlet tessellations[END_REF]Watson in 1981 (Bowyer 1981;[START_REF] Watson | Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes[END_REF].

Very fast open-source implementations who benefited from various optimizations through the decades are available: Tetgen [START_REF] Si | TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator[END_REF], CGAL [START_REF] Boissonnat | Triangulations in CGAL[END_REF] and Geogram [START_REF] Levy | Geogram[END_REF]. They are designed similarly and offer therefore similar performances, they triangulate a million points in about 6 seconds on one core of a high-end laptop.

The work of Célestin was to first accelerate, if possible, the sequential implementation of the incremental Delaunay insertion algorithm in 3D. He succeeded, to be 3 times faster, this gain in performance is to ascribe to the details of the specific data structures we have developed, to the optimization of geometric predicates evaluation, and to specialized adjacency computations.

Bowyer-Watson Algorithm

Given a point set S , let DT k be the Delaunay triangulation of the subset S k = {p 1 , . . . , p k } ⊂ S . The Delaunay kernel is the procedure to insert a new point p k+1 ∈ Ω(S k ) into DT k , and construct a new valid Delaunay triangulation DT k+1 of S k+1 = {p 1 , . . . , p k , p k+1 }. The Delaunay kernel can be written in the following abstract manner: The state-of-the-art incremental insertion algorithm (Algorithm 1) that we implemented has five main steps that are described in the following.

DT k+1 ← DT k -C(DT k , p k+1 ) + B(DT k , p k+1 ), ( 4
Init The triangulation is initialized with the tetrahedron formed by the first four non-coplanar vertices of the point set S . These vertices define a tetrahedron τ with a positive volume.

Sort The complexity of the Delaunay triangulation algorithm depends on the number of tetrahedra to traverse during the Walk and on the number of tetrahedra in the cavities. To minimize the former, the points are sorted before insertion so that two points that have close indices are close in space. For very large point sets, the Sort function becomes the bottleneck, we implemented a very fast sorting procedure (Section 4.1.3) that eliminates this limitation. HXT [START_REF] Levy | Geogram[END_REF], Tetgen [START_REF] Si | TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator[END_REF] and CGAL [START_REF] Boissonnat | Triangulations in CGAL[END_REF]. Timings in seconds are given for 5 million points (random uniform distribution). The ≈ prefix indicates that no accurate timing is available.

τ ← t ∈ B 12: return DT +1 (a) Walk ( , +1 ) ( , +1 ) +1 (b) Cavity +1 ( , +1 ) ( , +1 ) (c) DelaunayBall
Walk The goal of this step is to identify the tetrahedron τ k+1 enclosing the next point to insert p k+1 .

The search starts from a tetrahedron τ k in the last Delaunay ball B(DT k , p k ), and walks through the current triangulation DT k in direction of p k+1 (Figure 4.1a). We say that a point is visible from a facet when the tetrahedron defined by this facet and this point has negative volume. The Walk function thus iterates on the four facets of τ, selects one from which the point p k+1 is visible, and then walks across this facet to the adjacent tetrahedron. This new tetrahedron is called τ and the Walk process is repeated until none of the facets of τ sees p k+1 , which is equivalent to say that p k+1 is inside τ (see Figure 4.1a).

The visibility walk algorithm is guaranteed to terminate for Delaunay triangulations [START_REF] De Floriani | On sorting triangles in a delaunay tessellation[END_REF]. If the points have been sorted, the number of walking steps is essentially constant (J.-F. Remacle 2017). Our implementation of this robust and efficient walking algorithm is similar to the existing ones.

Cavity Once the tetrahedron τ ← τ k+1 has been identified, the function Cavity finds all tetrahedra whose circumsphere contain p k+1 and deletes them. The Delaunay cavity C(DT k , p k+1 ) is simply connected and contains τ (J. Shewchuk 1997), it is then built using a breadth-first search algorithm. The core and most expensive operation of the Cavity function is the inSphere predicate, which evaluates whether a point e is inside/on or outside the circumsphere of given tetrahedron. This Cavity function is thus an expensive function of the incremental insertion, which accounts for about 33% of the total computation time (Table 4.1). To accelerate this, sub-components of the inSphere predicate can be precomputed Section 4.1.4. DelaunayBall Once the cavity has been carved, the DelaunayBall function first generates a set of new tetrahedra adjacent to the newly inserted point p k+1 and filling up the cavity, and then updates the mesh structure. In particular, the mesh update consists in the computation of adjacencies between the newly created tetrahedra. This is the most expensive step of the algorithm, with about 40% of the total computation time (Table 4.1). To accelerate this step, the general purpose elementary operations like tetrahedron creation/deletion or adjacency computation can be replaced by batches of optimized operations making benefit from a cavity-specific data structure (Section 4.1.5).

Data Structure Dedicated to Triangulation

One key for enhancing performances of a Delaunay triangulation algorithm (and maybe any algorithm) is the optimization of the data structure used to store the object of interest. Various data structure designs have been proposed that are very flexible and allow representing hybrid meshes, high order meshes, add mesh elements of any type, or manage adjacencies around vertices, edges, etc. The versatility of such general purpose data structures has a huge cost, both in terms of storage and efficiency. Here, our aim is to have a structure as lightweight and fast as possible, dealing exclusively with 3D triangulations. The implementation proposed in HXT is coded in plain C language, with arrays of doubles, floats, and integers to store both the mesh topology and geometry. This seemingly old-style coding has important advantages in terms of optimization and parallelization because it compels us to use simple and straightforward algorithms. The mesh data only contains vertices and tetrahedra explicitly, and all topological and geometrical information can be deduced from it. One very useful and relatively expensive information is the adjacency between tetrahedra.

Vertices Vertices are stored in a single array of structures point3d_t (see Listing 4.1). For each vertex, in addition to the vertex coordinates, a padding variable is used to align the structure to 32 bytes (3 doubles of 8 bytes each, and an additional padding variable of 8 bytes sum up to a structure of 32 bytes) and conveniently store vertex related information. Memory alignment ensures that a vertex does not overlap two cache lines during memory transfer. Modern computers usually work with cache lines of 64 bytes. The padding variable in the vertex structure ensures that a vertex is always loaded in one single memory fetch. Moreover, aligned memory allows to take advantage of the vectorization capabilities of modern microprocessors.

Tetrahedra Each tetrahedron knows its 4 vertices and its 4 neighboring tetrahedra. These two types of adjacencies are stored in separate arrays. The main motivation for this storage is flexibility and,

4t 0 a 4t 1 + 4t 0 + 1 b 4t 2 + 4t 0 + 2 c 4t 3 + 4t 0 + 3 d - : 4t 1 b - 4t 1 + 1 c - 4t 1 + 2 d - 4t 1 + 3 e 4t 0 + : 4t 2 a - 4t 2 + 1 d - 4t 2 + 2 c - 4t 2 + 3 f 4t 0 + : 4t 3 a - 4t 3 + 1 b - 4t 3 + 2 d - 4t 3 + 3 g 4t 0 + Figure 4
.2: Four adjacent tetrahedra : t 0 , t 1 , t 2 , t 3 and one of their possible memory representations in the tetrahedra data structure given in Listing 4.1. tetrahedra.neighbor_ID[4t i + j]/4 gives the index of the adjacent tetrahedron opposite to tetrahedra.vertex_ID[4t i + j] in the tetrahedron t i and tetrahedra.neighbor_ID[4t i + j] gives the index where the inverse adjacency is stored.

once more, memory alignment. Keeping good memory alignment properties on a tetrahedron structure evolving with the implementation is cumbersome. In addition, it provides little to no performance gain in this case. On the other hand with parallel arrays, additional information per tetrahedron (e.g. a color for each tetrahedron, sub-determinants etc.) can be added easily without disrupting memory layout. Each tetrahedron is identified by the indices of its four vertices in the vertices structure. Vertex indices of tetrahedron t are read between positions 4*t and 4*t+3 in the global array tetrahedra.vertex_ID storing all tetrahedron vertices. Vertices are ordered so that the volume of the tetrahedron is positive. An array of double, sub_determinant, is used to store 4 values per tetrahedron. This space is used to speed up geometric predicate evaluation (see Section 4.1.4).

Adjacencies By convention, the i-th facet of a tetrahedron is the facet opposite the i-th vertex, and the i-th neighbor is the tetrahedron adjacent to that facet. In order to travel efficiently through the triangulation, facet indices are stored together with the indices of the corresponding adjacent tetrahedron, thanks to an integer division and its modulo. The scheme is simple. Each adjacency is represented by the integer obtained by multiplying by four the index of the tetrahedron and adding the internal index of the facet in the tetrahedron. Take, for instance, two tetrahedra t 1 and t 2 sharing facet f , whose index is respectively i 1 in t 1 and i 2 in t 2 . The adjacency is then represented as

tetrahedra.neighbor_ID[4t 1 +i 1 ]=4t 2 +i 2 and tetrahedra.neighbor_ID[4t 2 +i 2 ]=4t 1 +i 1 .
This multiplexing avoids a costly looping over the facets of a tetrahedron. The fact that it reduces by a factor 4 the maximum number of elements in a mesh is not a real concern since, element indices and adjacencies being stored as 64 bytes unsigned integers, the maximal number of element in a mesh is 2 62 4.6 10 18 , which is huge. Note also that division and modulo by 4 are very cheap bitwise operations for unsigned integers: i/4 = i 2 and i%4 = i&3.

Memory footprint

The key to an efficient data structure is the balance between its memory footprint and the computational cost of its modification. We do not expect to generate meshes of more than UINT32_MAX vertices, i.e. about 4 billion vertices on one single machine. Each vertex therefore occupies 32 bytes, 24 bytes for its coordinates and 8 bytes for the padding variable. On the other hand, the number of tetrahedra could itself be larger than 4 billion, so that a 64 bits integer is needed for element indexing. Each tetrahedron occupies 80 bytes, 4 × 4 = 16 bytes for the vertices, 4 × 8 = 32 bytes for the neighbors, 32 bytes again for the sub-determinants. In average a tetrahedral mesh of n vertices has a little more than 6n tetrahedra. Thus, our mesh data structure requires ≈ 6 × 80 + 32 = 512 bytes per vertex.

Fast Point Spatial Sorting

An efficient spatial sorting scheme has been proposed by [START_REF] Boissonnat | Incremental Construction of the Delaunay Triangulation and the Delaunay Graph in Medium Dimension[END_REF] and is used in state-of-the-art Delaunay triangulation implementations. The main idea is to first shuffle the vertices, and to group them in rounds of increasing size (a first round with, e.g., the first 1000 vertices, a second with the next 7000 vertices, etc.) as described by the Biased Randomized Insertion Order (BRIO) [START_REF] Amenta | Incremental constructions con BRIO[END_REF]. Then, each group is ordered along a space-filling curve. The space-filling curve should have the property that successive points on the curve are geometrically close to each other. With this spatial sorting, the number of walking steps between two successive cavities remains small and essentially constant (J.-F. Remacle 2017). Additionally, proceeding by successive rounds according to the BRIO algorithm tends to reduce the average size of cavities.

Space-filling curve coordinates

The Hilbert curve is a continuous self-similar (fractal) curve. It has the interesting property to be space-filling, i.e., to visit exactly once all the cells of a regular grid with 2 m × 2 m × 2 m cells, m ∈ N. A Moore curve is a closed Hilbert curve. in the sense that points close to each other on the curve are also close to each other in R 3 [START_REF] Haverkort | Locality and bounding-box quality of two-dimensional space-filling curves[END_REF][START_REF] Abel | A Comparative Analysis of some 2-Dimensional Orderings[END_REF]). Any 3D point set can be ordered by a Hilbert/Moore curve according to the order in which the curve visits the grid cells that contains the points. The Hilbert/Moore index is thus an integer value d ∈ {0, 1, 2, . . . , 2 3m -1}, and several points might have the same index. The bigger m is for a given point set, the smaller the grid cells are, and hence the lower the probability of having two points with the same Hilbert/Moore index. Given a 3D point set with n points, there are in average n/2 3m points per grid cell. Therefore, choosing m = k log 2 (n) with k constant ensures the average number of points in grid cells to be independent of n. If known in advance, the minimum mesh size can also be taken into account: m can be chosen such that a grid cell never contain more than a certain number of vertices, thus capturing non-uniform meshes more precisely. For Delaunay triangulation purposes, the objective is both to limit the number of points with the same indices (in the same grid cell) and to have indices within the smallest range as possible to accelerate subsequent sorting operations. There is thus a balance to find, so that Hilbert indices are neither too dense nor too sparse.

Computing the Hilbert/Moore index of one cell is a somewhat technical point. The Hilbert/Moore curve is indeed fractal, which means recursively composed of replicas of itself that have been scaled down by a factor two, rotated and optionally reflected. Those reflections and rotations can be efficiently computed with bitwise operations. Various transformations can be applied to point coordinates to emulate non-regular grids, a useful functionality we will resort to in Section 4.2.2. [START_REF] Hamilton | Compact Hilbert indices: Space-filling curves for domains with unequal side lengths[END_REF] give extensive details on how to compute Hilbert indices.

Radix sort Once the Hilbert/Moore indices have been computed, points can be sorted accordingly in a data structure where the Hilbert/Moore index is the key, and the point the associated value. An extremely well suited strategy to sort bounded integer keys is the radix sort [START_REF] Blelloch | Design and Implementation of a Practical Parallel Delaunay Algorithm[END_REF][START_REF] Zagha | Radix sort for vector multiprocessors[END_REF][START_REF] Sohn | Load balanced parallel radix sort[END_REF][START_REF] Satish | Designing efficient sorting algorithms for manycore GPUs[END_REF]), a non-comparative sorting algorithm working digit by digit. The base of the digits, the radix, can be freely chosen. Radix sort has a O(wn) computational complexity, where n is the number of {key, value} pairs and w the number of digits (or bits) of the keys. In our case, the radix sort has a complexity in O(mn) = O(n log(n)), where n is the number of points and m = k log 2 (n) the number of levels of the Hilbert/Moore grid. In general, m is small because a good resolution of the space-filling curve is not needed. Typically, the maximum value among keys is lower than the number of values to be sorted. We say that keys are short. Radix-sort is able to sort such keys extremely quickly.

Literature is abundant on parallel radix sorting and impressing performances are obtained on manycore CPUs and GPus [START_REF] Wassenberg | Engineering a Multi-core Radix Sort[END_REF][START_REF] Polychroniou | A comprehensive study of main-memory partitioning and its application to large-scale comparison-and radix-sort[END_REF][START_REF] Satish | Fast sort on CPUs and GPUs: a case for bandwidth oblivious SIMD sort[END_REF][START_REF] Bell | Thrust: A productivity-oriented library for CUDA[END_REF][START_REF] Sengupta | Scan primitives for GPU computing[END_REF]). However, implementations are seldom available and we were not able to find a parallel radix sort implementation properly designed for many-core CPUs. We implemented HXTSort, that is available independently as open source at https://www.hextreme.eu/hxtsort. Figure 4.3 compares the performances of HXTSort with qsort, std::sort and the most efficient implementations that we are aware of for sorting 64-bit key and value pairs. Our implementation is fully multi-threaded and takes advantage of the vectorization possibilities offered by modern computers such as the AVX512 extensions on the Xeon PHI. It has been developed primarily for sorting Hilbert indices, which are typically short. We use different strategies depending on the key bit size. These are the reason why HXTSort outperforms the Boost Sort Library, GNU's and Intel's parallel implementation of the standard library and Intel TBB when sorting Hilbert indices.

Improving Cavity: Spending Less Time in Geometric Predicates

The central operation of the Delaunay kernel is the construction of the cavity C(DT k , p k+1 ) formed by all the tetrahedra {a, b, c, d} whose circumscribed sphere encloses p k+1 (Figure 4.1). This step of the incremental insertion represents about one third of the total execution time in available implementations (Table 4.1). The cavity is initiated with a first tetrahedron τ containing p k+1 determined with the Walk function (Algorithm 1 and Figure 4.1a), and then completed by visiting the neighboring tetrahedra with a breadth-first search algorithm.

Optimization of the inSphere predicate The most expensive operation of the Cavity function is the fundamental geometrical evaluation of whether a given point e is inside, exactly on, or outside the circumsphere of a given tetrahedra {a, b, c, d} (Table 4.1). This is evaluated using the inSphere predicate that computes the sign of the following determinant:

inSphere(a, b, c, d, e) = a x a y a z a 2 1 b x b y b z b 2 1 c x c y c z c 2 1 d x d y d z d 2 1 e x e y e z e 2 1 = b x -a x b y -a y b z -a z b -a 2 c x -a x c y -a y c z -a z c -a 2 d x -a x d y -a y d z -a z d -a 2 e x -a x e y -a y e z -a z e -a 2
This is a very time consuming computation, and to make it more efficient, we propose to expand the 4 × 4 determinant into a linear combination of four 3 × 3 determinants independent of point e.

inSphere(a, b, c, d, e) = -(e x -a x ) b y -a y b z -a z b -a 2 c y -a y c z -a z c -a 2 d y -a y d z -a z d -a 2 + (e y -a y ) b x -a x b z -a z b -a 2 c x -a x c z -a z c -a 2 d x -a x d z -a z d -a 2 -(e z -a z ) b x -a x b y -a y b -a 2 c x -a x c y -a y c -a 2 d x -a x d y -a y d -a 2 + e -a 2 b x -a x b y -a y b z -a z c x -a x c y -a y c z -a z d x -a x d y -a y d z -a z
Being completely determined by the tetrahedron vertex coordinates, the four 3 × 3 determinants can be pre-computed and stored in the tetrahedron data structure when it is created. The cost of the inSphere predicate becomes then negligible. Notice also that the fourth sub-determinant is minus the tetrahedron volume. We can set it to a positive value to flag deleted tetrahedra during the breadth-first search, thereby saving memory space.

Maximal improvement

The maximal improvement of the Cavity function obtained with this optimization depends on the number of times the inSphere predicate is invoked per tetrahedron. First, in order to simplify further discussion, we will assume that the number of tetrahedra is about 6 times the number of vertices in the final triangulation. This means that each point insertion results in average in the creation of 6 new tetrahedra. On the other hand, we have seen in Section 4.1.3 that an appropriate point ordering ensures an approximately constant number of tetrahedra in the cavities. This number is close to 20 in a usual mesh generation context (J.-F. Remacle 2017). One point insertion thus results in the deletion of 20 tetrahedra, and the creation of 26 tetrahedra (all figures are approximations). This number is also the number of triangular faces of the cavity, since all tetrahedra created by the DelaunayBall function associate a cavity facet to the inserted vertex p k+1 . The inSphere predicate is therefore evaluated positively for the 20 tetrahedra forming the cavity and negatively for the 26 tetrahedra adjacent to the faces of the cavity, a total of 46 calls for each vertex insertion.

When n points have been inserted in the mesh, a total of 26n tetrahedra were created and the predicate has been evaluated 46n times. Thus, we may conclude from this analysis that the inSphere predicate is called approximately 46n/26n = 1.77 times per tetrahedron. In consequence, the maximal improvement that can be obtained from our optimization of the inSphere predicate is of 1 -1/1.77 = 43%. Storing the sub-determinants has a memory cost (4 double values per tetrahedron) and a bandwidth cost (loading and storing of sub-determinants). For instance, for n = 4. 10 6 , we observe a speedup of 32% in the inSphere predicate evaluations, taking into account the time spent to compute the stored sub-determinants.

Note that a second geometric predicate is extensively used in Delaunay triangulation. The orient3d predicate evaluates whether a point d is above/on/under the plane defined by three points {a, b, c} ∈ R 3 (J. Shewchuk 1997). It computes the triple product (ba) • ((ca) × (da)), i.e. the signed volume of tetrahedron {a, b, c, d}. This predicate is mostly used in the Walk function.

Implementation details Geometrical predicates evaluated with standard floating point arithmetics may lead to inaccurate or inconsistent results. To have a robust and efficient implementations of the inSphere and orient3d predicates, we have applied the strategy implemented in Tetgen [START_REF] Si | TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator[END_REF].

1. Evaluate first the predicate with floating point precision. This gives a correct value in the vast majority of the cases, and represents only 1% of the code.

2. Use a filter to check whether the obtained result is certain. A static filter is first used, then, if more precision is needed, a dynamic filter is evaluated. If the result obtained by standard arithmetics is not trustworthy, the predicate is computed with exact arithmetics (J. Shewchuk 1997).

3. To be uniquely defined, Delaunay triangulations requires each point to be inside/outside a sphere, under/above a plane. When a point is exactly on a sphere or a plane, the point is in a non-general position that is slightly perturbed to "escape" the singularity. We implemented the symbolic perturbations proposed by [START_REF] Edelsbrunner | Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms[END_REF].

Improving DelaunayBall: Spending Less Time Computing Adjacencies

The DelaunayBall function creates the new tetrahedra filling the cavity (Figure 4.1c), and updates the tetrahedron structures. In particular, tetrahedron adjacencies are recomputed. This is the most expensive step of the Delaunay triangulation algorithm as it typically takes about 40% of the total time (Table 4.1). In contrast to existing implementations, we strongly interconnect the cavity building and cavity retriangulation steps. Instead of relying on a set of elegant and independent elementary operations like tetrahedron creation/deletion or adjacency computation, a specific cavity data structure has been developed, which optimizes batches of operations for the specific requirements of the Delaunay triangulation (Listing 4.2).

Use cavity building information

The tetrahedra reached by the breadth-first search during the Cavity step (Section 4.1.4), are either inside or adjacent to the cavity. Each triangular facet of the cavity boundary is thus shared by a tetrahedron t 1 ∈ A(DT k , p k+1 ) outside the cavity, by a tetrahedron t 2 ∈ C(DT k , p k+1 ) inside the cavity, and by a newly created tetrahedron inside the cavity t 3 ∈ B(DT k , p k+1 ) (Figure 4.1). The facet of t 1 adjacent to the surface of the cavity defines, along with the point p k+1 , the tetrahedron t 3 . We thus know a priori that t 3 is adjacent to t 1 , and store this information in the cavityBoundaryFacet_t structure (Listing 4.2).

Fast adjacencies between new tetrahedra During the cavity construction procedure, the list of vertices on the cavity as well as adjacencies with the tetrahedra outside the cavity are computed. The last step is to compute the adjacencies between the new tetrahedra built inside the cavity, i.e. the tetrahedra of the Delaunay ball.

The first vertex of all created tetrahedra is set to be p k+1 , whereas the other three vertices {p 1 , p 2 , p 3 } are on the cavity boundary, and are ordered so that the volume of the tetrahedral element is positive, i.e., orient3d(p k+1 , p 1 , p 2 , p 3 ) > 0. As explained in the previous section, the adjacency stored at index 4t i + 0, which corresponds to the facet of tetrahedron t i opposite to the vertex p k+1 , is already known for every tetrahedron t i in B(DT k , p k+1 ). Three neighbors are thus still to be determined for each tetrahedron t i , which means practically that an adjacency index has to be attributed to 4t i + 1, 4t i + 2 and 4t i + 3. The internal facets across which these adjacencies have to be identified are made of the common vertex p k+1 and one oriented edge of the cavity boundary. Using an elaborated hash table with complex collision handling would be overkill. We prefer to use a double entry lookup table of dimension n × n, whose rows and columns are associated with the n vertices {p j } of the cavity boundary,to which auxiliary indices 0 ≤ i j < n are affected for convenience and stored in the padding variable of the vertex, With this, the unique index of an oriented edge p 1 p 2 is set to be n × i 1 + i 2 , which corresponds to one position in the n × n lookup table. So, for each tetrahedron t i in the Delaunay ball with vertices {p k+1 , p 1 , p 2 , p 3 }, the adjacency index 4t i + 1 is stored at position n × i 2 + i 3 in the lookup table, and similarly 4t i + 2 at position n × i 3 + i 1 and 4t i + 3 at position n × i 1 + i 2 . A square array with a zero diagonal is built proceeding this way, in which the sought adjacencies are the pairs of symmetric components.

Theoretically, there is no maximal value for the number n of vertices in the cavity but, in practice, we can take advantage of the fact that it remains relatively small. Indeed, the cavity boundary is the triangulation of a topological sphere, i.e. a planar graph in which the number of edges is e = 3 f /2, where f is the number of faces, and whose Euler characteristic is ne + f = 2. Hence, the number of vertices is linked to the number of triangular faces by n = f /2 + 2. As we have shown earlier that f is about 26, we deduce that there are about n = 15 vertices on a cavity boundary. Hence, a n max × n max lookup table, with maximum size n max = 32 is sufficient provided there are at most f max = 2(n max -2) = 60 tetrahedra created in the cavity. If the cavity exceptionally contains more elements, the algorithm smoothly switches to a simple linear search.

Once all adjacencies of the new tetrahedra have been properly identified, the DelaunayBall function is then in charge of updating the mesh data structure. This ends the insertion of point p k+1 . The space freed by deleted tetrahedra is reused, if needed additional tetrahedra are added. Note that almost all steps of adjacencies recovery are vectorizable.

About a Ghost

The Bowyer-Watson algorithm for Delaunay triangulation assumes that all newly inserted vertices are inside an element of the triangulation at the previous step (Algorithm 1). To insert a point p k+1 outside the current support of the triangulation, one possible strategy is to enclose the input vertices in a sufficiently large bounding box, and to remove the tetrahedra lying outside the convex hull of the point set at the end of the algorithm. A more efficient strategy, adopted in TetGen [START_REF] Si | TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator[END_REF], CGAL [START_REF] Boissonnat | Triangulations in CGAL[END_REF], Geogram [START_REF] Levy | Geogram[END_REF], is to work with the so-called ghost tetrahedra connecting the exterior faces of the triangulation with a virtual ghost vertex G. Using this elegant concept, vertices can be inserted outside the current triangulation support with almost no algorithmic change.

The ghost vertex G is the vertex "at infinity" shared by all ghost tetrahedra. The ghost tetrahedra cover the whole space outside the regular triangulation. Like regular tetrahedra, ghost tetrahedra are stored in the mesh data structure, and are deleted whenever a vertex is inserted inside their circumsphere. The accurate definition of the circumsphere of a ghost tetrahedron, in particular with respect to the requirements of the Delaunay condition, is however a more delicate question.

From a geometrical point of view, as the ghost vertex G moves away towards infinity from an exterior facet abc of the triangulation, the circumscribed sphere of tetrahedron abcG tends to the half space on the positive side of the plane determined by the points abc, i.e., the set of points x such that orient3d(a, b, c, d) is positive. The question is whether this inequality should be strict or not, and the robust answer is neither one nor the other. The circumcircle of a ghost triangle abG actually contains not only the open half-plane strictly above the line ← → ab, but also the line segment [ab] itself; the other parts of the line ← → ab being excluded. In 3D, the circumsphere of a ghost tetrahedron abcG contains the half-space L such that for any point d ∈ L, orient3d(a, b, c, d) is strictly positive, plus the disk defined by the circumcircle of the triangle abc. If the point d is in the plane abc, i.e. orient3d(a, b, c, d)=0, we thus additionally test if d is in the circumsphere of the adjacent regular tetrahedron sharing the triangle abc. This composite definition of the circumscribed circle makes this approach robust with minimal changes: only the inSphere predicate is modified in the implementation.

Because the 3 × 3 determinant of orient3d is used instead of the 4 × 4 determinant of inSphere for ghost tetrahedra, the approach with a ghost vertex is faster than other strategies (J. R. Shewchuk 1997). Note that the Walk cannot evaluate orient3d on the faces of ghost tetrahedra that connect edges of the convex hull to the ghost vertex. If the starting tetrahedron is a ghost tetrahedron, the Walk directly steps into the non-ghost adjacent tetrahedron. Then, if it steps in an other ghost tetrahedron abcG while walking toward p k+1 , we have orient3d(a, b, c, p k+1 ) > 0 and the ghost tetrahedron is inside of the cavity. The walk hence stops there, and the algorithm proceeds as usual.

Serial Implementation Performances

HXT serial Delaunay triangulation algorithm has about one thousand lines (without Shewchuk's geometric predicates (J. Shewchuk 1997)). It is open-source and available in Gmsh (www.gmsh.info). Overall, it is almost three time faster than other sequential implementations. Figure 4.4a and 4.4b show the performance gap between our implementation and concurrent software on a laptop with a maximum core frequency of 3.5Ghz, and on a many-core computer with a maximum core frequency of 1.5Ghz and wider SIMD instructions. Table 4.1 indicates that the main difference in speed comes from the more efficient adjacencies computation in the DelaunayBall function. Other software use a more general function that can also handle cavities with multiple interior vertices. But this situation does not happen with Delaunay insertion, where there is always one unique point in the cavity, p k+1 . Since our DelaunayBall function is optimized for Delaunay cavities, it is approximately three time faster, despite the additional computation of sub-determinants. The remaining performance gain is explained by our choice of simple but efficient memory aligned data structure.

Parallel Delaunay Triangulation

The availability of more and more massive point sets as well as the size of finite element meshes is a major challenge for Delaunay triangulation algorithms. It is now common to deal with meshes of over several hundred millions of tetrahedra [START_REF] Rasquin | Scalable fully implicit finite element flow solver with application to high-fidelity flow control simulations on a realistic wing design[END_REF][START_REF] Ibanez | PUMI: Parallel Unstructured Mesh Infrastructure[END_REF][START_REF] Vázquez | Alya: Multiphysics engineering simulation toward exascale[END_REF]. workflow where they are adapted to the solution.

Related work As we have seen in the previous section the algorithm to build the Delaunay triangulation is essentially sequentially and parallelizing 3D Delaunay triangulations and Delaunay-based meshing algorithms is a very challenging task.

The most part of today's clusters have two levels of parallelism. Distributed memory systems contain thousands of nodes, each node being itself a shared memory system with multiple cores. In recent years, nodes have seen their number of cores and the size of their memory increase, with some clusters already featuring 256-core processors and up to 12tb of RAM [START_REF] Nystrom | Bridges: A Uniquely Flexible HPC Resource for New Communities and Data Analytics[END_REF][START_REF] Fu | The Sunway TaihuLight supercomputer: system and applications[END_REF]. As many-core shared memory machines are becoming standard, Delaunay triangulation algorithms designed for shared memory should not only scale well up to 8 cores but to several hundred cores. To overcome memory and time limitations, Delaunay triangulations should be constructed in parallel making the most of both distributed and shared memory architectures. A triangulation can be subdivided into multiple parts, each constructed and stored on a node of a distributed memory cluster. Multiple methods have been proposed to merge independently generated Delaunay triangulations [START_REF] Cignoni | Parallel 3D Delaunay Triangulation[END_REF][START_REF] Chen | The Merge Phase of Parallel Divide-and-Conquer Scheme for 3D Delaunay Triangulation[END_REF][START_REF] Funke | Parallel d -D Delaunay Triangulations in Shared and Distributed Memory[END_REF][START_REF] Blelloch | Design and Implementation of a Practical Parallel Delaunay Algorithm[END_REF]). However, those methods feature complicated merge steps, which are often difficult to parallelize. To avoid merge operations, other distributed implementations maintain a single valid Delaunay triangulation and use synchronization between processors whenever a conflict may occur at inter-processor boundaries [START_REF] Okusanya | 3D PARALLEL UNSTRUCTURED MESH GENERATION[END_REF][START_REF] Chrisochoides | Parallel Delaunay mesh generation kernel[END_REF]. In finite-element mesh generation, merging two triangulations can be simpler because triangulations are not required to be fully Delaunay, allowing algorithms to focus primarily on load balancing [START_REF] Lachat | Parallel mesh adaptation using parallel graph partitioning[END_REF].

On shared memory machines, divide-and-conquer approaches remain efficient, but other approaches have been proposed since communication costs between different threads are not prohibitive to the contrary of distributed memory machines. To insert a point in a Delaunay triangulation, the kernel procedure operates on a cavity that is modified to accommodate the inserted point (Figure 4.1). Two points can therefore be inserted concurrently in a Delaunay triangulation if their respective cavities do not intersect, C(DT k , p k1 ) ∩ C(DT k , p k2 ) = ∅, otherwise there is a conflict. In practice, other types of conflicts and data-races should possibly be taken into account depending on the chosen data structures and the insertion point strategy. Conflict management strategies relying heavily on locks1 lead to relatively good speedups on small numbers of cores [START_REF] Kohout | Parallel Delaunay triangulation in E2 and E3 for computers with shared memory[END_REF][START_REF] Blandford | Engineering a compact parallel delaunay algorithm in 3D[END_REF][START_REF] Batista | Parallel geometric algorithms for multi-core computers[END_REF][START_REF] Foteinos | Dynamic Parallel 3D Delaunay Triangulation[END_REF] However, synchronization overheads prevent those strategies from scaling to an high number of cores. More promising approaches rely on a partitioning strategy [START_REF] Lo | Parallel Delaunay triangulation in three dimensions[END_REF][START_REF] Loseille | Unique cavity-based operator and hierarchical domain partitioning for fast parallel generation of anisotropic meshes[END_REF]. Contrarily to pure divide-and-conquer strategies for distributed memory machines, partitions can change and move regularly, and they do not need to cover the whole mesh at once.

In HXT, we propose to parallelize the Delaunay kernel using partitions based on a space-filling curve, similarly to [START_REF] Loseille | Unique cavity-based operator and hierarchical domain partitioning for fast parallel generation of anisotropic meshes[END_REF]. The main difference is that we significantly modify the partitions at each iteration level. Our code is designed for shared memory architecture only, we leave its integration into a distributed implementation for future work.

A Parallel Strategy Based on Partitions

Our second contribution is a multi-threaded version of the Delaunay kernel that is able to concurrently insert vertices. Moore curve coordinates are used to partition the point set, avoiding heavy synchronization overheads. Conflicts are managed by modifying the partitions with a simple rescaling of the space-filling curve. The performances of our implementation have been measured on three different processors, an Intel core-i7, an Intel Xeon Phi and an AMD EPYC, on which we have been able to compute 3 billion tetrahedra in 53 seconds. This corresponds to a generation rate of over 55 million tetrahedra per second.

There are multiple conditions a program should ensure to avoid data-races and conflicts between threads when concurrently inserting points in a DT. Consider thread t 1 is inserting point p k1 and thread t 2 is simultaneously inserting point p k2 .

1. Thread t 1 cannot access information about any tetrahedron in C(DT , p k2 ) and inversely. Hence:

(a) C(DT , p k1 ) ∩ C(DT , p k2 ) = ∅ (b) A(DT , p k1 ) ∩ C(DT , p k2 ) = ∅ and A(DT , p k2 ) ∩ C(DT , p k1 ) = ∅ (c) Thread t 1 cannot walk into C(DT , p k2
) and reciprocally, t 2 cannot walk into C(DT , p k1 ) 2. A tetrahedron in B(DT , p k1 ) and a tetrahedron in B(DT , p k2 ) cannot be created at the same memory index.

To ensure rule (1) it is sufficient to restrain each thread to work on an independent partition of the mesh (Figure 4.5). This lock-free strategy minimizes synchronization between threads and is very efficient. Each point of the Delaunay triangulation is assigned a partition corresponding to a unique thread. A tetrahedron belongs to a thread if at least three of its vertices are in that thread's partition.

To ensure (1a) and (1b), the insertion of a point belonging to thread is aborted if the thread accessed a tetrahedron that belongs to another thread or that is in the buffer zone. To ensure (1c), we forbid threads to walk in tetrahedra belonging to another thread. Consequently, a thread aborts the insertion of a vertex when (i) the Walk reaches another partition, or when (ii) the Cavity function reaches a tetrahedron in the buffer zone. To insert these points, the vertices are re-partitioned differently (Section 4.2.2), a procedure repeated until there is no more vertices to insert or until the rate of successful insertions has become too low. In that case, the number of threads is decreased (Section 4.2.3). When the number of points to insert has become small enough, the insertion runs in sequential mode to insert all remaining vertices. The first BRIO round is also inserted sequentially to generate a first base mesh. Nevertheless, the vast majority of points are inserted in parallel.

Rule ( 2) is obvious from a parallel memory management point of view. It is however difficult to ensure it without requiring an unreasonable amount of memory. As explained in Section 4.2.5, synchronization between threads is required punctually.

Partitioning and Re-partitioning with Moore Curve

We subdivide the n v points to insert such that each thread inserts the same number of points. Our partitioning method is based on the Moore curve, i.e. on the point insertion order implemented for the sequential Delaunay triangulation (Section 4.1.3). Space-filling curves are relatively fast to construct and have already been used successfully for partitioning meshes [START_REF] Aluru | Parallel domain decomposition and load balancing using spacefilling curves[END_REF][START_REF] Lieber | FD4: A Framework for Highly Scalable Load Balancing and Coupling of Multiphase Models[END_REF][START_REF] Devine | New challenges in dynamic load balancing[END_REF][START_REF] Loseille | Unique cavity-based operator and hierarchical domain partitioning for fast parallel generation of anisotropic meshes[END_REF]. Each partition of the n thread partitions is a set of grid cells that are consecutive along the Moore curve. To compute the partitions, i.e. its starting and ending Moore indices, we sort the points to insert according to their Moore indices (see Section 4.1.3). Then, we assign the first n v /n threads points to the first partition, the next n v /n threads to the second, etc (Figure 4.6c). The second step is to partition the current Delaunay triangulation in which the points will be inserted. We use once more the Moore indices to assign the mesh vertices to the different partitions. The ghost vertex is assigned a random index. The partition owning a tetrahedron is determined from the partitions of its vertices, if a tetrahedron has at least three vertices in a partition it belong to this partition, otherwise the tetrahedron is in the buffer zone. Each thread then owns a subset of the points to insert and a subset of the current triangulation.

Once all threads attempted to insert all their points, a large majority of vertices is generally inserted. To insert the vertices for which insertion failed because the point cavity spans multiple partitions (Figure 4.5), we modify significantly the partitions by modifying the Moore indices computation. We apply a coordinate transformation and a circular shift to move the zero index around the looping curve (Figure 4.6). Coordinates below a random threshold are linearly compressed, while coordinates above the threshold are linearly expanded.

Ensuring Termination

When the number of vertices of the current Delaunay triangulation is small, typically at the first steps of the algorithm, the probability that the mesh vertices belong to different partitions is very high and none of the point insertions may succeed. To avoid wasting precious milliseconds, the first BRIO round is always inserted sequentially.

Moreover, there is no guarantee that re-partitioning will be sufficient to insert all the points. And even when a large part of the triangulation has already been constructed, parallel insertion may enter an infinite failure loop. After a few rounds, the remaining points to insert are typically located in restricted volumes (intersection of previous buffer zones). Because re-partitioning is done on the not-yet-inserted points, the resulting partitions are also small and thin. This leads to inevitable conflicts as partitions get smaller than cavities. In practice, we observed that the ratio ρ of successful insertions decreases for a constant number of threads. If 80 out of 100 vertices are successfully inserted in a mesh (ρ k = 0.8), less that 16 of the 20 remaining vertices will be inserted at the following attempt (ρ k+1 < 0.8). Note that this effect is more important for small meshes, because the bigger the mesh, the relatively smaller the buffer zone and the higher the insertion success rate. This difference explains the growth of the number of tetrahedra created per second with the number of points in the triangulation (Figure 4.8).

To avoid losing a significant amount of time in re-partitioning and trying to insert the same vertex multiple times we gradually decrease the number of threads. Choosing the adequate number of threads is a question of balance between the potential parallelization gain and the partitioning cost that comes with each insertion attempt. When decreasing the number of threads, we decrease the number of Numbers of threads used to insert points in our parallel Delaunay triangulation implementation according to the number of points to insert, the mesh size and the insertion success at the previous step. 94.5% of points are inserted using 8 threads and 5% using 4 threads attempts needed. When the ratio of successful insertions is too low, ρ < 1/5, or when the number of points to insert per thread is under 3000, we divide the number of threads by two. Furthermore, if ρ k ≤ 1/n threads , the next insertion attempt will not benefit from multi-threading and we insert the points sequentially. In Table 4.2 are given the number of threads used at each step of the Delaunay triangulation of a million vertices depending on the number of points to insert, the size of the current mesh, and the success insertion ratio ρ. Note that the computation of Moore indices and the sorting of vertices to insert are always computed on the maximal number of threads available (8 threads in this case) even when the insertion runs sequentially.

Data Structures

The data structure for our parallel Delaunay triangulation algorithm is similar to the one used by our sequential implementation (Section 4.1.2). There are two small differences. First, the parallel implementation does not compute sub-determinants for each tetrahedron. Actually, the bandwidth usage with the parallel insertions is already near its maximum. Loading and storing sub-determinant becomes detrimental to the overall performance. Instead we store a 16-bit2 color flag to mark deleted tetrahedra. Second, each thread has its own Cavity_t structure (Listing 4.2) to which are added two integers identifying the starting and ending Moore indices of the thread's partition.

Memory footprint Because we do not store four sub-determinants per tetrahedra anymore but only a 2-bytes color, our mesh data structure is lighter. Still assuming that there is approximately 6n tetrahedra for n vertices, it requires a little more than 6 × 50 + 32 = 332 bytes per vertex.

Critical Operations

When creating new tetrahedra in the Delaunay ball of a point, a thread first recycles unused memory space by replacing deleted tetrahedra. The indices of deleted tetrahedra are stored in the cavity->deleted.tetrahedra_ID4.1 array of each thread. When the cavity->deleted.tetrahedra_ID array of a thread is empty, additional memory should be reserved by this thread to create new tetrahedra. This operation is a critical part of the program requiring synchronization between threads to respect the rule (2). We need to capture the current number of tetrahedra and increment it by the requested number of new tetrahedra in one single atomic operation. OpenMP provides the adequate mechanism, see Listing 4.3.

To reduce the number of time this operation is done, the number of tetrahedra is incremented atomically by at least 8192, and the cavity->deleted.tetrahedra_ID is filled with the new indices of tetrahedra. Those tetrahedra are conceptually deleted although they have never been in the mesh. The number 8192 was found experimentally among multiples of 5123 . Increasing it would reduce the number of time reservation of new tetrahedra is performed but it is not necessary. Indeed, since the critical operation occurs then in average every 1000 + insertions, the time wasted is very low. Therefore, increasing the default number of deleted tetrahedra would only wastes memory space for little to no gain.

When the space used by tetrahedra exceeds the initially allocated capacity, reallocation is implemented that doubles the capacity of the arrays of mesh tetrahedra. During that operation, memory is potentially moved at another location. No other operation can be performed at that time on the mesh. Therefore, the reallocation code is placed in between two OpenMP barriers (Listing 4.4). This synchronization event is very rare and does not impact performances. In general, a good estimation of the needed memory needed to reserve is possible and that critical section is never reached.

Parallel Implementation Performances

We are able to compute the Delaunay triangulation of over one billion tetrahedra in record-breaking time: 41.6 seconds on the Intel Xeon Phi and 17.8 seconds on the AMD EPYC. These timings do not include I/Os. We are able to generate three billion tetrahedra in 53 seconds on the EPYC. The scaling of our implementation regarding the number of threads is detailed in

Tetrahedral Mesh Generation

We finally show how this very efficient parallel Delaunay triangulation can be integrated in a Delaunay refinement mesh generator which takes as input the triangulated surface boundary of the volume to mesh.

Taking into account the boundary The parallel Delaunay triangulation algorithm that we presented in the previous section is integrated in a Delaunay refinement mesh generator. A tetrahedral mesh generator is a procedure that takes as input the boundary of a domain to mesh, defined by set of triangles t that defines the boundary of a closed volume, and that returns a finite element mesh, i.e. a set of tetrahedra T of controlled sizes and shapes which boundary is equal to the input triangulation: ∂T = t. From a coarse mesh that is conformal to t, a Delaunay-refinement based mesher inserts progressively vertices until element sizes and shapes follow the prescribed ranges. Generating a mesh is an intrinsically more difficult problem than constructing the Delaunay triangulation of given points because (1) the points are not given, (2) the boundary triangles must be facets of the generated tetrahedra, and (3) the tetrahedra shape and size should be optimized to maximize the mesh quality. Our objective in this section is to demonstrate that our parallel Delaunay point insertion may be integrated in a mesh generator. The interested reader is referred to the book by Frey and George (Frey and George 2007) for a complete review of finite-element mesh generation.

The mesh generation algorithm 2 proposed in this section follows the approach implemented for example in Tetgen [START_REF] Si | TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator[END_REF]. First, all the vertices of the boundary triangulation t are tetrahedralized to form the initial "empty mesh" T 0 . Then, T 0 is transformed into a conformal mesh T which boundary ∂T is equal to t : ∂T = t. The triangulation T is then progressively refined by (i) creating vertices S at the circumcenters of tetrahedra for which the circumsphere radius, r τ , is significantly larger than the desired mesh size, h, i.e. r τ /h > 1.4 (J. R. Shewchuk 1997). (ii) inserting the vertices in the mesh using our parallel Delaunay algorithm. The sets of points S to insert are filtered a priori in order not to generate short edges, i.e. edges of size smaller than 0.7h. We first use Hilbert coordinates to discard very close points on the curve, and implemented a slightly modified cavity algorithm that aborts if a point of the cavity is too close from the one to insert (Section 4.1.4). Points are thus discarded both in the filtering process and in the insertion process. Note that contrary to existing implementations, we insert as large as possible point sets S during the refinement step to take advantage of the efficiency of our parallel point insertion algorithm (Section 4.2).

We parallelized all steps of the mesh generator (Algorithm 2) except the boundary recovery procedure that is the one of Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in preand post-processing facilities[END_REF]) which is essentially based on Tetgen [START_REF] Si | TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator[END_REF]. The cavity algorithm has also been modified to accommodate possible constraints on the faces and edges of tetrahedra. With this modification, mesh refinement never breaks the surface mesh and constrained edges and faces can be included in the mesh interior. As a result, the mesh may not satisfy the Delaunay property anymore. Therefore, we must always ensure that cavities are star-shaped. This is achieved by a simple procedure that checks if boundary facets of the cavity are oriented such that they form a positive volume with the new point. If a boundary facet is not oriented correctly, indicating that the cavity is not star-shaped, the tetrahedron containing it is removed from the cavity and the procedure is repeated. In practice, these modifications do not affect the speed of point insertions significantly.

To generate high-quality meshes, an optimization step should be added to the algorithm to improve the general quality of the mesh elements and remove sliver tetrahedra. Mesh improvements are out of the scope of this part. However, first experiments have shown that optimization step should slow down the mesh generation by a factor two. For example, mesh refinement and improvement take approximately the same time in Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in preand post-processing facilities[END_REF]).

Small and Medium Size Test Cases on Standard Laptop

In order to verify the scalability of the whole meshing process, meshes of up to one hundred million tetrahedra were computed on a 4 core 3.5 GHz Intel Core i7-6700HQ with 1, 2, 4 and 8 threads. Those meshes easily fit within the 8Gb of RAM of this laptop.

Three benchmarks are considered in this section: (i) a cube filled with cylindrical fibers of random radii and lengths that are randomly oriented, (ii) a mechanical part and (iii) a truck tire. Surface meshes T ← RecoverBoundary(T 0 ) [START_REF] Si | TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator[END_REF] 4:

while T contains large tetrahedra do 5:

S ← SamplePoints(T )

6:

S ← FilterPoints(S ) 7:

T ← InsertPoints(T , S ) Modified Delaunay kernel 8:

return T are computed with Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in preand post-processing facilities[END_REF]. Mesh size on the surfaces is controlled by surface curvatures and mesh size inside the domain is simply interpolated from the surface mesh.

Illustrations of the meshes, as well as timings statistics are presented in Figure 4.10. Our mesher is able to generate between 40 and 100 million tetrahedra per minute. Using multiple threads allows some speedup, the mesh refinement process is accelerated by a factor ranging between 2 and 3 on this 4 core machine.

The last test case (truck tire) is defined by more than 7000 CAD surfaces. Recovering the 27 892 triangular facets missing from T 0 takes more than a third of the total meshing time with the maximal number of threads. Parallelizing the boundary recovery process is clearly a priority of our future developments. On this same example, the surface mesh was done with Gmsh using four threads. The surface mesher of Gmsh is not very fast and it took about the same time to generate the surface mesh of 6 881 921 triangles as to generate the volume mesh that contains over one hundred million tetrahedra using the same number of threads. The overall meshing time for the truck tire test case is thus about 6 minutes.

Large Mesh Generation on Many Core Machine

We further generated meshes containing over 300,000,000 elements on a AMD EPYC 64 core machine. Three benchmarks are considered: (i) two cubes filled with many randomly oriented cylindrical fibers of random radii and lengths, and (ii) the exterior of an aircraft. Surface meshes were also generated by Gmsh.

Our strategy reaches its maximum efficiency for large meshes. In the 500 thin fibers test case, over 700,000,000 tetrahedra were generated in 135 seconds. This represents a rate of 5.2 million tetrahedra per second. In the 500 thin fibers test case, boundary recovery cost was lower and a rate of 6.2 million tetrahedra per second was reached.

Remaining Challenges We optimized both the sequential and the parallel algorithm by specifying them exclusively for Delaunay triangulation purposes. Our parallel implementation of 3D Delaunay triangulation construction has one important drawback: partitioning is more costly for non-uniform point sets because Moore indices are not adaptive. Adaptive Hilbert or Moore curve computation is relatively simple [START_REF] Hamilton | Compact Hilbert indices: Space-filling curves for domains with unequal side lengths[END_REF][START_REF] Su | Rapid Delaunay triangulation for randomly distributed point cloud data using adaptive Hilbert curve[END_REF]) but may not be faster than our refinement and sorting in one batch approach. Another drawback of our implementation is that it might not be well suited for distributed memory architectures. However, nothing prevents our strategy from being included in any larger divide-and-conquer approach.

We additionally presented a lightweight mesh generator based on Delaunay refinement. This mesh generator is almost entirely parallel and is able to produce computational meshes with a rate above 5 million tetrahedra per second. One issue when parallelizing mesh generation is to ensure reproducibility. For mesh generation, it is usually admitted that mesh generation should be reproducible for a given number of threads. Our mesh generator has been made reproducible in that sense, with a loss of 20% in overall performance. This feature has been implemented as an option. On the other hand, Figures 4.10 and 4.11 show that, starting from the same input, the number of elements varies.

Making the code reproducible independently of the number of threads would dramatically harm its scalability.

Note that our mesh generator does not include the final optimization step of the mesh, mesh improvement, to obtain a mesh adequate for finite-element simulations. This crucial step of any industrial-grade meshing procedure is the major limitation of our mesh generator. The basic operations to perform: edge swap, edge collapse, edge split, vertex relocation and face swap operations will be done using the same parallel framework than the one we described to build the Delaunay triangulation. The final challenge is the parallelization of the boundary recovery procedure. A third limitation is known in the placement of new vertices. Our current strategy, placing new vertices at the circumcenter of existing tetrahedra, does not always offer a sufficient quality when compared with frontal approaches. The frontal point insertion strategy that has been proposed in [START_REF] Baudouin | A frontal approach to hex-dominant mesh generation[END_REF] should also be parallelizable using our partitioning strategy. Chapter 5

Hybrid Meshing à la Voronoi

When generating 3D meshes constrained by input surfaces, the shapes and sizes of the generated tetrahedra depend on the model geometry and on the input surface mesh quality. The temptation is great to remesh the surfaces and volumes simultaneously. The temptation to do this while simplifying the small features of the surface model is even greater. During my PhD thesis, unaware of the inherent extreme difficulties of mesh generation in 3D, we experimented the construction of a hybrid mesh by extending the work of surface remeshing presented in Section 3. model, restricted Voronoi cells to regions will intersect several times the model surfaces (Figure 3.9). The connected components sandwiched between two surfaces may then not contain any Voronoi vertex (for example cells A, B, and C on Figure 5.2 share only a restricted Voronoi edge e). By putting one vertex at the centroid of each restricted Voronoi cell to the surfaces, a prism corresponding to e can be built.

Putting in correspondance the cells, facets, edges, and vertices of the Voronoi diagram cut into several components by the model surfaces with the cells, facets, edges, and vertices of a hybrid mesh is a hazardous task. We proposed however to do so. One or two vertices of the final mesh are associated to each volumetric restricted Voronoi cell depending on the number of adjacent surface connected components.

Building the Mesh Cells

To build all the cells automatically whatever the configuration of their intersections with the model regions and surfaces, we propose an algorithm that consists of the following steps.

1. Identify the cell from the corresponding restricted Voronoi edges and vertices.

Build the cells one after another by adding successively their vertices, edges, and facets

The procedure is illustrated for a restricted Voronoi vertex to a thin section sandwiched between two boundary surfaces (Figure 5.3a).

Vertices For each of the four volumetric restricted Voronoi cell containing the restricted Voronoi vertex, there is one or two points to add to the cell. If the restricted cell does not intersect model surfaces, that point is the site of the cell. If it does, the centroids of these intersections (adjacent connected component of the restricted cell to the surfaces) are taken for the points. When there are more than two intersections, some points are merged at their centroid so that there is maximum two points per restricted cell. For example, on Figure 5.3a, each restricted cell that contains vertex v is adjacent to two surface connected components and corresponds to two points (Figure 5.3b). Edges There are two types of edges to add: those linking the points corresponding to the same Voronoi cell (Figure 5.3b), and those corresponding to Voronoi facets containing v (Figure 5.3c). Each of the facets 1, 2, 3, 5 corresponds to two edges; the small triangular facets 4 and 6 intersect once the model surfaces and correspond to a unique edge in the cell to build.

Facets The first facets to add are those corresponding to the four restricted edges containing v. They are built from the edges corresponding to the restricted facets containing the restricted edge. In our example, they are all triangles (Figure 5.3d), but in the general case they could be polygons with up to six vertices. The other facets to add correspond to the restricted Voronoi facets containing v that are adjacent to at least two surface connected components. The four gray quad facets on Figure 5.3d correspond to restricted facets 1, 2, 3 and 5 on Figure 5.3a.

Finalizing the current cell Once all vertices, edges, and facets have been added to the cell, the cell type is determined from the number of these elements. Theoretically, a cell created with the above algorithm may have four to eight vertices, six to sixteen edges, four to ten facets, each facet having up to six vertices. In practice, most of the cells are tetrahedra, prisms, or pyramids. The other cells are processed as follows. They are first sorted into three categories: (1) the cells to subdivide (cells with 7 or 8 vertices defining a valid volume, Figure 5.3d), (2) invalid cells (in which at least one edge is a diagonal of a facet, or in which two quad facets share more than two vertices, like 5-point cells which are relatively abundant, see Figure 5.4a) and (3) the cells that belong to none of the above categories. Cells to subdivide are subdivided by adding a vertex at their centroid and building pyramids and tetrahedra with their facets. To compute tetrahedra, pyramids, and prisms from invalid cells, a two step procedure is implemented. First, all facets that are inconsistent with some edges are split along these edges. Then the validity of cell facets must be ensured. Facets with more than four vertices are triangulated and quad facets that share three vertices with another facet are split (Figure 5.4b). This processing permits to reduce the number of invalid cells, but is not sufficient to obtain a final valid mesh in all cases.

Could it Work?

The results of the algorithm are quite impressive for simple models in which layer thicknesses do not vary significantly and that have smooth boundaries, such as three nested spheres. It is not only the number of points sampling the spheres that determines which layers are thin and meshed with prisms, but also the positions of these points and the way they were optimized (Figure 5.5). The heptoroid model is meshed from only 1010 sampling points in 11 seconds (Figure 5.7). However the model is pinched near its extremities because only one site samples both creases and each restricted Voronoi cell to the surface has one connected component. The surface and volume are remeshed simultaneously. The final mesh does not depend on the input triangulation of the model surfaces (the only requirement is that their meshes are conformal) and no preliminary identification of the thin sections, that are meshed with prisms, is required.

Results for synthetic geological models are less convincing (Figure 5.6). The thinnest layer of the simple fold model on the top of Figure 5.6 is completely meshed with prisms. In the second model, built from the first by adding two faults that cut and displace the layers (one of them not cutting through the model), the same thin layer is also meshed with prisms. The discontinuities induced by the faults are respected in the mesh. However the layer is pinched near the faults and near the model box, an approximation that may be undesirable. The third model is built from the second by adding an erosion surface, below which remaining layers can be very thin and which can intersect other surfaces tangentially. Computational time (between several seconds and several minutes on a laptop with a 8 cores 1.73GHz Inter Core i7 processor and 8 GB of RAM) are quite good, it depends on the number of sampling points, on the input model number of triangles, and on the number of invalid cells to process.

The method has the same drawbacks than the surface remeshing method on which it is based Section 3.2.3. The outpush mesh is not valid because some output cell do not define a valid volume, and some other cells may intersect neighboring cells. Like for octree meshing methods, the difficulty is that all possible configurations for the intersections between subdivision cells and the model can and will occur. Having a control on these intersections is one of the reasons why the points are optimized so that their (restricted) Voronoi diagram is centroidal, and why many (re)meshing methods based on such an optimization refine the model sampling so that cell intersections with the model have a unique connected component (and verify the topological ball property [START_REF] Edelsbrunner | Triangulating topological spaces[END_REF]).

Missing pieces to make it work However, we did not think (then) that the fact that the algorithm fails for some cells is a show stopper, because, before a fully satisfactory solution is found, whenever the algorithm does not successfully generate a polyhedral cell, it is still possible to generate tetrahedra in cavities corresponding to the cells where the algorithm failed. The second possibility is to completely change the cell building strategy, and instead of building a mesh from scratch, progressively modify the Delaunay triangulation of the sampling points by local valid operations (vertex displacement, vertex/edge/facet duplication) that would permit to obtain the desired conformal mixed element mesh. For instance, as shown on Figure 5.4, whenever a restricted Voronoi cell intersects the set of surfaces twice, it may generate a cell with 5 vertices, due to the considered vertex that will be duplicated. In addition, note that when duplicating a vertex and moving the two instances at different locations (i.e., at the centroid of the restricted Voronoi cell connected components), special care is needed to avoid generating inverted elements (i.e. with negative orientation/Jacobian).

Food for thought Mesh generation in 3D is hard, because geometry in 3D is hard. Hoping to easily extend a meshing method from 2D to 3D is a dream that never comes true. The third dimension means trouble, and we often cannot imagine how bad it can get. So, if 2D is not a piece of cake, do not even try to think about 3D.

Chapter 6

Hybrid Meshing: Indirect Method Motivations Hexahedral meshes are considered by many of the finite element practitioners to be superior to tetrahedral meshes. Compared to tetrahedral meshes, hexahedral meshes have important numerical properties: faster assembly (J.-F. [START_REF] Remacle | GPU accelerated spectral finite elements on all-hex meshes[END_REF], high accuracy in solid mechanics (E. [START_REF] Wang | Back to elements-tetrahedra vs. hexahedra[END_REF]), or for quasi-incompressible materials [START_REF] Benzley | A comparison of all hexagonal and all tetrahedral finite element meshes for elastic and elasto-plastic analysis[END_REF]). Yet, no robust meshing technique is able to process general 3D domains. And generating hexahedral meshes in an automatic manner is still considered as the biggest challenge in mesh generation [START_REF] Shepherd | Hexahedral mesh generation constraints[END_REF]. The generation of ustructured hexahedral meshing is an active topic, e.g. [START_REF] Gao | Robust hex-dominant mesh generation using field-guided polyhedral agglomeration[END_REF][START_REF] Lyon | HexEx: robust hexahedral mesh extraction[END_REF][START_REF] Solomon | Boundary Element Octahedral Fields in Volumes[END_REF][START_REF] Fang | All-hex meshing using closed-form induced polycube[END_REF]) but some of the underlying theoretical challenges remain to explore as we will see in Chapter 7.

An alternative to full-hexahedral meshes are hybrid meshes that contains hexahedra (ideally a majority), prisms, pyramids and tetrahedra. To build those meshes promising, several author proposed to start from a tetrahedral mesh and merge its tetrahedra to build the output mesh [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF]Yamakawa et al. 2003;Levy and Y. Liu 2010b;[START_REF] Baudouin | A frontal approach to hex-dominant mesh generation[END_REF][START_REF] Huang | Boundary aligned smooth 3D cross-frame field[END_REF][START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF][START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF]. These methods take advantage of the existence of robust algorithms to generate tetrahedral meshes (see Chapter 4). Combining tetrahedra permit to produce meshes composed of hexahedra associated to prisms, pyramids and tetrahedra. The four steps of these indirect hex-dominant meshing methods can be summarized as follows (see also Figure 6.1):

1. A set of mesh vertices V is initially sampled in the domain.

2.

A tetrahedral mesh T is built by connecting V, e.g. using a Delaunay kernel like [START_REF] Si | TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator[END_REF].

3. The set of potential cells H (hexahedra, prisms, pyramids) that can be defined by combining tetrahedra of T is built.

4. A maximal subset H c ⊂ H constituted of cells that can be part of the same final mesh is determined. The final hex-dominant mesh is obtained adding the remaining not selected tetrahedra T .

Contribution The ultimate goal is to combine all tetrahedra into hexahedra, i.e. obtain a final full-hexahedral mesh. All steps of indirect method are then crucial. Previous works primarily focus on the first step of placing the final mesh vertices. In this chapter, we focus on the third step. Our input is a tetrahedral mesh T of a given point set V Our contribution is an efficient algorithm that performs the identification of combination of tetrahedra into hexahedra, prisms and pyramids, and builds the set H of all possible combinations of tetrahedral elements of an input mesh T . All identified cells are valid for engineering analysis. Each potential hexahedron/prism/pyramid is computed only once. Around 3 millions potential hexahedra are computed in 10 seconds on a laptop. The method is about 10 times faster than previous methods while identifying all possible combinations to the contrary of method relying on a predefined set of patterns of the decomposition of a hexahedron into tetrahedra [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF][START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF][START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF], or those relying on patterns of edge connections in a hexahedron (Yamakawa et al. 2003;[START_REF] Baudouin | A frontal approach to hex-dominant mesh generation[END_REF]).

(2) Tetrahedral mesh (3) All possible combinations of tets in hexahedra, prisms, pyramids (4) Hex-dominant mesh (1) Point set The C++ code implementing the methods is open-source and available at https://www.hextreme.eu/download/.

State of the Art

Before giving details on subdivisions of hexahedra , prisms or pyramids into tetrahedra (Section 6.1.1) and on methods used to identify these in an existing mesh (Section 6.1.2), we define the terms and notations used throughout this chapter. Definitions We have to be very clear that all the cells we are considering are finite element cells hexahedra /prisms/pyramids valid for finite element simulations. A very important point is that their quadrilateral facets are not planar, but are bilinear surfaces. We require the Jacobian determinant to be strictly positive at any point inside the cell. The following conventions will be used throughout this chapter (see Figure 6.3).

A triangulation of a hexahedron/prism/pyramid is a triangulation of the vertices of the cell that respect the cell boundary, in other words it is a subdivision of the interior of the hexahedron/prism/pyramid into a set of conformal tetrahedra without any additional vertex (Figure 6.2). The tetrahedra induce a subdivision of each quadrilateral facet into two triangles by a diagonal boundary edge. We further define the boundary tetrahedra as the tetrahedra connecting the four vertices of a cell quadrilateral facet (Figure 6.4.2). In previous works [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF][START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF][START_REF] Baudouin | A frontal approach to hex-dominant mesh generation[END_REF][START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF], boundary tetrahedra are called slivers. We do not use that term which refers to a geometrical property (degeneracy) of tetrahedra. The triangulation is determined by interior tetrahedra (Figure 6.4.1). Indeed, the addition or removal of one or several boundary tetrahedra does not modify the cell.

Decomposing Hexahedra into Tetrahedra

There are as many subdivisions of the general pyramid as there are subdivisions of its planar base [START_REF] De Loera | Triangulations: structures for algorithms and applications[END_REF]. We are considering square pyramids, there are then exactly two triangulations of the pyramid. The ordinary triangular prism is the result of the product of a triangle with an edge: prism(D 3 ) = D 3 × D 2 . A prism(D n ) has exactly n! triangulations that are all equivalent to one another by affine symmetries [START_REF] De Loera | Triangulations: structures for algorithms and applications[END_REF]. The prism has then 6 triangulations that are all equivalent.

74 triangulations of the 3-cube I 3 = [0, 1] 3 (De Loera et al. 2010) :

1. Every triangulation of the 3-cube contains either a regular tetrahedron (i.e. a tetrahedron whose 6 edges are of equal lengths) or a diameter, i.e. an interior edge joining two opposite vertices (red edges on Figure 6.5). 2. There are 2 triangulations with a regular tetrahedron, symmetric to one another. The triangulations containing an interior edge are completely classified modulo symmetries by their dual complex which can be one of the last five shown on Figure 6.5. There are respectively 8, 24, 12, 24, 4 triangulations in each class.

A dual complex (Figure 6.5) is a practical way to visualize the 6 different possible decompositions (tetrahedrizations) of the 3-cube. In the dual complex, one vertex corresponds to one tetrahedron and two vertices are connected by an edge if the corresponding tetrahedra are adjacent through a triangular facet. A 2-cell of the dual complex (cycle in the dual graph) corresponds to an interior edge of the tetrahedrization (red on Figure 6.5). In a meshing context, these different possible decompositions of the 3-cube were identified by [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF] who enumerate the feasible dual complex graphs, called RF-graph in their paper. In the RF-graph, additional dashed edges connect tetrahedra that are adjacent to the same quadrilateral facet (Figure 6.5).

Triangulations of the real cube

Recently, the work of [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF] was extended by [START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF] decomposition patterns into seven tetrahedra (Figure 6.6). The hexahedron is split into two prisms by a tetrahedron without any facet on the hexahedron boundary and containing two interior edges. For this tetrahedron to have a strictly positive volume, it is sufficient to work in finite precision, i.e. move slightly one of its vertices.

Combining Tetrahedra into Hexahedra

Two main approaches have been proposed. The first relies on a predefined set of patterns of the decomposition of a hexahedron into tetrahedra [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF][START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF][START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF], the second on patterns of edge connections in a hexahedron (Yamakawa et al. 2003;[START_REF] Baudouin | A frontal approach to hex-dominant mesh generation[END_REF]). Both strategies do not build the largest set of potential hexahedra H, and would miss the hexahedra on Figure 6.7.

To compute the set H of potential hexahedra and other cells that may be built by combining the elements of a tetrahedral mesh T without modifying its connectivity there are two known approaches. [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF] propose to find combinations of tetrahedra into hexahedra by searching the adjacency graph of T for all occurrences of the cube decomposition dual complexes (Figure 6.5). The problem of matching subgraphs in large sparse graphs is solved using standard data mining algorithms that operate on graphs. The same technique is used by (Levy and Y. Liu 2010a; J. [START_REF] Huang | Boundary aligned smooth 3D cross-frame field[END_REF][START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF]) and (Sokolov, Ray, Untereiner, and Levy 2016) who consider four decompositions into seven tetrahedra (Figure 6.6).

The second approach proposed by (Yamakawa et al. 2003) relies on the vertices and edges of the tetrahedral mesh T . Local searches are performed into the vertex-edge graph of T using two patterns. These vertex connectivity patterns generalize those proposed by [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF] and relax partially the dependency on the tetrahedral mesh. This method has been implemented by [START_REF] Baudouin | A frontal approach to hex-dominant mesh generation[END_REF] where a third pattern taking into account configurations with an interior flat tetrahedron was added. Some other approaches like H-Morph [START_REF] Owen | H-Morph: an indirect approach to advancing front hex meshing[END_REF] combine tetrahedra into hexahedra, while allowing for modifications of the connectivity and geometry of the input tetrahedral mesh (tetrahedron flips, node insertions, and node displacement). This great flexibility can make the algorithm intractable, but one advantage is that the mesh is valid throughout the procedure. (2) A decomposition into eight tetrahedra. This is a counter example to (Sokolov, Ray, Untereiner, and Levy 2016)'s claim that there is no decomposition of the hexahedron into more than seven interior tetrahedra. Figure 6.8: Two hexahedra with different edges and faces may be defined from the same set of tetrahedra.

Motivations for a New Approach

Important observations led us to work on improving these existing techniques. First, they do not identify the largest set H of potential hexahedra. On Figure 6.7 we gave two valid hexahedra that would neither be found by [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF]'s method nor by (Yamakawa et al. 2003)'s method. The first is a decomposition that encompasses one internal vertex, a configuration that may occur when a Steiner point is added when generating the tetrahedra. The second is a decomposition that has 8 interior tetrahedra. It is a counter example to (Sokolov, Ray, Untereiner, and Levy 2016)'s claim that there is no hexahedron decomposition with more than 7 interior tetrahedra. Both decompositions are not identified when searching for hexahedra made of 5, 6, or 7 interior tetrahedra. Neither are they by (Yamakawa et al. 2003)'s algorithm since none of their constitutive tetrahedra has three facets on the hexahedron boundary.

Second, as mentioned by (Yamakawa et al. 2003), several hexahedra may be defined using the same decomposition pattern by modifying the ordering of the vertices (Figure 6.8). The hexahedra have different edges and different faces while having the same tetrahedral decomposition. The hexahedron on the left being a perfect cube, the one on the right is undoubtedly invalid (zero Jacobian determinant), but were the vertices in a more general position, both could be valid. Third, the existing methods identify the same hexahedron several times. That number is as high as the number of corner tetrahedra in the decomposition in (Yamakawa et al. 2003)'s approach and depends on dual complex symmetries in [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF])'s approach.

Algorithm to Combine Tetrahedra into Hexahedra

2D Algorithm

Combinatorics Given a set V of n vertices labeled from 1 to n, let us first compute all possible quadrilaterals that can be defined from these vertices. We define a numbered quadrilateral abcd by the order of its 4 vertices and we define a non-oriented quadrilateral abcd by its edges ab, bc, cd and da ignoring orientation. Generating all numbered quadrilaterals that can be built from V is a combinatorial problem solved by a permutation generation algorithm that oup of 4 vertices of the n labels. When V = {1, 2, 3, 4}, the output is the set of the 24 permutations of 4 values. These 24 permutations define 3 different non-oriented quadrilaterals, each of them corresponding to 8 equivalent permutations. Adding constraints on the relative order of the vertices of one quadrilateral a < b, a < c, b < d, we obtain Algorithm 3 which fulfill our first objective and compute all possible non-oriented quadrilaterals that can be built from V. The output for V = {1, 2, 3, 4} is now the 3 non-oriented quadrilaterals that may be generated from 4 vertices: 1234, 1243, 1324.

Geometry Now, associate to each labeled vertex a point of R 2 . Among the 3 possible quadrilaterals that can be defined from 4 points of R 2 , only one is valid, i.e. has non-intersecting edges (Figure 6.9). We further associate to the set of vertices V a triangulation T and modify Algorithm 3 such that it generates all quadrilaterals whose edges are edges of the triangulation. The search for b and d is then restricted to the set of vertices connected to a. Similarly c should be connected through an edge to both b and d. The last step of the procedure is to identify the triangles subdividing each quadrilateral. Vertex selection order is now a, b, d, c instead of a, b, c, d since the choice of c depends on both b and d. All steps of the identification of quadrilaterals in a simple mesh are detailed on Figure 6.10b.

The advantage of this approach over the classical algorithms pairing adjacent triangles, e.g. (J.-F. [START_REF] Remacle | Blossom-Quad: A non-uniform quadrilateral mesh generator using a minimum-cost perfect-matching algorithm[END_REF]) is that it identifies quadrilaterals which encompass one (or more) vertex. An example is the quadrilateral {1245} on Figure 6.10a that encompasses vertex {3}. The other advantages of the algorithm are that it is easy to add geometrical quality tests (edge lengths, angles of the quadrilateral under construction) and that its parallelization is trivial. Its complexity may seem prohibitive but a vertex of a 2D triangulation is connected to an average of 6 other vertices. valid if the Jacobian determinant is strictly positive at any point of the element. All cells that have a negative Jacobian determinant are discarded.

The quality of a finite element cell is defined as the minimal value taken by the scaled Jacobian determinant over the element. If this value is inferior or equal to zero, the cell is invalid. In the first-order finite element cells we are considering (hexahedra, prisms, and pyramids) the maximum quality of the element is bounded by the quality at the corners, itself bounded by the quality of the facets sharing this vertex Q cell < min(Q corners ) < min(Q f acets ).

• The quality of a quadrilateral face corner abd is evaluated as the sinus of the angle made by the incident edges:

sin( ab, ad)

• The quality of a triangle facet corner abc [START_REF] Johnen | Robust and efficient validation of the linear hexahedral element[END_REF]) is evaluated as: • The quality of a pyramid base corner abde is evaluated as:

3 det(J) 2 3 J 2 F where J = J p J -1 I with J I =              1 0 1 2 0 1 1 2 0 0 √ 2 2              J p = ab ad ae
where . F denotes the Frobenius norm.

These quality values give an upper bound of the overall cell quality that we use to accelerate dramatically the cell identification of Algorithm 5. Indeed a new upper quality bound for the cell under construction can be computed when a vertex completing a face or a corner is added. When this bound becomes smaller than the required minimum quality, the cell construction is terminated. Additional quality tests on the planarity of quadrilateral facets, or on edge lengths could be added.

Guaranteeing that the Jacobian determinant of the final cells is strictly positive, is more challenging and time consuming and is the last test performed using the method of (Johnen, J.-F. Remacle, et al. 2013;[START_REF] Johnen | Robust and efficient validation of the linear hexahedral element[END_REF]. We modified the implementation so that it is exact and robust to floating point errors.

Hex-Dominant Meshing

To demonstrate that hex-dominant meshes may be generated from the set of potential cells identified by the previously described algorithm, we choose a subset of all compatible hexahedra, prisms, and pyramids.

Cell compatibility To have a valid final mesh, chosen cells should be mutually compatible. Two cells (hexahedron, prism, pyramid) are compatible if all following conditions are satisfied:

1. They share no interior tetrahedron. The incompatibilities between the remaining tetrahedra (not selected to build any cell) and the cells should also be accounted for. There are at least two possible strategies to manage them: [START_REF] Owen | Pyramid elements for maintaining tetrahedra to hexahedra conformability[END_REF] propose to raise pyramids on each non-conformal quadrilateral face and [START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF] propose to subdivide the pyramid or hexahedron incident to a non conformal contact into pyramids and tetrahedra. Both methods insert a new vertex, the apex of the pyramid or a point inside the neighboring cell algorithm. Their major drawback is that they increase the proportion of tetrahedra and pyramids compared to the proportion of hexahedra. In the meshes produced for this chapter, the compatibility condition between tetrahedra and cells is relaxed. Some quadrilateral faces will be adjacent to one or two triangles. This mesh should then be used with finite element solvers capable of handling these type of non-conformities.

Graph formalization

The selection of the cells of the final mesh can be reformulated as a Maximum Weight Independent Set (MWIS) problem [START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF]. Let us consider the graph G that has one vertex for each of the cell that may be built by combining tetrahedra of the input mesh T , and one edge linking each pair of compatible cells. The objective is then to find the largest subgraph in which all vertices are linked to one another. This is the Maximum Clique Problem (MCP), which is in general NP-hard. We may further associate a weight to each vertex depending on the cell quality. Since the compatibility graph G is usually very dense, it is advantageous to replace it with its complement, the incompatibility graph G * . The goal is then to find the solution to the Maximum Weighted Independent Set problem (MWIS).

When the graph G * contains up to a few hundreds of vertices, the optimal solution may be found by enumerating all independent sets and comparing their total weights [START_REF] Wu | A review on algorithms for maximum clique problems[END_REF]. However such an algorithm cannot be expected to terminate in a reasonable amount of time for graphs with a few thousands of vertices, let alone graphs with a few millions of vertices like the one we obtain. Reviewing the methods to solve the relevant MWIS problem is out of the scope. The interested reader is referred to [START_REF] Verhetsel | Solving the Maximum Weight Independent Set Problem: Application to Indirect Hexahedral Mesh Generation[END_REF]) and references therein for applicable methods in the specific case of indirect hex-dominant meshing.

Greedy solution The strategy we develop to obtain a hex-dominant mesh is the one used by previous works: greedily compute an approximate solution to the MWIS problem [START_REF] Baudouin | A frontal approach to hex-dominant mesh generation[END_REF][START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF][START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF]). The vertices (potential cells) are sorted by decreasing weight (quality), and independent vertices (compatible cells) are iteratively added to the solution in decreasing order of weight. This solution could be improved by taking advantage of the locality of the problem and optimizing small disjoint subgraphs [START_REF] Verhetsel | Solving the Maximum Weight Independent Set Problem: Application to Indirect Hexahedral Mesh Generation[END_REF]. Our non-optimized sequential implementation runs typically in a few seconds. The resulting hex-dominant meshes for three of the input tetrahedral meshes are shown on Figure 6.11.

Results and Discussion

We have applied our algorithm to 12 different tetrahedral meshes. Those were generated using the point placement strategy described in [START_REF] Baudouin | A frontal approach to hex-dominant mesh generation[END_REF]) and implemented in Gmsh (www.gmsh.info). They have between 127 and more than 3 million vertices (Table 6.1). The results of the tetrahedra combination algorithm in terms of numbers of detected potential hexahedra, prisms, pyramids and computational times are given in Table 6.2. The number of potential hexahedra, prisms, pyramids mostly depends on the number of vertices of the input tetrahedral mesh and on the minimal required quality.

Performances As expected, for a given input mesh, the higher the minimal quality, the faster the algorithm. For example, the running time on dataset Knuckle decrease from 214s to about 9s when quality increases (Table 6.2). Discarding as soon as possible cells that have a low quality is key to the efficiency of the algorithm. The multi-threaded version of our algorithm is very fast with about 300,000 potential hexahedra built per second. The algorithm for prisms generates about 900,000 prisms per second and the one for pyramids about more than 1.5 millions per second. All timings and performances are given for a laptop with 16Go RAM and an Intel Core TM i7-6700HQ CPU @2.60 GHz processor. For all cells, the running time clearly depends almost only on the number of potential cells detected. Note that for models Knuckle and Los1 the computed potential cells are only counted since they do not fit in the 16 Go RAM. We estimate that this method is more than ten times faster than the state of the art pattern matching method to identify combinations of tetrahedra into cells [START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF]. However, providing a valuable comparison is delicate because quality criteria have a strong impact on performances and are not explicitly stated in previous works, preventing comparison. Moreover no timings are provided for this one step of the hex-dominant meshing workflow.

Comparison of identified hexahedra with pattern based methods

We compare the number of potential hexahedra identified by our algorithm with the number of potential hexahedra that would be identified by pattern-matching methods [START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF][START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF] or vertex combination (Yamakawa et al. 2003) in Table 6.3. To count the potential hexahedra matching one of the six cube decompositions or one of the four decompositions into seven tetrahedra we compare the dual complex graphs. To count the potential hexahedra that would be detected by (Yamakawa et al. 2003), we count those containing a tetrahedron that has three facets on the hexahedron facets. On small models, our algorithm detects 4 to 5% more potential hexahedra than the existing methods. That number depends on the input tetrahedral mesh. This demonstrates that using a predefined set of templates for the hexahedron triangulation does not permit to detect all the potential hexahedra in a given tetrahedral mesh.

On larger meshes, the small difference between methods can be explained by the point placement strategy of the input tetrahedral mesh. The points are generated by propagation from the boundary of ,909,206 25,424,213 12,802,385 42.21 30.00 6.30 0.2 2,341,875 3,706,438 1,745,660 5.07 4.48 0.81 0.2 16,152,752 20,505,357 7,698,238 28.56 23.60 3.61 0.4 1,262,784 3,032,209 1,490,139 3.15 3.82 0.75 0.4 10,300,168 18,629,804 6,821,354 19.39 21.42 ,688,872 122,469,704 64,495,784 220.75 161.08 28.29 0.2 2,983,490 4,376,413 1,874,339 5.84 5.11 0.86 0.2 51,459,040 93,888,327 36,416,856 132.79 122.33 19.81 0.4 1,683,753 3,777,653 1,598,601 3.63 4.41 0.78 0.4 34,537,750 87,271,830 34,099,717 96.46 115.21 3,536,997 3,175,201 3,513,863 3,522,483 CrankShaft_2 4,405,892 3,919,768 4,383,763 4,392,081 Fusee_2 4,585,236 4,110,253 4,568,332 4,574,594 FT47_b 8,374,930 7,384,695 8,343,306 8,355,231 FT47 13,846,837 12,133,218 13,806,781 13,821,507 Fusee_3 17,048,021 14,983,008 17,010,173 17,023,768 Los1 21,908,307 19,308,740 21,865,212 21,880,467 Knuckle 86,553,836 79,544,742 86,346,178 86,433,942 (Yamakawa et al. 2003) the model. Where the fronts collide, a roughly 2-dimensional surface, point placement is not optimal. It is in this area that out method makes a difference, and the bigger the mesh is, the relatively smaller this area. Note that the higher the required minimum scaled Jacobian, the smaller the difference between the number of potential hexahedra detected by our method and the number of hexahedra detected by the existing methods. This is no surprise since the best quality is obtained for hexahedra that are close to the perfect cube which has a limited number of decompositions.

Challenges The algorithm presented in this chapter solves one step of the indirect hex-dominant meshing workflow that is the identification in a given tetrahedral mesh, of all the possible combinations of tetrahedral elements into hexahedra, prisms, pyramids. The algorithm identifies all the valid cells elements since no assumption is made on subdivisions into tetrahedra, each cell detected only once, bad quality cells are discarded early, and the algorithm is easy to implement. The C++ code is open-source and available at https://www.hextreme.eu/download/ and in www.gmsh.info. This is however only one step of the workflow to build hex-dominant meshes and all steps have a crucial impact on the final output. Indeed, the quality and number of potential hexahedra we can generate highly depend on the input tetrahedral mesh and then on the placement of its vertices. For example is the model Caliper shown on Figure 6.11 has complex geometrical features and points are non optimally placed around these preventing good quality hex-meshing. One key question is the optimization of such a mesh if that is possible, how can we increase the quality of such hybrid meshes (Chapter 9). The placement of the vertices is also one very important point for the generation of good hex-dominant mesh and is a very active research subject. The choice of the cells of the final mesh is the second key. We have shown that there may be up 30 times more cells than they are vertices in the input mesh, making the construction of the incompatibility graph very costly and the exact resolution of the MWIS problem intractable.

Part III

Theoretical Mesh Generation Chapter 7 Theoretical Tetrahedral Meshing Motivations The unexpected observation that there actually are much more patterns of subdivisions of the hexahedron into tetrahedra than previously thought when I implemented the algorithm that combines tetrahedra into hexahedra presented in Chapter 6 raised the (simple) question: How many more exactly? It is well known that the cube has five subdivisions into 6 tetrahedra and one subdivision into 5 tetrahedra. However, all hexahedra are not cubes and moving the vertex positions increases the number of subdivisions. Recent hexahedral dominant meshing methods try to take these configurations into account for combining tetrahedra into hexahedra, but fail to enumerate them all. The objective of this chapter is to answer this important theoretical question: How many different triangulations of the hexahedron into tetrahedra are there? And, more generally, is it possible to generate all the triangulations of a given polyhedron? and if it is, how do we do it? How many triangulations of the hexahedron are there? A finite element hexahedron is combinatorially a 3-dimensional cube whose square facets have been triangulated and whose geometry is defined as the image of a reference cube by a trilinear mapping denoted ϕ on Figure 7.1. Our hexahedra, in the physical space, are not strictly cubes, moreover they are only valid as long as the mapping from the reference cube is injective. Validity is difficult to evaluate because the Jacobian determinant is not constant over the hexahedron e.g. [START_REF] Chandrupatla | Introduction to finite elements in engineering[END_REF].

A set of tetrahedra is a valid triangulation of the hexahedron {12345678} (Figure 7.1) if it is a valid combinatorial 3-manifold, has 8 vertices and if its boundary contains the 8 vertices, 12 edges {12}, {14}, {15}, {23}, {26}, {34}, {36}, {48}, {56}, {58}, {67}, {78} and 12 triangular facets that can be paired into 6 quadrilaterals {1234}, {1265}, {1485}, {2376}, {3487}, {5678}.

Bounds on the number of tetrahedra that can be build with 8 vertices to mesh a 3-ball may determined from the Euler characteristic. Each triangulation is a 3-ball with Euler characteristic χ = 1, where χ = ve + ft. Then ve inte bd + f int + f bdt = 1. Since there are 4 triangular faces per tetrahedron and 2 tetrahedra per interior triangular face, we have 4t = 2 f int + f bd . Since the number of boundary edges e bd , and boundary triangular facets f bd are fixed the number of tetrahedra t in a hexahedron (without internal vertices) depends only on the number of interior edges e int . Moreover, there are at most n 2e bd edges in a cell with n vertices we then have trivial bounds on the number of tetrahedra in the triangulation of the cells. t = 5 + e int and 5 ≤ t hex ≤ 15. See also [START_REF] Edelsbrunner | Tetrahedrizing point sets in three dimensions[END_REF]) for additional combinatorial results. To enumerate the triangulations of the hexahedron we need discrete mathematics where enumeration is part of the job and where powerful concepts unheard of in the engineering world: matroids, chirotopes (no mythology, real math) permit to actually solve the problem. We solve the purely combinatorial problem of the hexahedron triangulations and show that there are 174 triangulations of the hexahedron up to isomorphism that have between 5 and 15 tetrahedra (Table 7.1). We prove this result by reviewing all the combinatorial triangulations of the 3-ball with eight vertices. We then solve the geometrical problem. The result that exactly 171 hexahedron triangulations have a geometrical realization in R 3 . To solve the realization problem we implemented the method proposed in [START_REF] Firsching | Realizability and inscribability for simplicial polytopes via nonlinear optimization[END_REF]. We exhibit triangulations with 14 and 15 tetrahedra of the hexahedron, contradicting the belief that there might only be subdivisions with up to 13 tetrahedra [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF][START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF]. This is an important theoretical result that links hexahedral meshes to tetrahedral meshes.

• Pellerin, J., K. Verhetsel, and J.-F. Remacle (2018). "There are 174 subdivisions of the hexahedron into tetrahedra". en. In: ACM Transactions on Graphics 37.6, pp. 1-9

The code and all result triangulations are available both in the supplementary materials of the paper and on the site of the project at https://www.hextreme.eu. This research was supported by the European Research Council (project HEXTREME, ERC-2015-AdG-694020).

Triangulations of the Cube

The triangulations of the hexahedron, limited to the specific geometrical configuration of the cube, were studied in discrete geometry [START_REF] De Loera | Triangulations: structures for algorithms and applications[END_REF] and in mesh generation methods that combine tetrahedra into hexahedra [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF][START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF][START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF].

The cube Following [START_REF] De Loera | Triangulations: structures for algorithms and applications[END_REF]), let us first look at the different types of tetrahedra that can be built from the eight vertices of a cube I 3 = [0, 1] 3 . Suppose that one facet of the tetrahedron is contained in one of the six cube facets, there are four possible vertices to built a non-degenerate tetrahedra that represent three cases due to symmetry (Figure 7.2):

• Corner tetrahedra have three facets on the cube faces.

• Staircase tetrahedra have two facets on the cube faces, they have two facets on the hexahedron boundary.

• Slanted tetrahedra have a unique facet on the cube faces.

The last type of tetrahedra are Core tetrahedra that have no facet on the cube faces. For the 3-cube, this list is complete and irredundant [START_REF] De Loera | Triangulations: structures for algorithms and applications[END_REF]. From this classification of the tetrahedra in a cube triangulation follows the classification of triangulations of the 3-cube. The 3-cube has exactly 74 triangulations that are classified in 6 classes (Figure 7.4) (see [START_REF] De Loera | Triangulations: structures for algorithms and applications[END_REF]) and references therein).

• Every triangulation of the 3-cube contains either a regular tetrahedron (i.e. a tetrahedron whose 6 edges are of equal lengths) or a diameter, i.e. an interior edge joining two opposite vertices.

• There are two triangulations of the first type, symmetric to one another. The triangulations of the second type are completely classified, modulo symmetries by their dual complexes (Figure 7.4).

The dual complex of a triangulation is a graph in which there is a vertex for each tetrahedron and one edge if corresponding tetrahedra share a triangle (Figure 7.4). By construction, a 2-cell of the dual complex corresponds to an interior edge of the triangulation.

There are at least two ways to prove that these six triangulations are indeed the six possible types of triangulations of the cube. Discrete geometers rely on the determination of the cell containing the barycenter of the cube and on symmetry considerations about the vertex arrangement around that cell [START_REF] De Loera | Triangulations: structures for algorithms and applications[END_REF]. [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF] enumerate the possible dual complexes.

The 3-cube in finite precision Floating point numbers have a finite precision and represent the cube I 3 = [0, 1] 3 up to some tolerance. Moreover, sets of points that are exactly cospherical and coplanar are often specifically processed because these non-general positions are problematic for many geometrical algorithms [START_REF] Edelsbrunner | Simulation of simplicity: a technique to cope with degenerate cases in geometric algorithms[END_REF]. Therefore, to combine elements of tetrahedral meshes, [START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF]) and [START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF] consider subdivisions of the cube into two prisms linked by a sliver, i.e. a very flat tetrahedra connecting four points that are coplanar up to a given tolerance. Note that [START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF] aim at enumerating the triangulations of more general configurations, but the geometrical argument used in the proof limits its validity to the cube represented in finite precision.

In a finite precision 3-cube, two new types of non-degenerate tetrahedra appear (Figure 7.3):

• Core diagonal tetrahedra have no facet on the hexahedron boundary. They subdivide the hexahedron into two prisms, creating four additional types of triangulations into seven tetrahedra (Figure 7.5).

• Boundary tetrahedra connect the four vertices of a facet of the cube. We exclude these tetrahedra from the triangulations we consider because they are not strictly inside the hexahedron. Indeed, the hexahedron facets are bilinear patches strictly inside boundary tetrahedra.

Together with the previous four types of tetrahedra in the cube they constitute the only possible configurations for a tetrahedron in a hexahedron triangulation.

Cube trianglations

Dual graphs Figure 7.5: The four triangulations of the finite precision cube with seven tetrahedra. Dashed edges connect pairs of tetrahedra incident to the same hexahedron facet.

Combinatorial Triangulations of the Hexahedron

The geometrical problem of the triangulation of the hexahedron is extremely difficult and we first consider the purely combinatorial side of the question. A set of tetrahedra is a valid combinatorial triangulation of the hexahedron {12345678} (Figure 7.1) if it is a valid combinatorial 3-manifold with 8 vertices and if its boundary contains the 8 vertices, 12 edges {12}, {14}, {15}, {23}, {26}, {34}, {36}, {48}, {56}, {58}, {67}, {78} and 12 triangular facets that can be paired into 6 quadrilaterals {1234}, {1265}, {1485}, {2376}, {3487}, {5678}. Theorem 7.2.1. The hexahedron {12345678} has 174 combinatorial triangulations up to isomorphism that do not contain any boundary tetrahedra.

To demonstrate Theorem 7.2.1 we first use a result about the triangulations of the 3-sphere with 9 vertices to compute all triangulations of the 3-ball with 8 vertices. At the next step, we select triangulations whose boundary is the boundary of a triangulated hexahedron. Finally the valid combinations are classified into isomorphism classes.

Triangulations of the 3-Ball

The sphere, or 2-sphere, is the 2-dimensional surface boundary of a 3-dimensional ball. The 3-sphere generalizes this, and is defined as the 3-dimensional boundary of a 4-dimensional ball. There is a direct link between the triangulations of the 2-sphere (3-sphere) and the triangulations of the 2-ball (3-sphere). Indeed, a triangulation of the 2-sphere can be constructed from the triangulation of a 2-ball by building a cone as illustrated on Figure 7.6. The cone is built by adding point v over a polygon and connecting it to the five boundary edges. The inverse transformation, the removal of one point v of the sphere triangulation as well as all triangles incident to v, permits to obtain the triangulation of a ball. (The interested reader is referred to [START_REF] Ziegler | Lectures on Polytopes[END_REF] for more details on this construction).

We enumerate all triangulations of the hexahedron, i.e. 3-dimensional ball with 8 vertices. They can be obtained from the triangulations of the 3-sphere with 9 vertices. The 1296 triangulations of the 3-sphere have been enumerated by [START_REF] Altshuler | Neighborly 4-polytopes with 9 vertices[END_REF]Altshuler and Steinberg 1974;[START_REF] Altshuler | An enumeration of combinatorial 3-manifolds with nine vertices[END_REF][START_REF] Altshuler | The classification of simplicial 3-spheres with nine vertices into polytopes and nonpolytopes[END_REF]) and are available online [START_REF] Lutz | The Manifold Page[END_REF]. Nine triangulations of the 3-ball with eight vertices can be built from each of the 1296 triangulations by removing one of the vertices v i , i = 1, . . . , 9 and its link, i.e. all tetrahedra incident to v i .

Could the naive brute force work? Let t be the number of tetrahedra, f the number of facets, e the number of edges and v = 8 the number of vertices. The hexahedron boundary is fixed to have 8 vertices, 12 facets and 18 edges. Each tetrahedron has 4 facets, each interior facet is incident to 2 tetrahedra, then 4t = 2 f + 12. Injecting this in the Euler-Poincaré characteristic of the triangulation χ = ve + ft = 1, we have t = e -13. There are at least 18 edges and at most 8 2 = 28 edges in the triangulation of a hexahedron, so 18 ≤ e ≤ 28, then 5 ≤ t ≤ 15. It is then possible to write a brute-force naive algorithm to enumerate the isomorphism-free triangulations of the hexahedron {12345678}. The exponential complexity of such an algorithm make it very expensive, but it ran in several hours in parallel thanks to the implementation of aggressive optimizations. .9: The correspondence between a hexahedron triangulations and its decomposition graph for triangulation 8_E.

Triangulations of the Hexahedron

The next step is to test if the obtained 3-ball triangulations are valid hexahedron triangulations. They are if their boundary matches the triangulated boundary of a hexahedron. The triangulation of the boundary of a hexahedron has 8 vertices and 18 edges. Among these, 12 are fixed and there are 2 possibilities to place the remaining 6 diagonals of the quadrilateral facets. We have then 2 6 = 64 possible triangulations. These triangulations can be classified into 7 equivalence classes, i.e. there are 7 triangulations of the hexahedron boundary up to isomorphism (Figure 7.7).

Comparison of Hexahedron Triangulations

The final step is to compare the obtained triangulations of the hexahedron and determine those that are the same up to isomorphism. We use a graph formalism proposed in [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF] to represent and compare the hexahedron triangulations.

The decomposition graph (Figure 7.8) of the triangulation of a hexahedron is an edge-colored graph where there is one vertex per tetrahedron, one black (plain) edge between vertices if the corresponding tetrahedra are adjacent, and one grey (dashed) edge between vertices if the corresponding tetrahedra have triangle faces on the same hexahedron facet. By construction, each vertex is of degree 4 and simple chordless cycles of plain edges correspond to interior edges of the triangulation.

There is a one-to-one correspondence between triangulations without boundary tetrahedra and decomposition graphs (Figure 7.9). Suppose that we fix the vertex labels of the triangulation such that {12345678} is a hexahedron. Then 12 of its edges are fixed ({12}, {14}, {15}, {23}, {26}, {34}, {36}, {48}, {56}, {58}, {67}, {78}). There is then one choice to subdivide each of the 6 facets into two triangles. Once the 6 boundary diagonal edges are chosen, the remaining degrees of freedom to triangulate the hexahedron are controlled by the interior facets. The plain edges of the decomposition graph represent the interior facets, while the 6 dashed edges represent the diagonal edges.

Two combinatorial triangulations are isomorphic, if and only if their decomposition graphs are isomorphic, i.e. are the same with respect to the relabeling of their vertices. More formally, an isomorphism of graphs G and H is a bijection between the vertex sets of G and H, f : V(G) → V(H), such that any two vertices u and v of G are adjacent if and only if f (u) and f (v) are adjacent in H.

The decomposition graphs of the triangulations of the hexahedron obtained in the previous section are therefore classified into isomorphism classes to obtain our main combinatorial result. Our reference implementation, attached in the supplementary materials, uses the nauty library [START_REF] Mckay | Practical graph isomorphism, II[END_REF]. This concludes the demonstration of Theorem 7.2.1.

Geometrical Realization

We now tackle the geometrical side of the problem of the triangulations of the hexahedron in R 3 . There is indeed no guarantee that a valid geometrical triangulation exists for each of the 174 combinatorial triangulations.

More formally, we are solving a geometrical realization problem for each combinatorial triangulation T , i.e. we are searching positions of the 8 vertices, A in R 3 , such that the geometrical triangulation defined by T is valid. A triangulation is valid if (1) the union of the tetrahedra is the convex hull of the points σ∈T conv σ = conv A and (2) there is no unwanted intersection between tetrahedra. conv σ ∩ conv σ = conv(σ ∩ σ ) ∀σ, σ ∈ T . Theorem 7.3.1. 171 of the 174 combinatorial triangulations of the hexahedron have a geometrical realization in R 3 .

We provide the realizations along with the code at in the supplemental material of (Pellerin, Verhetsel, et al. 2018). The 3 triangulations that do not have a realization have 15 tetrahedra and are illustrated on Figure 7.12. This is a particular case of the general problem of finding an embedding of an abstract triangulations of dimension m in R d . For our case, m = 3 and d = 3, the problem complexity is itself an open-problem [START_REF] Matoušek | Hardness of embedding simplicial complexes in Rd[END_REF]. For some values it is solvable in polynomial time (e.g. m = 2, d = 2), for others it is a NP-hard problem (e.g. m = 2, d = 4), for others it is algorithmically undecidable (e.g. m = 4, d = 5).

Our proof follows the strategy described in [START_REF] Pfeifle | Computing Triangulations Using Oriented Matroids[END_REF][START_REF] Schewe | Nonrealizable minimal vertex triangulations of surfaces: showing nonrealizability using oriented matroids and satisfiability solvers[END_REF][START_REF] Firsching | Realizability and inscribability for simplicial polytopes via nonlinear optimization[END_REF]) which consists in (1) solving the underlying combinatorial problem of the point configurations that can realize a given triangulation T , before (2) finding point coordinates in R 3 corresponding to one configuration. The main steps are described in Algorithm 6.

Useful Definitions

To encode the combinatorial properties of the point sets realizing a given combinatorial triangulation, we use oriented matroids. This is an important topic of discrete mathematics (Björner 1999), here we restrict definitions to the studied problem.

Let E be a finite set, r ∈ N with r ≤ |E|, and a mapping χ : E r → {-1, +1}. M = (E, χ) is a uniform oriented matroid of rank r if the following two conditions are satisfied:

1. The mapping χ is alternating, i.e. for all permutations π of {1, . . . , r} χ(e π( 1) , e π(2) , . . . , e π(r) ) = sgn(π)χ(e 1 , e 2 , . . . , e r ) where sgn(π) is the permutation sign, i.e. 1 for even permutations and -1 for odd permutations. A = solveSystem(inequalities) SCIP solver [START_REF] Gleixner | The SCIP Optimization Suite 5[END_REF] 8:

if A ∅ then return A Realization is found 9: undecidedRealizability = false 10: for all χ ∈ S do Check non-realizability 11:

if !findFinalPolynomial(χ) then 12: undecidedRealizability = true 2. For all σ ∈ E r-2 and all subsets {e 1 , e 2 , e 3 , e 4 } ⊆ E \ σ, we have that:

{-1, +1} ⊆ {χ(σ, e 1 , e 2 )χ(σ, e 3 , e 4 ),
-χ(σ, e 1 , e 3 )χ(σ, e 2 , e 4 ), χ(σ, e 1 , e 4 )χ(σ, e 2 , e 3 )}

The mapping χ is called the chirotope of the oriented matroid. Given a set of points A in R d , we can build a chirotope for a matroid of rank r = d + 1 can be constructed by considering the determinant of the points in homogeneous coordinates using v → 1 v . Set for any (a 1 , a 2 , . . . , a d+1 ) ∈ A

χ (1, 2, . . . , d + 1) = sign det 1 1 . . . 1 a 1 a 2 . . . a d+1
The chirotope is a sign function which gives the orientation of all possible subsets of r = d + 1 of points in A. The chirotope of a point set A ∈ R 2 (respectively R 3 ) is the orientation of all the triangles (respectively tetrahedra) that can be built from A (Figure 7.10). The chirotope encodes the point set structure, and, as we will see in Section 7.3.2, all necessary information to check if A is a realization of the combinatorial triangulation T .

Note that we only consider point sets A in general position (no 4-coplanar of 5-copsherical points) restricting ourselves to uniform oriented matroids, χ cannot take the value zero.

1 2 3 4 1 2 3 4 χ(1,2,3) = + χ(1,2,4) = + χ(1,3,4) = + χ(2,3,4) = + χ(1,2,3) = + χ(1,2,4) = + χ(1,3,4) = + χ(2,3,4) = -

Combinatorial Formulation of Geometrical Realizations

We first review the five constraints that the oriented matroids of point sets realizing a combinatorial triangulation T must enforce, before seeing how the oriented matroids are determined.

Constraints The point set is in general position therefore the corresponding oriented matroid should:

1. be uniform, χ 0 (no 4-coplanar or 5-cospherical points)

2. be acyclic, a property shared by the oriented matroids of all point sets.

3. be valid, an oriented matroid mapping. All r -2 subsets should satisfy condition (2) of the definition in Section 7.3.1.

The tetrahedra of the geometrical triangulation are valid, therefore the oriented matroid should:

4. define valid simplices, ∀σ ∈ T, χ(σ 1 , . . . , σ d+1 ) = +1. We assume that they are consistently oriented and are all positive in T .

The geometrical triangulation is valid, i.e. all tetrahedra of T intersect properly, i.e. for all pairs of tetrahedra σ and σ of T we have conv σ ∩ conv σ = conv(σ ∩ σ ). The intersection property can also be translated to constraints on the oriented matroids, let us see how.

If a point set A is in general position, any subset of d + 1 points is affinely independent while any subset of d + 2 points is affinely dependent, i.e. there exists λ ∈ R d+2 such that:

d+2 i=1 λ i a i = 0 and d+2 i=1 λ i = 0
We can then build two subsets of A, A + = {a i | λ i > 0} and A -= {p i | λ i < 0}, whose convex hulls intersect at a unique point (Radon's theorem). If two different simplices σ and σ of T intersect, then there is a subset of d + 2 points of vert(σ) ∪ vert(σ ) that are affinely dependent and corresponds to a pair (X + , X -). In oriented matroid jargon, this pair (X + , X -) is called a circuit of A, a fundamental notion related to another definition of oriented matroids (Björner 1999). All n d+2 subsets of (d + 2) points of A define a circuit (X + , X -) that can be computed from

X + = {i | (-1) i χ(1, 2, . . . , i -1, i + 1, . . . , d + 2) = +1} X -= {i | (-1) i χ(1, 2, . . . , i -1, i + 1, . . . , d + 2) = -1}
This give us the final constraint on oriented matroids of points sets realizing a combinatorial triangulation T , they should: 5. not admit any circuit (X + , X -) such that X + ⊆ σ and X -⊆ σ where σ and σ are two tetrahedra of T .

Solving The five constraints described above fix partially the values taken by χ for all n d+1 subsets of point labels. To completely determine the admissible oriented matroids for a triangulation T , we follow [START_REF] Schewe | Nonrealizable minimal vertex triangulations of surfaces: showing nonrealizability using oriented matroids and satisfiability solvers[END_REF], and solve an instance of a boolean satisfiability problem SAT that is satisfiable, if and only if, an admissible oriented matroid exists for a triangulation T . For each combination c = (c 1 , c 2 , . . . , c d+1 ) (where the elements are ordered by increasing index), we define a Boolean variable [c]. Setting [c] to true (resp. false) corresponds to setting χ(c) to 1 (resp. -1). We express the remaining axioms and conditions using those n d+1 Boolean variables, logical negations (¬ The chirotope must satisfy :

[c]), logical conjunctions ([c 1 ] ∧ [c 2 ] is true iff both [c 1 ] and [c 2 ] are true) and logical disjunctions ([c 1 ] ∨ [c 2 ] is true iff at least one of [c 1 ] and [c 2 ] is true).
1. the oriented matroid definition, that gives a constraints on all r = 4 subsets of elements. We use the formulation as a logical formula given by (Schewe 2010).

2. positivity. For any (d + 1)-simplex t ∈ T , [t] must be true.

3. the intersection property. For any signed circuit (X + , X -) with X = {p 1 , p 2 , . . . , p d+2 } that the oriented matroid cannot accept, we write

        p i ∈I [p 1 , p 2 , . . . , p i-1 , p i+1 , . . . , p d+2 ]         ∨          p j ∈J ¬[p 1 , p 2 , . . . , p j-1 , p j+1 , . . . , p d+2 ]         
where

I = {p i ∈ X + | (-1) i = -1} ∪ {p i ∈ X -| (-1) i = +1} J = {p i ∈ X + | (-1) i = +1} ∪ {p i ∈ X -| (-1) i = -1}
From the values of the n d+1 Boolean variables, we can determine the chirotope of an oriented matroid admissible for the abstract triangulation T . We must now find a point set, the chirotope of which matches.

We used the MiniSAT [START_REF] Sorensson | Minisat v1. 13-a sat solver with conflict-clause minimization[END_REF] solver to find all admissible chirotopes for 174 combinatorial triangulations of the hexahedron. Our implementation is provided in the supplemental material.

Result Three combinatorial triangulations 15_A, 15_B, 15_J do not admit an oriented matroid that enforce the five constraints and are therefore not realizable.

Geometric Realizations

Once the oriented matroids fulfilling all the constraints related to a given combinatorial triangulation are computed, the final step is to find real coordinates for the 8 vertices realizing one of these oriented matroids. We follow the approach of (Firsching 2017) who proposes to solve the system of n d+1 polynomial inequalities of the form:

χ(1, 2, . . . , d + 1) det 1 1 . . . 1 p 1 p 2 . . . p d+1 > 0
where the possible values of χ(1, 2, . . . , d + 1) have been computed in Section 7.3.2.

Solving In order to find feasible solutions we used SCIP, a solver that relies on branch and bound techniques [START_REF] Gleixner | The SCIP Optimization Suite 5[END_REF]. To encode the strict inequalities numerically, we rewrite the strictly positive constraint as superior or equal to an value (we used = 10 -4 ). When a solution exists, SCIP quickly converges and we obtain a realization of the combinatorial triangulation T .

Result All the 171 triangulations that admit an oriented matroid have a realization. All realizations are provided along our implementation in supplemental material. This completes the demonstration of Theorem 7.3.1. Figure 7.11 gives realizations for the configurations with 8 tetrahedra.

Convex Geometric Realization

We furthermore studied the realization problem with an additional constraint on the 8 point positions to ensure that they are in convex position, i.e. all vertices are on their convex hull (none is in its interior). This is a reasonable condition for hexahedra to be valid for finite element computations.

The convexity condition can be directly taken into account when testing admissible chirotopes by excluding all configurations in which a vertex is inside a tetrahedron defined by four other vertices.

Result Over the 171 realizable triangulations of the hexahedron, 13 do not have a realization such that their vertices are in convex position. 12_U, 12_Y, 13_T, 13_W, 14_A, 14_N, 14_O, 14_P, 14_R, 15_F, 15_H, 15_I do not admit a uniform oriented matroid and are therefore not realizable. Triangulation 15_G admits an oriented matroid, but SCIP does not terminate when searching for a realization (Figure 7.12). To prove that there is none, a different approach is necessary. A property of non-realizable matroids is that they admit a final polynomial, which can be found, if bi-quadratic, using the practical method of [START_REF] Bokowski | On the finding of final polynomials[END_REF]). It consists in proving the infeasibility of a linear program that encodes the second property of oriented matroid definition. 15_G oriented matroid has a final polynomial and is therefore not realizable with points in convex position.

Discussion

We demonstrated that there are 174 combinatorial triangulations of the hexahedron up to isomorphism and that 171 of them have a geometric realization. This result is consistent with the previous results on the triangulations of the 3-cube presented in [START_REF] De Loera | Triangulations: structures for algorithms and applications[END_REF]. However, Theorem 4 in [START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF], that states that there are only 10 triangulations of the hexahedron with 5, 6, or 7 tetrahedra, is incorrect. We described 161 additional geometrical triangulations of the hexahedron. The smallest triangulations for which no valid hexahedron may exist have 12 tetrahedra. The largest triangulations for which a valid hexahedron exists have 15 tetrahedra, by valid we mean a a strictly positive Jacobian (Figure 7.12). The existence of these triangulations contradicts the belief that there may only be subdivisions of the hexahedron with up to 13 tetrahedra. The practical consequence is that the methods solely based on the 10 triangulations of the cube to identify groups of tetrahedra that can be merged into a hexahedron do not compute all the possible hexahedra as claimed [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF]Yamakawa et al. 2003;[START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF][START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF]. Adapting these methods to take into account dozens more patterns is theoretically possible, but consequences on performances would be catastrophic.

Occurrences in real meshes

The theoretical problem of the triangulation of the hexahedron is an important step for a better understanding of the methods combining elements of a tetrahedral mesh to build hexahedra. Indeed, in practice, when combining the element of a tetrahedral mesh, many of the additional patterns we identified occur Figure 7.13.

To identify combinations of tetrahedra into hexahedra in a input tetrahedral mesh, we used the algorithm presented in Chapter 6. This algorithm provably computes all the possible combinations of eight vertices whose connectivity matches the one of a hexahedron. It further guarantees that all the identified hexahedra are valid for finite element computations using the exact test of (Johnen, J.-F. Remacle, et al. 2013). We compute all valid hexahedra for the input tetrahedral meshes provided Model Nb V 5 6 7 8 9 10 11 12 13 14 15 Total Cube 127 1 5 2 1 0 0 0 0 0 0 0 9 Fusee 11,975 1 5 5 7 13 9 4 0 0 0 0 44 CShaft 23,245 1 5 5 7 13 17 6 0 0 0 0 54 Fusee_1 71,947 1 5 5 7 7 1 0 0 0 0 0 26 Caliper 130, 572 1 5 5 7 7 1 0 0 0 0 0 26 CShaft2 140,985 1 5 5 7 8 0 1 0 0 0 0 27 Fusee_2 161,888 1 5 5 7 7 1 0 0 0 0 0 26 FT47_b 221,780 1 5 5 7 8 2 0 0 0 0 0 28 FT47 370,401 1 5 5 7 7 2 0 0 0 0 0 27 Fusee_3 501,021 1 5 5 7 8 4 0 0 0 0 0 30 Los1 583,561 1 5 5 7 7 2 0 0 0 0 0 27 Knuckle 3,058,481 1 5 5 7 8 2 0 0 0 0 0 28 Table 7.2: Number of triangulations patterns per number of tetrahedra counted in the available input data of [START_REF] Pellerin | Identifying combinations of tetrahedra into hexahedra: A vertex based strategy[END_REF]).

#points 5 6 7 8 9 10 11 12 13 14 15 Total 3,000 1 5 2 7 13 16 4 0 0 0 0 51 10,000 1 5 5 7 13 19 10 2 0 0 0 62 20,000 1 5 5 7 13 19 15 2 0 0 0 67 100,000 1 5 5 7 13 20 24 5 0 0 0 80 500,000 1 5 5 7 13 20 28 12 1 0 0 92 1,000,000 1 5 5 7 13 20 30 14 0 0 0 95 2,000,000 1 5 5 7 13 20 30 15 4 1 0 101 5,000,000 1 5 5 7 13 20 30 16 4 1 0 102 10,000,000 1 5 5 7 13 20 31 16 6 1 0 105 with that paper (Table 7.2) as well as for the Delaunay triangulations of random point sets (Table 7.3).

The decomposition graphs of all identified hexahedra are compared to the 171 possible triangulations and classified.

There is a clear link between the quality of the point set for hexahedral mesh generation and the number of patterns of combinations of tetrahedra into hexahedra. In random point sets, all existing realizable patterns up to 9 tetrahedra are found with as few as 3,000 randomly distributed vertices. In point sets generated for hexahedral meshing with the placement strategy described in [START_REF] Baudouin | A frontal approach to hex-dominant mesh generation[END_REF], the number of patterns is smaller and we note that the number of combination patterns is relatively independent of both the model and the number of vertices. Two models stand out and the high number of patterns is correlated to a lower number of hexahedra of worse quality in the final mesh. We further observed that occurrences of triangulations with more than 11 tetrahedra are rare. Since only the meshes are based on Delaunay triangulations, it is also probable that some combinatorial triangulations only have realizations that are not Delaunay.

Challenges Using our results, we would be able to decide on the existence of a triangulation of eight vertices whose boundary matches a hexahedron. We could also produce all existing triangulations of 8 given points. Since the principle of our demonstration is general, it could be extended to enumerate the triangulations of any polyhedron with quadrilateral and triangular faces and could be used to characterize non-triangulable polyhedra, a major difficulty in many geometric and combinatorial problems about which little is known theoretically [START_REF] Rambau | On a generalization of Schönhardt's polyhedron[END_REF]. We expect the theoretical advance of this chapter to be an important step toward a more complete study of the geometrical properties of hexahedral cells and will help develop robust indirect hexahedral mesh generation algorithms. Chapter 8

Theoretical Hexahedral Meshing

Motivations One fundamental problem in hexahedral mesh generation is that there is no algorithm to generate a hexahedral mesh conformal to a given quadrilateral boundary. The existence of hexahedral meshes for all even quadrangulations of the topological sphere has been proven by [START_REF] Mitchell | A characterization of the quadrilateral meshes of a surface which admit a compatible hexahedral mesh of the enclosed volume[END_REF]), yet seemingly innocuous quadrangulations such as the 16-quadrangle pyramid (Schneider's pyramid) or the 8-quadrangle tetragonal trapezohedron, also called octogonal spindle, (Figure 8.1) are notorious failure cases of general purpose meshing methods. [START_REF] Eppstein | Linear complexity hexahedral mesh generation[END_REF] shows how the interior of a quadrangulated sphere can be meshed with a linear number of hexahedra; this construction was later generalized to all domains which admit hexahedral meshes [START_REF] Erickson | Efficiently hex-meshing things with topology[END_REF]. Both methods reduce the problem to meshing a few quadrangulated spheres, but neither provide explicit hex meshes for these cases. One method only is able to generate hexahedral meshes for these templates [START_REF] Carbonera | A constructive approach to constrained hexahedral mesh generation[END_REF]. The huge drawback is that it requires 5396 n hexahedra to construct a non-degenerate hexahedral mesh of a ball bounded by n quadrangles, and is not directly applicable to arbitrary domains.

How do we generate combinatorial hexahedral meshes? In this chapter we are interested in the combinatorial constrained hexahedral meshing problem. The algorithms described are dedicated to building combinatorial hexahedral meshes and the geometrical problem is ignored for the most part of this work. Given a quadrangulation of the topological sphere Q, the objective is to determines a set of combinatorial cubes H such that:

1. the intersection of any two hexahedra h 1 , h 2 ∈ H is a combinatorial face shared by h 1 and h 2 (i.e.

the empty set, a vertex, an edge, or a quadrangle); 2. all quadrangular faces are shared by at most two hexahedra;

3. the set of boundary facets (adjacent to exactly one hexahedron) exactly matches Q.

This is an extremely challenging problem, even when the subsequent problem of finding a geometrical embedding is ignored, and for which no practical method exists.

Contributions This chapter presents results of my collaboration with Kilian Verhetsel, a brilliant PhD student I co-advised during my post-doc in Belgium. By applying constrained programming techniques we successfully solved problems that we only dreamt solving, such as the Schneider's pyramid hexahedral meshsing challenge. This project was funded by European Research Council project HEXTREME, ERC-2015-AdG-69402 awarded to Jean-François Remacle.

In this work with Kilian we first focused on the the hexahedral meshing of Schneider's pyramid. The goal was to determine the smallest hexahedral mesh of the pyramid. Following the idea explored to enumerate triangulations Chapter 7, we considered all possible groups of 8 vertices, all hexahedra, that can be built without creating an invalid mesh. The algorithm of this vertex-based enumeration of hexahedral meshes permit to compute lower bounds on the number of hexahedra and vertices required to mesh a given boundary Section 8.2. For small boundaries, that same algorithm can compute the smallest hexahedral mesh and we propose to locally remesh groups of adjacent hexahedra Section 8.3.

We compute a 44 hexahedra mesh of Schneider's pyramid by remeshing the 88 hex mesh of (Yamakawa et al. 2010).

The major contribution is the first practical algorithm to build constrained combinatorial hexahedral meshes (Section 8.4). The key idea of the method is to exploit the equivalence between quad flips in the boundary and the insertion of hexahedra glued to this boundary. The tree of all sequences of flipping operations is explored, searching for a path that transforms the input quadrangulation Q to the boundary of a cube. For larger meshes the search may stop when a sequence transforming Q to a set of pre-computed hexahedral meshes is found Section 8.5. Combined with an efficient backtracking search, it allows small shellable hexahedral meshes to be found for all even quadrangulations with up to 20 quadrangles. The 54, 943 such quadrangulations were meshed using no more than 72 hexahedra. This algorithm is able to compute small hexahedral meshes of quadrangulations for which the previously known best solutions could only be built by hand or contained thousands of hexahedra.

We finally prove that an arbitrary ball bounded by n quadrangles can be meshed using only 78 n hexahedra Section 8.6.1. This very significantly lowers the previous upper bound of 5396 n. The construction of Erickson is made fully explicit by computing hexahedral meshes for its two quadrangulated templates.

• Verhetsel, K., J. Pellerin, and J.-F. Remacle (2018). "A 44-element mesh of Schneiders' pyramid:

Bounding the difficulty of hex-meshing problems". In: Proceedings of the 27th International Meshing Roundtable. Springer Best paper award

• Verhetsel, K., J. Pellerin, and J.-F. Remacle (2019a). "A 44-element mesh of Schneiders' pyramid: Bounding the difficulty of hex-meshing problems". en. In: Computer-Aided Design 116, p. 102735

• Verhetsel, K., J. Pellerin, and J.-F. Remacle (2019b). "Finding hexahedrizations for small quadrangulations of the sphere". en. In: ACM Transactions on Graphics 38.4, pp. 1-13

Reference implementation of all algorithms presented in this chapter, as well as all result meshes are available at https://www.hextreme.eu.

State of the Art

Hexahedral mesh generation is a thriving field of research, with a variety of proposed methods. These include multi-block decomposition methods using frame-field parametrizations [START_REF] Kowalski | Block-structured Hexahedral Meshes for CAD Models Using 3D Frame Fields[END_REF][START_REF] Nieser | CubeCover-parameterization of 3D volumes[END_REF][START_REF] Liu | Singularity-constrained octahedral fields for hexahedral meshing[END_REF][START_REF] Lyon | HexEx: robust hexahedral mesh extraction[END_REF], hex-dominant meshing methods (Yamakawa et al. 2003;[START_REF] Gao | Robust hex-dominant mesh generation using field-guided polyhedral agglomeration[END_REF][START_REF] Baudouin | A frontal approach to hex-dominant mesh generation[END_REF][START_REF] Sokolov | Hexahedral-dominant meshing[END_REF][START_REF] Pellerin | Identifying combinations of tetrahedra into hexahedra: A vertex based strategy[END_REF], octree-based methods [START_REF] Maréchal | Advances in octree-based all-hexahedral mesh generation: Handling sharp features[END_REF][START_REF] Ito | Octree-based reasonable-quality hexahedral mesh generation using a new set of refinement templates[END_REF][START_REF] Zhang | A robust 2-refinement algorithm in octree and rhombic dodecahedral tree based all-hexahedral mesh generation[END_REF][START_REF] Qian | Sharp feature preservation in octree-based hexahedral mesh generation for CAD assembly models[END_REF], and polycube-based methods [START_REF] Gregson | All-hex mesh generation via volumetric polycube deformation[END_REF][START_REF] Han | Constructing hexahedral shell meshes via volumetric polycube maps[END_REF][START_REF] Yu | Optimizing polycube domain construction for hexahedral remeshing[END_REF][START_REF] Fang | All-hex meshing using closed-form induced polycube[END_REF]. Most methods do not address the problem of generating meshes with a given boundary quadrangular mesh. The rest of this section focuses on methods tackling the constrained problem.

Existence Proofs

Existence theorem for ball inputs [START_REF] Thurston | Hexahedral decomposition of polyhedra[END_REF] and [START_REF] Mitchell | A characterization of the quadrilateral meshes of a surface which admit a compatible hexahedral mesh of the enclosed volume[END_REF] independently showed that a ball bounded by a quadrangulated sphere can be meshed with hexahedra, if and only if, the number of quadrangles on the boundary, n, is even. The proof is based on the dual cell complex of quadrangular and hexahedral meshes. The dual complex of a quadrangular mesh is obtained by placing a vertex at the center of each quadrangle and adding edges between the vertices corresponding to adjacent quadrangles. Grouping edges traversing opposite edges of a same quadrangle, the dual complex is interpreted as an arrangement of curves (Figure 8.2). Similarly, the dual of a hexahedral mesh can be interpreted as an arrangement of surfaces [START_REF] Murdoch | The spatial twist continuum: A connectivity based method for representing all-hexahedral finite element meshes[END_REF]). In Mitchell proof, an arrangement of surfaces bounded by the dual arrangement of curves of the input quadrangulation is first constructed. For curves with an even number of self-intersections (including curves with no self-intersections), a disk is constructed inside the domain and a regular homotopy between a circle and the curve can be used to create a manifold bounded by that curve. Curves with an odd number of self-intersections are paired up arbitrarily. For each pair, a manifold bounded by the two curves is constructed by computing a regular homotopy between the two of them. This arrangement is not in general the dual of a hexahedral mesh, so the next step of the construction is to add new surfaces completely inside the ball until all connectivity requirements of a hexahedral mesh are met.

Linear-complexity meshing The construction of Mitchell can necessitate up to Ω(n 2 ) hexahedra where n is the number quadrangles. Eppstein showed this was the case and proposed a different construction which guarantees the use of O(n) hexahedra [START_REF] Eppstein | Linear complexity hexahedral mesh generation[END_REF]. His algorithm first subdivides each quadrangle into two triangles, so that a tetrahedral mesh of the interior can be computed.

After subdividing each tetrahedron into four hexahedra, a hexahedral mesh is obtained. However, its boundary does not match the initial input quadrangulation. This is solved by inserting buffer cells: for each quadrangle, add a cube, and glue one of its face to the original quadrangle; then, subdivide the opposite face into six quadrangles. The six new quadrangles are matched with those obtained from subdividing the original quadrangles during the previous step. The four remaining sides of the buffer cells are carefully subdivided into either two or three quadrangles, so that each buffer cell is bounded by an even number of quadrangles. Mitchell proof can then be invoked to show that each buffer cell can be subdivided into a finite number of hexahedra.

Generalization to other inputs Generalizing the previous results, [START_REF] Erickson | Efficiently hex-meshing things with topology[END_REF] gives necessary and sufficient conditions for the existence of a hexahedral mesh of a domain Ω bounded by a quadrangulation Q. The requirement is that every null-homologous subgraph of the input quadrangulation (i.e. every subgraph which bounds an embedded surface of Ω) contain an even number of edges. The construction of the hexahedral mesh is similar to the one proposed by Eppstein, and also starts by computing a tetrahedral mesh of the domain, subdividing it into a hexahedral mesh, and inserting buffer cells to get a complete mesh with the correct boundary. The last step is to subdivide the buffer cells of two different types (Figure 8.3) into hexahedra, which is again shown to be possible Mitchell 1996. Like Eppstein, Erickson does not give an explicit construction of the hexahedral meshes of the buffer cells that are the base of its proofs.

A constructive method [START_REF] Carbonera | A constructive approach to constrained hexahedral mesh generation[END_REF] give the first completely explicit construction. Their algorithm first adds hexahedra inside the domain, guaranteeing that the dual arrangement of the boundary of the remaining region contain no self-intersecting curve. Buffer cells are then inserted to transition to a mesh where each quadrangle has been subdivided into four quadrangles. The rest of the domain is then filled using pyramids. A complete hexahedral mesh is obtained after subdividing the pyramids into hexahedra. Given a topological ball bounded by n quadrangles, their construction produces a mesh of 76 n hexahedra. This mesh is degenerate: it contains quadrangles sharing multiple edges and hexahedra sharing multiple faces. A combinatorially valid mesh can be obtained by further refining the mesh [START_REF] Mitchell | Pillowing doublets: refining a mesh to ensure that faces share at most one edge[END_REF]. This method, however, requires as many as 5396 n hexahedra to build a hexhahedral mesh bounded by quadrangulation of size n.

Constrained Hexahedral Meshing in Practice

The methods presented in the previous section are impractical, they create far more hexahedra than necessary. In practice, less general methods have been used to obtain smaller meshes.

Whisker Weaving Whisker Weaving has been proposed by [START_REF] Tautges | The whisker weaving algorithm: a connectivitybased method for constructing all-hexahedral finite element meshes[END_REF]. The idea is to use a topological advancing front to construct the dual of a hexahedral mesh. The algorithm initially assumes that the final mesh will contain one dual surface for each dual curve of the input quadrangulation. Hexahedra are created inside the domain by creating intersections between three of these sheets, until the entire domain is filled. To choose between the multiple possible operations, heuristics based on geometric information such as the dihedral angle of faces are used. These heuristics are often not enough to completely fill the domain.

Dual cycle elimination (Müller-Hannemann 1999) proposed a method based on dual cycle eliminations. At each step of the algorithm, one of the curves of the dual mesh is removed, matching this elimination as the insertion of a layer of hexahedra. The new boundary after removing this cycle bounds the part of the input domain which has not been meshed yet. This process is repeated until the boundary matches that of a single cube. This method succeeds for certain classes of input quadrangulations, but fails for the common cases where the dual contains self-intersecting curves (e.g. Figure 8.2).

Searching for hexahedral meshes Some specific cases have particularly attracted the attention of the research community: Schneiders' pyramid [START_REF] Schneiders | Open problem[END_REF] and the octogonal spindle (Figure 8.1). In 2002, [START_REF] Yamakawa | Hexhoop: modular templates for converting a hex-dominant mesh to an all-hex mesh[END_REF] introduced the HEXHOOP template family and constructed by hand a hexahedral mesh of Schneiders' pyramid with 118 hexahedra. Later on, they improved their solution, building an 88-element mesh (Yamakawa et al. 2010). These ad-hoc constructions have long been the smallest known solutions for these two polyhedra.

Vertex-based Enumeration of Hexahedral Meshes

Following the work on triangulations Chapter 7, the idea we explored first was to consider all possible groups of 8 vertices that can be built without creating an invalid mesh. This means using computer search to exhaustively explore the space of all possible hexahedral meshes matching an input quadrangulation up to a given number of vertices. The search space is usually too large for a solution to be found except in small configurations. This allows to define cavity operations to locally modify the mesh and reduce the number of hexahedra Section 8.3. The algorithm is a backtracking algorithm that lists all possible combinatorial hexahedral meshes with a prescribed boundary. It is used to prove that there is no hexahedral mesh of Schneiders' pyramid with strictly fewer than 12 interior vertices. Using the same approach, we also prove that there is no hexahedral mesh of the octagonal spindle with strictly fewer than 21 interior vertices.

Backtrack Search Algorithm

Given ∂H, a combinatorial quad-mesh of a sphere or a handlebody, H max a maximum number of hexahedra, and V max a maximum number of vertices, our algorithm lists all combinatorial hexahedral meshes H such that:

• the boundary of H is ∂H,

• the number of hexahedra |H| is at most H max ,

• the total number of vertices in H is at most V max .

This problem has similarities with problems commonly encountered in constraint programming: (i) efficiently filtering a large set of potential solutions and (ii) managing solutions having multiple equivalent representations. Our implementation adopts concepts and strategies from this field. For a more general study of these problems, we refer the reader to (Rossi et al. 2006).

In Algorithm 7, the hexahedra are built one at a time by choosing a sequence of 8 vertices. At each step, all possible candidates for one of the 8 vertices are considered and the algorithm branches for each possibility. Each branch corresponds to the addition of a vertex to the current hexahedron. When a complete solution is determined, or when the search fails (no available candidates to complete a hexahedron), the algorithm backtracks to the previous choice. This process is repeated until all possibilities have been explored. This algorithm corresponds to the exploration of a search tree (Figure 8.4) where each branching node represents the choice of a vertex, and the leaves represent either solutions or failure points where the algorithm backtracks. The search tree has an exponential size in the maximum number of hexahedra in a solution. This high complexity is managed by pruning branches that cannot contain a solution and by using efficient implementations of all performed operations.

Search Space Reduction Strategies

In this section, we describe the key points of our implementation of Algorithm 7, all of which aim at reducing the search space explored by the algorithm:

• the order in which the hexahedra are constructed is crucial -we use an advancing-front strategy and start the construction of hexahedra from the boundary;

• an efficient filtering algorithm that eliminates candidate vertices that would create incompatible combinatorial hexahedra in the solution;

• a method to manage the high number of symmetries of this problem;

• the order in which the current hexahedron vertices are selected;

• a strategy to use topological invariants in order to filter out branches that do not contain any solutions. C ← Filter-Candidates(∂H, S , C)

9: if |C 1 | = • • • = |C 8 | = 1 then 10: S ← S ∪ {(v 1 , v 2 , v 3 , v 4 , v 5 , v 6 , v 7 , v 8 )}
C 4+i ← {0, 1, . . . , V max -1} \ {v 1 , v 2 , v 3 , v 4 } \ Interior-Diagonals(v i ) 7:
for each j ∈ {1, . . . , 4} do 8:

if i j then 9: C 4+i ← C 4+i \ Interior-Diagonals(v j ) \ Edges(v j ) 10: if i = j + 2 mod 4 then 11: C 4+i ← C 4+i \ Quad-Diagonals(v j ) return (C 1 , . . . , C 8 )
Advancing Front Construction While the hexahedra of a combinatorial mesh can be arbitrarily reordered, constructing them in a specific order makes the algorithm significantly faster. We use a classical advancing front generation strategy and require the hexahedron under construction to share a face with a front of quadrangles. There are then only four vertices needed to complete a hexahedron. The quadrangle front is constituted of the interior facets that are in only one hexahedron, or of boundary facets that are in no hexahedra. At the root of the search tree, it is set to be the boundary ∂H. An interior facet is added to the front after its first appearance in the mesh. The facet is removed from the front when it is added to the partial solution. When the front becomes empty, the boundary of the solution matches the input (Figure 8.5).

Filtering Candidate Vertices For each of the eight vertices of the hexahedron under construction, we store a set of candidate vertices that could be part of the solution. Some of these vertices would make the current hexahedron incompatible with some already existing hexahedra. Therefore when initiating the construction of a hexahedron, or when adding a vertex to a hexahedron, vertices that cannot be added without creating incompatibilities between the current hexahedron and the already built hexahedra are filtered out. The following rules are used to eliminate candidates:

1. the sets of edges, quadrangle diagonals, and interior diagonal of hexahedra are disjoint; 2. no two hexahedra may share an interior diagonal; 3. if one facet diagonal matches an existing quadrangle diagonal, so must the second one; 4. all eight vertices must be different.

To enforce these rules, our implementation tracks three sets of vertices for each vertex v: the sets of vertices u such that (u, v) is an edge, the diagonal of a quadrangle, or the diagonal of a hexahedron. These sets are updated whenever a new hexahedron is created. Because the execution time of the search algorithm blows up as the number of vertices increases, the number of vertices each set contains is always small, making them good candidates for being represented as bit-sets.

Symmetry breaking Combinatorial meshes are characterized by their large number of symmetries, a major challenge when operating on combinatorial hexahedral meshes. Indeed, a combinatorial hexahedral mesh has many equivalent representations:

1. n interior vertices can be labelled in n! different ways (boundary vertices labels are in input) 2. m hexahedra can be constructed therefore labelled in m! different ways; 3. each hexahedron, written as an ordered sequence of 8 vertices, has 1680 = 8!/24 equivalent vertex labelling.

The advancing front strategy defines the order in which the solution hexahedra are constructed (symmetry 2). This also uniquely determines the order of vertices in a hexahedron (symmetry 3). To prevent the relabelling of interior vertices (symmetry 1), we add value precedence constraints to our problem [START_REF] Law | Global constraints for integer and set value precedence[END_REF]. A solution H found by the algorithm can be written as an array of 8|H| integers, writing down the vertices of each hexahedron in the order in which they were constructed by the algorithm. In an array, x precedes y when the first occurrence of x is before the first occurrence of y. Enforcing a total precedence order on interior vertices, we guarantee that only one of their permutations is a solution.

Optimization of hexahedron construction

The efficiency of the backtracking search depends on the size of the search tree needed to explore all possibilities. A cheap approach to reduce the number of nodes in the search tree is to choose the vertex with the smallest set of candidate vertices when deciding which vertex to branch on [START_REF] Beck | Trying again to fail-first[END_REF]). This does not affect the correctness of the algorithm, as long as a vertex with more than one candidate is selected.

A similar heuristic is used to select the quadrangle on top of which the next hexahedron should be built: for each quadrangle on top of which a hexahedron may be built, the sets of candidates for the four remaining vertices are computed. The quadrangle which minimizes the product of the sizes of these four sets is then selected. This allows the search to detect failures early when it is not possible to build any hexahedron on top of an existing quadrangle.

Using topological properties during the search The meshes that we are searching for, in addition to being valid combinatorial meshes, need to be meshes of the interior of the input surface and some some topological requirements must be met by any solution. These requirements are used not only to filter out branches of the search tree that do not contain any valid meshes, but also to reject combinatorial meshes of different 3-manifolds with the same boundary. Only topological invariants that can efficiently be computed are considered during the search, since no efficient algorithms to recognize 3-manifolds are known, even for the 3-sphere [START_REF] Schleimer | Sphere recognition lies in NP[END_REF].

At any point, the partial mesh is maintained to be oriented: every quadrangle that appears in two hexahedra must have an opposite orientation in each of them. Whenever one of the faces of the hexahedron under construction is identified with an existing quadrangle because they share a diagonal, the other two vertices are selected so that the two quadrangles have opposite orientation.

Meshes of a topological ball, or any hexahedral mesh that could be used as a subset of a mesh of the ball, have a bipartite graph [START_REF] Eppstein | Linear complexity hexahedral mesh generation[END_REF]. Let A and B be the two parts of the partition. The sets of candidate vertices are reduced so that edges between two elements of A or of B are never created. The partition is computed initially based on the input quadrangulation, and updated every time a hexahedron is added to the partial mesh.

The last topological property that we use is related to homology groups, computed over Z 2 . A detailed introduction to homology computations can be found in [START_REF] Hatcher | Algebraic topology[END_REF]. We define a k-chain as a set of elements of dimension k -a 1-chain is a set of edges, a 2-chain a set of quadrangles, and a 3-chain a set of hexahedra. The boundary of a k-chain C is the (k -1)-chain ∂C whose elements are the faces contained in an odd number of elements of C. A k-chain whose boundary is empty is referred to as a cycle. A k-cycle which is the boundary of a (k + 1)-chain is called null-homologous (Figure 8.6).

For the domains that we consider, all 2-cycles are null-homologous. Since hexahedra have an even number of faces, any set of hexahedra is bounded by an even number of quadrangles. Therefore, if a 2-cycle containing an odd number of quadrangles is ever created, the search can backtrack, as it will never be possible to create a set of hexahedra bounded by this cycle. We therefore compute a basis for the space of 2-cycles using Gaussian elimination, and verify that every element of the basis contains an even number of quadrangles.

Parallel search implemntation

The exploration a search tree can be parallelized in a natural way by exploring different subtrees in parallel, making the algorithm much faster on parallel architectures (Figure 8.7). We use an approach similar to the embarrassingly parallel search of [START_REF] Régin | Embarrassingly parallel search[END_REF].

The main challenge to overcome is that some subtrees are multiple orders of magnitude larger than other ones without any possibility to determine it ahead of time.

We solve this issue by attributing many subtrees to each worker thread, so that all threads must on average perform the same amount of work (we used 4096 subtrees per thread). At the start of the search, the tree is explored in a breadth-first manner until a layer with enough subproblems is reached. The nodes of this layer are then explored in parallel by independent worker threads using Algorithm 7.

Lower Bounds for Schneiders Pyramid and the Octagonal Spindle

Using Algorithm 7, we computed lower bounds for the number of vertices and hexahedra required to mesh Schneiders' pyramid and the octagonal spindle (Figure 8.1). The algorithm is run multiple times, and we increment either V max or H max between each run. At each step, we verify that no solution was found by the algorithm. The time required to compute these bounds increases exponentially as the bounds become tighter (Figure 8 C ← C ∪ {h} return C

Cavity Based Hexahedral Remeshing Algorithm

The algorithm described in Section 8.2 can be used to find the smallest hexahedral mesh with a given boundary. In this section, we use this algorithm to compute upper bounds for the number of hexahedra required to mesh Schneiders' pyramid. From the 88-element solution of (Yamakawa et al. 2010), we locally simplify the mesh. By simplification we mean decreasing the number of hexahedra (Figure 8.9). The realized operations may be viewed as a generalized form of cube flips [START_REF] Bern | Flipping cubical meshes[END_REF] that substitute a set of hexahedra by another set without changing their boundary. However, instead of having a predefined set of operations, as e.g [START_REF] Tautges | Topology modification of hexahedral meshes using atomic dual-based operations[END_REF], the algorithm proposed in this section is able to operate on any group of hexahedra and automatically determines them at execution time.

Globally minimizing the number of hexahedra in the mesh is a computationally demanding task. Our algorithm therefore selects a small subset of the mesh, or cavity, and focuses on modifying the connectivity of the mesh only within this cavity. Our hexahedral mesh simplification algorithm is based on Algorithm 7. From a geometric hexahedral mesh it outputs a geometric hexahedral mesh whose boundary is strictly identical and which has fewer elements.

The mesh simplification procedure has three main steps:

1. the selection of a cavity, the group of hexahedra to simplify, C;

2. finding the smallest hexahedral mesh C min compatible with the cavity boundary ∂C and replacing the cavity with this smaller mesh;

3. untangling the hexahedra to determine valid coordinates for the mesh vertices.

Figure 8.9: The number of elements in a mesh can be reduced by operating locally on a cavity.

Cavity selection

The cavity selection algorithm is a greedy algorithm that starts from a random element of the input hexahedral mesh (Algorithm 9). When the target size, in terms of number of hexahedra, is reached, this process stops. The choice of a target cavity size is a trade-off between the cost of finding the hexahedral meshes of the cavity and the likelihood that the mesh can be simplified by remeshing the cavity. Cavities with many hexahedra are more likely to accept smaller meshes, but the cost of finding the smallest hexahedral mesh C min increases exponentially with the number of hexahedra in the cavity. In practice, we start by considering relatively small cavities containing up to 10 hexahedra, and increase this limit when no improvement is possible. We require the cavity to contain at least 4 interior vertices. Indeed, when there are no interior vertices (e.g. with a stack of hexahedra), it is not possible to remove any hexahedra. As the number of interior vertices increases, so does the likelihood that the cavity can be simplified.

Cavity remeshing To find a smaller mesh of the boundary of a cavity C, we first solve the combinatorial problem, i.e. we find the smallest combinatorial hexahedral mesh of ∂C, and then solve the geometric problem of finding valid coordinates for the modified mesh vertices.

The combinatorial problem of finding the smallest mesh of ∂C is a direct application of Algorithm 7 which enumerates all combinatorial meshes of a given surface. The maximum number of hexahedra H max of the solution is set to a smaller value than |C|. Changing the parity of a hexahedral mesh is known to be a difficult operation [START_REF] Schwartz | Construction techniques for cubical complexes, odd cubical 4-polytopes, and prescribed dual manifolds[END_REF]), so we set H max to |C| -2. We also set the limit to the number of interior vertices V max to one less than the number of interior vertices in C to accelerate the search.

There is a subtle but important difference between meshing a cavity in an existing mesh and meshing a stand-alone polyhedron: the hexahedra inside the cavity must be compatible with the other elements of the input mesh. An example where new elements from a cavity are not compatible with elements adjacent to the cavity is given in Figure 8.10. A 3-element cavity is replaced by 2 quadrangles, but one of these two quadrangles shares two edges with an existing element, which is an invalid configuration. To guarantee that the algorithm does not break the mesh validity, the data structures used to filter out inadequate vertex candidates (Section 8.2.2) are modified to take into account the hexahedra that are not part of the cavity.

Geometrical untangling

The previous step of the algorithm found a new connectivity for the mesh. The simplified mesh obtained by using this result is not valid in general because the interiors of hexahedra may intersect (Figure 8.11). To obtain a valid geometric mesh we use the untangling algorithm described in [START_REF] Toulorge | Robust untangling of curvilinear meshes[END_REF]. The vertices are iteratively moved until all hexahedra in the mesh are valid. If the untangling fails, connectivity changes are undone. The validity of the final mesh is evaluated with the method proposed by [START_REF] Johnen | Robust and efficient validation of the linear hexahedral element[END_REF].

Remeshing Schneiders' Pyramid and the Octagonal Spindle

A 66-element mesh of Schneiders' pyramid In this section, we build a hexahedral mesh of Schneiders' pyramid and the smallest known mesh of the octagonal spindle (Figure 8.13) using the remeshing algorithm. Starting from the 88-element solution of (Yamakawa et al. 2010), the number of hexahedra is reduced to 66 by locally simplifying groups of hexahedra. Table 8.1 shows the sizes of the different cavities simplified by our algorithm. It takes a few minutes for our algorithm to reduce the number of hexahedra in the mesh from 88 down to 66 in our final mesh. Figure 8.12 shows the changes to the connectivity of the mesh performed in two different iterations of the algorithm. The vertices had to be moved to obtain a valid mesh, but the combinatorial boundary remains the same. For example, for the second pair of cavities in the figure, the same 30 facets can be seen before and after the remeshing operation: there is a central facet, surrounded by a ring of five quadrangles, followed by three rings of six quadrangles, followed by one more ring of five quadrangles surrounding a single face. We also determined a combinatorial mesh with 64 hexahedra and 59 interior vertices, on which the untangling failed. Shortly after the first version of our conference paper was published on ArXiv an alternative algorithm was proposed by [START_REF] Xiang | A 36-Element Solution To Schneiders' Pyramid Hex-Meshing Problem And A Parity-Changing Template For Hex-Mesh Revision[END_REF]. Starting from a single cube, the algorithm of considers all possible ways to add one hexahedron while maintaining a mesh which is both valid and combinatorially equivalent to a topological ball. The process is stopped when the boundary of the mesh matches the target quadrangulation. The search space considered is smaller, and they however build a mesh of the pyramid with 36 hexahedra (Figure 8.13).

A 40-element mesh of the octogonal spindle We used the 72-element mesh constructed in one of the intermediate steps detailed in Table 8.1 to create a 40-element mesh of the octagonal spindle. Indeed, Schneiders' pyramid and the octagonal spindle are related by their dual graphs: replace each quadrangle by a vertex, and create an edge between each pair of adjacent quadrangles. If the dual edges that traverse opposite edges of a quadrangle in the primal are grouped together, a simple arrangement of curve is obtained. The dual of Schneiders' pyramid is obtained by adding one curve to the dual of the octagonal spindle (Figure 8.2).

This relationship determines a method to create hexahedral meshes of the octagonal spindle from certain meshes of Schneiders' pyramid. The dual of a hexahedral mesh is a simple arrangement of surfaces [START_REF] Murdoch | The spatial twist continuum: A connectivity based method for representing all-hexahedral finite element meshes[END_REF], where each surface is bounded by zero, one, or multiple curves of the dual arrangement of the boundary quad mesh (Figure 8.2). If the dual of a mesh of Schneiders' pyramid contains a surface which is bounded only by the curve present in Schneiders' pyramid but not in the octagonal spindle, that surface can be removed by collapsing all the edges that it traverses [START_REF] Borden | Hexahedral Sheet Extraction[END_REF]. The resulting mesh is a hexahedral mesh of the octagonal spindle.

In general, this operation may produce a degenerate mesh, with hexahedra sharing multiple quadrangles or quadrangles sharing multiple edges. This is what prevented the construction of a mesh of the octagonal spindle from the 88-element mesh of (Yamakawa et al. 2010). Applying this construction to our 72-element mesh, however, we obtain a new mesh the octagonal spindle, with 40 hexahedra and 42 interior vertices (Figure 8.14). This is the smallest known mesh of the octagonal spindle.

A 44-element mesh of Schneiders pyramid Our 44-element mesh of the pyramid was obtained by adding 4 hexahedra to our mesh of the spindle (Figure 8.13). One limitation of the hex-mesh simplification algorithm described in this section is that its execution becomes expensive, because finding the smallest hexahedral mesh of a cavity becomes exponentially more time consuming as its size increases. A cheaper algorithm could be designed by finding a small set of local operations and an algorithm to choose which of them to perform in order to reduce the size of the mesh.

Quad Flip Based Enumeration of Hexahedral Meshes

Enumerating hexahedral meshes from their vertices explore the space of all hexahedral meshes, but is extremely expensive, because of the very strong topological constraints. In this section we restrict the search space to a subset of all hexahedral meshes, the space of so-called shellable meshes. These hexahedral meshes can be explored efficiently using quad flips (Section 8.5.1). These quad flips are a set of operations to modify quadrilateral meshes and whose application can be interpreted as the construction of a hexahedron. They are the base of the practical algorithm we propose. Given a quadrangulated sphere Q, a hexahedral mesh bounded by Q is built by exploring the space of flipping operations that can be applied to Q. A solution is then obtained by finding a sequence of operations that transforms Q into the boundary of a cube (Section 8.4). When this search space is too large, the algorithm instead searches for a sequence of operations transforming Q into the boundary of any mesh within a library of pre-computed hexahedral meshes (Section 8.5). The shellable meshes space is explored in its entirety by considering all possible sequences of quad flips that correspond to valid hexahedral meshes (Section 8.4.1). Because many different sequences of flipping operations represent the same mesh, most of this work focuses on how to account for the symmetries of the input quadrangulation in order to avoid generating different sequences of quad flips corresponding to isomorphic hexahedral meshes (Section 8.4.2).

Algorithm 10 Search: Enumerate shellable hexahedral meshes Input: Q: A quadrangulation of the sphere; H: A partial mesh; H max : maximum number of hexahedra in a solution; V max : maximum number of vertices in a solution.

Construction of Shellable Hexahedral Meshes Using Quad Flips

Shellable hexahedral mesh The method only considers a specific class of meshes: shellable hexahedral meshes. Shellability is an important and useful combinatorial concept in the study of polytopes and cell complexes [START_REF] Ziegler | Lectures on Polytopes[END_REF]. Slightly different notions of shellability are found in the literature. We use that of pseudo-shellings [START_REF] Bern | Flipping cubical meshes[END_REF] or topology-preserving shellings (Müller-Hannemann 1999). This type of shelling is an ordering of the hexahedra (H 1 , H 2 , . . . , H n ) of a hexahedral mesh such that any prefix 0≤i<k H i is homeomorphic to a ball (Figure 8.15). This definition implies that any hexahedron H k must intersect the union of the previous hexahedra in one of six possible patterns.

The hex mesh -quad flips correspondence Gluing a hexahedron to one of these patterns modifies the boundary of the mesh locally (Figure 8.16). The transitions between these patterns are known as quad flips or bubble moves [START_REF] Funar | Cubulations mod bubble moves[END_REF]. These flipping operations are therefore a valuable building block to explore the space of shellable meshes.

Note that not all hexahedral meshes admit a shelling order -see for example Furch's ball [START_REF] Furch | Zur grundlegung der kombinatorischen topologie[END_REF]. Hence, by relying on these flipping operations to build hexahedral meshes, our method is inherently unable to construct certain meshes. Nonetheless, we can guarantee that a solution still exists: all quadrangulations of the sphere with an even number of quadrangles admit a shellable hexahedral mesh [START_REF] Bern | Flipping cubical meshes[END_REF]. Identifying possible flips For each quadrangulation Q visited during the search, all possible quad flips need to be identified. Each flip corresponds to a different hexahedron that can be inserted in the mesh. The algorithm successively tries adding all of them to the current mesh. Because flips are performed by starting from the target boundary, the hexahedra that are constructed during this process form the reverse of a shelling order (this is one difference with the search method of [START_REF] Xiang | A 36-Element Solution To Schneiders' Pyramid Hex-Meshing Problem And A Parity-Changing Template For Hex-Mesh Revision[END_REF]). (Müller-Hannemann 1999) construct the hexahedra in the same order, but our method, instead of only considering one mesh, explores the entire tree of possible sequences of quad flips.

+ + + +

The identification of all possible flips is split into two steps: first, the boundary Q is inspected to identify all occurrences of the 6 patterns from Figure 8.16. Second, those flips that correspond to the insertion of hexahedra that would make the mesh invalid are filtered out.

Checking the validity of the result mesh The hexahedron inserted by performing a flip is obtained by computing the union of the pattern before and after the flip. To determine whether or not this hexahedron is compatible with the mesh constructed so far, an efficient test is devised by considering three relations between the vertices of the mesh:

1. E, the edges of the mesh; 2. D Q , the diagonals of the quadrangles in the mesh; 3. D H , the interior diagonals of the hexahedra in the mesh.

These relations are disjoint in any combinatorial hexahedral mesh. For example, if a pair (u, v) is an edge, it is not the diagonal of any quadrangles or hexahedra. This leads to an efficient implementation of the test: simply maintain the three sets E, D Q , and D H , and verify that, after adding a new hexahedron:

1. the sets E, D Q , and D H remain disjoint; 2. new quadrangle facets in the hexahedron share no diagonals with any other quadrangle in the mesh;

3. none of the four interior diagonals of the new hexahedron are an interior diagonal of angleother hexahedron.

It is easy to verify that when any two hexahedra share only a vertex, an edge, or a quadrangle, these conditions are met. To verify their sufficiency, consider two hexahedra with an invalid intersection pattern. If an interior diagonal of one hexahedron is contained in the other hexahedron, one of the rules is always violated: rule 3 is violated if it is also an interior diagonal of the second hexahedron, and rule 1 is violated if it is an edge or the diagonal of a quadrangle. The only remaining cases to consider are those where the shared vertices are part of two distinct quadrangles. In all of those cases, the diagonal of one of those quadrangles appears in the other one. If it appears as an edge, rule 1 is violated; if not, both quadrangles have a shared diagonal, violating rule 2.

Symmetry Breaking

There are many distinct sequences of quad flips which represent identical hexahedral meshes. It is thus important to only consider a single representation for each hexahedral mesh constructed during the search, otherwise most of the computation time would be spent generating different representations of equivalent solutions.

A technique commonly used to deal with this type of issue is to define a canonical representation for objects under constructions, so that all those that belong to a given isomorphism class are transformed into the same representative element [START_REF] Burton | The pachner graph and the simplification of 3-sphere triangulations[END_REF][START_REF] Brinkmann | Fast generation of planar graphs[END_REF]. A significant portion of the execution time is then spent computing the canonical representations of partial solutions, which may completely change after every operation [START_REF] Jordan | Parallel Enumeration of Triangulations[END_REF]. The symmetry breaking method used within our algorithm instead compares partial solutions directly, and exploits the tree-shaped structure of the search in order to reuse results from previous computations.

The strategy described in this section is based on Symmetry Breaking via Dominance Detection (SBDD) [START_REF] Fahle | Symmetry Breaking[END_REF]. Consider the search tree explored by the algorithm: its nodes are partial meshes constructed during the search, and edges correspond to the insertion of new hexahedra through quad flips. The objective is to prune from this search tree nodes that correspond to meshes that have already been explored (up to symmetry). This is accomplished using the following steps:

1. The automorphism group of the input quadrangulation is pre-computed; 2. Traverse the search tree and encode fully explored subtrees into a sequence S ; 3. Determine if each new node should be pruned or notion by comparing it against the nodes stored in S .

Computing the automorphism group Given a quadrangulation Q, we compute the set of its symmetries, known as its automorphism group. A permutation σ of the vertices of Q is a symmetry if it preserves the set of quadrangles: the image of any quadrangle (a, b, c, d), (σ(a), σ(b), σ(c), σ(d)) is a quadrangle of Q, and any quadrangle (a, b, c, d) is the image of a quadrangle (σ -1 (a), σ -1 (b), σ -1 (c), σ -1 (d)). Note that the orientations of the quadrangles may be reversed by σ.

Symmetries are computed one after another, by fixing some quadrangle q A ∈ Q and assuming that its image under a symmetry σ is known to be q B ∈ Q. There are 8 different ways to map the vertices of q A onto the vertices of q B , corresponding to the 8 symmetries of a quadrangle. The entire permutation σ is uniquely determined by this part of the map: the quadrangles adjacent to q A must be the images of the quadrangles adjacent to q B under σ, and the quadrangles adjacent to those must also be images of each other, and so on, until the whole quadrangulation has been traversed. This process is well-defined because each edge is in at most two quadrangles.

The entire set of symmetries is computed by considering all 8|Q| possible ways to map an arbitrary quadrangle q A to any other quadrangle of Q. If an assumption is correct, a symmetry σ is obtained; if not, a contradiction will be reached when trying to construct the symmetry (two vertices mapping onto the same target vertex, or a single vertex with two images under σ). Because q A must be the image of some quadrangle under any symmetry σ, this process yields the entire automorphism group.

In the worst case, the entire automorphism group is determined in O(|Q| 2 ) operations. In practice, this quadratic time algorithm outperforms more complex linear time algorithms designed for planar graph isomorphism [START_REF] Eppstein | Linear complexity hexahedral mesh generation[END_REF][START_REF] Colbourn | Linear time automorphism algorithms for trees, interval graphs, and planar graphs[END_REF]) when applied to small quadrangulations, thanks to well-tuned heuristics. In particular, our implementation stops the algorithm as soon as

H 1 H 2 H 3 H 4 H 5 H 6 H 7 H 8 H 1 H 2 H 3 H 4 é H 5 H 6 é H 7
Figure 8.17: A partially explored search tree and the sequence used to compare the current node (in white) against the previously explored part of the tree. No-goods are shown in red, and partially explored subtrees in blue.

two vertices of different degree are mapped onto one another by the permutation under construction [START_REF] Brinkmann | Fast generation of planar graphs[END_REF]. Moreover, because this method does not use the planarity of the graph, it is also more general. The only requirement is that the input be a pseudomanifold: a combinatorial cell complex in which every facet is contained in at most two distinct cells. Indeed, a variant of this method will be used to compare hexahedral meshes in Section 8.4.2, by having quadrangles take over the role of edges.

Encoding the already explored partial meshes An efficient traversal of the search tree requires the search to stop as soon as the mesh under construction is the symmetric counterpart of a mesh that has previously been constructed. In the previous section, the set of symmetries that need to be considered was determined. This section now focuses on efficiently encoding the set of hexahedral meshes that have been constructed during the search.

Of course, the search tree is exponentially large, making it impossible to store every single mesh that is constructed during the search. Instead, SBDD only stores information about the roots of maximal fully explored subtrees, known as no-goods [START_REF] Gent | Symmetry in Constraint Programming[END_REF]. The current node should then be pruned if and only if the mesh under construction is the symmetric counterpart of one of the children of one of the no-goods. Note that a no-good is referred as such even if some of its children are solutions, since it is not desirable to compute the symmetric counterparts of those solutions.

No-goods can be stored efficiently thanks to the structure of the search tree. Recall that each node within the search tree corresponds to a partial hexahedral mesh, and each edge corresponds to the insertion of a hexahedron. Nodes with a common ancestor in the tree then share a common set of hexahedra as a prefix, and this prefix only needs to be stored once Figure 8.17. Upon visiting a new node, the most recently added hexahedron is inserted in the sequence, followed by a special branching symbol, indicating that the rest of the sequence will encode the children of this node. Upon backtracking, everything up to and including the last branching symbol of the sequence is removed.

Dominance detection and pruning

The last part of our symmetry breaking method is the test used to prune nodes of the search tree that do not need to be explored because any solution that could be found by doing so has already been found. These nodes are said to be dominated by one of the no-goods, i.e. they are the symmetric counterpart of one of the children of one of the nodes that have been previously explored and stored in the sequence Figure 8.17.

Since the search involves exploring exponentially many nodes, this dominance test must be implemented without explicitly comparing the current node against all previously explored nodes. Instead, this test is broken down into two steps:

1. Find a no-good such that all its hexahedra are contained in the current partial solution 2. Determine if the hexahedra that are in the partial solution but missing from the no-good could be inserted using flipping operations.

The sequence S constructed in the previous section is very valuable for this: not only does it save space by factoring out a common prefix, but it also saves time by allowing this prefix to be processed only once.

Let H be the current partial solution. The first step is to search within S for a partial mesh whose hexahedra are a subset of H. The process to find such a partial mesh is similar to the algorithm used to compute the automorphism group initially (Section 8.4.2). The goal is to construct σ, which maps the vertices of some partial mesh encoded in S to vertices of the current solution H, such that all hexahedra in the no-good are preserved by the map σ. The construction of σ again begins from an initial assumption, namely that the images of all boundary vertices through σ are known. Because boundary quadrangles must be preserved by σ, the set of possible initial assumptions is precisely the automorphism group that was previously computed.

The dominance test algorithm is executed once for each element of the automorphism group and consists in a traversal of S during which the map σ is extended. The process ends either upon finding a partial mesh contained in H or upon reaching a contradiction. For each hexahedron h found in S , we attempt to extend σ such that h maps to some hexahedron of the current solution H. Each hexahedron created by a quad flip shares at least one quadrangle with the boundary or with a previously created hexahedron. Because of this, each hexahedron in S has at least one quadrangle whose symmetric counterpart is known. It is therefore possible to search for the hexahedron h within the current solution H that contains this quadrangle (and has not already been determined to be the symmetric counterpart of another hexahedron).

If such a hexahedron h exists, it must be the symmetric counterpart of the hexahedron h encoded in S , and the map σ is extended accordingly. If h corresponds to a fully explored node (shown in red in Figure 8.17), the current partial solution H contains the symmetric counterparts of all hexahedra of the corresponding no-good. If, however, h corresponds to a partially explored node (shown in blue in Figure 8.17), and its symmetric counterpart cannot be found, the traversal ends early because all subsequent partial meshes encoded in S contain h, which is not in the current solution. In all other cases, the traversal of the sequence continues.

Clearly, containing all hexahedra from some no-good is a requirement for a node being dominated -all children of the no-good share this common prefix. There could still be cases where none of the children of this no-good contain all the hexahedra that are in the current node. In other words, it may be impossible to find a sequence of quad flips which inserts the missing hexahedra when starting from the no-good. Testing for the existence of such a sequence may appear intractable at first, because shellability is an NP-Complete property [START_REF] Goaoc | Shellability is NP-Complete[END_REF]. Thankfully, a correct test only needs not to produce any false positives, since false positives are the only reason a part of the search tree would incorrectly get pruned, causing solutions to be missed. Furthermore, because shellable meshes tend to accept many different shelling orders, there is a straightforward algorithm meeting this requirement and which very often computes the correct result: try a small number of permutations (say 10), then give up if no reverse shelling order was found.

Creation of hexahedra on the boundary

Precomputed mesh Output

Figure 8.18: Given a quadrangulation of the topological sphere, our algorithm creates hexahedra on the boundary until the unmeshed cavity matches the boundary of a pre-computed hex mesh that is merged to obtain the final combinatorial hexahedral mesh.

A Practical Algorithm to Compute Hexahedrization of Small Quadrangulations

The exhaustive search described in the previous section can only be used within small limits on the maximum number of hexahedra, because of its exponential execution time. In many cases, finding a complete shelling by searching exhaustively is too difficult: the sequence of flips to construct the smallest solution is too long, and the search tree contains many paths which transform the initial boundary into one which is more difficult to mesh, instead of being closer to a solution. Instead of searching for a sequence of quad flips that transforms the initial boundary Q into a cube (see line 5 in Algorithm 10), the key idea for solving larger cases is to stop the algorithm when a known configuration is found. For that purpose, we compute all boundaries that can be shelled with at most n hexahedra (say n ≤ 11). Using a list of all such boundaries and one of their shellings (Section 8.5.1), this variant of the algorithm can efficiently look up boundaries in the list during the search. This allows complete solutions to be constructed from any sequence of flips leading to any of the boundaries in the pre-computed set.

Pre-computing Small Shellable Meshes

Consider the flip graph for quadrangulations of the sphere: its nodes represent quadrangulations of the sphere, and arcs between these nodes represent a flip between two quadrangulations. A breadth-first traversal of this graph starting from the cube and stopped at depth n generates all quadrangulations that can be obtained using a sequence of up to n flips. To deal with cycles in this graph, previously visited quadrangulations are stored in a hash table. The hash value for quadrangulations is constructed from a signature based on the valence of vertices, and the isomorphism test two quadrangulations is performed using a variation on the symmetry breaking algorithm where the two starting quadrangles are part of different quadrangulations.

The signature used by our algorithm is a histogram of the valences of each vertex, followed by the number of edges connecting vertices of valence v a and valence v b , for any v a and v b where this number is non-zero. While this choice of signature causes collisions, it can be computed quickly and new entries can be inserted without necessarily needing a slower computation to find a unique canonical representation.

Not all quadrangulations generated in this breadth first search admit a shelling with up to n hexahedra: interpreting the flips performed during the traversal of the graph as the insertion of hexahedra, these hexahedra may not all be compatible. By explicitly testing for compatibility while performing the breadth first search (Algorithm 11), we obtain a greedy construction similar to the procedure outlined by [START_REF] Xiang | A 36-Element Solution To Schneiders' Pyramid Hex-Meshing Problem And A Parity-Changing Template For Hex-Mesh Revision[END_REF]: the hexahedra that are found are those which admit a shelling such that any prefix is the smallest shellable hexahedral mesh for the corresponding boundary. For large values of n, it is not clear that such a shelling should always exist, but we can verify this property for small values of n. For every quadrangulation found during the breadth-first search but without a hexahedral mesh found by Algorithm 11, Algorithm 10 is used to verify that there is indeed no shellable hexahedral mesh with at most n hexahedra. This test was performed for n ≤ 10, and no counter-examples were found.

From Algorithm 11, a table of 69, 043, 690 boundaries that can be meshed with up to 11 hexahedra is constructed (Table 8.2).

Using the Pre-computed Table

If at any point during the search, the boundary of the unmeshed region matches one of the pre-computed quadrangulations, the shelling of that quadrangulation is used to finish the meshing of that region.

The idea is to use the shelling computed in the previous section to fill the unmeshed region. Simply combining the two solutions is not always possible: this may produce an invalid mesh where, for example, two hexahedra share multiple quadrangles (Figure 8.19). Algorithm 10 could be used to compute all shellings of the unmeshed region with up to n hexahedra. If this search finds a shelling compatible with the partial solution constructed so far, a solution can be generated, at the cost of an additional computation.

Algorithm 11 Generate-Shellings: Generate small shellable hexahedral meshes Input n: maximum size for the generated hexahedral meshes Output H: a set of hexahedral shellings with up to n hexahedra in each mesh. Even if no such shelling was found, a solution can be constructed from any shelling of the unmeshed region, without performing an additional search or storing multiple hexahedrizations for each boundary: first construct a copy of the boundary of the unmeshed region, then, for each quadrangle of this boundary, create a hexahedron to connect each quadrangle to its copy. The hexahedra that have been inserted in this manner are guaranteed to be compatible with any hexahedrization of the unmeshed region, allowing a complete mesh to be constructed. When this approach is used, the first solution found by the algorithm is not in general the smallest. However, when the smallest solution contains a large number of hexahedra, this approach can construct solutions in many cases where methods with stronger guarantees fail to find any, because it adds several hexahedra without branching.

Efficient access to the pre-computed table is performed using a binary search. We create an array of all the quadrangulations we found, sorted by their signatures. To find the hexahedral mesh corresponding to a given quadrangulation, its signature is computed and an isomorphism test is performed on all quadrangulations in the table that have the same signature.

Results

A New Bound on Hexahedral Mesh Size

Only one previous solution to the constrained hexahedral meshing problem gives a completely explicit construction [START_REF] Carbonera | A constructive approach to constrained hexahedral mesh generation[END_REF]. This method requires 5396 n hexahedra to construct a valid mesh bounded by n quadrangles. In the following, we prove that this bound can be lowered to 78 n using the construction proposed by [START_REF] Erickson | Efficiently hex-meshing things with topology[END_REF].

Using our search algorithm (Section 8.5), we found hexahedral meshes for both types of buffer cells that Erickson's construction needs, along with geometric realizations using linear hexahedra (Figure 8.20) obtained by applying existing mesh untangling techniques [START_REF] Toulorge | Robust untangling of curvilinear meshes[END_REF][START_REF] Livesu | Practical hex-mesh optimization via edgecone rectification[END_REF], although the solutions have a very low minimum scaled Jacobian (Table 8.3). The meshes that we found contain 37 and 40 hexahedra. Because gluing multiple buffer cells together as needed by the construction would create a degenerate mesh, a hexahedron is added on each boundary quadrangle. The resulting meshes of the buffer cells have 57 and 62 hexahedra respectively, giving the following result:

Theorem 8.6.1. Let Ω be a compact and connected subset of R 3 bounded by a 2-manifold ∂Ω. Given a quadrangulation Q of ∂Ω, each component of Q containing an even number of quadrangles, and a triangulation T of Ω (splitting each quadrangle of Q into two triangles), if there is a combinatorial hexahedral mesh of Ω bounded by Q, then there is one with no more than 62|Q| + 8|T | hexahedra. In particular, if Ω is a ball (hence ∂Ω is a sphere) and |Q| is even, there is a combinatorial hexahedral mesh bounded by Q with no more than 78|Q| hexahedra.

Proof. Follow the construction of [START_REF] Erickson | Efficiently hex-meshing things with topology[END_REF] using the templates that we computed. There is one buffer cell for each boundary quadrangle, and each tetrahedron of the triangulation T is split into 4, 7, or 8 hexahedra. In the worst case, each buffer cell will be meshed with 62 hexahedra, and each tetrahedron will be split into 8 hexahedra.

If Ω is a ball, there is always a triangulation T with 2|Q| tetrahedra, obtained by arbitrarily splitting each quadrangle into two triangles, adding a vertex inside the domain, and joining each triangle to this new vertex by a tetrahedron. The bound for this special case is therefore 62|Q| + 8 × 2|Q| = 78|Q|.

A similar bound can be obtained for quadrangulations with an odd-number of quadrangles in some of their components. In that case, hexahedra are added to connect pairs of odd components, and 8.6.1 is used to compute the number of hexahedra to mesh the rest of the domain. 8.6.2 One Hexahedrization for each Quandrangulation of the Sphere We used the algorithm described in Section 8.5 to compute hexahedrizations for all even quadrangulations of the sphere containing up to 20 quadrangles (Table 8.4). The 54, 943 input quadrangulations were generated using plantri [START_REF] Brinkmann | Generation of simple quadrangulations of the sphere[END_REF]. We pre-computed shellable hexahedral meshes with up to 11 hexahedra. Of the 69, 043, 690 boundaries that were pre-computed, only 130 are included in the list of inputs. Nonetheless, in about 20% of all instances, the search for a solution terminates almost immediately after loading the set of pre-computed solutions (Figure 8.21). Only a few additional seconds are enough to find hexahedrizations bounded by most quadrangulations of the sphere.

Which are the most challenging boundaries? There are however some more difficult cases, requiring over an hour of computation time (Figure 8.23). The trapezohedron bounded by n faces, obtained by generalizing the tetragonal trapezohedron of Figure 8.1, is usually among the most difficult cases of a given size, requiring meshes with an intricate internal structure in order to be filled. For example, the smallest solution found for the 20-face decagonal trapezohedron contained 72 hexahedra, strictly more than any of the other boundaries (Figure 8.22). Similarly, the 16-face octagonal trapezohedron required 67 hexahedra, with the decagonal trapezohedron being the only boundary for which all solutions found were larger. The trapezohedra are also among the boundaries that require the most time before any solution could be found. The 14-face heptagonal trapezohedron is the second most time consuming input, requiring 2h 50min, and the 20-face decagonal trapezohedron is the third, requiring 2h 43min.

In the worst case, shown on Figure 8.23, it took 6h 15min before a 58-element mesh was found. Schneiders' pyramid and the octogonal spindle The quadrangulations of Figure 8.1 are not particularly difficult to mesh using our method. Indeed, most quadrangulations of up to 18 quadrangles require more computation time than these two cases -and even finding the smallest known solutions is orders of magnitude easier than finding any solutions for the cases shown on Figure 8.23. On a 4-core Intel R Core TM i7-7700HQ CPU, after pre-computing a list of shellable meshes with up to 10 hexahedra, the 36-element mesh of Schneiders' pyramid originally found by [START_REF] Xiang | A 36-Element Solution To Schneiders' Pyramid Hex-Meshing Problem And A Parity-Changing Template For Hex-Mesh Revision[END_REF]) is found within 3 seconds of search. A 44-element mesh of the tetragonal trapezohedron is found within 4 seconds and it takes 31 seconds to find the smallest known 40-element mesh constructed in (Verhetsel, Pellerin, and J.-F. Remacle 2019a). Statistics for the geometric realizations found for both meshes are shown on Table 8.3. This research, by allowing hexahedral meshes to be computed for any small quadrangulation of the sphere, opens up a wide array of new possibilities. It is an important step for hex-dominant meshing (Yamakawa et al. 2003), as this solves the combinatorial aspect of the problem of filling the cavities that those methods leave. This method also offers new insights into the structure of block decompositions for configurations for which state-of-the-art techniques such as (H. [START_REF] Liu | Singularity-constrained octahedral fields for hexahedral meshing[END_REF] fail to generate a valid block structure, by allowing combinatorial meshes to be computed in some of those cases.

Chapter 9

Perspectives

The best meshes are invisible. When you are watching the latest 3D animation movie or playing your favorite video game, you are not thinking (at all) about the meshes defining the geometry of all the 3D shapes around. If you are running numerical simulations, you wish you should not have to see the mesh, but you have to. If the mesh is not good enough the simulation will not run or will give incorrect results. You may dream about a one-click meshing method that builds the mesh you exactly want (without even knowing the criteria yourself) completely automatically at the right level of detail, fully adapted to the simulations you will be running next. Such a general automatic method might never exist and mesh generation will probably remain an expensive step in engineering analysis for many years.

The research contributions presented in this manuscript open promising perspectives to improve our understanding of the underlying combinatorial problems of mesh generation, and more particularly to develop algorithms for the generation and optimization of hexahedral and hybrid meshes. Chapter 4 presented a parallel algorithm to generate automatically tetrahedral meshes for general domains. Up to this date, there is no equivalent hexahedral meshing algorithm. My point of view is that the main difference is theoretical. Major theoretical results on triangulations lead to robust provable algorithms on one hand and to more practical algorithms used for numerical simulations. For hexahedral meshes, theoretical results are scarce and mostly topological. The contributions presented in Chapter 8 are important but remain topological. One main perspective is to study the geometrical counterparts of these results.

Subdividing and Combining Polyhedra

Studying the underlying theoretical problems is one way forward to new algorithms in mesh generation. My theoretical contributions Chapter 7 and Chapter 8 show that an adequate strategy is to first solve the combinatorial problems before taking into account geometry. Both the combination of tetrahedra (or other polyhedra) and subdivision into tetrahedra (or other polyhedra) are two basic problems on which we have but a limited knowledge and theoretical understanding. Exploring and understanding the subdivision (that is but the meshing) of simple polyhedra is the key to the development of new, potentially disruptive meshing algorithms. The objective would be to constitute a database of subdividing patterns of polyhedra into tetrahedra, hexahedra, or other polyhedra. In meshes, to the contrary of the usually admitted definition of polyhedra, the cell faces are not necessarily planar. Quadrilateral facets of finite elements are bilinear. The range of admissible geometries is extremely large and the only constraint to enforce is that the mapping must be injective, therefore the determinant of the Jacobian cannot take the zero value. The bilinear faces can be approximated by two triangles, defining this way, polyhedral with triangular facets.

Subdivision in tetrahedra Enumerating the triangulations of polyhedral cells with triangular of quadrilateral facets is related to two fundamental problems: the existence of a triangulation of a given polyhedron and the computation of the minimal and maximal subdivision of a polyhedron. The first one is an NP-complete problem [START_REF] Ruppert | On the difficulty of triangulating three-dimensional Nonconvex Polyhedra[END_REF]) and the second is NP-hard [START_REF] Below | The complexity of finding small triangulations of convex 3-polytopes[END_REF].

The solution to this problem is related to the vanishing not well understood, nor well characterized boundary between triangulable and non-triangulable polyhedra. This is a central problem for the boundary recovery step in tetrahedral meshing. Discrete geometry results and methods can be used to enumerate combinatorial triangulation as well as discriminate between realizable and non-realizable triangulations for a small number of vertices. Additional work is necessary to extend these results to other polyhedra than the hexahedron.

Subdivision in hexahedra

The most important perspective is to explore the subdivision of general polyhedra into hexahedra. Prisms, hexahedra themselves, and tetrahedra can very easily be subdivided into respectively 4, 6, and 8 hexahedra. As we have seen, there is no simple solution for the square based pyramid. This problem was described initially by [START_REF] Schneiders | A grid-based algorithm for the generation of hexahedral element meshes[END_REF]. We know that finding a good quality solution with a small number of hexahedra is not possible, since we proved that any solution must at least contain 17 hexahedra and we suspect that the smallest solution has 36 hexahedra (Chapter 8). The method developed to prove this result is particularly promising. The combinatorial meshing problem (in short, geometry is ignored) is extremely constrained and particularly difficult for hexahedral meshes. To ever be able to take the results we presented into account in meshing algorithms, we need to better understand which are the base operations and to take into account geometry. The method could be used for remeshing hexahedral meshes while ensuring the elements to be valid for computations.

Remark The meshing problem is an ill-posed problem, in most cases there are an infinite number of meshes of a given domain because there is an infinite number of positions for the vertices of the mesh. When the vertices are fixed, we know that the set of all the meshes that can be built from the points is finite, but we also know that we cannot compute that set in the general case in 3D. There is however a potentially huge impact in the exploration of all the possible meshes that can be generated from a given set of points.

Practical Meshing

Significant improvements of the main steps of the indirect meshing method (Section 6.3) could permit to build an algorithm that can be adapted to the needs of multiple applications.

Mesh vertices sampling Geometrical features (low angles, thin layers, high curvature) must be taken into account when sampling the mesh vertices because these geometrical characteristics are often in conflict with the carateristics desired for the cells (type of cells, minimal angle, edge size). 3D models, whatever the application, geological modeling, aeronautics, etc., include thin layers that are a challenges for meshing software because either a high number of points are added, or the points should be aligned to build prisms and hexahedra. No existing method is able to build these cells robustly. The measures of the difficulty to locally mesh a given geometrical model in Section 2.3 could be extended to locally determine the type of cell to build, i.e. to design measures that caraterise the size, the orientation, and the type of cells to generate.

Tetrahedral mesh generation Each step of mesh generation raises questions, and one important for indirect mesh generation is which tetrahedral mesh is adequate to combine tetrahedra. This is an open question and the geometrical and topological properties of the triangulations from a point of view of their ultimate combination in hexahedra are unknown. The answer is most probably not the Delaunay triangulation that is the base of most algorithms. It maximizes the smallest angle in the mesh, but shows not particular qualities to build tetrahedra with angles closer to 90 degrees.

Combine tetrahedra in hexahedra, prisms and pyramids Building all possible cells by combining tetrahedra into other cells is a now a solved problem (Section 6.3). But choosing the cells of the final mesh is much more difficult problem and one main obstacle to the usability of the method. From my point of view, the abstract formulation of the problem as a graph problem is not fully adapted to the meshing problem because meshing is a local geometrical problem.

Mesh improvement

The indirect method is very robust but will probably never be able to generate full-hexahedral meshes. Let us suppose that the obtained hex-dominant mesh could be modified so that all pyramids, prisms and tetrahedra are replaced by hexahedra. The real question is then: Do local operations that improve a hybrid mesh in a provable way exist? If they do exist, they are the key to adapt a hybrid mesh to particular application needs. Local operations, also called cavity operations, are the operations that allow the automatization of tetrahedral meshing algorithms, as we have seen in Section 4.1.1 and Section 8.3. A cavity operation takes as input a set of neighboring mesh cells that are deleted (the cavity) and replaced by another set of cells that have exactly the same boundary. For hexahedral meshes the base operations are not all known [START_REF] Ledoux | Topological modifications of hexahedral meshes via sheet operations: a theoretical study[END_REF]) and the status for hybrid meshes is worse. There might not even be a tentative list of the operations where pyramids, tetrahedra, hexahedra, prisms could be involved. Such operations would be extremely useful for robust and efficient meshing of complex geometries and mesh adaptation of hybrid meshes.

Final Words

Can you mesh my model? Having been working on mesh generation algorithm for some years, I have been asked this question many times. Answering right away is often not a good idea because tight definitions of the input and output are required to give a reliable answer. What is the input (exact file format, geometry, characteristics)? Is the input a valid solid model (watertight, tagged components, no holes, no gaps, no intersections, no hidden defects)? What is the desired output mesh (cell type, size, quality measure, level of detail, angles)? What is the mesh purpose? Depending on the exact answers, generating the desired mesh can be trivial, easy, difficult, extremely difficult, or simply cannot be done. Relatively easy in two dimensions, the difficulty of mesh generation is widely underestimated in three dimensions and mesh generation is a fascinating problem where one can always be surprised.

Could a software answer that question? In other words, can we predict the success or failure of a meshing algorithm on a given input geometry? This is one of the key of automatic mesh generation. The main problem are the geometrical features of 3D models that are (always) in conflict with the desired mesh properties that have been for decades an incredible challenge in mesh generation. The increasing size and complexity of models emphasizes this challenge. Whatever the type of mesh, whether tetrahedral, hexahedral, or hybrid, evaluating efficiently what can be done and cannot be done would be a great help to all users. Different constraints should be accounted for, constraints on cells (points, segments, faces, cells) that must be part of the final mesh, constraints on the number of points or cells inside a volume, on the distance between the mesh and the input model and so on and so forth.

A meshing dream Imagine we can compute all the subdivisions of a 3D polyhedron in tetrahedra, hexahedra, prisms, pyramids, or combinations of these. Imagine that these computations are extremely fast and that the cell qualities can be evaluated as fast. I believe then than most of the algorithmic challenges in mesh generation would be over. This is however a dream. For hexahedra, we have been stuck for a long time at trying to compute one subdivision. For tetrahedra, we can in most cases compute one subdivision in tetrahedra, but computing all of them to choose the best remains out of reach. I believe that studying those theoretical questions could be one of the key to find really new algorithms for 3D meshing, even if they do not appear at first relevant in practice. I believe that understanding in depth these problems that are encountered by our current algorithms is one of the key to more robust and more efficient algorithms.

Figure 2 . 1 :

 21 Figure 2.1: A structured mesh and a surface model representation of the same synthetic geological model.

Figure 2 Figure 2 . 3 :

 223 Figure 2.2: A 3D subsurface surface delimited by surfaces that represent boundaries between rock layershorizons -and by a network of discontinuities that displace these layers

Figure 2 . 4 :

 24 Figure 2.4: Some mesh invalidity configurations.

Figure 2

 2 Figure 2.5: A good looking model can be invalid. This subsurface model of the Annot sandstones (Salles et al. 2011) has degenerated mesh elements 5 colocated vertices associated to 1 duplicated triangle and 8 degenerate triangles, result in the invalidity of one surface and its two adjacent regions. 2 degenerate lines containing a unique degenerate edge, (2 pairs of colocated vertices) resulting in the invalidity of 2 surfaces each containing a degenerate facet built on this edge.

Figure 2 . 7 :

 27 Figure 2.7: Synthetic geological models leading to typical challenging configurations for meshing algorithms. (dimensions 16km × 9.3km × 5km).

Figure 2 . 8 :

 28 Figure 2.8: Cross-section views of some sources of complexity in geological structural models.

( 7 )

 7 Volume sampling measures Instead of being computed globally, the number of model components can be computed in the neighborhoods of points sampling the model. We propose to compute the numbers of regions, surfaces, lines, and corners in the cells of a model subdivision determined by a Centroidal Voronoi diagram of points distributed inside the model (see definition in Section 3.2.2).

Figure 2 . 10 :

 210 Figure 2.10: Classification of the 3D synthetic models using the number of model entities in 10 000 restricted centroidal Voronoi cells.

Figure 3 . 1 :Figure 3 . 2 :

 3132 Figure 3.1: General workflow of the the method to automatically correct 2D models.

Figure 3 . 3 :

 33 Figure 3.3: Geometrical and topological operations performed to automatically correct geological 2D models.

Figure 3

 3 Figure 3.4: (a) Input non-watertight cross-section (note the gaps and the overlaps on detailed views). (b) Sealed output cross-section. Gaps and overlaps have been removed, resolving non-watertightness features.

  Figure 3.5: Comparison of the reference simulation with the automatically simplified model M6 in terms of amplification and velocities.

Figure 3 . 6 :

 36 Figure 3.6: 3D Voronoi diagram (a) 200 sites are distributed in a box. (b) Solid slice in the Voronoi diagram of the sites cut by the box. (c) One Voronoi cell.

Figure 3 . 7 :

 37 Figure 3.7: 3D Voronoi-Delaunay dual relationship (a) Voronoi vertex V is equidistant from the sites of cells A, B, C and D. (b) To each Voronoi facet containing V (numbers 1 to 6) corresponds one segment linking the sites of the 2 cells sharing the facet. (c) To each Voronoi edge containing V (numbers 1 to 4) corresponds a triangle linking the sites of the 3 cells sharing this edge.

Figure 3 . 8 :

 38 Figure 3.8: Restricted Voronoi diagrams (RVD) and their dual restricted Delaunay triangulations (RDT) to a B-rep model surfaces or regions.

Figure 3

 3 Figure 3.9: (a) A 3D Voronoi cell. (b) Its restriction to the model regions has 4 connected components. (c) Its restriction to the model surfaces has 3 connected components.

Figure 3 .

 3 Figure 3.10: Optimization of 100 sampling points on a sphere. After optimization the restricted Voronoi diagram has compact well-shaped cells and the restricted Delaunay triangulation has almost equilateral triangles.

Figure 3 . 11 :

 311 Figure 3.11: Nested spheres remeshing. (a) 100 optimized points are sandwiched between two spheres (b) Each restricted Voronoi cell has 2 connected components (c) Dual of the connected components of the restricted Voronoi diagram (see Figure 3.12).

Figure 3 . 12 :

 312 Figure 3.12: Remeshing two close surfaces. (a) 3 points (A, B and C) are sandwiched between two close surfaces, their restricted Voronoi cell s have two connected components, but in the dual restricted Delaunay triangulation there is only one triangle ABC. (b) Each point is replaced by two points. (c) Two triangles A 1 B 1 C 1 and A 2 B 2 C 2 are built corresponding to the points v 1 and v 2 shared by three restricted cell connected components.

Figure 3 . 13 :

 313 Figure 3.13: Remeshing a surface with a boundary. (a) 21 points sample the star; Voronoi cells of the white points intersect twice the boundary, those of the black points intersect it once or not at all. (b) Each white point is replaced by two points, one per intersection of the cell with the boundary (c) The elements dual of the polygons of the final mesh are: regular restricted Voronoi vertices (white), restricted Voronoi vertices neighboring at least a restricted Voronoi cell with two points (gray), and the segments on the restricted Voronoi edges intersecting twice the boundary. (d) Final mesh is made of regular Delaunay triangles and quads.

Figure 3 .

 3 Figure 3.14: Configurations leading to modifications of the model. (a) The central restricted Voronoi cell component intersect 3 Corners (A, B, C), the polygons to build with these points ABGF and ADEC intersect. We merge the 3 corners. (b) Contact lines (black lines) cut the cell into 6 connected components. The 4 triple Corners (A, B C, D) are merged. (c & d) If d min is inferior to the given resolution, points A and B are merged.

Figure 3 . 15 :

 315 Figure 3.15: Coal veins remeshing. 29 sub-vertical surfaces delimit thin coal veins. 1000 sampling points are sufficient to remesh the model decreasing the number of triangles from nearly one million to 35 thousands.

Figure 3 . 16 :

 316 Figure 3.16: Forward model remeshing with varying resolutions. Input model shows 3 three challenges for remeshing, very thin layers, major layer thickness variations and low-angle contacts between horizons due to onlapping geometries.

Figure 3 . 17 :

 317 Figure 3.17: Remeshing the Detachment, Leipzig, and Lambda models

Figure 3 . 19 :

 319 Figure 3.19: Cloudspin model remeshing. (a) Input model with very small throw near fault ends. (b) Output surfaces, remeshing was done with 5000 sampling points, contact lines are locally merged.

Figure 3 . 20 :

 320 Figure 3.20: Nancy model remeshing. (a) Input model (b) Model remeshed with 10000 sampling points.

Figure 3 .

 3 Figure 3.21: Challenging model remeshing: Sandbox and Annot.

Figure

  Figure 3.22: Remeshing the Clyde model with 30000 and 10000 sampling points

. 1 )

 1 where the Delaunay cavity C(DT k , p k+1 ) is the set of all tetrahedra whose circumsphere contains p k+1 (Figure 4.1b), whereas the Delaunay ball B(DT k , p k+1 ) is a set of tetrahedra filling up the polyhedral hole obtained by removing the Delaunay cavity C(DT k , p k+1 ) from DT k (Figure 4.1c).

Algorithm 1

 1 Sequential computation of the Delaunay triangulation DT of a set of vertices S Input: S Output: DT (

Figure 4 . 1 :

 41 Figure 4.1: Insertion of a vertex p k+1 in the Delaunay triangulation DT k . (a) The triangle containing p k+1 is obtained by walking toward p k+1 . The Walk starts from τ ∈ B p k . (b)The Cavity function finds all cavity triangles (orange) whose circumcircle contains the vertex p k+1 . They are deleted, while cavity adjacent triangles (green) are kept. (c) The DelaunayBall function creates new triangles (blue) by connecting p k+1 to the edges of the cavity boundary.

1 :

 1 The aligned mesh data structure mesh_t we use to store the vertices and tetrahedra of a Delaunay triangulation in 3D.

Figure 4 . 3 :

 43 Figure 4.3: Performances of HXTSort for sorting {key, value} pairs on an Intel Xeon Phi TM 7210 CPU and comparison with widely used implementations.

Figure 4 . 4 :

 44 Figure 4.4: Performances of our sequential Delaunay triangulation implementation (Algorithm 1) on a laptop (a) and on a slow CPU having AVX-512 vectorized instructions (b). Timings are in seconds and exclude the initial spatial sort.

Figure 4 . 6 :

 46 Figure 4.6: Partitioning of 20 points in 2D using the Moore indices, on the right the supporting grid of the Moore curve is transformed and the curve is shifted. In both cases, each partition contains 5 points. Indeed, the starting and ending Moore index of each partition are defined in a way that balances the point insertions between threads.

Figure 4 . 7 :

 47 Figure 4.7: Strong scaling of our parallel Delaunay for a random uniform distribution of 15 million points, resulting in over 100 million tetrahedra on 3 machines: a quad-core laptop, an Intel Xeon Phi with 64 cores and a dual-socket AMD EPYC 2 × 32 cores.

Figure 4 . 8 :

 48 Figure 4.8: Number of tetrahedra created per second by our parallel implementation for different number of points. Tetrahedra are created more quickly when there is a lot of points because the proportion of conflicts is lower. An average rate of 65 million tetrahedra created per second is obtained on the EPYC.

Figure 4 . 7 .

 47 We obtain a good scaling until the number of threads reach the number of cores, i.e. 4 cores for the Intel i7-6700HQ, 64 cores for the Intel Xeon Phi 7210 and the AMD EPYC 7551. Comparisons of implementation are given on on a laptop (Figure 4.9a) and on the Intel Xeon Phi (Figure 4.9b).

Figure 4 . 9 :

 49 Figure 4.9: Comparison of our parallel implementation with the parallel implementation in CGAL (The CGAL Project 2016) and Geogram (Levy 2016) with on a high-end laptop (a) and a many-core computer (b). All timings are in seconds.

Algorithm 2

 2 Mesh generation algorithmInput: A set of triangles t Output: A tetrahedral mesh T 1: function Parallel Mesher(t)

Figure 4 . 10 :

 410 Figure 4.10: Performances of our parallel mesh generator on a Intel Core TM i7-6700HQ 4-core CPU. Wall clock times are given for the whole meshing process for 1 to 8 threads. They include IOs (sequential), initial mesh generation (parallel), as well as sequential boundary recovery (BR), and parallel Delaunay refinement for which detailed timings are given.

Figure 4 . 11 :

 411 Figure 4.11: Performances of our parallel mesh generator on a AMD EPYC 64-core machine. Wall clock times are given for the whole meshing process for 1 to 127 threads. They include IOs (sequential), initial mesh generation (parallel), as well as sequential boundary recovery (BR), and parallel Delaunay refinement for which detailed timings are given.

  2.3. The method we propose generates a hybrid mesh (tetrahedra, prisms, pyramids) from the intersection of the Voronoi diagram of given points and a surface model. The Voronoi diagram determines a global subdivision of the model while the mesh building strategy is local: the type of the generated elements depends on the local thickness of the region considered and on the positions of the points sampling the model. The vertices, edges, facets, and cells of the final volumetric mesh are determined from the combinatorial analysis of the intersections between the model elements and a Voronoi diagram of sampling points. Where the intersections of the Voronoi cells with the model surfaces have a unique connected component, tetrahedra are modified to fit the input triangulated surfaces. Where these intersections are more complicated, a correspondence between the elements of the Voronoi diagram and the elements of the mixed-element mesh is used to build the final volumetric mesh.Preliminary results were presented at the International Meshing Roundtable in 2012 and a more complete version in 2014.• Pellerin, J., B.Lévy, and G. Caumon (2012). "A Voronoi-Based Hybrid Meshing Method". In: 21st International Meshing Roundtable, Research Notes. San Jose, CA, USA • Pellerin, J., B. Lévy, and G. Caumon (2014). "Toward Mixed-element Meshing based on Restricted Voronoi Diagrams". en. In: Procedia Engineering 82. Proceedings of the 23rd International Meshing Roundtable, London, UK, pp. 279-290 5.1 Twisting the Restricted Delaunay Triangulation5.1.1 Main Idea The idea we developed is modifying the restricted version of the Voronoi-Delaunay duality (Section 3.2.2) and propose a correspondence that considers the connected components of both the intersection of the Voronoi diagram with the regions and the intersection of the Voronoi diagram with the surfaces. From a centroidal Voronoi diagram Section 3.2.3.1, we propose to take into account the connected components of the restricted Voronoi diagram to the regions, surfaces to build the final mesh. In the regular restricted Delaunay triangulation the different connected components are ignored Figure 3.8 When restricted Voronoi cells to the regions may intersect one surface (Figure 5.1a). In the classical Voronoi-Delaunay relationship, the dual of a point shared by such cells is a tetrahedron linking the four sites. We propose to modify the tetrahedron vertices and instead of the sampling points, take the centroids of the adjacent restricted Voronoi diagram to the surfaces (Figure 5.1b). In thin sections of the

Figure 5 . 1 :

 51 Figure 5.1: Delaunay tetrahedron modification near one surface.

Figure 5 . 2 :

 52 Figure5.2: Prism built to mesh a thin section, it corresponds to the restricted Voronoi edge to this thin section e.

Figure 5 . 3 :

 53 Figure 5.3: Cell building. (a) Restricted Voronoi vertex v to a thin section is shared by four restricted Voronoi cells (A, B, C, D), 6 facets, and 4 edges. (b) To each restricted cell correspond two vertices and one edge. (c) To each facet correspond one (facets 4 and 6) or two edges (facets 1, 2, 3 and 5). (d) Final cell has eight vertices, 4 triangle facets, and 4 quad facets.

Figure 5 . 4 :

 54 Figure 5.4: Invalid cells with 5 vertices. (a) Vertices D 1 and D 2 correspond to the same Voronoi cell. Facets ABD 1 D 2 and BCD 1 D 2 share three vertices. (b) Cutting these facets with BD 2 results in a tetrahedron plus a triangular facet.

Figure 5 . 5 :

 55 Figure 5.5: Two meshes of 4 nested spheres computed with two different samplings

Figure 5 . 7 :

 57 Figure 5.7: Mixed element mesh inside a smooth heptoroid surface, 516 tetrahedra (gray), 665 prisms (white), and 507 pyramids (black).

Figure 6 . 1 :Figure 6 . 2 :

 6162 Figure 6.1: Indirect hexahedral dominant meshing principle.

Figure 6 . 3 :Figure 6 . 4 :

 6364 Figure 6.3: 3D cell templates: tetrahedron, pyramid, prism, and hexahedron. These linear finite element cells have straight edges, planar triangulat facets, and bilinear quadrangular facets.

Figure 6 . 5 :

 65 Figure 6.5: The six types of triangulations of the 3-cube and their dual complex representations.

Figure 6 . 6 :

 66 Figure 6.6: The four types of triangulations of an almost perfect cube into 7 tetrahedra proposed by[START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF] and their dual complex representation.

Figure 6 . 7 :

 67 Figure 6.7: Two hexahedra not identified by existing combination methods. (1) A decomposition with an interior vertex v.(2) A decomposition into eight tetrahedra. This is a counter example to[START_REF] Sokolov | Hexahedral-Dominant Meshing[END_REF])'s claim that there is no decomposition of the hexahedron into more than seven interior tetrahedra.

Algorithm 3 Figure 6 . 9 :

 369 Figure 6.9: For 4 points in R 2 , only 1 over 3 possible combinatorial quadrilateral is valid. Note that in R 3 , all 3 quadrilaterals define valid bilinear quadrilateral facets.

Figure 6 .

 6 Figure 6.10: Search tree for quadraliterals in a triangulation. Left: Input 2D mesh. Right: Each branch reaching a depth of 4 defines a quadrilateral. 5 quadrilaterals are identified: {1243}, {1235}, {1345}, {1245}, {2354}.

•

  The quality of a hexahedron corner abde is evaluated as the scaled Jacobian:|( ab × ad) • ae|ab ad ae • The quality of a prism corner abcd is evaluated as: 2 ( ab × ac) • ad 3 √ 3 ab + ac + bc ab ac bc ad

2.

  If they share 4 vertices, they share a quadrilateral face connecting these 4 vertices; 3. If they share 3 vertices, they share a triangle face connecting these 3 vertices; 4. If they share 2 vertices, they share an edge connecting these 2 vertices;

Figure 6 .

 6 Figure 6.11: Examples of three hexahedral dominant meshes generated by a greedy selection (white: hexahedra, red: tetrahedra, yellow: prisms, black: pyramids). The complete workflow, from the loading of the tetrahedral mesh till the writing of the mixed-cell mesh, typically runs in less than a minute.

Figure 7 . 1 :

 71 Figure 7.1: The geometry and connectivity of the hexahedra we study. A Finite Element hexahedron is the image of a reference cube by a trilinear mapping, its connectivity is the one of the cube.

Figure 7 . 2 :Figure 7 . 3 :

 7273 Figure 7.2: The four types of tetrahedra in a cube I 3 = [0, 1] 3 .

Figure 7 . 4 :

 74 Figure 7.4: The six triangulations of the 3-cube I 3 = [0, 1] 3 up to isomorphism and their dual complexes. Each class corresponds to 2 to 24 equivalent triangulations.

Figure 7 . 6 :

 76 Figure 7.6: Building a cone over a triangulation of a 2-dimensional pentagon into 3 triangles gives a triangulation of the 2-sphere into 8 triangles.

Figure 7 . 7 :Figure 7 . 8 :

 7778 Figure 7.7: The 7 triangulations of a hexahedron boundary up to isomorphism.

  Figure 7.9: The correspondence between a hexahedron triangulations and its decomposition graph for triangulation 8_E.

Algorithm 6

 6 Geometrical realization Input: T combinatorial triangulation Output: A a realization 1: clauses = encodeSAT(T ) Section 7.3.2 2: S = solveSAT(clauses) MiniSAT solver (Sorensson et al. 2005) 3: if S == ∅ then

Figure 7 .

 7 Figure 7.10: The two possible 4-point configurations in R 2 have different chirotopes. The chirotope encodes the orientation of all 4 triangles, only the one whose orientation changes is drawn.

Figure 7 .

 7 Figure 7.11: Six triangulations of the hexahedon into 8 tetrahedra and their geometrical realizations. Colors correspond to tetrahedron indices. Hexahedron edges are grey, while tetrahedron edges are white. The seventh triangulation with 8 tetrahedra, 8_E, is presented on Figure 7.9.

Figure 7 .

 7 Figure 7.12: The 11 combinatorial triangulations of the hexahedron into 15 tetrahedra. Four are realizable with points in convex position. The hexahedra are valid, their Jacobian is strictly positive (left). Three are not realizable (center). Four are realizable but with points in non-convex position (right).

Figure 7 .

 7 Figure 7.13: The link between tetrahedral meshes (a) and hexahedral meshes (b) depends on hexahedron triangulations. (c) One of the 174 triangulations of hexahedron {12345678} has 8 tetrahedra {1258}, {5286}, {1826}, {1246}, 1486}, {2467}, {4867}, {2347} .

Figure 8 . 1 :

 81 Figure 8.1: Two famous examples of boundary quadrangulation for which no hexahedral mesh was known in 1996. Left is Schneiders'pyramid (Schneiders 1996) and right the octogonal Spindle

Figure 8 . 3 :

 83 Figure 8.3: If hexahedrizations of these two quadrangulated spheres exists, then it is possible to construct hexahedrizations for all other quadrangulated surface[START_REF] Erickson | Efficiently hex-meshing things with topology[END_REF]. But this was left to the reader to build by hand as an exercise.

  Figure 8.4: Searching all quadrilateral meshes of a polygon with up to one interior point. The search tree leaves are either valid solutions, or correspond to detected failure points where Algorithm 7 backtracks.

Figure 8 . 5 :

 85 Figure 8.5: Each new element must share a face with the front of boundary faces (red).

Figure 8 . 6 :

 86 Figure 8.6: (left) A null-homologous 1-cycle in a quad-mesh, bounding 3 quadrangles; (right) a cycle which is not null-homologous, because it surrounds a hole within the mesh.

Figure 8 . 7 :

 87 Figure8.7: Time to explore a search tree in parallel on a machine with two AMD EPYC 7551 CPUs (32 cores each, 2 threads per core). Using 64 threads, the speed-up is of 48 for Schneiders' pyramid and 52 for the octagonal spindle.

Figure 8 . 10 :

 810 Figure 8.10: Replacing the cavity with a valid mesh sharing the same boundary still produces an invalid mesh by creating two quadrangles sharing two edges.

Figure 8 .

 8 Figure 8.11: A valid change to the connectivity of the mesh can create a geometrically invalid mesh, fixed by moving the vertices.

Figure 8 .

 8 Figure 8.12: (top) Removal of two hexahedra from Schneiders' pyramid; (bottom) removal of six hexahedra. The initial cavity (left) and the remeshed cavity (right) have the same combinatorial boundary (top: 18 facets; bottom: 30 facets). Colors highlight the correspondence between faces.

Figure 8

 8 Figure 8.13: Comparison of our 44-element mesh of Schneiders' pyramid (left) with the smallest known 36-element solution (right). Both admit two planar symmetries.

Figure 8 .

 8 Figure 8.14: A layer of hexahedra in a 72-element mesh of Schneiders' pyramid (left) is removed (center). After merging the endpoints of each removed edge, a mesh of the octagonal spindle is obtained (right).

Figure 8 .

 8 Figure 8.16: Equivalence between quad flips and hex creation. Adding one hexahedron to the the top quadrangles, correspond to a flip in the quadrangulation resulting in the bottom grey quadrangles. Six patterns are obtained by taking the symmetric counterpart of the two flips on the left.

Figure 8 .

 8 Figure 8.20: Hexahedrizations of the two types of buffer cubes used to mesh arbitrary domains in the algorithm of (Erickson 2014). (top) 37 hexahedra to mesh the 20-quadrangle cell; (bottom) 40 hexahedra to mesh the 22-quadrangle cell. Colors correspond to the different sides of the original cubes (shown on the left).

Figure 8 . 21 :

 821 Figure 8.21: Time to compute hexahedrizations for the quadrangulations of the sphere with up to 20 quadrangles. Run on a machine with two AMD EPYC 7551 CPUs (32 cores per CPU).

Figure 8 . 22 :

 822 Figure 8.22: Sizes of the smallest hexahedrizations found for the quadrangulations of the sphere with up to 20 quadrangles.

Figure 8 . 23 :

 823 Figure 8.23: The four most time-consuming quadrangulated spheres to mesh using our method. Each required over an hour of computation time on 64 cores.

Table 2 .

 2

1: Number of entities of the open models

Table 3 . 1 :

 31 Comparison of seismic simulation on reference model and on two simplified models. Errors on Vx and Vz are in m/s.

Table 3 .

 3 2: Main features and challenges of the 12 remeshed models.

  3.22: Remeshing the Clyde model with 30000 and 10000 sampling points

	Dataset	#Pts	#ver	#tri	#S	#L	#C	Angles (deg)	Avg qual.	Timing (s)
								Min	< 30 (%)		Sampling	Mesh
	Veins	Input	471337	923286	29	0	0	0.01	22.07	0.55
		1000	19928	35484	29	0	0	1.04	3.46	0.81	153	33
	Forward	Input	10151	13588	46	76	40	0.01	20.21	0.53
		1000	2457	3400	47	71	35	0.56	3.25	0.82	12	1
		5000	7961	12849	46	74	36	1.85	1.00	0.88	28	1
		10000	14096	23845	46	77	41	1.30	0.58	0.89	46	2
	Detachment	Input	61480	109098	50	84	46	7.86	2.19	0.79
		15000	25165	44599	50	79	41	2.60	0.50	0.90	61	3
		30000	42219	76514	50	82	44	3.85	0.24	0.91	85	3
	Leipzig	Input	9286	11344	188	320	166	0.98	7.54	0.73
		5000	8578	11694	186	287	135	4.05	1.61	0.84	13	1
		10000	14911	22281	186	295	143	6.08	0.96	0.86	23	2
	Lambda	Input	24528	37553	132	256	177	0.09	15.52	0.62
		1000	3416	3711	144	236	147	0.98	8.35	0.73	20	2
		10000	16113	24223	134	242	155	0.57	1.46	0.86	33	3
	DFN	Input	7876	7723	435	307	481	0.00	40.70	0.33
		30000	38081	62070	435	300	480	1.00	0.88	0.88	48	3
	HC	Input	39919	70684	80	140	80	0.12	19.54	0.60
		30000	37255	65198	80	141	81	0.78	0.25	0.91	72	5
	Cloudspin	Input	18313	30049	97	124	112	0.00	25.59	0.52
		5000	10778	16494	91	117	103	0.45	3.52	0.82	21	3
		10000	17725	28725	94	134	124	0.30	2.25	0.85	29	3
	Clyde	Input	41355	69343	227	387	303	0.01	20.39	0.56
		10000	15551	23367	206	318	244	0.72	2.51	0.85	38	4
		30000	38884	64850	220	354	282	0.15	1.30	0.88	66	9
	Nancy	Input	59115	85775	753	1307	774	0.00	25.98	0.50
		10000	24840	30445	719	1096	626	0.13	6.88	0.75	43	21
		50000	79087	119309	741	1259	774	0.05	2.53	0.84	106	49
	Annot	Input	76204	130403	332	590	300	0.00	20.61	0.56
		3000	12737	18253	301	455	212	0.51	5.69	0.77	40	7
		20000	41761	68240	311	522	264	0.55	2.31	0.85	80	10
	Sandbox	Input	72927	109267	500	688	713	0.00	21.98	0.53
		30000	52498	77688	503	890	897	0.00	1.11	0.85	93	28

Table 3

 3

.3: Remeshing result statistics. For each model, input model and produced results are compared in terms of mesh sizes, numbers of components (surfaces, lines, and corners), and quality (minimum angle, percentage of triangles with an angle inferior to 30 degrees, average quality)

Table 4 .

 4 Geogram TetGen CGAL

	Sequential_Delaunay	12.7	34.6	32.9	33.8
	Init + Sort	0.5	4.2	2.1	1.3
	Incremental insertion	12.2	30.4	30.8	32.5
	Walk	1.0	2.1	1.6	1.4
	orient3d	0.7	1.4	1.1	≈ 0.5
	Cavity	6.2	11.4	≈ 10	14.9
	inSphere	3.2	6.2	5.6	10.5
	DelaunayBall	4.5	12.4	≈ 15	15.3
	Computing sub-determinants	1.3	/	/	/
	Other operations	0.5	4.5	≈ 4	≈ 1

1: Timings for the different steps of the Delaunay incremental insertion (Algorithm 1) for four implementations: HXT (Marot, Pellerin, and J. Remacle 2018), Geogram

  .Remacle et al. presented an interesting strategy that checks if insertions can be done in parallel by synchronizing threads with barriers (J.-F.[START_REF] Remacle | A two-level multithreaded Delaunay kernel[END_REF]. Vertices are partitioned such that each vertex belongs to a single thread. A triangle can only be modified by a thread that owns all of its three vertices. Triangles that cannot be modified by any thread form a buffer zone.

				mesh vertices
				points to insert
	t 1	bu er zone	2 t	con icts
	Figure 4.5:

Table 4

 4

.2:

  mb 43.8 mb 404.8 mb 3.8 gb Geogram 6.7 mb 30.5 mb 268.6 mb 2.7 gb CGAL 14.1 mb 66.7 mb 578.8 mb 5.7 gb Table 4.3: Comparison of the maximum memory usage of our parallel implementation, CGAL (The CGAL Project 2016) and Geogram (Levy 2016) when using 8 threads. low memory footprint, we were able to compute the tetrahedralization of N = 10 9 vertices (about 6 billion of tetrahedra) on an AMD EPYC machine that has 512 GB of RAM. The experimental memory usage shown in Table 4.3 differ slightly from the theoretical formula because (i) more memory is allocated than what is used (ii) measurements represent the maximum memory used by the whole program, including the stack, the text and data segment etc.

	# vertices	10 4	10 5	10 6	10 7
	Ours	6.9
					Thanks to this

  When there are less deleted tetrahedra than there are tetrahedra in the Delaunay ball, 8192 new "deleted tetrahedra" indices are reserved by the thread. As the mesh data structure is shared by all threads, mesh->tetrahedra.num must be increased in one single atomic operation. Memory allocation for new tetrahedra is synchronized with OpenMP barriers.

	if(cavity -> to_create .num > cavity -> deleted .num)
	{
			uint64_t nTetNeeded = MAX (8192 , cavity -> to_create .num) -cavity ->
				deleted .num;
			uint64_t nTet;
			# pragma omp atomic capture
			{
				nTet = mesh -> tetrahedra .num;
				mesh -> tetrahedra .num += nTetNeeded ;
			}
			reallocTetrahedraIfNeeded (mesh);
			reallocDeletedIfNeeded (state , cavity -> deleted .num + nTetNeeded ); Hxt
			10	Geogram 1.5.4 (8 threads)
			for ( uint64_t i=0; i< nTetNeeded ; i++){ CGAL 4.12
				cavity -> deleted . tetrahedra_ID [cavity -> deleted .num+i] = 4*( nTet+i);
	}	Time [s]	mesh -> tetrahedra .color[nTet+i] = DELETED_COLOR ; cavity -> deleted .num += nTetNeeded ; } 1 # vertices 10 4 Ours 0.032 Geogram 0.041 CGAL 0.037	10 5 0.13 0.19 0.24	10 6 0.85 1.73 2.20	10 7 7.40 17.11 23.37
	0.1 Listing 4.3: void reallocTetrahedraIfNeeded ( mesh_t * mesh) { if(mesh -> tetrahedra .num > mesh -> tetrahedra . allocated_num ) { # pragma omp barrier // all threads are blocked except the one doing the reallocation # pragma omp single { uint64_t nTet = mesh -> tetrahedra .num; alignedRealloc (&mesh -> tetrahedra . neighbor_ID , nTet *8* sizeof ( uint64_t )); alignedRealloc (&mesh -> tetrahedra .vertex_ID , nTet *8* sizeof ( uint32_t )); alignedRealloc (&mesh -> tetrahedra .color , nTet *2* sizeof ( uint16_t )); mesh -> tetrahedra . allocated_num = 2* nTet; } // implicit OpenMP barrier here } } 10 5 10 6 10 7 Number of points (random uniform distribution) (a) 4-core Intel Core TM i7-6700HQ CPU. 10 4 10 5 10 6 10 7 10 8 0.1 10 100 Number of points (random uniform distribution) Time [s] Hxt Geogram 1.5.4 (64 threads) CGAL 4.12 # vertices 10 4 10 5 10 6 10 7 10 8 Ours 0.11 0.43 1.17 4.48 28.95 Geogram 0.10 0.54 4.58 43.70 / Listing 4.4: 10 4 1 CGAL 0.27 0.48 2.44 20.15 /

(b) 64-core Intel Xeon Phi TM 7210 CPU.

Table 6 . 1 :

 61 Characteristics of the input tetrahedral meshes on which tests are performed. For each model are given the number of vertices, the number of tetrahedra, the number of triangles defining the model boundary, the sequential timings to load the mesh and build the data structures used by our algorithm in seconds. The input meshes are available at: www.hextreme.eu.

			Cells			Timings (s)			Cells		Timings (s)
	Q min	#Hexes	#Prisms	#Pyr.	Hex	Pri.	Pyr.	Q min	#Hexes	#Prisms	#Pyr.	Hex	Pri.	Pyr.
	Cube							Fusee 2
	0	710	1,596	1,172 0.03 0.01 0.01	0	4,586,779	5,866,663 3,122,057	9.46	6.88 1.21
	0.2	349	1,052	661 <0.01 <0.01 <0.01	0.2	3,147,148	4,614,298 2,013,799	6.17	5.37 0.92
	0.4	308	1,010	639 <0.01 <0.01 <0.01	0.4	1,940,213	3,978,611 1,746,731	4.13	4.79 0.85
	0.6	129	661	615 <0.01 <0.01 <0.01	0.6	455,155	1,849,333 1,489,604	1.27	2.44 0.78
	0.8	64	218	138 <0.01 <0.01 <0.01	0.8	114,224	520,350	333,250	0.40	0.86 0.46
	Fusee							FT47_b
	0	149,731	251,131	158,671 0.50 0.37 0.07	0	8,374,128	9,845,963 4,954,902 16.34 11.28 2.06
	0.2	95,259	193,440	120,496 0.30 0.28 0.06	0.2	6,111,946	8,046,676 3,177,657 11.11	9.11 1.45
	0.4	50,931	148,390	99,949 0.19 0.23 0.05	0.4	3,413,741	6,963,922 2,798,967	6.88	8.09 1.34
	0.6	14,978	73,397	74,831 0.07 0.13 0.04	0.6	869,890	3,257,991 2,349,396	2.33	4.13 1.22
	0.8	4,187	20,323	18,312 0.02 0.05 0.03	0.8	184,921	808,736	745,059	0.65	1.40 0.78
	CrankShaft						FT47
	0	309,938	516,217	327,071 1.06 0.78 0.14	0	13,842,934 15,994,194 8,098,013 26.79 19.32 3.42
	0.2	196,441	400,698	255,842 0.65 0.60 0.12	0.2 10,225,710 13,177,806 4,968,253 18.28 14.90 2.33
	0.4	97,194	291,381	205,388 0.38 0.47 0.11	0.4	5,978,752 11,638,277 4,441,345 11.65 13.54 2.18
	0.6	27,549	138,308	143,964 0.15 0.26 0.09	0.6	1,481,030	5,392,542 3,886,408	3.76	6.76 1.99
	0.8	6,381	35,157	33,804 0.05 0.10 0.06	0.8	325,062	1,428,695 1,080,358	1.05	2.32 1.24
	Fusee 1						Fusee 3
	0	1,692,873 2,323,521 1,283,513 3.90 2.85 0.52	0	17,052,534 20,251,770 10,377,769 33.06 24.08 4.43
	0.2 1,098,993 1,792,440	866,043 2.40 2.17 0.40	0.2 12,581,618 16,270,311 6,218,061 22.75 18.75 2.94
	0.4	655,375 1,496,380	739,789 1.58 1.88 0.37	0.4	8,214,177 14,639,094 5,562,577 15.55 16.63 2.75
	0.6	167,752	714,225	610,352 0.53 0.99 0.33	0.6	1,770,306	6,701,699 4,969,021	4.41	8.39 2.59
	0.8	46,215	207,852	139,521 0.17 0.36 0.21	0.8	414,528	1,859,542 1,047,155	1.34	2.92 1.54
	Caliper						Los1
	0	3,536,954 4,675,695 2,508,509 7.96 5.73 1.03	0	21

Table 6 .

 6 2: Number of valid cells (hexahedra, prisms and pyramids) identified in 12 tetrahedral meshes (Table6.1) with our algorithm for different minimal quality values. Running times are given for a laptop with 16Go RAM and an Intel Core TM i7-6700HQ CPU @2.60 GHz processor. Input tetrahedral meshes are available at: www.hextreme.eu

	19.37

Table 6 .

 6 3: Number of valid potential hexahedra detected by our algorithm and comparison with the number of hexahedra that correspond to patterns used by previous methods;:[1] [START_REF] Meshkat | Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh[END_REF], [2][START_REF] Botella | Indirect unstructured hex-dominant mesh generation using tetrahedra recombination[END_REF] Sokolov, Ray, Untereiner, and Lévy 2017), [3]

Table 7 . 1 :

 71 Numbers of combinatorial triangulations of the hexahedron with 5 to 15 tetrahedra up to isomorphism.

Table 7 .

 7 3: Number of triangulations patterns per number of tetrahedra counted in the Delaunay triangulations of random point sets.

  Algorithm 7 Recursive enumeration of the hexahedral meshes of the interior of ∂H Input: ∂H, the boundary 1: S , a partial solution 2: C = (C 1 , . . . , C 8 ), the sets of candidate vertices for the current hexahedron 3: if boundary of S == ∂H then

	4:	Print solution S
	5: else if |S | = H max then
	6:	Backtrack
	7: else
	8:

  Algorithm 8 Initialize-candidates(S ): Compute the sets candidate vertices Input: S , a set of hexahedra 1: return C = (C 1 , . . . , C 8 ), the sets of candidate vertices for the next hexahedron. 2: (v 1 , v 2 , v 3 , v 4 ) is a quadrangle that is an output mesh facet. 3: for each i ∈ {1, . . . , 4} do

	4:	C i ← {v i }
	5: for each i ∈ {1, . . . , 4} do
	6:

  .8). Any hexahedral mesh of Schneiders' pyramid has at least 18 interior vertices and 17 hexahedra. Any hexahedral mesh of the octagonal spindle has at least 29 interior vertices and 21 hexahedra. The time to prove lower bounds for the number of interior vertices V int and the number of hexahedra H required to mesh a polyhedron increases exponentially. This is due to the exponential size of the search tree explored by the algorithm.

	10 0 10 1 10 2 10 3 Bounding V int for Schneiders' pyramid 20 22 24 26 28 30 32 10 2 10 4 10 6 10 8 10 10 10 -1 10 0 10 1 10 2 10 3 4 6 8 10 12 14 10 2 10 4 10 6 10 8 10 10 V max Time # Branches H max Bounding |H| for Schneiders' pyramid Time # Branches Figure 8.8: Algorithm 9 Cavity selection algorithm Time [s] Time [s] Input: H, the mesh; n, the size of the cavity Output: A cavity C of n elements 1: h ← a random element of H 2: C ← {h} 3: while |C| n do 4: h ← a random element of H \ C sharing a facet with a hexahedron in C 10 -1 10 0 10 1 10 2 10 3 20 22 24 26 28 Bounding V int for the octagonal spindle 30 32 34 10 2 10 4 10 6 10 8 10 10 10 -1 10 0 10 1 10 2 8 10 12 14 16 18 10 2 10 4 10 6 10 8 10 10 Time # Branches V max H max Bounding |H| for the octogonal spindle Time # Branches Theorem 8.2.1. Theorem 8.2.2. 10 -1 5:	# Branches # Branches

Table 8 .

 8

1: Cavity remeshing operations performed by our hex-mesh simplification algorithm on Yamakawa's 88-element mesh of Schneiders' pyramid

(Yamakawa et al. 2010

).

1 :

 1 if H max = |H| then

	2:	return
	3: else if Visited-Symmetric-Counterpart(H) then
	4:	return	Section 8.4.2
	5: else if Q ≈ a cube then
	6:	h ← the hexahedron bounded by Q
	7:	if Is-Compatible(H, h) then
	8:	Output-Solution(H ∪ {h})
	9: for all quad flip F do
	10:	(Q , h) ← Perform-Flip (F, Q)	Section 8.4.1
	11:	if Num-Vertices(H ∪ {h}) ≥ V max then
	12:	continue
	13:	if Is-Compatible(H, h) then
	14:

Search(Q , H ∪ {h}, H max , V max )

Figure 8.15: A shelling of a quadrangulation on the top. Not a shelling on the bottom because a hole appears after inserting the first 4 quadrangles.

Table 8 .

 8 2: Number of combinatorial quadrangulated boundaries that can be shelled with up to H max hexahedra. Timings are given for a single thread on an Intel R Core TM i7-7700HQ CPU.

Figure

8

.19: Insertion of a buffer layer (light gray) to guarantee final mesh validity when using precomputed cavity meshes.

  Table 8.3: Statistics for the geometric meshes computed for the test cases of Figures 8.1 and 8.3

	Template	|Q| |V bnd | |H| |V total |	|E| per valence sc. jacobian 3 4 5 min max
	Tetragonal trapezohedron	8	10	40	52	40 75	4	0.35 0.42
	Schneiders' pyramid	16	18	36	51	32 62	4	0.12 0.49
	Erickson's buffer cell (1)	20	22	37	53	43 49	4	0.31 0.63
	Erickson's buffer cell (2)	22	24	40	55	44 48	9	0.031 0.45

Table 8 .

 8 4: Statistics on the number of hexahedra, number of vertices, and edge valences for the combinatorial meshes we computed for the even quadrangulations of the sphere with up to 20 quadrangles.

http://www.mmgtools.org/

A lock is a synchronization mechanism enforcing that multiple threads do not access a resource at the same time. When a thread cannot acquire a lock, it usually waits.

An 8-bit char would also work (not less or it would create data races) but the color flag is also used to distinguish volumes in our mesh generator described in Section 4.3

It is common to choose a multiple of the page size (usually

bytes) to minimize TLB misses.

Acknowledgments

project, and the R&D managers of my current employer who supported this project.

I would like want to thank all those who helped me

3D Algorithm

Combinatorics The 3D algorithm is direct extension of the 2D algorithm to identify combinations of tetrahedra into hexahedra. We modify an algorithm generating all possible 8-subset of the set of labeled vertices V to generate exactly once all oriented hexahedra (Algorithm 4). The first corner of the hexahedron is built by choosing vertices {a, b, d, e} such that b > a, d > b and e > b. This corner sets the orientation of the hexahedron (Figure 6.3). Orientation can be ignored by setting e > d. The four other vertices are chosen to be greater than a. The output of Algorithm 4 for V = {1, 2, 3, 4, 5, 6, 7, 8} is a set of 1,680 oriented hexahedra. This is consistent with the well-known fact that there are 24 permutations of the labeled vertices that do not modify the orientation, the edges, or the faces of a hexahedron since 8! = 1, 680 × 24. When orientation is ignored, there are 48 such permutations and then 840 different hexahedra. The complexity of Algorithm 4 is of course catastrophic and cannot be used as is for large point sets.

Algorithm 4 Generates exactly once each potential oriented hexahedron in a vertex set V.

Geometry Let us now associate to each labeled vertex of V a point of R 3 . Among the n 8 × 1680 possible hexahedra, only a small subset will be geometrically valid. To restrict the search space for hexahedra, we build a triangulation (tetrahedrization) of the points. The key idea is to consider only hexahedra whose edges are triangulation edges and whose quadrilateral facets are a combination of two triangle facets. Algorithm 5 outputs all the hexahedra that may be built using the edges and facets of the triangulation. Note that the existence of the triangle facets must be checked explicitly, since, in a 3D triangulation, the existence of edges {ab}, {bc}, {bd}, {da} does not guarantee that any of the triangles {abc}, {acd}, {abd}, or {adc} exist in the triangulation. These tests are performed when computing the possible cells with Algorithm 5 in order to skip invalid configurations and accelerate the procedure. We also ensure that two merged triangles belong to the same parts of the input model, or model faces.

To compute all the hexahedra that can be generated by combining tetrahedra of an input mesh T the last step is to compute the tetrahedra of T that are inside the hexahedron. A tetrahedron is inside a cell if it either (i) have its four vertices be vertices of the cell (and the cell is convex) or (ii) have one facet on the boundary and a volume of the same sign than the cell, or (iii) be adjacent, through a facet that is not on the cell boundary, to a tetrahedron that respect (i) or (ii). At that at this step the real cell boundary is not yet determined, since there are two choices to triangulate each quadrilateral facet and 4 triangle facets are considered. We further ensure that all tetrahedra should belong to the same part of the model, otherwise the cell is discarded. Note that the boundary tetrahedra, as defined in Section 6.1, are not considered to be inside the cell and are ignored.

The algorithm to identify prisms is similar. To build pyramids it is however much faster to iterate on all triangular facets of T and to create three pyramids for each of them by changing the apex to be one of the three vertices of the facet.

Bounding the quality of generated cells

To have an efficient detetection of adequate cells by Algorithm 5, it is crucial to discard invalid or bad quality cells as soon as possible. The quality and validity of a cell depend only on the coordinates of its vertices. We recall that we consider that a cell is