Combinatorial Techniques for Hexahedral Mesh Generation

The main objects discussed throughout this thesis are meshes, volumetric and hexahedral meshes in particular. The properties of meshes are studied through two interrelated branches of mathematics: topology and geometry.

The topological properties of a mesh can be derived from its combinatorial structure, which we call a topological or combinatorial mesh (Figure 4). Given a finite set of vertices, we represent surface meshes as a set of polygons, only storing the incidence relations between them, e.g. by labeling the vertices and representing polygons as ordered lists of vertex labels. Similarly, we represent volumetric meshes by a set of cells (also known as elements). Each cell is a combinatorial polyhedron, bounded by polygonal facets. Such cells are encoded, for example, by picking a convention to order their vertices (Figure 5). The intersection of any two cells must be exactly one face of both cells: the empty set, a vertex, an edge, a polygonal facet or the cell itself.

This combinatorial object corresponds to a topological space \(X\) using the following construction: each k-dimensional face is represented by a copy of the unit ball \(\left\{\mathbf{x} \in \mathbb{R}^k \mid ||\mathbf{x}||^2 \le 1\right\}\). The k-skeleton \(X^{(k)}\) is the disjoint union of these balls. The boundary of each k-dimensional face (for \(k > 1\)) is the union of balls of dimension \((k - 1)\), corresponding to the boundary facets, and attached to the skeleton \(X^{(k-1)}\) using a gluing map from the unit sphere \(\left\{\mathbf{x} \in \mathbb{R}^k \mid ||\mathbf{x}||^2 = 1\right\}\) to \(X^{(k-1)}\). The topological space \(X\) is the union of its skeletons \(X = \bigcup_{k = 0}^{n} X^{(k)}\). Cell complexes are used extensively in algebraic and computational topology [Munkres 1984; Edelsbrunner 2001; Hatcher 2002; Kaczynski et al. 2003].

A geometric mesh is one where each cell is a convex polyhedron, defined as the convex hull of its vertices \(\mathbf{x}_1, \dots, \mathbf{x}_n\):

\[\left\{\sum_{i = 1}^{n} \lambda_i\mathbf{x}_i \;\middle|\; \sum_{i = 1}^{n} \lambda_i = 1 \mbox{ and } \forall i. 0 \le \lambda_i \le 1 \right\}\]

Equivalently, convex polyhedra may be defined as the intersection of half-spaces delimited by the planes that contain their facets, i.e. they must have planar faces. An extensive discussion on the properties of convex polyhedra and higher-dimensional polytopes is found in [Ziegler 1995].

A geometric realization of a topological mesh assigns coordinates to each vertex such that each cell is a convex polyhedron. A volumetric mesh is realizable if it has a geometric realization in \(\mathbb{R}^3\).

For finite elements applications, a less restrictive type of geometric mesh is used: trilinear hexahedral meshes [Knupp 1990]. Given coordinates in \(\mathbb{R}^3\) for each vertex, a hexahedron whose vertices are \((\mathbf{x}_{000}, \mathbf{x}_{001}, \mathbf{x}_{011}, \mathbf{x}_{010}, \mathbf{x}_{100}, \mathbf{x}_{101}, \mathbf{x}_{111}, \mathbf{x}_{110})\) is defined as the image of a unit cube under a trilinear map \(f : [0, 1]^3 \rightarrow \mathbb{R}^3\), obtained by taking the products of linear functions (denoted \(F_0\) and \(F_1\)), weighted by the vertices:

\[f(u, v, w) = \sum_{i = 0}^{1} \sum_{j = 0}^{1} \sum_{k = 0}^{1} F_i(u)F_j(v)F_{k}(w)\mathbf{x}_{ijk}\]

where

\[\begin{align*}F_0(t) &= 1 - t \\ F_1(t) &= t \end{align*}\]

A trilinear mesh is valid if every point of the physical domain is the image of precisely one point under the trilinear maps that define the mesh. This corresponds to the Jacobian determinant of \(f\) being strictly positive everywhere for every hexahedron (assuming a convention on the orientation of hexahedral elements).

Cell complexes can describe a wide range of spaces. For meshing purposes, the domains that we seek to mesh are usually manifolds, or manifolds with boundary. An n-manifold is a space where each vertex \(v\) has a neighborhood homeomorphic to an open subset of \(\mathbb{R}^n\).

In combinatorial terms, manifolds are characterized in terms of the links and stars of faces in a cell complex. The star of a face \(f\), \(\mathrm{St}(f)\) is the set of all cells that contain \(f\), and all of their faces. The link of \(f\), \(\mathrm{Lk}(f)\), is the set of faces in \(\mathrm{St}(f)\) that do not contain \(f\). A combinatorial n-manifold with boundary (for \(n \le 3\)) is a topological mesh such that the link of every face of dimension \(k\) is homeomorphic to a sphere \(S^{n-k}\) (for interior faces) or a ball \(B^{n-k}\) (for boundary faces) [Edelsbrunner 2001; Rourke and Sanderson 2012]. In particular, the link of each vertex is homeomorphic to \(S^{n-1}\) for interior vertices or \(B^{n-1}\) for boundary vertices (Figure 6).

The faces of a topological mesh are partially ordered by a relation \(a \preceq b \Leftrightarrow a \mbox{ is a face of } b\). Reinterpreting a mesh by inverting this relation yields its dual mesh [Basak 2010]. Given a surface mesh \(\mathcal{M}\), its dual has one vertex for each polygon in \(\mathcal{M}\), and one dual edge connecting any two dual vertices corresponding to primal polygons that share an edge. Each vertex of the primal mesh \(\mathcal{M}\) then corresponds to a polygonal cell of the new dual mesh. The dual of a volumetric mesh \(\mathcal{M}\) is also a volumetric mesh, each vertex of which corresponds to a polyhedral cell of the original mesh. Dual edges connect the vertices corresponding to primal cells that share a facet and dual facets correspond to primal edges.

For quadrangular and hexahedral mesh, we add structure to the dual mesh by grouping the dual edges or polygonal facets into curves or sheets respectively. In a quadrangular mesh, a dual curve groups all edges corresponding to opposite edges of a combinatorial quadrangle. In a hexahedral mesh, any two dual polygonal facets that correspond to parallel primal edges of a combinatorial cube are grouped together in a dual sheet (Figure 7). The boundary of a dual sheet is either empty, or one or more dual curves of the boundary quadrangular mesh.

Whisker Weaving. Whisker Weaving has been proposed by [Tautges et al. 1996]. The idea is to use a topological advancing front to construct the dual of a hexahedral mesh with a prescribed boundary. 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. [Kremer et al. 2014] extended this algorithm with heuristics to determine the elimination order of the dual cycles, taking into accounts geometric properties to handle concave objects.

Octrees-based methods are a type of automatic hex-meshing algorithm that convert the combinatorial dual of an octree into a hexahedral mesh [Schneiders 2000; Maréchal 2009; Gao et al. 2019]. The octree structure is first refined according to the local feature size. A set of fixed templates is used to transform size-transitions present in the octree into a valid hexahedral mesh. The choice of templates directly impacts the quality of the resulting mesh, so they must be designed carefully to avoid poorly shaped hexahedra [Tong et al. 2024].

A major challenge for this family of algorithms is to accurately capture features that are not aligned with the octree. Given a set of sharp features, octree-based methods can project vertices onto the sharp features so that they are preserved in the final mesh. This often requires a high level of refinement in order to capture small features, limiting the practical applicability of those techniques for complex inputs.

Frame fields assign three mutually orthogonal directions to each point in space. Lines where this orientation field vanishes are known as singularities. In meshing terms, these lines correspond to edges that are surrounded by a number of hexahedra different from four [Beben 2020].

Methods based on frame field start by generating a smooth, boundary-aligned orientation field and attempt to extract a block-structured mesh whose singularities match that of the frame field [Nieser et al. 2011; Li et al. 2012; Liu et al. 2018]. Currently, however, there is no guarantee that the singularities in the generaed frame-fields correspond to that of a valid hexahedral mesh. Invalid configurations are generated for most non-trivial models. For example, there can be vertices where only two singularities meet, one of valence three and the second of valence five — which is not a valid configuration for an all-hexahedral mesh. [Liu and Bommes 2023] formulate the problem so as to avoid such invalid configurations, but there is still no guarantee the global orientation field corresponds to a valid hexahedral mesh.

Polycubes are solids made up of a finite number of cubes glued face-to-face. This structure makes it easy to subdivide them into a regular arrangement of hexahedra. Polycube-based algorithms use such an arrangement to mesh an object by computing a mapping between it and a polycube [Gregson et al. 2011; Livesu et al. 2013]. If the distortion of this polycube mapping is low enough, then a hexahedral mesh of the polycube can be deformed into a valid mesh of the target model (Figure 10).

[Lin et al. 2008] create such polycubes by first segmenting the input model into different parts corresponding to features of the input domain. Each part is then approximated by a coarse polycube chosen from a set of basic primitives. This coarse representation of the model often results in a high distortion mapping, making it unsuitable for meshing applications.

[Gregson et al. 2011] instead deform the input model, so as to align all boundary faces to one of the six axis-aligned directions (\(\pm x\), \(\pm y\), \(\pm z\)). The resulting polycube is then meshed with a structured hexhadedral mesh which is mapped back onto the input shape. Polycut [Livesu et al. 2013] also uses a deformation-based approach, but improves the segmentation of the input into axis-aligned patches using a hill climbing algorithm. This leads to a mapping with less distortion while requiring fewer irregular vertices. To reduce artifacts due to the high sensitivity of parametrization methods to the mesh embedding, [Mandad et al. 2022] propose an intrinsic formulation. [Protais et al. 2022] focus on the robust construction of very coarse meshes from such polycube mappings while still preserving relevant geometric features.

Even though these method allow the construction of high-quality, block-structured hexahedral meshes, the deformation process often results in inverted hexahedra. The regular connectivity of polycube meshes that does not contain singularities also makes it difficult to correctly represent certain boundary features.

For many models, the all-hexahedral methods mentioned so far fall short of producing a satisfactory hexahedral mesh. Hex-dominant are a more robust alternative thanks to the inclusion of non-hexahedral elements, while still filling most of the volume with well-shaped hexahedra.

H-Morph [Owen 2001] computes a hex-dominant with a prescribed quadrangular mesh as its boundary. The input is a tetrahedral mesh constrained so that its boundary is a subdivision of the target quad mesh. The algorithm traverses this tetrahedral mesh in an advancing front fashion, applying local transformations to rearrange tetrahedra until they can be combined into hexahedra.

More recent hex-dominant meshing techniques directly create a tetrahedral mesh with its vertices placed so as to facilitate recombination. Such a distribution of vertices is obtained by packing rectangular cells [Yamakawa and Shimada 2003], using the \(L_p\) Centroidal Voronoi Tesselation [Botella et al. 2016; Lévy and Liu 2010], or frontal point insertion [Baudouin et al. 2014]. The tetrahedra are combined into hexahedra by searching for occurrences of the triangulations of a hexahedron [Pellerin et al. 2018b]. This search may rely on an explicit set of templates [Botella et al. 2016; Yamakawa and Shimada 2003], or by traversing the mesh to look for groups of 12 edges whose connectivity match that of a cube [Pellerin et al. 2018a]. [Ray et al. 2018] use an orientation field to guide the construction of a hex-dominant mesh, extracting the hexahedra using standard methods [Nieser et al. 2011]. A quad-dominant mesh of the boundary is computed ahead of time, and hexahedra that self-intersect or intersect the boundary are excluded. The final hex-dominant mesh is obtained by filling the leftover space using a Delaunay triangulation.

The hex-dominant meshes produced by these methods are in general non-conforming, as they may contain hexahedra that share a quadrangular face with two tetrahedra. Non-conformities are usually resolved as a post-processing step. Alternatively, [Gao et al. 2017b] generate a conforming mesh directly by allowing a small number of cells with arbitrary polyhedra, in addition to hexahedra, tetrahedra, prisms and pyramids.

Each variable is represented by stacking up layers of the gadget shown in Figure 14. Each layer can be combined into exactly three hexahedra optimally in two different ways, one to represent true being assigned to a variable, one to represent it being false. When two layers are stacked on top of one another, the optimal combination corresponds to consistently picking one value for the variable. Creating a hex-dominant mesh that mixes hexahedra corresponding to the two different variable states is only possible by using fewer hexahedra per layer or including low-quality hexahedra (Figure 15). A threshold of 0.6 on the minimum scaled Jacobian is used to preclude the latter. This same threshold will be used for all other gadgets we describe.

Each clause \(x \lor y \lor z\) is represented by 10 tetrahedra that can be combined to form one hexahedron in three different ways (Figure 16). Each of the three hexahedra represents selecting a literal (\(x\), \(y\), or \(z\)) whose value is true.

To ensure the selected literal is consistent with the variable gadget, the clause gadget must be connected to three variable gadgets, so that the only way to optimally combine both the clause gadget and the variable gadget is if the clause is satisfied.

These constraints are enforced by connecting three of the vertices from the last layer in the variable gadget to the clause. These three vertices are all part of a single quadrangular face in the hexahedra corresponding to the negation of the corresponding literal (i.e. the state where \(x\) is false if \(x\) is included in the clause, and the state where \(x\) is true if \(\lnot x\) is included instead). They are merged with three vertices of a single quadrangular face of the hexahedron used to represent the corresponding literal, while still ensuring that this hexahedron remains compatible with the two hexahedra corresponding to the other two literals.

The gadgets used to represent clauses need to be connected to the last layer of a variable gadget, but each variable needs to be connected to multiple clauses. To allow this, we use a third gadget to create a T-junction in the middle of the stack of layers used to represent variable.

The T-junction gadget is constructed by adding three tetrahedra to a stack of two layers of the variable gadget and attaching another layer in a perpendicular orientation (Figure 18). The resulting tetrahedral mesh admits two optimal recombination into 10-hexahedra, corresponding to the two states a variable can take.

Any 3-SAT instance \(\mathcal{C}\) can be represented by a tetrahedral mesh \(T\) using a polynomial number of the gadgets we described (Figures 14,16, and 18): each variable \(x_i\) is encoded by stacking \(O(|\mathcal{C}|)\) layers of the variable gadget. One instance of the clause gadget is used for each clause. For each occurrence of \(x_i\) in a clause, two consecutive layers are replaced with a T-junction, and more layers are inserted to connect \(x_i\) to the clause. \(\mathcal{C}\) is satisfiable if and only if \(T\) can be recombined into a hex-dominant mesh containing \(n\) hexahedra with a minimum scaled Jacobian of at least \(0.6\), where \(n\) is the sum of the number of hexahedra that each gadget can optimally be combined into.

By construction, any solution to the original 3-SAT instance corresponds to such a hex-dominant mesh. To verify that there is no way to combine \(T\) into \(n\) hexahedra of sufficient quality when the formula is not satisfiable, consider the incompatibility graph \(G\). The nodes of \(G\) are the candidate hexahedra that can be obtained by combining tetrahedra from \(T\). Any two nodes \(h_i\) and \(h_j\) of \(G\) are connected if and only if they are incompatible, either because the corresponding hexahedra geometrically overlap or because they only share part of face. The graph \(G\) can be partitioned into one sub-graph for each gadget and one sub-graph for each pair of directly connected gadgets, since no candidate hexahedra span across more than two gadgets.

The maximum number of hexahedra in a hex-dominant mesh obtained from \(T\) is at most the sum of the sizes of the maximum independent sets of each sub-graph, an independent set of a graph being a set of nodes none of which are connected by an arc. This corresponds to ignoring the incompatibility constraints across sub-graphs. To verify that the optimal number of hexahedra can only be achieved if the 3-SAT instance \(\mathcal{C}\) admits a solution, we consider the sub-graphs that correspond to all pairs of directly connected gadget \(G_1\) and \(G_2\). There are only a small number of cases used in the construction: the connection between consecutive layers of the variable gadget, the connection between layers and T-junctions, and the 6 different ways a clause can be connected to a variable gadget — corresponding to the 3 literals of each clause, and whether or not any of each of them appears with a negation.

Let \(G_{12}\) be the sub-graph of \(G\) that contains all candidate hexahedra that cross both gadgets. In all cases, we verify that any independent set that contains a hexahedron \(h \in G_{12}\) is sub-optimal, so that a better local solution can always be obtained by instead recombining \(G_1\) and \(G_2\) optimally. Since the optimal number of hexahedra \(n\) assumes that all gadgets were combined optimally and the optimal hex-dominant mesh for each gadget encodes an assignment for the variables in \(\mathcal{C}\), a hex-dominant mesh containing \(n\) hexahedra can only be obtained if \(\mathcal{C}\) is satisfiable.

This shows our recombination problem is NP-Hard. Because it is in NP (e.g. one can construct the incompatibility graph \(G\) and compute its maximum independent set), it also follows that it is NP-Complete.

Given \(\partial 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 \(\partial 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 we are solving 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].

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. Algorithm 1 corresponds to the exploration of a search tree (Figure 20) 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.

In this section, we describe the key points of our implementation of Algorithm 1, 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, i.e. different representations of the same meshes;
• 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 suitable solutions.

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 \(\partial{}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 21).

Filtering out 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:

• the sets of edges, quadrangle diagonals, and interior diagonal of hexahedra are disjoint;
• no two hexahedra may share an interior diagonal;
• if one facet diagonal matches an existing quadrangle diagonal, so must the second one;
• 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. interior vertices can be relabelled (Figure 22) — for boundary vertices, the algorithm uses the same labels as the input;
2. the hexahedra of the solution can be constructed in a different order (Figure 23);
3. for a given hexahedron, written as an ordered sequence of 8 vertices, there are \(1,680 = 8! / 24\) ways to reorder these vertices while leaving then hexahedron unchanged (Figure 24).

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. 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.

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. This implies topological requirements that 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 [Schleimer 2011].

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 [Eppstein 1999a]. 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, in particular those computed over \(\mathbb{Z}_2\). A detailed introduction to homology computations can be found in [Hatcher 2002]. We define a \(k\text{-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\text{-chain}\) \(C\) is the \((k-1)\text{-chain}\) \(\partial{}C\) whose elements are the faces contained in an odd number of elements of \(C\). A \(k\text{-chain}\) whose boundary is empty is referred to as a cycle. A \(k\text{-cycle}\) which is the boundary of a \((k+1)\text{-chain}\) is called null-homologous (Figure 25).

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.

The utility and correctness of these tests were verified by testing the implementation using different filtering rules, and analyzing the topology of the output meshes using Regina [Burton et al. 1999], a standard piece of software to analyze low-dimensional manifolds.

The exploration of 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 26). We use an approach similar to the embarrassingly parallel search of [Régin et al. 2013]. The main challenge to overcome is that some subtrees are multiple orders of magnitude larger than other ones without any possibility to determine which ones 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. At the start of the search, the tree is explored in a breadth-first manner until a layer with at least \(kT\) subproblems is reached, where \(T\) is the number of threads and \(k\) a constant to control the amount of work per thread (we used \(k = 4096\)). The nodes of this layer are then explored in parallel by independent worker threads using Algorithm 1.

Using Algorithm 1, we computed lower bounds for the number of vertices and hexahedra required to mesh Schneiders' pyramid and the octagonal spindle (Figure 19). 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 27).

The following bounds were proven in this manner:

1. Any hexahedral mesh of Schneiders' pyramid has at least 18 interior vertices and 17 hexahedra.
2. Any hexahedral mesh of the octagonal spindle has at least 29 interior vertices and 21 hexahedra.

The cavity selection algorithm is a greedy algorithm that starts from a random element of the input hexahedral mesh (Algorithm 3). 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 \(\mathcal{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.

To find a smaller mesh of the boundary of a cavity \(\mathcal{C}\), we first solve the combinatorial problem, i.e. we find the smallest combinatorial hexahedral mesh of \(\partial \mathcal{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 \(\partial \mathcal{C}\) is a direct application of Algorithm 1 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 \(|\mathcal{C}|\). Changing the parity of a hexahedral mesh is known to be a difficult operation [Schwartz and Ziegler 2004], so we set \(H_{max}\) to \(|\mathcal{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 \(\mathcal{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 29. 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 4.1.2) are modified to take into account the hexahedra that are not part of the cavity.

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 30). To obtain a valid geometric mesh we use the untangling algorithm described by [Toulorge et al. 2013]. 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 [Johnen et al. 2017].

Given a quadrangulation of the sphere \(Q\), we describe an algorithm to enumerate all hexahedral meshes bounded by \(Q\) and which can be constructed using quad flips (Algorithm 4). To force the algorithm to terminate, the search is limited to meshes with a maximum of \(H_{\max}\) hexahedra and a maximum \(V_{\max}\) of vertices. In Section 5.2, we then extend the approach to search for hexahedral meshes that are prohibitively large for an exhaustive search of this kind.

The algorithm detailed in this section does not search the entire space of hexahedral meshes, but only the space of so-called shellable meshes, which can be explored efficiently using quad flips (Section 5.1.2). This space is explored in its entirety by considering all possible sequences of quad flips that correspond to valid hexahedral meshes (Section 5.1.3). Because many different sequences of flipping operations represent the same mesh, most of this section 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 5.1.4).

Our 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 [Ziegler 1995]. Slightly different notions of shellability are found in the literature. We use that of pseudo-shellings [Bern et al. 2002] or topology-preserving shellings [Müller-Hannemann 1999]. This type of shelling is an ordering of the hexahedra \((H_1, H_2, \dots, H_n)\) of a hexahedral mesh such that any prefix \(\bigcup_{0 \le i < k} H_i\) is homeomorphic to a ball (Figure 31).

This definition implies that any hexahedron \(H_k\) must intersect the union of the previous hexahedra in one of six possible patterns. Gluing a hexahedron to one of these patterns modifies the boundary of the mesh locally (Figure 32). The transitions between these patterns are known as quad flips or bubble moves [Funar 1999]. These flipping operations are therefore a valuable building block to explore the space of shellable meshes.

Note that not all combinatorial meshes admit a shelling order (Figure 33) [Lutz 2003]. 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 [Bern et al. 2002].

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. [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 32. Second, those flips that correspond to the insertion of hexahedra that would make the mesh invalid are filtered out.

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}\), represented as bitsets, and verify that, after adding a new hexahedron:

1. the three sets \(E\), \(D_{Q}\), and \(D_{H}\) remain disjoint;
2. the new quadrangles 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 some other 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.

The insertion of the last hexahedron requires special treatment. This step does not correspond to a quad flip: when the boundary of the unmeshed region is isomorphic to the boundary of a cube, a hexahedron is inserted to finish the mesh. Detecting whether or not the current boundary corresponds to that of the cube is straightforward: simply verify that the boundary has exactly 6 faces.

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, lest most of the computation time be spent generating different representations of equivalent solutions.

One 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 [Burton 2011; Brinkmann and McKay 2007]. A significant portion of the execution time is then spent computing the canonical representations of partial solutions, which may completely change after every operation [Jordan et al. 2018]. 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) [Fahle et al. 2001]. 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. first, the automorphism group of the input quadrangulation is pre-computed;
2. then, as the search tree is traversed, fully explored subtrees are encoded into a sequence \(S\);
3. for each new node, we determine whether or not it should be pruned by comparing it against the nodes stored in \(S\).

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 for two quadrangulations is performed using a variation on Algorithm 5 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 8), we obtain a greedy construction similar to the procedure outlined by [Xiang and Liu 2018]: 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 8, Algorithm 4 is used to verify that there is indeed no shellable hexahedral mesh with at most \(n\) hexahedra. This test was performed for \(n \le 10\), and no counter-examples were found.

From Algorithm 8, a table of \(69,043,690\) boundaries that can be meshed with up to 11 hexahedra is constructed (Table 1). These results also coincide with those obtained by [Xiang and Liu 2018].

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 37). Algorithm 4 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.

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.

The trilinear embedding of the meshes available to us were constructed to maximize quality measures used to evaluate meshes in finite element applications. The faces are non-planar because these measures do not seek to minimize the extent to which faces are warped.

Using SCIP [Bolusani et al. 2024], we compute a numerical mesh ensuring the volume of the tetrahedron spanned by the vertices \(\mathbf{a}\), \(\mathbf{b}\), \(\mathbf{c}\), \(\mathbf{d}\) of each face is within a small tolerance \(\epsilon\) of zero:

\[-\epsilon \le \frac{1}{6} [(\mathbf{b} - \mathbf{a}) \times (\mathbf{c} -\ \mathbf{a})] \cdot (\mathbf{d} - \mathbf{a}) \le \epsilon\]

All boundary vertices have the coordinates fixed using an embedding of the boundary polyhedron. This process informs the choice for a combinatorial mesh to extend to a geometric mesh: whereas this process finds hexahedral mesh with volumes for the faces with \(\epsilon \le 10^{-8}\) from our 44-element mesh, it is not able to find such solutions for the 36-element mesh of the pyramid originally found by [Xiang and Liu 2018]. This suggests it may not be possible to extend the 36-element mesh to a geometric mesh.

From this approximate numerical solution, we construct a new mesh where every face will be planar. We keep track of the set \(F\) of vertices whose coordinates in the geometric mesh have been determined. These exact coordinates are computed by applying a sequence of the following three operations:

1. If there is a hexahedron \(h\) containing a vertex \(v\) such that the 6 vertices that share a quadrangle of \(h\) with \(v\) are in \(F\), compute \(v\) by intersecting the 3 planes that contain it and add it to \(F\).
2. If there is a quadrangle \(q\) containing a vertex \(v\) such that the other 3 vertices of \(q\) are in \(F\), compute \(v\) by projecting it onto the plane of the remaining vertices and add it to \(F\).
3. If none of those are possible, directly use a rational approximation of the coordinates of \(v\) in the numerical solution and add it to \(F\). For this operator, the vertices located on the symmetry planes (i.e. the fixed points of the symmetries) are prioritized, since they have fewer degrees of freedom. In particular the vertices at the intersection of all symmetry planes are of the form \((0, y, 0)\) and only have one degree of freedom.

After adding a vertex \(v\) to \(F\), its symmetric counterparts under each symmetry \(\sigma\) are also added to \(F\), with their coordinates computed from the image of \(v\) under the corresponding geometric symmetry. Due to the high number of planarity constraints across the mesh, the location of all vertices become completely constrained after defining only a small number of them.

Each operation tends to increase the size of the fractions used to represent rational vertices. To keep them from growing too large, making each consecutive operation slower, the coordinates in the numerical solution are initially quantized to the nearest multiple of \(\frac{1}{256}\) (except for the boundary vertices so that the boundary faces remain planar).

We manually find a sequence of operations that yields a mesh with all faces planar by construction, using 10 vertices from the quantized mesh (operation 3), 3 projections of a vertex onto a quadrangle (operation 2) and computing the other vertices with the first operation whenever possible (Figure 50). The resulting mesh may still include flipped hexahedra depending on the initial vertices, despite all faces being planar. Changing the location of one vertex by a small amount in the initial numerical solution and performing these operations again results in the first geometric mesh of Schneiders' pyramid, as well as the octagonal spindle (Figure 51).