Identifying combinations of tetrahedra into hexahedra: a vertex based strategy
Jeanne Pellerin, Amaury Johnen, Kilian Verhetsel, Jean-Francois Remacle
IIdentifying combinations of tetrahedra into hexahedra: a vertex based strategy
Jeanne Pellerin a, ∗ , Amaury Johnen a , Kilian Verhetsel a , Jean-Franc¸ois Remacle a a Universit´e catholique de Louvain, iMMC, Avenue Georges Lemaitre 4, bte L4.05.02, 1348 Louvain-la-Neuve, Belgium
Abstract
Indirect hex-dominant meshing methods rely on the detection of adjacent tetrahedra that may be combined to formhexahedra, prisms and pyramids. In this paper we introduce an algorithm that performs this identification and buildsthe set H of all possible combinations of tetrahedral elements of an input mesh T into hexahedra, prisms, or pyra-mids. All identified cells are valid for engineering analysis. First, all combinations of eight / six / five vertices whoseconnectivity in T matches the connectivity of a hexahedron / prism / pyramid are computed. The subset of tetrahedra of T triangulating each potential cell is then determined. Quality checks allow to early discard poor quality cells and todramatically improve the e ffi ciency of the method. Each potential hexahedron / prism / pyramid is computed only once.Around 3 millions potential hexahedra are computed in 10 seconds on a laptop. We finally demonstrate that the set ofpotential hexes built by our algorithm is significantly larger than those built using predefined patterns of subdivisionof a hexahedron in tetrahedral elements. Keywords:
Finite element; Hex-dominant mesh; Indirect meshing; Triangulation; Prism; Pyramid
1. Introduction
In this paper we propose a new algorithm to identify all the hexahedra that may be built by combining tetrahedraof a given mesh. Hexahedral meshes are considered by most of the finite element practitioners to be superior totetrahedral meshes (see e.g. [1]). Yet, no robust meshing technique is able to process general 3D domains. Andgenerating hexahedral meshes in an automatic manner is still considered as the ultimate goal in mesh generation[2]. Recently, promising techniques producing meshes composed of a majority of hexahedra have been proposed[3, 4, 5, 6, 7]. These methods take advantage of the existence of robust algorithms to generate tetrahedral meshesand combine tetrahedra to produce meshes composed of a majority of hexahedra associated to prisms, pyramids andtetrahedra . The four steps of these indirect hex-dominant meshing methods can be summarized as follows (see alsoFigure 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 [8].3. The set of potential cells H (hexahedra, prisms, pyramids) that can be defined by combining tetrahedra of T isbuilt.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 (cid:48) .To reach the ultimate goal and combine all tetrahedra into hexahedra, i.e. obtain a final full-hexahedral mesh, all stepsare crucial. Previous works primarily focus on the first step of placing the final mesh vertices. In this paper, we focuson the third step. Our input is a tetrahedral mesh T of a given point set V and we output the set of all hexahedra,prisms, and pyramids that may be built by combining tetrahedra of T . ∗ Corresponding author: Universit´e catholique de Louvain, iMMC, Avenue Georges Lemaitre 4, bte L4.05.02, 1348 Louvain-la-Neuve, Belgium.Tel.: +
32 10 47 23 55. E-mail address: [email protected]
Email addresses:
[email protected] (Jeanne Pellerin),
[email protected] (Amaury Johnen),
[email protected] (Kilian Verhetsel),
[email protected] (Jean-Franc¸ois Remacle)
Preprint submitted to Elsevier January 8, 2018 a r X i v : . [ c s . C G ] J a n
2) Tetrahedral mesh (3) All possible combinationsof tets in hexahedra, prisms, pyramids (4) Hex-dominant mesh(1) Point set
Figure 1: Indirect hexahedral dominant meshing principle. (1) v (2) Figure 2: Two hexahedra not identified by existing combination methods. (1) A decomposition with an interior vertex v. (2) A decomposition intoeight tetrahedra. This is a counter example to [7]’s claim that there is no decomposition of the hexahedron into more than seven interior tetrahedra.
Identifying combinations of tetrahedra into pyramids is trivial. There are three possible pyramids for each facetshared by two tetrahedra of T . Identifying combinations of tetrahedra into prisms is a bit more challenging sincethree tetrahedra properly connected should be identified. However, when identifying combinations of tetrahedra intohexahedra, there are at least ten di ff erent subdivisions of a hexahedron into five, six, or seven tetrahedra ( § H . For example, they donot detect the two hexahedra of Figure 2 (more details in § ffi cient. Thealgorithm does not rely on any pattern. An algorithm variation permits to compute all the possible potential prisms.After reviewing the main methods to combine tetrahedra into hexahedra , prisms or pyramids ( § § § § ++ codeimplementing the methods of this paper is open-source and available at https: // / download / .
2. Background
Before giving details on subdivisions of hexahedra , prisms or pyramids into tetrahedra ( § § b cde f gha bcd efa b cd ea bc d Figure 3: 3D cell templates: tetrahedron, pyramid, prism, and hexahedron.
We have to be very clear that all the cells we are considering are finite element cells hexahedra / prisms / pyramidsvalid for finite element simulations. A very important point is that their quadrilateral facets are not planar, but arebilinear surfaces. We require the Jacobian determinant to be strictly positive at any point inside the cell. The followingconventions will be used throughout this paper (see Figure 3).The pyramid { abcde } has: • { a } , { b } , { c } , { d } , { e }• { ab } , { bc } , { cd } , { ad } , { ae } , { be } , { ce } , { de }• { abe } , { bce } , { cde } , { ade }• { abcd } The prism { abcdef } has: • { a } , { b } , { c } , { d } , { e } , { f }• { ab } , { bc } , { ca } , { ad } , { be } , { cf } , { de } , { ef } , { fd }• { abc } , { def }• { abed } , { efcb } , { acfd } The hexahedron { abcdefgh } has: • { a } , { b } , { c } , { d } , { e } , { f } , { g } , { h } , •
12 straight line edges: { ab } , { bc } , { cd } , { ad } , { ae } , { bf } , { cg } , { dh } , { ef } , { fg } , { gh } , { eh }• { abcd } , { efgh } , { abfe } , { dcgh } , { bcgf } , { adhe } ,A triangulation (tetrahedrization) of a hexahedron / prism / pyramid is a triangulation of the vertices of the cell thatrespect the cell boundary, in other words it is a subdivision of the interior of the hexahedron / prism / pyramid into a setof conformal tetrahedra without any additional vertex. The tetrahedra induce a subdivision of each quadrilateral facetinto two triangles by a diagonal boundary edge. We further define the boundary tetrahedra as the tetrahedra connectingthe four vertices of a cell quadrilateral facet (Figure 4.2). In previous works [3, 6, 5, 7], boundary tetrahedra are calledslivers. We do not use that term which refers to a geometrical property (degeneracy) of tetrahedra. The triangulation isdetermined by interior tetrahedra (Figure 4.1). Indeed, the addition or removal of one or several boundary tetrahedradoes not modify the cell. We first review the subdivisions of the pyramid, of the prism and of the 3-cube into tetrahedra [9].3 b cde f gh
1. Interior tetrahedra a b cde f gh
2. Boundary tetrahedra
Figure 4: Interior tetrahedra and boundary tetrahedra of a hexahedron triangulation. Boundary tetrahedra connect the four vertices of the samequadrilateral facet.
Figure 5: The six types of triangulations of the 3-cube and their dual complex representations. In the dual complexes, there is on vertex pertetrahedron , one plain edge linking adjacent tetrahedra and one dashed edge linking tetrahedra incident to the same quadrilateral facet. The interior(red) edge is the cube diameter and corresponds to a cell in the dual complex.
Triangulation of the pyramid.
There are as many subdivisions of the general pyramid as there are subdivisions of itsplanar base [9]. We are considering square pyramids, there are then exactly two triangulations of the pyramid.
Triangulation of the prism.
The ordinary triangular prism is the result of the product of a triangle with an edge:prism( D ) = D × D . A prism( D n ) has exactly n ! triangulations that are all equivalent to one another by a ffi nesymmetries [9]. The prism has then 6 triangulations that are all equivalent. Triangulation of the 3-cube I = [0 , . The 3-cube has exactly 74 triangulations [9] :1. Every triangulation of the 3-cube contains either a regular tetrahedron (i.e. a tetrahedron whose 6 edges are ofequal lengths) or a diameter, i.e. an interior edge joining two opposite vertices (red edges on Figure 5).2. There are 2 triangulations with a regular tetrahedron, symmetric to one another. The triangulations containingan interior edge are completely classified modulo symmetries by their dual complex which can be one of thelast five shown on Figure 5. There are respectively 8, 24, 12, 24, 4 triangulations in each class.A dual complex (Figure 5) is a practical way to visualize the 6 di ff erent possible decompositions (tetrahedrizations)of the 3-cube. In the dual complex, also called dual graph, one vertex corresponds to one tetrahedron and two verticesare connected by an edge if the corresponding tetrahedra are adjacent through a triangular facet. A 2-cell of thedual complex (cycle in the dual graph) corresponds to an interior edge of the tetrahedrization (red on Figure 5). Ina meshing context, these di ff erent possible decompositions of the 3-cube were identified by [3] who enumerate thefeasible dual complex graphs, called RF-graph in their paper. In the RF-graph, additional dashed edges connecttetrahedra that are adjacent to the same quadrilateral facet (Figure 5). Triangulations of the real cube.
Recently, the work of [3] was extended by [6] and [7] who proposed four additionaldecomposition patterns into seven tetrahedra (Figure 6). The hexahedron is split into two prisms by a tetrahedronwithout any facet on the hexahedron boundary and containing two interior edges. For this tetrahedron to have astrictly positive volume, it is su ffi cient to work in finite precision, i.e. move slightly one of its vertices.4 igure 6: The four types of triangulations of an almost perfect cube into 7 tetrahedra proposed by [6] and their dual complex representation. Bounds on the number of tetrahedra.
The Euler characteristic gives a relationship between the number of tetrahedra t in a cell decomposition and its number of interior edges e int . Each triangulation is a 3-ball with Euler characteristic χ =
1, where χ = v − e + f − t . Then v − e int − e bd + f int + f bd − t =
1. Since there are 4 triangular faces per tetrahedronand 2 tetrahedra per interior triangular face, we have 4 t = f int + f bd . Since the number of boundary edges e bd , andboundary triangular facets f bd are fixed for each type of cell, the number of tetrahedra t in a hexahedron , prism orpyramid decomposition (without internal vertices) depends only on the number of interior edges e int . Moreover, thereare at most (cid:16) n (cid:17) − e bd edges in a cell with n vertices we then have trivial bounds on the number of tetrahedra in thetriangulation of the cells. • For the hexahedron: t hex = + e int and 5 ≤ t hex ≤ • For the prism: t prism = + e int and 3 ≤ t prism ≤ • For the pyramid: t prism = + e int and 2 ≤ t pyramid ≤ To compute the set H of potential hexahedra and other cells that may be built by combining the elements of a tetra-hedral mesh T without modifying its connectivity there are two known approaches. [3] propose to find combinationsof tetrahedra into hexahedra by searching the adjacency graph of T for all occurrences of the cube decompositiondual complexes (Figure 5). The problem of matching subgraphs in large sparse graphs is solved using standard datamining algorithms that operate on graphs. The same technique is used by [11, 12, 6] and [7] who consider four de-compositions into seven tetrahedra (Figure 6). The second approach proposed by [4] relies on the vertices and edgesof the tetrahedral mesh T . Local searches are performed into the vertex-edge graph of T using two patterns. Thesevertex connectivity patterns generalize those proposed by [3] and relax partially the dependency on the tetrahedralmesh. This method has been implemented by [5] where a third pattern taking into account configurations with aninterior flat tetrahedron was added. Some other approaches like H-Morph [13] 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 advantageis that it maintains a valid mixed mesh throughout the procedure. Important observations led us to work on improving these existing techniques. First, they do not identify thelargest set H of potential hexahedra. On Figure 2 we gave two valid hexahedra that would neither be found by [3]’smethod nor by [4]’s method. The first is a decomposition that encompasses one internal vertex, a configuration thatmay occur when a Steiner point is added when generating the tetrahedra. The second is a decomposition that has8 interior tetrahedra. It is a counter example to [7]’s claim that there is no hexahedron decomposition with morethan 7 interior tetrahedra. Both decompositions are not identified when searching for hexahedra made of 5, 6, or 75 igure 7: Two hexahedra with di ff erent edges and faces may be defined from the same set of tetrahedra. interior tetrahedra. Neither are they by [4]’s algorithm since none of their constitutive tetrahedra has three facets onthe hexahedron boundary.Second, as mentioned by [4], several hexahedra may be defined using the same decomposition pattern by modi-fying the ordering of the vertices (Figure 7). The hexahedra have di ff erent edges and di ff erent faces while having thesame tetrahedral decomposition. The hexahedron on the left being a perfect cube, the one on the right is undoubtedlyinvalid (zero Jacobian determinant), but were the vertices in a more general position, both could be valid. Third, theexisting methods identify the same hexahedron several times. That number is as high as the number of corner tetra-hedra in the decomposition in [4]’s approach and depends on dual complex symmetries in [3]’s approach. With thoseobservations in mind, we believe that an algorithm that finds all potential hexahedra in a tetrahedral mesh should notbe based on a predefined set of patterns.
3. An algorithm to combine tetrahedra into hexahedra
In this section, we detail our algorithm to detect combinations of tetrahedra into hexahedra or prisms and findall cells that may be generated by combining elements of a given input tetrahedral mesh T . We first explain the 2dversion of the algorithm that combines triangles into quadrilaterals and show its relationship to algorithms generatingpermutations and combinations. The 2D version of the algorithm is built by modifying an algorithm generating the 4-permutations of a n -set.Given a set V of n vertices labeled from 1 to n , let us first compute all possible quadrilaterals that can be defined fromthese vertices. We define a numbered quadrilateral abcd by the order of its 4 vertices and we define a non-orientedquadrilateral abcd by its edges ab , bc , cd and da ignoring orientation. Generating all numbered quadrilaterals thatcan be built from V is a combinatorial problem solved by Algorithm 1 which generates all possible permutations of4 vertices of the n labels. When V = { , , , } , the output is the set of the 24 permutations of 4 values (Table 1).These 24 permutations define 3 di ff erent non-oriented quadrilaterals, each of them corresponding to 8 equivalentpermutations. Adding constraints on the relative order of the vertices of one quadrilateral a < b , a < c , b < d , weobtain Algorithm 2 which fulfill our first objective and compute all possible non-oriented quadrilaterals that can bebuilt from V . The output for V = { , , , } is now the 3 non-oriented quadrilaterals that may be generated from 4vertices (Table 2). Algorithm 1: V . Data: V vertex set Result: Q set of potential quads foreach a in V do foreach b in V, b (cid:44) a do foreach c in V, c (cid:60) { a , b } do foreach d in V, d (cid:60) { a , b , c } do Q ← Q ∪ { a , b , c , d } ; Algorithm 2:
Unique quadrilaterals in V . Data: V vertex set Result: Q set of potential quads foreach a in V do foreach b in V, b > a do foreach c in V, c > a, c (cid:44) b do foreach d in V, d > b, d (cid:44) c do Q ← Q ∪ { a , b , c , d } ; Figure 8: For 4 points in R , only 1 over 3 possible combinatorial quadrilateral is valid. Note that in R , all 3 quadrilaterals define valid bilinearquadrilateral facets. Table 1: Output of Algorithm 1 for V = { , , , } is the set of the24 permutations of 4 values define 3 quadrilaterals. Table 2: Output of Algorithm 2 for V = { , , , } is the set of the3 possible quadrilaterals. Let us now associate to each labeled vertex a point of R . Among the 3 possible quadrilaterals that can be definedfrom 4 points of R , only one is valid, i.e. has non-intersecting edges (Figure 8). We further associate to the set ofvertices V a triangulation T and modify Algorithm 2 such that it generates all quadrilaterals whose edges are edgesof the triangulation (Algorithm 3). 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 thetriangles 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 9b.The advantage of this approach over the classical algorithms pairing adjacent triangles (e.g. [14]) is that it identifiesquadrilaterals which encompass one (or more) vertex. An example is the quadrilateral { } on Figure 9a thatencompasses vertex { } . The other advantages of the algorithm are that it is easy to add geometrical quality tests (edgelengths, angles of the quadrilateral under construction) and that its parallelization is trivial. Its complexity may seemprohibitive but a vertex of a 2d triangulation is connected to an average of 6 other vertices. Algorithm 3:
Vertex based search algorithm of all potential quadrilaterals in a triangulation.
Data: T triangulation of vertex set V Result: Q set of potential quads foreach a in V do foreach b in neighbors( a ) , b > a do foreach d in neighbors( a ) , d > b do foreach c in neighbors( { b , d } ) , c > a do Q ← Q ∪ { abcd } In this section we extend the 2D algorithm in 3D where the goal is to identify combinations of tetrahedra intohexahedra. Similarly to what we did in 2D, we can modify an algorithm generating all possible 8-subset of the set oflabeled vertices V to generate exactly once all oriented hexahedra (Algorithm 4). The first corner of the hexahedronis built by choosing vertices { a, b, d, e } such that b > a , d > b and e > b . This corner sets the orientation of thehexahedron (Figure 3). Orientation can be ignored by setting e > d . The four other vertices are chosen to be greaterthan a . The output of Algorithm 4 for V = { , , , , , , , } is a set of 1,680 oriented hexahedra. This is consistentwith the well-known fact that there are 24 permutations of the labeled vertices that do not modify the orientation, the7 (a) ∅ (b) Figure 9: Search tree of Algorithm 3 on a example. Left: Input 2D mesh. Right: Each branch reaching a depth of 4 defines a quadrilateral. Fivequadrilaterals are identified: { } , { } , { } , { } , { } . Algorithm 4:
Generates exactly once each potential oriented hexahedron in a vertex set V . Data: V vertex set Result: H set of potential quads foreach a in V do foreach b in V, b > a do foreach d in V, d > b do foreach e in V, e > b do foreach c in V, c > a , e (cid:60) { b , d , e } do foreach f in V, f > a , f (cid:60) { b , d , e , c } do foreach h in V, h > a , h (cid:60) { b , d , e , c , f } do foreach g in V, g > a , g (cid:60) { b , d , e , c , f , h } do H ← H ∪ { abcde f gh } ; edges, or the faces of a hexahedron since 8! = , ×
24. When orientation is ignored, there are 48 such permutationsand then 840 di ff erent hexahedra.The complexity of Algorithm 4 is catastrophic, but it may be nonetheless useful for relatively small point sets. Forexample, we managed to recompute the full-hex mesh of the Schneider’s pyramid subdivision of [15] from its 104points taking into account strong quality constraints on the hexahedra generated. Our algorithm.
Let us now associate to each labeled vertex of V a point of R . Among the (cid:16) n (cid:17) × T .The algorithm to identify prisms is built with the same principles. To build pyramids it is however much faster toiterate on all triangular facets of T and to create three pyramids for each of them by changing the apex to be one ofthe three vertices of the facet. In addition to the vertices and edges of the hexahedron / prism / pyramid, the cell boundary facets and the tetrahedrameshing its interior should be computed to fully define the cell. For the decomposition to be valid, each quadrilateral8 lgorithm 5: Vertex based search algorithm of the set of all potential hexahedra H in a tetrahedral mesh T . Data: V vertex set, T tetrahedrization of V Result: H set of potential hexahedra foreach a in V do foreach b in neighbors( a ) , b > a do foreach d in neighbors( a ) , d > b do foreach e in neighbors( a ) , e > b , e (cid:44) d do foreach c in neighbors( { b , d } ) , c > a , c (cid:44) e do if ! is quad face( a , b , c , d ) then continue; foreach f in neighbors( { b , e } ) , f > a , f (cid:60) { b , d , e , c } do if ! is quad face( a , b , f , e ) then continue; foreach h in neighbors( { d , e } ) , h > a , h (cid:60) { b , c , f } do if ! is quad face( a , d , h , e ) then continue; foreach g ∈ neighbors( { c , f , h } ) , g > a , g (cid:60) { b , d , e } do if ! is quad face( d , c , g , h ) then continue; if ! is quad face( e , f , g , h ) then continue; if ! is quad face( b , c , g , f ) then continue; HexahedronTets = compute tets( { abcde f gh } ) ; H ← H ∪ { abcde f gh , HexahedronT ets } ; facet should be subdivided into two triangles that are facets of the input tetrahedral mesh. The existence of thesetriangle faces must be checked explicitly. Indeed, 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 } do exist in the triangulation. These tests areperformed when computing the possible cells with Algorithm 5 in order to skip invalid configurations and acceleratethe procedure. It is then guaranteed, that the boundaries of the cells defined by each set of ordered vertices output bythe combination algorithm correspond to a set of triangular facets of the input tetrahedral mesh. We also ensure thattwo merged triangles belong to the same parts of the input model, or model faces.For each cell, the last step is to determine the interior tetrahedra. Starting from a tetrahedron that is inside the cellwe propagate to the adjacent tetrahedra and determine if they are inside the cell too. For a tetrahedron to be insidethe cell, it should either (i) have its four vertices be vertices of the cell, or (ii) have one facet on the boundary and avolume of the same sign than the cell, or (iii) be adjacent, through a facet that is not on the theoretical cell boundary,to a tetrahedron that respect (i) or (ii). The di ffi culty is that at this step the real cell boundary is not yet determinedsince there are two choices to triangulate each quadrilateral facet. The four triangle facets are then part of what wepreviously called the theoretical boundary. All the tetrahedra should belong to the same part of the model, when theydo not the cell is discarded. Note that the boundary tetrahedra, as defined in §
2, are not considered to be inside thecell and are ignored. ffi ciency and flexibility of the algorithm To improve the e ffi ciency of Algorithm 5 or of its prism variation, it is crucial to discard invalid or bad qualitycells as soon as possible. The quality and validity of a cell depend only on the coordinates of its vertices. We recallthat we consider that a cell is valid if the Jacobian determinant is strictly positive at any point of the element. All cellsthat 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 overthe element. If this value is inferior or equal to zero, the cell is invalid. In the first-order finite element cells we areconsidering (hexahedra, prisms, and pyramids) the maximum quality of the element is bounded by the quality at thecorners, 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 ( (cid:126) ab , (cid:126) ad )9 • The quality of a triangle facet corner abc [16] is evaluated as:2 (cid:107) (cid:126) ab × (cid:126) ac (cid:107) √ (cid:107) (cid:126) ab (cid:107) + (cid:107) (cid:126) ac (cid:107) + (cid:107) (cid:126) bc (cid:107)(cid:107) (cid:126) ab (cid:107) (cid:107) (cid:126) ac (cid:107) (cid:107) (cid:126) bc (cid:107)• The quality of a hexahedron corner abde is evaluated as the scaled Jacobian: | ( (cid:126) ab × (cid:126) ad ) · (cid:126) ae |(cid:107) (cid:126) ab (cid:107) (cid:107) (cid:126) ad (cid:107) (cid:107) (cid:126) ae (cid:107)• The quality of a prism corner abcd is evaluated as:2 ( (cid:126) ab × (cid:126) ac ) · (cid:126) ad √ (cid:107) (cid:126) ab (cid:107) + (cid:107) (cid:126) ac (cid:107) + (cid:107) (cid:126) bc (cid:107)(cid:107) (cid:126) ab (cid:107) (cid:107) (cid:126) ac (cid:107) (cid:107) (cid:126) bc (cid:107) (cid:107) (cid:126) ad (cid:107)• The quality of a pyramid base corner abde is evaluated as:3 det ( J ) (cid:107) J (cid:107) F where J = J p J − I with J I = √ J p = (cid:18) (cid:126) ab (cid:126) ad (cid:126) ae (cid:19) where (cid:107) . (cid:107) F denotes the Frobenius norm.These quality values give an upper bound of the overall cell quality that we use to accelerate dramatically thecell identification of Algorithm 5. Indeed a new upper quality bound for the cell under construction can be computedwhen a vertex completing a face or a corner is added. When this bound becomes smaller than the required minimumquality, the cell construction is terminated. Additional quality tests on the planarity of quadrilateral facets, or on edgelengths could be added.To guarantee that the Jacobian determinant is strictly positive, a lower bound of the quality is needed. Thatcomputation is more challenging and time consuming and is done when the upper bound test passes. We implementedthe mathematically exact tests proposed by [17, 18]. We additionally modified it so that it is exact and robust tofloating point errors.
4. Results
We have applied our algorithm to 12 di ff ffi ciencyof the algorithm. The multi-threaded version of our algorithm is very fast with about 300,000 potential hexahedrabuilt per second. The algorithm for prisms generates about 900,000 prisms per second and the one for pyramids aboutmore than 1.5 millions per second. All timings and performances are given for a laptop with 16Go RAM and anIntel (cid:114) Core TM i7-6700HQ CPU @2.60 GHz processor. For all cells, the running time clearly depends almost only onthe number of potential cells detected. Note that for models Knuckle and Los1 the computed potential cells are onlycounted since they do not fit in the 16 Go RAM. 10e estimate that our method is more than ten times faster than the state of the art pattern matching method toidentify combinations of tetrahedra into cells [7]. However, providing a valuable comparison is delicate becausequality criteria have a strong impact on performances and are not explicitly stated in previous works, preventingcomparison. Moreover no timings are provided for this one step of the hex-dominant meshing workflow. We compare the number of potential hexahedra identified by our algorithm with the number of potential hexahedrathat would be identified by pattern-matching methods [6, 7] or vertex combination [4] in Table 5. To count the potentialhexahedra matching one of the six cube decompositions or one of the four decompositions into seven tetrahedra wecompare the dual complex graphs. To count the potential hexahedra that would be detected by [4], we count thosecontaining a tetrahedron that has three facets on the hexahedron facets. On small models, our algorithm detects4 to 5% more potential hexahedra than the existing methods. That number does depend on the input tetrahedralmesh.On larger meshes, the small di ff erence between methods can be explained by the point placement strategy of theinput tetrahedral mesh. The points are generated by propagation from the boundary of the model. Where the frontscollide, a roughly 2-dimensional surface, point placement is not optimal. It is in this area that out method makes adi ff erence, and the bigger the mesh is, the relatively smaller this area. Note that the higher the required minimumscaled Jacobian, the smaller the di ff erence between the number of potential hexahedra detected by our method andthe number of hexahedra detected by the existing methods. This is no surprise since the best quality is obtained forhexahedra that are close to the perfect cube which has a limited number of decompositions. To demonstrate that hex-dominant meshes may be generated from the set of potential cells identified by ouralgorithm, 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.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;The incompatibilities between the remaining tetrahedra (not selected to build any cell) and the cells should also beaccounted for. There are at least two possible strategies to manage them: [19] propose to raise pyramids on eachnon-conformal quadrilateral face and [6] propose to subdivide the pyramid or hexahedron incident to a non conformalcontact into pyramids and tetrahedra. Both methods insert a new vertex, the apex of the pyramid or a point inside < < able 4: Number of valid cells (hexahedra, prisms and pyramids) identified in 12 tetrahedral meshes (Table 3) with our algorithm for di ff erentminimal quality values. Running times are given for a laptop with 16Go RAM and an Intel (cid:114) Core TM min min < < < < < < < < < < < < able 5: Number of valid potential hexahedra detected by our algorithm and comparison with the number of hexahedra that correspond to patternsused by previous methods. Model Ours [3] [6, 7] [4]Cube 710 672 696 706Fusee 149,627 123,696 142,783 146,022CrankShaft 310,181 251,806 294,182 302,131Fusee 1 1,692,188 1,532,747 1,683,543 1,686,709Caliper 3,536,997 3,175,201 3,513,863 3,522,483CrankShaft 2 4,405,892 3,919,768 4,383,763 4,392,081Fusee 2 4,585,236 4,110,253 4,568,332 4,574,594FT47 b 8,374,930 7,384,695 8,343,306 8,355,231FT47 13,846,837 12,133,218 13,806,781 13,821,507Fusee 3 17,048,021 14,983,008 17,010,173 17,023,768Los1 21,908,307 19,308,740 21,865,212 21,880,467Knuckle 86,553,836 79,544,742 86,346,178 86,433,942 the neighboring cell algorithm. Their major drawback is that they increase the proportion of tetrahedra and pyramidscompared to the proportion of hexahedra. In the meshes produced for this paper, the compatibility condition betweentetrahedra and cells is relaxed. Some quadrilateral faces will be adjacent to one or two triangles. This mesh shouldthen 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 Indepen-dent Set (MWIS) problem [6]. Let us consider the graph G that has one vertex for each of the cell that may be built bycombining tetrahedra of the input mesh T , and one edge linking each pair of compatible cells. The objective is then tofind 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 thecompatibility graph G is usually very dense, it is advantageous to replace it with its complement, the incompatibilitygraph 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 allindependent sets and comparing their total weights [20]. However such an algorithm cannot be expected to terminatein a reasonable amount of time for graphs with a few thousands of vertices, let alone graphs with a few millions ofvertices like the one we obtain. Reviewing the methods to solve the relevant MWIS problem is out of the scope ofthis paper. The interested reader is referred to [21] and references therein for applicable methods in the specific caseof indirect hex-dominant meshing. Greedy solution.
The strategy we develop to obtain a hex-dominant mesh is the one used by previous works: greedilycompute an approximate solution to the MWIS problem [5, 6, 7]. The vertices (potential cells) are sorted by decreasingweight (quality), and independent vertices (compatible cells) are iteratively added to the solution in decreasing orderof weight. This solution could be improved by taking advantage of the locality of the problem and optimizing smalldisjoint subgraphs [22]. Our non-optimized sequential implementation runs typically in a few seconds. The resultinghex-dominant meshes for three of the input tetrahedral meshes are shown on Figure 10.
5. Conclusion
In this paper we solved one step of the indirect hex-dominant meshing workflow that is the identification in agiven tetrahedral mesh, of all the possible combinations of tetrahedral elements into hexahedra, prisms, pyramids.Our algorithm identifies all the valid cells elements since no assumption is made on subdivisions into tetrahedra, eachcell detected only once, bad quality cells are discarded early, and the algorithm is easy to implement. The C ++ codeis open-source and available at https: // / download / Caliper Q min min min
Figure 10: 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 lessthan a minute.
Acknowledgments
This research is supported by the European Research Council (project HEXTREME, ERC-2015-AdG-694020).
References [1] S. E. Benzley, E. Perry, K. Merkley, B. Clark, G. Sjaardema, A comparison of all hexagonal and all tetrahedral finite element meshes forelastic and elasto-plastic analysis, in: In Proceedings, 4th International Meshing Roundtable, 1995, pp. 179–191.[2] J. F. Shepherd, C. R. Johnson, Hexahedral mesh generation constraints, Engineering with Computers 24 (3) (2008) 195–213.[3] S. Meshkat, D. Talmor, Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh, InternationalJournal for Numerical Methods in Engineering 49 (1-2) (2000) 17–30.[4] S. Yamakawa, K. Shimada, Fully-automated hex-dominant mesh generation with directionality control via packing rectangular solid cells,International journal for numerical methods in engineering 57 (15) (2003) 2099–2129.[5] T. C. Baudouin, J.-F. Remacle, E. Marchandise, F. Henrotte, C. Geuzaine, A frontal approach to hex-dominant mesh generation, AdvancedModeling and Simulation in Engineering Sciences 1 (1) (2014) 1.[6] A. Botella, B. Levy, G. Caumon, Indirect unstructured hex-dominant mesh generation using tetrahedra recombination, Computational Geo-sciences 20 (3) (2016) 437–451.[7] D. Sokolov, N. Ray, L. Untereiner, B. Levy, Hexahedral-Dominant Meshing, ACM Transactions on Graphics 35 (5) (2016) 1–23.[8] H. Si, TetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator, ACM Transactions on Mathematical Software 41 (2) (2015) 1–36.[9] J. A. De Loera, J. Rambau, F. Santos, Triangulations: structures for algorithms and applications, no. v. 25 in Algorithms and computation inmathematics, Springer, Berlin ; New York, 2010, oCLC: ocn646114288.[10] H. Edelsbrunner, F. Preparata, D. West, Tetrahedrizing point sets in three dimensions, Journal of Symbolic Computation 10 (3-4) (1990)335–347.[11] B. Levy, Y. Liu, L p Centroidal Voronoi Tessellation and its applications, ACM Transactions on Graphics 29 (4) (2010) 1.[12] J. Huang, Y. Tong, H. Wei, H. Bao, Boundary aligned smooth 3d cross-frame field, ACM Press, 2011, p. 1.[13] S. J. Owen, S. Saigal, H-Morph: an indirect approach to advancing front hex meshing, International Journal for Numerical Methods inEngineering 49 (1-2) (2000) 289–312.[14] J.-F. Remacle, J. Lambrechts, B. Seny, E. Marchandise, A. Johnen, C. Geuzaine, Blossom-Quad: A non-uniform quadrilateral mesh generatorusing a minimum-cost perfect-matching algorithm, International Journal for Numerical Methods in Engineering 89 (9) (2012) 1102–1119.[15] S. Yamakawa, K. Shimada, 88-Element solution to Schneiders’ pyramid hex-meshing problem, International Journal for Numerical Methodsin Biomedical Engineering 26 (2010) 1700–1712.[16] A. Johnen, C. Geuzaine, T. Toulorge, J.-F. Remacle, E ffi cient computation of the minimum of shape quality measures on curvilinear finiteelements (submitted), Computer-Aided Design.[17] A. Johnen, J.-F. Remacle, C. Geuzaine, Geometrical validity of curvilinear finite elements, Journal of Computational Physics 233 (2013)359–372.[18] A. Johnen, J.-C. Weil, J.-F. Remacle, Robust and E ffi cient Validation of Linear Hexahedral Elements, in: 26th International MeshingRoundtable, Barcelona, Spain, 2017.[19] S. J. Owen, S. A. Canann, S. Saigal, Pyramid elements for maintaining tetrahedra to hexahedra conformability, ASME Applied MechanicsDivision-Publications-AMD 220 (1997) 123–130.[20] Q. Wu, J.-K. Hao, A review on algorithms for maximum clique problems, European Journal of Operational Research 242 (3) (2015) 693–709.[21] K. Verhetsel, Solving the Maximum Weight Independent Set Problem: Application to Indirect Hex-Mesh Generation, Master’s thesis, Uni-versit catholique de Louvain (2017).[22] K. Verhetsel, J. Pellerin, A. Johnen, J.-F. Remacle, Solving the Maximum Weight Independent Set Problem: Application to Indirect Hexahe-dral Mesh Generation, in: 26th International Meshing Roundtable, Research Notes, Barcelona, Spain, 2017.[23] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, International Journalfor Numerical Methods in Engineering 79 (11) (2009) 1309–1331.cient Validation of Linear Hexahedral Elements, in: 26th International MeshingRoundtable, Barcelona, Spain, 2017.[19] S. J. Owen, S. A. Canann, S. Saigal, Pyramid elements for maintaining tetrahedra to hexahedra conformability, ASME Applied MechanicsDivision-Publications-AMD 220 (1997) 123–130.[20] Q. Wu, J.-K. Hao, A review on algorithms for maximum clique problems, European Journal of Operational Research 242 (3) (2015) 693–709.[21] K. Verhetsel, Solving the Maximum Weight Independent Set Problem: Application to Indirect Hex-Mesh Generation, Master’s thesis, Uni-versit catholique de Louvain (2017).[22] K. Verhetsel, J. Pellerin, A. Johnen, J.-F. Remacle, Solving the Maximum Weight Independent Set Problem: Application to Indirect Hexahe-dral Mesh Generation, in: 26th International Meshing Roundtable, Research Notes, Barcelona, Spain, 2017.[23] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, International Journalfor Numerical Methods in Engineering 79 (11) (2009) 1309–1331.