Figure 1. Hex-meshes can be categorized according to their structural
regularity, which depends on the amount of irregular edges and vertices
present in the mesh (red dots) and on how they are connected to each other
(thick lines). The first three meshes were computed with a polycube method
(Livesu et al., 2013), the latter with an octree method (Gao et al., 2019).

# Hex-Mesh Generation and Processing: a Survey

Nico Pietroni<EMAIL_ADDRESS>University of Technology Sydney,
Australia  ,  Marcel Campen<EMAIL_ADDRESS>Osnabrück University, Germany  ,
Alla Sheffer<EMAIL_ADDRESS>University of British Columbia, Canada  ,
Gianmarco Cherchi<EMAIL_ADDRESS>University of Cagliari, Italy  ,  David
Bommes<EMAIL_ADDRESS>University of Bern, Switzerland  ,  Xifeng
Gao<EMAIL_ADDRESS>Tencent America, USA  ,  Riccardo Scateni
<EMAIL_ADDRESS>University of Cagliari, Italy  ,  Franck Ledoux
<EMAIL_ADDRESS>CEA, France  ,  Jean-François Remacle  jean-
<EMAIL_ADDRESS>Université catholique de Louvain, Belgique  and
Marco Livesu<EMAIL_ADDRESS>CNR IMATI, Italy

(2022)

###### Abstract.

In this article, we provide a detailed survey of techniques for hexahedral
mesh generation. We cover the whole spectrum of alternative approaches to mesh
generation, as well as post processing algorithms for connectivity editing and
mesh optimization. For each technique, we highlight capabilities and
limitations, also pointing out the associated unsolved challenges. Recent
relaxed approaches, aiming to generate not pure-hex but hex-dominant meshes,
are also discussed. The required background, pertaining to geometrical as well
as combinatorial aspects, is introduced along the way.

hexahedral mesh, integer-grid map, polycube, block decomposition, dual sheets,
frame field

††copyright: acmlicensed††journalyear: 2022††doi:
10.1145/1122445.1122456††journal: TOG††ccs: Computing methodologies Shape
analysis††ccs: Computing methodologies Mesh models††ccs: Computing
methodologies Volumetric models

##  1\. Introduction

Volume meshes explicitly encode both the surface and the interior of an
object, thus offering a richer representation than surface meshes. They are
primarily used in industrial and biomedical applications, where volume
elements are exploited to encode various information, such as structural and
material properties, permitting to simulate and precisely estimate the
physical behavior of an object subject to external or internal forces, or the
dynamics involving multiple objects interacting in the same environment.
Alongside tetrahedra, hexahedral elements are the most prominent solid
elements used to represent discrete volumes in computational environments.
Meshes entirely or partially made of hexahedra have been used for many years
as the computational domain to solve partial differential equations (PDEs)
that are relevant for the automobile, naval, aerospace, medical and geological
industries to name a few, and are at the core of prominent software tools used
by such industries, such as (Distene SAS, 2020; ANSYS, 2021; Altair, 2021;
CUBIT, 2021; CoreForm, 2021a; Tessaels, 2021).

In academic research, the generation and processing of hexahedral meshes have
been studied for more than 30 years now. Despite the huge effort that various
scientific and industrial communities have spent so far, the computation of a
high-quality hexahedral mesh conforming to (or approximating) a target
geometry remains a challenge with various open aspects for which no fully
satisfactory solutions have been provided yet. Some of the known methods are
extremely robust and scale well on complex geometries; some others produce
high-quality meshes; some others are fully automatic. But no known method
successfully combines all these properties into a single product. The hex-
meshing problem had been so elusive that it was even once termed the holy
grail of mesh generation (Blacker, 2000). Ever since, many advancements in the
field have been made, while major challenges still remain.

In the last decade, the Computer Graphics community has contributed
significantly to the hex-meshing problem, proposing both seminal ideas and
practical algorithms. In this survey, we wish to summarize this work, also
reporting on previous methods developed by other scientific communities.

The engineering community has already produced a few surveys on this topic,
but they are either no longer up to date (Schneiders, 2000; Owen, 1998;
Tautges, 2001; Blacker, 2000) or focus just on a particular narrow subset of
the available techniques (Shepherd and Johnson, 2008; Armstrong et al., 2015;
Sarrate Ramos et al., 2014). We wish to create a comprehensive entry point for
researchers and practitioners dealing with hexahedral meshing. We therefore
embrace the whole field, revisiting and structuring a vast amount of
literature, and covering basic topological (Sec. 2) and geometrical (Sec. 3)
concepts, all kinds of approaches to hexahedral mesh generation (Sec. 4),
operators to edit mesh connectivity and to perform refinement or coarsening
(Sec. 5), mesh optimization and untangling (Sec. 6), visual exploration (Sec.
7), and also addressing the recent trend of methods for hex-dominant meshing
(Sec. 4.9). Last but not least, in the final part of the survey, we highlight
the current challenges the field is facing and indicate interesting directions
for future work (Sec. 8).

##  2\. Hex-Mesh Structure

A hexahedral mesh has structural aspects (concerning the connectivity of mesh
elements) and geometric aspects (concerning the elements’ shape and their
embedding or immersion in space). In this section we focus on the diverse set
of structural aspects, and consider geometry in Sec. 3.

###  2.1. Primal structure

In terms of connectivity, a hexahedral mesh is a 3-dimensional _cell complex_
, $\mathcal{H}=(V,E,F,C)$, consisting of vertices $V$ (0-cells), edges $E$
(1-cells), facets $F$ (2-cells), and cells $C$ (3-cells). The facets $F$ are
also often referred to as faces, and the 3-cells $C$ are, given the context,
often referred to as hexahedra or hexes. In a _pure hexahedral mesh_ , each
facet is a topological quadrilateral (i.e., incident to four edges) and each
cell is a topological cube (i.e., incident to six such facets). If a
relatively small number of facets and cells are of different type (e.g.,
tetrahedra, prisms, or pyramids) a mesh is called _hexahedral dominant_.

On top of this connectivity structure, a mesh is equipped with a geometric
structure, typically an embedding (or immersion) in $\mathbb{R}^{3}$ (Sec. 3).

Often, instead of assuming fully generic CW or $\Delta$ cell complexes
(Hatcher, 2000), more restricted connectivity definitions are used for
practical purposes (Erickson, 2013). A very common one is to assume that each
cell has eight _distinct_ vertices, i.e., no hexahedron is self-adjacent at a
vertex, edge, or facet. Similarly, pairs of edges, facets, or hexes being
adjacent via more than one vertex, edge, or facet, respectively, may be ruled
out. This simplifies data structures and algorithms; furthermore, many
applications assume each hex to be embedded in a geometrically simple way
(e.g. straight edges, ruled facets, cf. Sec. 3) which rules out such self-
adjacency and multi-adjacency anyway. Sec. 4.2 discusses further application-
dependent structural assumptions and requirements.

####  2.1.1. Singularities

The most _regular_ hexahedral mesh is an (infinite) Cartesian grid, where each
vertex, edge, and facet is incident to 8, 4, and 2 hexahedra, respectively.
General hexahedral meshes contain elements of different local connectivity,
which are accordingly called irregular or singular. Irregular facets simply
correspond to the boundary of the mesh, i.e., all facets that are incident to
a single hexahedron.

Since interior facets cannot be irregular and vertex singularities are never
isolated (Liu et al., 2018), structurally most interesting is the set of
irregular edges, which forms the so-called singularity graph.

Figure 2. Left: Singular edges of valence $k=\\{3,4,5\\}$. Right: There is a
1:1-correspondence between configurations of vertices in a hexahedral mesh and
triangulations of the sphere. Image from (Liu et al., 2018).

##### Singularity Graph

The fundamental building block of the singularity graph is a singular edge of
valence $k$, i.e., an interior edge incident to $k\neq 4$ hexahedra, or a
boundary edge incident to $k\neq 2$ hexahedra (cf. Fig. 2). While a single
integer $k$ is sufficient to characterize the structural type of an edge, the
specification of vertex types is more complex. As observed by Nieser et al.
(2011) there is a 1:1-correspondence between vertex configurations in a
hexahedral mesh, and triangulations of the 2-sphere. This can be understood by
observing that the intersection of a cube with a sphere centered at one of
it’s corners results in a triangular patch (Fig. 2). Hence, different vertex
types can be specified by enumerating all triangulations of the sphere.
Restricting to (the practically most relevant) edge valences $\\{3,4,5\\}$, it
turns out that only $11$ different configurations of interior vertex types
exist (Sabin, 1997; Liu et al., 2018). Specifically, for an interior vertex it
is impossible to be incident to a single singular edge, and in case of two
incident singular edges they can only be of identical type. Consequently, the
singularity graph is formed by singular arcs, which are chains of singular
edges with identical type. These singular arcs either terminate at the
boundary, or connect to other singular arcs at singular vertices, cf. Fig. 3.

Figure 3. The singularity graph is formed by interior singular edges of
valence $k=3$ (green) and $k=5$ (blue) and boundary edges of valence $k=1$
(green). Chains of edges with homogeneous type connect or terminate at
singular vertices (red). Image from (Liu et al., 2018).

###  2.2. Dual structure

In a hexahedral mesh, regardless of its level of structural regularity (Sec.
2.4), each cell has a constant number of 6 facets and each facet has a
constant number of 4 edges. Conversely, however, each vertex may have a
varying number of incident edges, and each edge a varying number of incident
facets.

One may consider the (polyhedral) cell complex that is _dual_ to a hexahedral
mesh: For each $k$-cell of the primal mesh there is a $(3-k)$-cell in
1:1-correspondence in the dual mesh and incidence relationships are adopted.
The above regularity of cells and facets in the primal mesh translates into
regularity of vertices and edges in the dual. Concretely, except at the mesh’s
boundary, each dual vertex has a constant number of 6 incident dual edges, and
each dual edge has a constant number of 4 incident dual facets. Conversely,
dual facets and dual cells are polygons and polyhedra of varying structure.
Further details and facts about the dual complex can be found in (Tautges and
Knoop, 2003).

Depending on the algorithmic context, it may be advantageous to consider the
primal or this dual view of a hexahedral mesh. A key reason is the following:
While vertices and edges of the primal mesh may be regular or singular, the
vertices and edges of the dual mesh are all regular; this is due to the fact
that all primal facets are quadrilaterals and all primal cells are hexahedra.
The following useful definition of _opposite_ edges at a regular vertex and
_opposite_ facets at a regular edge therefore applies everywhere in the dual
mesh.

##### Opposite Elements

At each regular interior vertex $v$, there are 6 incident edges. For each edge
$e_{1}$ of these, there is exactly one edge $e_{2}$ among these 6 that does
not share a facet with $e_{1}$; the edges $e_{1}$ and $e_{2}$ are called
_opposite at_ $v$. At each regular interior edge $e$, there are 4 incident
facets. For each facet $f_{1}$ of these, there is exactly one facet $f_{2}$
among these 4 that does not share a cell with $f_{1}$; the facets $f_{1}$ and
$f_{2}$ are called _opposite at_ $e$. In the primal setting, this concept of
opposite edge is relevant for algorithms that trace internal arcs in the mesh,
e.g., connecting pairs of singular vertices. Similarly, opposite facets are
useful to flood internal facet sheets bounded by singular arcs, e.g., to
perform a coarse block decomposition of a given mesh (Sec. 2.3).

Figure 4. Two examples (top and bottom) of the relation between dual and
primal mesh representation. Left: dual cell complex, formed by an arrangement
of 2-manifold sheets. Right: corresponding hexahedral mesh. Image from
(Murdoch et al., 1997).

In the dual setting, this opposite relation can be used to define the
following:

##### Sheets

Consider the transitive closure of the dual facets’ opposite-relation. Its
equivalence classes are called _sheets_ (also referred to as twist planes
(Murdoch et al., 1997) or (pseudo-)hyperplanes (Tautges and Knoop, 2003)).
These sheets are 2-manifold surfaces (with boundary, possibly self-
intersecting) formed by dual facets.

##### Chords

Analogously, the equivalence classes of the dual edges’ opposite-relation’s
transitive closure are referred to as _chords_ (Murdoch et al., 1997; Borden
et al., 2002a) (or _polychords_ (Daniels et al., 2008)).

The combinatorial continuity of opposite dual facets across dual edges has
inspired the early name _spatial twist continuum_ for this dual sheet based
perspective.

It follows that the entire dual complex can be viewed as an arrangement of
intersecting manifold surfaces (sheets): dual vertices are formed by three
intersecting sheets, chords are formed by two intersecting sheets and split
into dual edges by transversely crossing sheets, sheets are split into dual
facets by crossing sheets, and dual cells are the spatial compartments
enclosed by sheets. Conceptually, a sheet corresponds to one _layer_ of
hexahedra in the primal mesh; how this layer is composed of individual
hexahedra, however, is not defined by this sheet itself but by sheets that
cross this sheet transversely. Fig. 4 illustrates this primal-dual
relationship.

All this is in close analogy to dual complexes in the case of quadrilateral
meshes. These can be viewed as arrangements of intersecting 1-manifolds
(Campen et al., 2012; Campen and Kobbelt, 2014). Quite differently, though,
sheets can be topologically quite complex, they may have arbitrary genus and
an arbitrary number of boundary loops, whereas in the quadrilateral mesh case
each 1-manifold may only be either a closed loop curve, or an open-ended curve
starting and ending at the mesh boundary.

###  2.3. Block structure

Each hexahedral mesh can be decomposed into disjoint blocks, where each block
is a regular grid of hexahedra. Conversely, the mesh can be viewed as disjoint
union of such blocks. As an extreme example, each hexahedron could be
considered an individual block (of size $1\times 1\times 1$). We can
distinguish conforming and non-conforming block decompositions: a block
decomposition is conforming iff each side of each block coincides with one
other block side (except at the mesh boundary).

Of particular practical relevance are decompositions that are conforming, and
among these those that are coarse, i.e., that consist of a small number of
blocks. Meshes rarely have a unique conforming block decomposition. The
_coarsest_ conforming block decomposition is sometimes referred to as the
mesh’s _base complex_ (Bommes et al., 2011; Gao et al., 2015; Razafindrazaka
and Polthier, 2017).

It is worth pointing out that the term base complex is sometimes used with
alternative meanings (Livesu et al., 2013; Eck and Hoppe, 1996; Dong et al.,
2006; Hormann et al., 2008), for instance to refer to a coarse cell complex
that is used as a domain for (cross)-parametrization. Note that this is not
entirely unrelated though: a common use case of these parametrizations is
structured remeshing; the resulting meshes typically exhibit a block structure
induced by the underlying domain complex.

The base complex has the following defining property: a facet is part of a
block side if and only if it is transitively incident to a singular edge via
opposite facets (as defined in Sec. 2.2). This suggests a simple way to
extract the base complex of a given hexahedral mesh: starting from all facets
incident to any singular edge, iteratively expanding through opposite facets
across regular edges until termination. Due to the practical relevance of
semi-regular hexahedral meshes (Sec. 2.4) mesh generation algorithms that take
the coarseness of the implied base complex into account are of particular
interest.

###  2.4. Structure regularity

Similarly to quad-meshes (Bommes et al., 2013b), hex-meshes can be roughly
organized into four classes depending on the degree of regularity of their
topological structure. The concept of mesh regularity is closely related with
the relative amount of irregular vertices present in the mesh and with how
these vertices are connected to each other.

  * •

regular (or structured) meshes have the topology of a gridded cube (Fig. 1
left). These meshes are extremely convenient for storing and processing
because of their simple connectivity: each internal vertex has exactly the
same number of neighbors, with a consistent ordering. This allows for
efficient storage and optimal query time, and also makes the computation of
local quantities (e.g., finite differences) straightforward. There are,
however, severe limits in the class of shapes they can represent: mapping an
object containing long protrusions or deep cavities to a cube requires to
dramatically distort the grid, likely resulting in a mesh with no practical
usefulness due to the poor shape of its elements. Moreover, the rigid global
structure does not allow for localized refinement: if more vertices are
necessary around a specific area, the entire grid must be refined in order to
maintain the _pure hex_ property;

  * •

semi-regular (or semi-structured, also block-structured) hex-meshes are
obtained by gluing in a conforming way several regular grids (also called
_patches_ or _blocks_). In a semi-regular hex-mesh all vertices that are
internal to a patch are regular. Only vertices that lie at the edges or
corners of a block may possibly be irregular. Semi-regular meshes represent
the most important class in terms of applications, and are often the result of
a manual or semi-manual meshing process. Differently from regular meshes they
allow for higher flexibility and can be used to represent shapes of arbitrary
complexity. At the same time, they contain a limited amount of irregular
vertices, connected to each other so as to define a coarse block layout (Fig.
1, middle left) which can be exploited by dedicated data structures for
cheaper storage and fast querying (Tautges, 2004), and is also useful in a
variety of applications that exploit the tensor product structure of its
elements (e.g., IGA (Hughes et al., 2005));

  * •

valence semi-regular meshes also contain a limited amount of irregular
vertices, but they are not connected in a way that induces a coarse block
decomposition into few regular grids (Fig. 1, middle right). Meshes of this
kind are often produced by modern hex-meshing algorithms such as frame field
based methods, which introduce few singularities, but do not specifically
address their connectivity pattern;

  * •

irregular (or unstructured) hex-meshes contain a large fraction of irregular
vertices (Fig. 1, right). Meshes of this kind are often produced via
voxelization or other grid-based methods: portions of the object that do not
align with the ambient Cartesian grid exhibit a typical staircase effect,
triggering a proliferation of irregular vertices on the surface. Irregular
meshes are not suited for applications that exploit the coarse block structure
induced by the mesh connectivity, because the number of blocks is close to the
number of hexahedra in the mesh.

As for the quad-mesh case (Bommes et al., 2013b) the boundaries between semi-
regular, valence semi-regular, and irregular meshes are blurred. Nevertheless,
from an applicative perspective there is a substantial difference between
these three classes and there exists a variety of structure enhancement
algorithms that are specifically designed to improve mesh regularity (Sec.
5.5). The whole taxonomy can be understood in terms of the ratio between the
number of irregular vertices and the total amount mesh vertices ($r_{V}$), and
the ratio between the number of blocks and the total number of mesh elements
($r_{B}$). If both $r_{V}$ and $r_{B}$ are high, the mesh is irregular; if
$r_{V}$ is low but $r_{B}$ is high, the mesh is valence semi regular; if both
$r_{V}$ and $r_{B}$ are low the mesh is semi-regular. Finally, if the number
of blocks is exactly 1, the mesh is regular. Providing actual thresholds to
precisely define what high and low mean, is an application dependent matter.

Figure 5. An integer-grid map $f$ deforms the shape (left) such that its
boundary aligns with a Cartesian grid of integer isolines (right).
Consequently, its inverse $f^{-1}$ deforms the covered grid cells into a
structure-aligned mesh. Grid-preserving cuts (dashed green) enable irregular
vertices (red) in the mesh.

###  2.5. Integer-Grid Maps

Integer-grid maps (IGM) are a class of maps that exist in arbitrary dimensions
and which by construction induce structured meshes. The central idea, as
illustrated in Fig. 5, is to embed an $n$-dimensional shape into an
$n$-dimensional voxel grid such that the inverse map deforms the set of
covered voxels into a shape-aligned structured mesh. So far, integer-grid maps
have been studied for 2-manifolds to generate quadrilateral meshes (Kälberer
et al., 2007; Bommes et al., 2013a) and for 3-manifolds to generate hexahedral
meshes (Nieser et al., 2011; Liu et al., 2018). Similar to the parametrization
of a general manifold, an integer-grid map can be decomposed into multiple
charts. However, in order to guarantee that the inversely mapped voxels stitch
conformingly, it is necessary to require specific transition functions that
preserve the voxel grid. Assuming that the vertices of the voxel grid are
given by integer coordinates $\mathbb{Z}^{n}$, the grid-preserving transition
functions are exactly (i) integer translations and (ii) symmetry
transformations of an $n$-cube. Such transition functions are essential to
generate meshes with interior singularities, as for instance, the singular
vertex (red) in Fig. 5.

Mathematically, a map requires three properties to be an _integer-grid map_ :
(i) grid-preserving transition functions, (ii) local injectivity, and (iii)
singularities and boundaries mapping to integer-grid entities. A thorough
definition can be found in (Liu et al., 2018).

Integer-grid maps are sufficiently expressive to describe all potential
hexahedral meshes. We can trivially generate a chart for each hexahedron that
maps it to the voxel $[0,1]^{3}$. In this sense, integer-grid maps can be seen
as an alternative representation of hexahedral meshes that has proven highly
valuable for designing powerful generation algorithms.

Reformulating the hexahedral mesh generation task as a map optimization
problem offers many advantages. First of all, the optimization of low-
distortion maps is a well-studied topic with a rich body of theory and
algorithms that serve as a strong foundation. Moreover, the map optimization
perspective enables multiple geometrically motivated continuous relaxations
that are crucial for efficiently finding good approximate solutions of the
hard underlying mixed-integer problem, e.g. frame-fields to find suitable
singularities (cf. Sec. 4.8), or seamless maps to estimate the required
integer translations (cf. (Nieser et al., 2011), (Brückler et al., 2021)).

While a naive direct optimization formulation for a hexahedral mesh needs to
explicitly encode and deal with the full set of (inherently discrete) elements
and their connectivity, most of that becomes implicit in the map formulation,
enabling not only straightforward continuous relaxations but moreover a
reduced set of discrete variables. A simple but instructive example consists
of a regular block covering $n\times m\times o$ voxels in the IGM image.
Stretching the image along the first coordinate axes corresponds to a
continuous relaxation of the discrete action of changing the integer dimension
$n$. Note that from the map perspective, the block is indeed fully
characterized by only three integers $n$, $m$, and $o$, while a direct mesh
optimization would need to deal with $(n+1)\times(m+1)\times(o+1)$ (discrete)
vertices and their nontrivially-constrained connectivity. The number of
integer degrees of freedom of a general integer-grid map is proportional to
the number of singularities and topological handles. Consequently, in case of
pre-determined singularities the resulting discrete search space is
comparatively small since typically highly regular meshes with only few
singularities are desired.

Figure 6. Integer-grid map optimization algorithms turn the highly discrete
mesh optimization task into a more tractable mixed-integer map optimization.
All mesh quality criteria are directly related to properties of the map, as
indicated by blue dashed arrows.

Optimizing for a low-distortion map has two positive effects, (i) it directly
promotes well-shaped elements of high quality in the output hex-mesh, and (ii)
it demotes the occurrence of spurious singularities.

A conceptual overview of interpreting mesh optimization as map optimization is
shown in Fig. 6. The advantages related to superior continuous relaxations and
compact discrete search spaces explain the popularity and success of integer-
grid map based approaches in the automatic generation of structured meshes.

A special case of integer-grid maps with a single chart and no interior
singularities are polycube maps, which are discussed in more detail in Sec.
4.7. The optimization of polycubes, therefore, targets a deformation of the
input shape such that its surface aligns with the surface of the voxel grid.
Despite the continuous nature of the deformation, the resulting optimization
problem is nonetheless of mixed-integer type. The discrete degrees of freedom
are the choices of how to map surface normals since valid solutions are only
the coordinate axes $\\{\pm x,\pm y,\pm z\\}$ and a low-distortion assignment
is not known a priori.

Frame-field based methods, which are discussed in more detail in Sec. 4.8,
target the generation of hexahedral meshes in two stages. From a high-level
perspective, the first stage estimates the rotational part of the Jacobian of
an integer-grid map, i.e., a frame-field, while the second stage constructs
the map by inheriting the frame-field singularities. The decomposition is
beneficial because the first stage can be formulated in a representation that
automatically deals with the symmetry of the hexahedron. As a consequence,
frame field singularities can be optimized solely with continuous degrees of
freedom, converting the extremely difficult direct optimization of integer-
grid maps into a tractable form.

##  3\. Hex-Mesh Geometry

Besides its combinatorial and topological structure (Sec. 2), a hexahedral
mesh’s geometry, i.e., its embedding or immersion, typically in
$\mathbb{R}^{3}$, plays an essential role in most applications.

This concerns the question of geometric fidelity (to what extent the mesh
conforms to the target shape) and the question of element quality. This latter
question is concerned with the _shape_ of a mesh’s individual hexahedra or the
distortion of _maps_ defining these hexahedra as deformations of an ideal
(_reference_ or _master_) element.

Depending on the application context, various geometric requirements may be in
place: the mesh may be required to conform to a given boundary mesh or to
interpolate it within some prescribed tolerance; facets may be required to be
planar or to be convex; the above maps may be required to be locally injective
or even to have bounded distortion in some particular sense. In the context of
mesh generation (Sec. 4) the concrete requirements can have a significant
influence on the hardness of the meshing problem. Many methods so far are
unable to provide strict guarantees regarding such requirements, especially
when they are asked for in combination.

Also the relevant notion of element quality, and the effect of low or high-
quality elements, are application dependent. In the context of simulations by
means of finite element methods (FEM), element quality can have a crucial
impact on error estimates and convergence rates, thus simulation speed and
accuracy (Ciarlet, 2002; Zlámal, 1968), and relevant quality measures depend
on the type of simulation. In Sec. 4.2 these varying requirements are
discussed further.

##### The Trilinear Element

The geometry or embedding of hexahedral meshes is often represented by means
of coordinates assigned to their vertices. This alone is sufficient only for
simple applications. More commonly, the geometry of edges, faces, and cells
has to be defined as well. A particularly simple (and common) scenario is the
assumption of _trilinear_ elements (linear edges, bilinear faces, trilinear
cells), as this does not require the specification of any further
information—all other mesh elements’ geometric embedding in $\mathbb{R}^{3}$
are derived from the vertex positions via multilinear interpolation.
Precisely, a hexahedron’s embedding (with vertex positions $p_{ijk}$) is
defined via a _geometric map_ $\tau:[0,1]^{3}\to\mathbb{R}^{3}$ (also called
_isoparametric map_) as follows:

$\tau(u,v,w)=\sum_{i=0}^{1}\sum_{j=0}^{1}\sum_{k=0}^{1}p_{ijk}B_{k}(w)B_{j}(v)B_{i}(u),$

where

$B_{0}(t)=1-t\text{ and }B_{1}(t)=t.$

The hexahedral element effectively is the image of an ideal cube $[0,1]^{3}$
under this map. Note that the edges are straight line segments under this map;
the faces are ruled surfaces (planar iff the four corner vertices are
coplanar). More generally, this definition can be extended to higher-order
elements using higher-order basis functions $B_{i}^{n}$ (e.g., Bernstein
polynomials of degree $n$, giving rise to tensor-product Bézier elements
(Prautzsch et al., 2002)). In these higher-order cases, additional control
points (besides the vertex points) come into play as coefficients for a higher
number of basis functions.

The assessment of these elements’ quality (or even just validity) is an
application-dependent matter. In some cases it may be just the _shape_ of the
region $\tau([0,1]^{3})\subset\mathbb{R}^{3}$ that is of relevance, in others
its concrete _parametrization_ , given by the map $\tau$, is crucial.

###  3.1. Geometric map

A particularly common measure of quality is the determinant of the geometric
map’s Jacobian $J_{\tau}$. It quantifies to what extent the hexahedron,
defined through $\tau$, deviates (in terms of volume distortion) from the cube
$[0,1]^{3}$. Note that $\det J_{\tau}$ depends on parameters $(u,v,w)$. Due to
this dependence, the quality of an element (in contrast to the quality at a
particular point) rather needs to be assessed by the extremal value
$\min_{[0,1]^{3}}\det J_{\tau}$.

Note that $\det J$, while measuring volume distortion, is blind to angle
distortion; it cannot distinguish sheared cubes from cubes. Additional angle-
aware measures are thus often taken into account (Sec. 3.2).

####  3.1.1. Element Validity

If $\min_{[0,1]^{3}}\det J_{\tau}\leq 0$, the geometric map $\tau$ is non-
injective and the implied element is said to be irregular. Sometimes a
distinction is made between degeneration ($\det J_{\tau}=0$) and inversion or
fold-over ($\det J_{\tau}<0$).

In the context of the finite element method, irregular elements must be
considered invalid (Mitchell et al., 1971; Knupp, 2000); with such elements,
depending on the concrete setting, one may yield inaccurate solutions or no
solutions at all (Barrett, 1996), solutions are invalidated (Roca et al.,
2012), or calculations cannot be continued (Salagame and Belegundu, 1994). Due
to this crucial importance, specialized _untangling_ methods for the purpose
of irregular element removal in hexahedral meshes have been proposed (Sec. 6),
that attempt to achieve $\min_{[0,1]^{3}}\det J_{\tau}>0$.

####  3.1.2. Computation

The evaluation of $\det J_{\tau}$ at a concrete parameter point $(u,v,w)$ is
quite easy. For the computation of the extrema $\min/\max_{[0,1]^{3}}\det
J_{\tau}$, however, there is no closed-form expression. As this is
particularly relevant to _certify_ regularity, simply probing at a number of
well-distributed parameter points is a risky approach.

##### Determinant bounds:

Like $\tau$, the Jacobian determinant is a polynomial in $(u,v,w)$. It can
thus be expressed in the Bernstein basis as well:

$\det J_{\tau}(u,v,w)=\sum_{ijk}d_{ijk}B_{k}(w)B_{j}(v)B_{i}(u).$

Due to this basis’ implied convex hull property (due to $0\leq B_{i}(t)\leq 1$
for $t\in[0,1]$) the function value is bounded from below by the smallest
coefficient $\min_{ijk}d_{ijk}$ and from above by the largest coefficient
$\max_{ijk}d_{ijk}$. The coefficients $d_{ijk}$ are easily computed from the
vertex points $p_{ijk}$. For the particular case of trilinear hexahedral
elements, this is discussed in (Johnen et al., 2017). The same principle
applies to higher-order elements as well as to simplicial (rather than tensor-
product) elements (Johnen et al., 2013; Dey, 1999; Luo et al., 2002; Gravesen
et al., 2014; Mandad and Campen, 2020).

These bounds can be quite loose. They can, however, be tightened arbitrarily
by re-expressing $\det J_{\tau}$ piecewise over subdomains of $[0,1]^{3}$
(Hernandez-Mederos et al., 2006). This is accomplished (via affine
reparametrization) using Bézier subdivision (Prautzsch et al., 2002). Under
repeated subdivision, the coefficients (and thus the derived bounds) converge
to the actual function value (Prautzsch and Kobbelt, 1994; Leroy, 2008).

For use cases where precise knowledge of the Jacobian determinant’s value
range is not relevant but only injectivity is to be certified, simpler
(possibly loose) conservative tests can be employed (Zhang, 2005). Various
even simpler hypothetical tests (trying to derive bounds from determinant
values at vertices or along edges) were shown to be false (Knupp, 1990; Zhang,
2005).

##### Relaxation:

Through sum-of-squares (SOS) relaxation, the non-convex problem of finding the
Jacobian determinant polynomial’s global minimum (i.e., $\min_{[0,1]^{3}}\det
J_{\tau}$) can be replaced by a convex problem (Marschner et al., 2020). If a
sufficiently high degree is chosen for the formulation of this replacement
problem, the global minima coincide. A sufficient degree was determined
empirically; a formal guarantee is outstanding.

Table 1. List of alternative metrics for hexahedral elements, from the Verdict library (Stimpson et al., 2007). Normal ranges are intended for elements not having degeneracies. Metric | Overall | Acceptable | Value for
---|---|---|---
range | range | unit cube
Diagonal | $[0,1]$ | $[0.65,1]$ | 1
Dimension | $[0,\infty)$ | app. dep. | 1
Distortion | $[0,1]$ | $[0.5,1]$ | 1
Edge Ratio | $[1,\infty)$ | — | 1
Jacobian | $(-\infty,\infty)$ | $[0,\infty)$ | 1
Max. Edge Ratio | $[1,\infty)$ | $[1,1.3]$ | 1
Max. Asp. Frobenius | $[1,\infty)$ | $[1,3]$ | 1
Mean Asp. Frobenius | $[1,\infty)$ | $[1,3]$ | 1
Oddy | $[0,\infty)$ | $[0,0.5]$ | 0
Relative Size Squared | $[0,1]$ | $[0.5,1]$ | —
Scaled Jacobian | $[-1,1]$ | $[0.5,1]$ | 1
Shape | $[0,1]$ | $[0.3,1]$ | 1
Shape and Size | $[0,1]$ | $[0.2,1]$ | —
Shear | $[0,1]$ | $[0.3,1]$ | 1
Shear and Size | $[0,1]$ | $[0.2,1]$ | —
Skew | $[0,1]$ | $[0,0.5]$ | 0
Stretch | $[0,1]$ | $[0.25,1]$ | 1
Taper | $[0,\infty)$ | $[0,0.5]$ | 0
Volume (signed) | $(-\infty,\infty)$ | $[0,\infty)$ | 1

###  3.2. Shape quality

Besides metrics based on the pointwise assessment of the geometric map, there
exist a variety of metrics based simply on the vertex positions that have been
proposed in the literature to assess the quality of hexahedral elements or
have been exploited in specific applications. The documentation of the Verdict
library (Stimpson et al., 2007) – a de facto standard for finite element mesh
quality assessment – exhaustively reports per-hex metrics, as well as
associated bounds and commonly acceptable ranges. We succinctly report these
metrics in Tab. 1. For more details on how each metric is formulated, we point
the reader directly to the original source. It must be noted, though, that the
question whether an element is good or at least acceptable can be highly
application dependent; in FEM, for instance, elements far from being cube-
shaped (in particular anisotropically stretched elements) can be ideal – if
they are aligned suitably, in a PDE-guided or even solution-adaptive manner
(Knupp, 2007).

##  4\. Hex-Mesh generation

In this section, we survey all mesh generation techniques present in the
literature to date. We firstly provide a general introduction about input and
output requirements. Then, algorithms will be organized according to the
meshing paradigm they implement. The generation of hybrid, in particular hex-
dominant, meshes containing spurious non-hexahedral elements is also discussed
(Sec. 4.9). Finally, Tab. 2 summarizes the main properties of each class of
hex-meshing algorithms reported in this survey.

###  4.1. Input

Input data can be either a surface or a volume mesh describing the target
geometry. Methods that take a surface mesh or other surface description and
produce a conforming hexahedralization are often called _direct_ (Shepherd and
Johnson, 2008), as opposed to _indirect_ methods, which typically operate on a
supporting tetrahedral mesh and produce hexahedra by modifying this mesh
(through splitting, clustering, etc.) or by computing some volumetric mapping
encoded on the vertices of this supporting mesh.

The most trivial form of indirect hex-meshing consists of splitting each
tetrahedron into four hexahedra via midpoint refinement (Li et al., 1995).
This technique is trivial to implement and always guarantees a correct result.
However, it produces an unstructured mesh with an overly dense singular
structure, also containing four times more elements than the input mesh.
Therefore, this approach is unsuitable for real applications. As will come
clear in the remainder of this section, indirect hex-meshing has evolved
significantly since these early days and now comprises highly advanced tools
to convert a tet-mesh into a much coarser hex-mesh with cleaner singularity
structure. Notably, indirect approaches that cluster tetrahedra to form
hexahedra are quite predominant in hex-dominant meshing (Sec. 4.9).

Most of the techniques discussed in this section make assumptions on the
topology and geometry of the input mesh and are not able to operate on meshes
containing topological (e.g., open boundaries, holes, or non-manifold
elements) or geometric (e.g., intersecting or degenerate elements) defects.
Methods that operate on a supporting tetrahedral mesh may leverage robust
tetrahedralization techniques such as (Hu et al., 2018; Hu et al., 2020;
Diazzi and Attene, 2021). Methods that operate on surface meshes can sanitize
their inputs with known robust surface processing algorithms, such as (Attene
et al., 2013; Cherchi et al., 2020; Zhou et al., 2016; Attene, 2010).

In addition to the target geometry, algorithms may optionally take as input a
variety of other desiderata, such as target edge lengths or density fields to
control local element size, or a list of features that the output mesh should
conform to. Typical features are geometric curves on the outer surface (i.e.,
sharp creases), but there may also be additional ones – both internal and
external – such as separation membranes between different materials, or other
forms of semantic attributes. Finally, methods based on guiding fields (see
Sec. 4.8) may also take as input some additional parameters that control the
field generation, or may even assume the whole guiding field as an input by
itself.

Figure 7. Topology and geometry are tightly coupled: the number of hexahedra
incident to a singular edge directly bounds the inner angles, thus affecting
the geometric quality of the elements (numbers below each configuration refer
to the Scaled Jacobian of the geometric map). Image from (Livesu et al.,
2021).

###  4.2. Output

Output meshes must satisfy a variety of requirements, some of them strictly,
some others loosely. In the following we list the most important topological
and geometric requirements, also connecting them with specific applications
that demand their fulfillment. The main topological desiderata are:

  * •

element type:  methods that strive for pure hexahedral meshing must ensure
that all their cells are topological cuboids made of 8 vertices, 12 edges, and
6 quadrilateral faces. This requirement is loosened for hex-dominant methods,
where spurious non-hex elements may be present in the output mesh. This
topological freedom is not unlimited, and may be bounded by the specific
application. In fact, methods for the numerical solution of PDEs often require
non-hex elements to belong to a restricted class of polyhedra. For example,
the Poly-Spline Finite Element Method (Schneider et al., 2019a) demands that
all mesh elements (non-hexahedra included) have quadrilateral faces, and
enforces this property through mesh subdivision if the input mesh does not
fulfill this requirement. Similar restrictions are also imposed by alternative
methods;

  * •

local structure:  topological limitations may apply not only at a local (per
element) level, but also involve clusters of adjacent cells. For instance, the
Poly-Spline Finite Element Method (Schneider et al., 2019a) requires that two
non-hex cells are not face-, edge-, or vertex-adjacent, and also that non-hex
cells are not exposed on the boundary. More generally, many methods that
employ higher order basis functions can handle just a few local
configurations, and put constraints on the local mesh patterns. This holds for
both hex and hex-dominant meshes. For example, the blended spline method for
unstructured hexahedral meshes proposed in (Wei et al., 2018) embraces only a
small fraction of the possible singularities that are created by the meshing
methods surveyed in this section. To this end, the intricate mesh connectivity
generated by grid-based methods can be extremely challenging (Livesu et al.,
2021);

  * •

global structure:  depending on how the singular elements align, the mesh may
or may not have a coarse block structure (Sec. 2.4). While basic numerical
schemes like the Finite Element Method operate at a local (per element) level
and may not exploit this property, block-structured meshes may be highly
important for methods that employ tensor product constructions per block, for
multi-grid solvers that rely on a hierarchy of nested meshes, and also for
mesh compression (Tautges, 2004);

  * •

conformity:  some hex-dominant methods restrict their output to a narrow class
of alternative polyhedra (e.g., permitting only tetrahedra and hexahedra). On
the positive side, this restricts the alternative types of cells that
applications must handle. On the negative side, the resulting meshes may be
_non-conforming_ , meaning that structural discontinuities arise between
elements that are geometrically but not topologically adjacent (due to
T-junctions). Topological continuity can be restored using special
_connectors_. For example, a mesh containing two tetrahedra that are jointly
adjacent to a hexahedron can be made conforming by adding a zero volume
element containing one quadrilateral (at the hex side) and two triangular (at
the tet side) facets. Nevertheless, the resulting meshes (with or without
connectors) are not supported by all numerical solvers, and dedicated
numerical schemes (e.g., Discontinuous Galerkin (Chan et al., 2016)) must be
used.

From the geometric point of view, the output meshes should faithfully
represent the target shape, preserve its prescribed features (if any), and be
composed of well-shaped elements. More precisely, the main geometric
desiderata are:

  * •

fidelity:  geometric fidelity is achieved by construction by methods that
conform to an input quadrilateral mesh. Conversely, many other methods
typically deviate from the target geometry and may only produce a geometric
approximation of it. Just a handful of methods provide strict guarantees on
the maximum (Hausdorff) distance from the reference geometry (e.g., (Gao et
al., 2019)), whereas in the majority of the cases an unbounded approximation
of the input geometry is produced. Depending on the complexity of the input
shape, significant deviation from the target geometry may be present;

  * •

features:  special care must be paid for input features such as sharp creases.
While the general requirement is the same as for geometric fidelity,
imprecision in the geometric approximation of features is both aesthetically
much more evident and may also have a significant impact in the solution of
the PDE (e.g., when studying the aerodynamic flow around creased objects).
Feature alignment requires that sequences of edges of the hex-mesh conform to
feature curves, otherwise some deviation is inevitable, regardless of
resolution (Fig. 8);

  * •

quality:  the assessment of the quality of a mesh is a major topic by itself
(Knupp, 2007) that is only touched upon in this survey (Sec. 3). It is
important to note that the relation between mesh quality and, e.g., the
quality of a numerical solution of a PDE may heavily depend on the concrete
PDE as well as on the solver at hand. While a common requirement is that all
mesh elements are _valid_ (everywhere positive Jacobian determinant of the
geometric map), different numerical schemes may demand the fulfillment of
additional requirements. Shape regularity criteria for the Finite Element
Method (FEM) are mostly concerned with star-shapedeness and avoidance of large
angles (Ciarlet, 2002; Zlámal, 1968; Shewchuk, 2002). As recently shown, these
methods can be modified in order to even cope with badly shaped elements,
locally selecting higher order basis that compensate for the lack of geometric
quality (Schneider et al., 2018). In Computational Fluid Dynamics (CFD) it can
be beneficial to use meshes that are _orthogonal_ , meaning that the interface
between two shared elements and the line connecting their centroids form a
right angle (Moraes et al., 2013; Aqilah et al., 2018). The Virtual Element
Method (Beirão da Veiga et al., 2014) assumes that all mesh faces are planar.
Considering this jungle of metrics that are relevant for one numerical method
or the other, general purpose algorithms are often not suited to address these
specific criteria at the mesh generation stage, but mainly strive to create
meshes with valid elements, possibly addressing further quality concerns in
post processing (Sec. 6).

The methods surveyed in the following typically aim to create good meshes
according to a subset of the criteria above. Fully and equally embracing both
topological and geometric requirements at once can be a huge challenge, and
many methods put a stronger focus on one aspect over the other. Some methods
focus more on the topological aspects and may produce well structured meshes
containing (near-)degenerate or even invalid elements. Some others may
guarantee valid elements or even lower bounds on certain geometric quality
measures but produce meshes with a highly irregular topological structure.
Both flaws can potentially be alleviated to some extent in post-processing,
using dedicated algorithms for structure simplification (Sec. 5.5) or
geometric enhancement (Sec. 6). Certainly, topology and geometry are coupled
to some extent. For instance, a mesh with poor topological structure often
inevitably contains poorly shaped elements as well (Fig. 7).

Figure 8. Incorporating an input feature network (left, bold lines) into the
output hex-mesh is not possible if the connectivity does not align to it
(middle left), even refining the mesh (middle right). Key to feature
preservation is the ability to align surface edges to the input network,
carefully positioning mesh singularities (right). Image from (Livesu et al.,
2020).

###  4.3. Advancing/Receding front

First attempts to algorithmically generate hexahedral meshes were made by
extending 2D advancing-front algorithms that generated full quadrilateral
meshes. Starting from a quad-meshed boundary, algorithms like (Blacker and
Meyers, 1993; Blacker, 1996) incrementally insert hexahedra starting from the
boundary. The volume is progressively filled until final small voids are
solved with simple patterns made of a few hexahedral cells. Such an approach
is challenging on two main points. First, fronts can collide during their
generation, and geometrical intersection must be performed. Owen and Sunil
(2000) solve this issue by preserving a hybrid mesh during the whole process.
Every created hex is inserted into this mesh, and front collisions are easily
detected. The second point is much more problematic: there is no guarantee
that the process will eventually generate a usable full hexahedral mesh.
Starting from an even number of quads on the boundary of a remaining void, a
structural decomposition into a set of hexahedral elements is guaranteed to
exist (Mitchell, 1996), but the geometrical quality of hexes can be very low.
And if one ends up with an odd number of quads surrounding a remaining void,
one cannot fill it up with hexahedral elements at all, necessitating a
backtracking of the front propagation (with no general guarantee to perform
better the next time). The main reason for this inflexibility lies rooted in
the fact that one cannot easily perform structural modifications on a 3D
hexahedral mesh in a local manner (cf. Sec. 5).

Considering the problem as being over-constrained, the next generation of
advancing-front algorithms do not start from a quadrilateral boundary mesh,
but rather from the geometric surfaces (Staten et al., 2005, 2006; Staten et
al., 2010a). Complete layers of hexahedral cells are inserted in the domain
until they collide. Final cavities are easier to fill, but this process can
fail, too. In (Ruiz-Gironés et al., 2012), the authors adopt the advancing-
front technique. Considering that the final cavities that remain may be
difficult to mesh, they use an inside-outside mesh generation approach that
requires as an extra input an inner seed, which is a hexahedral mesh of a
possible final cavity. Two solutions of the Eikonal equation are then
computed: one going inward from the boundary of the geometric domain; another
one going outward starting from the surface mesh of the inner seed. Both
solutions are then combined to define a smooth distance function, and an
advancing-front algorithm is performed to expand the quadrilateral surface
mesh of the inner seed towards the unmeshed external boundary using the
distance function to locate points of each layer of cells. This process is
used in practice to mesh the outside of objects like aircraft (for
aerodynamics problems, for instance). But it remains limited to geometric
domains that are homeomorphic with the sphere, and the domain must not have
sharp features.

In general, advancing-front approaches are not reliable enough to generate a
good quality hexahedral mesh for general domains. They strongly depend on the
boundary mesh structure and the compatibility of this structure with the
restrictive structure of hexahedral meshes. Often, this compatibility is not
given, since the boundary mesh generation process is unaware of the structural
and geometric constraints imposed by the to-be-created hexahedral mesh. As a
consequence, they, e.g., fail to connect fronts when they collide (see Figure
9). Moreover, most of the proposed works deal with the extra constraint of
starting from a pre-meshed boundary. This constraint is strongly related to
the meshing process, which consists of meshing a complex assembly of parts
where meshes must be conforming along part interfaces.

Figure 9.  Example of advancing-front progression to fill in a geometric 3D
shape starting from different pre-meshed boundary. On the left, it succeeds in
getting a valid hex-mesh, while it fails on the right. Image from (Ledoux and
Weill, 2008).

###  4.4. Dual approaches

Taking a dual perspective in the context of mesh generation, i.e., focusing on
the dual representation of a hexahedral mesh (Sec. 2.2), has proven to offer
certain benefits.

##### Dual advancing front

For one, the interpretation of the advancing front approach (discussed in Sec.
4.3) in the dual domain can reveal interesting structures, simplify the
formulation of constraints and rules, and provide additional intuition. This
dual view is taken in the so-called Whisker Weaving (Tautges et al., 1996)
method and its variants (Folwell and Mitchell, 1999; Ledoux and Weill, 2008;
Kawamura et al., 2008). These start from a prescribed surface quad-mesh that
is to be matched by the hexahedral mesh to be constructed. Accordingly, the
quad-mesh’s dual loops form the prescribed boundaries of the hex-mesh’s dual
sheets. The algorithms’ objective thus is to determine the dual sheets – in
particular their mutual intersection combinatorics – inside these prescribed
boundary curves.

The addition of a next hexahedron in the course of an advancing front approach
can be interpreted in the dual as (combinatorially) fixing the intersection of
three dual sheets (Tautges et al., 1996) or as locally (combinatorially)
contracting one of the sheet boundary loops, conceptually fixing part of the
dual sheet and leaving the loop as the boundary of that part of the sheet that
is yet to be determined (Folwell and Mitchell, 1999). The dual view enables
the formulation of local and semi-local rules and heuristics to more favorably
steer the incremental mesh construction process (Ledoux and Weill, 2008;
Folwell and Mitchell, 1999). Nevertheless, issues such as poorly shaped
elements, inverted elements, or high valence vertices in the result are not
easy to avoid in general, even with this dual perspective.

A particular challenge for this approach is posed by the (very common)
existence of self-intersecting dual loops in the prescribed boundary quad-
mesh. While there is no general theoretical obstacle to the successful meshing
of these, such loops need to be brought into pairwise or manifold
correspondence and be filled by common sheets of non-trivial topology. It is
unclear how the process can be steered to naturally establish this required
structure in general; therefore, degenerate elements (so-called knives, Fig.
32) and inverted elements are common in the result in these cases. Various
strategies (with more or less severe negative side effects on quality) have
been proposed to modify the quad-mesh to get rid of such self-intersections in
advance (Folwell and Mitchell, 1999; Kawamura et al., 2008; Müller-Hannemann,
2001, 2002).

##### Dual sheet-by-sheet

Besides these alternative interpretations of advancing front methods, the dual
perspective gives rise to a further class of methods, less local and
incremental. A general challenge faced by algorithms that attempt to construct
hex-meshes in an incremental fashion (like those discussed in Sec. 4.3) is to
ensure that things work out in the end. Without careful look-ahead, one may
easily end up in intermediate configurations that cannot be completed in
either a valid or a qualitatively reasonable manner. Algorithms that, by
contrast, approach the problem of mesh generation in a global manner, e.g.,
via global optimization formulations (cf. Sec. 4.8), on the other hand, can be
computationally much more intensive.

The dual perspective permits an interesting incremental approach on a _semi-
local_ level. Instead of individual cells, entire dual sheets can be
considered as the atomic entities for incremental mesh generation in the dual
domain. For the case of quadrilateral mesh generation, which is in close
analogy to the hexahedral mesh generation scenario, the advantages of this
semi-local dual view for the purpose of incremental construction have been
discussed in depth (Campen et al., 2012; Campen and Kobbelt, 2014). Similar
properties hold in the hexahedral case, as is exploited by a number of
algorithmic approaches. However, while in the quad case the dual is formed by
chords, which are 1-manifolds (i.e., either a loop or a curve with two
endpoints, possibly self-intersecting in points), in the hex case the dual
consists of sheets, which are 2-manifolds of arbitrary genus and with an
arbitrary number of holes, possibly self-intersecting in curves. Therefore,
the problem is of significantly higher complexity and algorithms often
restrict to sub-classes of problem instances for simplicity, such as objects
of genus 0, sheets with a single boundary loop, or dual loops without self-
intersections. An idea of inserting dual sheets in a divide-and-conquer manner
was outlined by (Calvo and Idelsohn, 2000). A concrete algorithm for
incremental hex-mesh construction based on sequential dual sheet generation is
described by (Müller-Hannemann, 2001). The boundary geometry along an entire
candidate sheet is assessed in the decision-making process. In contrast to
related methods that can be interpreted as operating in a sheet-by-sheet
manner (Folwell and Mitchell, 1999; Ledoux and Weill, 2008), this algorithm
preserves an invariant through all intermediate stages that strictly avoids
combinatorially invalid configurations. This obviates the need for
intermediate repair operations and guarantees the absence of degenerate
elements such as knives or wedges (Fig. 32). On the downside, the more
restrictive sheet selection rules that are in place to ensure the invariant
can bring the algorithm to an early halt. Rather expensive back-tracking
strategies can be used as a remedy to some extent. By the introduction of
additional rules for the selection of sheet operations (Kremer et al., 2014)
in particular non-convex shapes can be handled in a more geometry-aware
manner, commonly leading to less distorted (or less inverted) mesh elements.

##### Free boundary

The above methods assume that a quadrilateral mesh of the domain boundary is
given, effectively as a starting point for the incremental construction. The
ability to prescribe a boundary mesh can be seen as an advantageous feature in
some scenarios (e.g., when adjacent domains are to be meshed separately but
compatibly). In others, it rather is a limitation: it restricts the meshing
algorithm from the set of all hex-meshes suitable for the domain to a (small)
subset. Recently there have been first attempts to construct hex-meshes on a
sheet-by-sheet basis _without_ predetermined boundary structure. Instead, they
exploit interactive user guidance along the domain boundary (Takayama, 2019)
(Fig. 10), or loosely follow principal curvature directions (Livesu et al.,
2020) to construct loops which then serve as sheet boundaries. It is worth
remarking that the latter method essentially outputs a subdivided version of
the primal mesh that is implied by the sheets; this has the effect that the
sheets appear as primal facet sheets in the output mesh. Nevertheless,
conceptually both methods are to be viewed as dual approaches.

Figure 10. Interactive sheet-based hex-mesh modeling. Image from (Takayama,
2019).

Due to the larger search space compared to methods with prescribed boundary
mesh, they conceptually have the potential to achieve results of better
quality – but at the same time are computationally more expensive and require
a user in the loop (Takayama, 2019) or make simplifications sometimes leading
to meshes that contain some non-hex elements (Livesu et al., 2020).

In this context, the interesting question is that of efficient geometric sheet
representation – while in the above methods assuming a prescribed boundary-
mesh, a non-geometric combinatorial representation was employed for
simplicity. An implicit representation by means of a level set formulation has
proven efficient (Takayama, 2019). It, however, does not support self-
intersecting sheets, which would grant higher flexibility and enable better
mesh quality in various cases. Another, discrete sheet representation space is
described by (Roca and Sarrate, 2008), embedding sheets in the facets of a
particular tessellation of the domain; a concrete algorithm that operates in
this space has not been addressed yet.

##### Dual validity

Generally, when constructing hex-meshes out of dual sheets, it needs to be
considered that not any arrangement of intersecting sheets implies a primal
hex-mesh. A number of conditions need to be satisfied so as to avoid non-
manifold configurations and self-adjacent elements, as detailed by (Mitchell,
1996). Violating sheet arrangements can be modified, often through the
insertion of additional sheets, to ensure these conditions are met (Folwell
and Mitchell, 1999). As these modifications not rarely have a negative impact
on (geometrical and structural) mesh quality, a relevant challenge is to avoid
the need for them right from the start.

###  4.5. Domain decomposition

Early proposals for automatic domain decomposition relied on simple
topological operations like submapping and sweeping (White et al., 1995), that
were mainly trying to incorporate the knowledge of the users upon the two-
dimensional domain to expand the decomposition to the third dimension with a
sweeping step.

##### Sweeping.

Given a volume represented by a closed surface, by identifying two patches
where one serves as the _source_ and the other one as _target_ , a hexahedral
mesh can be generated through sweeping the quad-meshed _source_ over the
volume to the _target_ (Shih and Sakurai, 1996). This simple idea is very
suitable for CAD models since many shapes are formed by extrusion. The first
batch methods using such an idea focus on shapes that can be easily meshed by
identifying one source and one target, which are called one-to-one methods
(Blacker, 1996; Liu and Gadh, 1997; Liu et al., 1999). However, for slightly
complex CAD models, more source or target patches have to be involved in
decomposing the extrusion geometry into simpler one-to-one sub-volumes for
easy processing.

A step ahead towards automatic decomposition is presented by Lu and colleagues
(2001) that suggest recognizing in a CAD model the characteristics of portions
that can be treated as submappable. The pipeline uses first a feature
recognition, then a cutting plane identification, and, finally, a
decomposition to mesh each portion with predetermined schemes. Along this
direction, a set of many-to-one and many-to-many approaches are developed (Lai
et al., 2000; White et al., 2004; Scott et al., 2006; Wu and Gao, 2014; Wu et
al., 2018). These methods often rely on specific rules to detect line and
planar features, such as various angle thresholds, so that the 3D model can be
decomposed into sub-volumes having the same sweeping direction. If the
decomposition is successful, various node insertion tricks for the sweeping
can be employed to ensure the high quality of the generated hex-mesh (Knupp,
1998; Staten et al., 1999; Ruiz-Gironés et al., 2011). There are also
approaches that allow multiple sweeping directions by computing a hierarchical
sub-geometry structure (Miyoshi and Blacker, 2000).

Kowalski and colleagues (2012) introduce the notion of fundamental sheets
(fun-sheets), noticing that a hexahedral mesh is layered, in opposition to the
lack of reference surfaces typical of tetrahedral meshes. Starting from a tet-
mesh, converted in a hex mesh and identifying these fun-sheets, using topology
and geometry of the shape, they obtain a better decomposition that catches the
intrinsic characteristics of the shape. This approach is further enhanced in
(Wang et al., 2017).

An interesting approach to the problem is the one presented by Lu and
colleagues (2017). They design and implement a sketch-based decomposition tool
and evaluate its performance on a group of beginners and experienced users.
They conclude that visual assistance and a geometric reasoning engine can help
to obtain excellent results from a semi-automatic decomposition.

##### Medial descriptors

Medial descriptors are a valid proxy in helping to realize a domain
decomposition. Both mechanical objects and free-forms are possible to identify
characteristics, mainly the skeleton catching the crucial elements of the
shape’s mutual relations. A three-dimensional shape’s skeleton is, in fact, a
topological representation of the shape capable of providing information
regarding the various boundary entities’ relative positions. The skeleton has
been used in multiple methods to help in generating a hexahedral mesh inside
the shape (see e.g., Fig. 11). When dealing with mechanical objects, usually
containing boxes, it is vital to use the general skeleton (or medial object),
including surfaces. Shapes more related to biology approximable with a
collection of generalized cones can be easily represented by their curve-
skeleton (or medial axis). Both these proxies have been used to guide the hex-
meshing.

Figure 11. Skeleton driven hex-meshing starts from an input surface mesh and
line skeleton (left), around which a tubular structure composed of hexahedral
boxes is initialized (middle left). Refining this structure and projecting it
on the target surface yields a hexahedral mesh (middle right) where the
distribution of the mesh elements aligns with the skeleton guiding curves
(right closeup). Image from (Livesu et al., 2016).

Price and colleagues introduced the possibility to use the topological
skeleton of the shape to produce a hex-mesh. They apply it first on convex
shapes (Price et al., 1995), and then on solids with flat and concave edges
(Price and Armstrong, 1997). The idea is to decompose the domain so that each
sub-domain can be hex-meshed using a midpoint subdivision scheme (Li et al.,
1995). Each sub-domain is meshed using basic primitives that can be placed
using the skeleton and used as elementary blocks to mesh the original domain.
The topological information guides the choice of the correct primitive. There
are limitations in the approach since high-valence boundary vertices do not
have elementary schemes placing them.

Instead of using the skeleton, Sheffer and colleagues (1999) start from the
embedded Voronoi graph of the domain, which is simpler to create. Using a set
of configurations that include the Voronoi graph’s local topology, it can
decompose the domain in sweepable subdomains that can be combined and smoothed
to yield the final decomposition of the whole domain. Through the computation
of a harmonic field, a general 3D model can be decomposed into 2D curved
slices where quad-mesh templates can be used to form a large structure
decomposition of the 3D model (Gao et al., 2016).

Zhang and colleagues (2007) exploit the particular shape of the vascular
structure to devise a method that uses the curve skeleton as a basis for the
meshing. It is the first proposal in which there is decomposition in tubular
subdomains that are quite simple to mesh via sweeping. The uniform diameter of
the typical vases treated in the application does not pose the problem of
resolution in the elements. Usai and colleagues (2015) use the curve-skeleton
to derive a quadrilateral base complex given the triangular mesh of shape. The
surface decomposition can be expanded to the domain’s interior and lead to a
method for hex-meshing (2016). In this work, a scheme for keeping the mesh
elements uniform while the diameter of the subdomains changes is introduced
and applied. Another similar approach (Livesu et al., 2017) employs solid
cylindrical parameterizations to map from the curve-skeleton to the
cylindrical subdomains. This choice allows a simple but effective way to use
the topological information to generate the hex-mesh.

All the methods described in the previous paragraph work fine only for models
resembling collections of generalized cones.

Quadros (2014) also uses the skeleton as a starting point for meshing and,
combining it with an advancing front approach, can create hex-dominant meshes.
The surface and the skeleton jointly contribute to form what the author calls
corridors that are the basis for meshing the domain with an advancing front
method.

Cai and Tautges (2015) propose an approach that heavily relies on integer
programming due to the classification of the edges for their parameterization.
It is in line with the topological methods since it introduces a new set of
templates that, once applied to the class of objects they use in their
experiments: mechanical parts.

Another interesting approach (Liu et al., 2015) mixes skeletal representation
of the shape and polycubes to guide the creation of the hex-mesh. The
resulting meshes are non-conforming, including T-junctions. Another type of
non-conforming decomposition, the so-called motorcycle complex (Brückler et
al., 2021), can be constructed guided by a seamless parametrization (cf. Sec.
2.5). This decomposition has hexahedral subdomains only, and can be refined
into a conforming hexahedral mesh.

Once a suitable decomposition is computed, submapping or sweeping approaches
can likely generate hex-meshes with satisfactory quality. For example, (Wu et
al., 2017) can be employed to first generate a quadrilateral mesh for
interfacing surfaces while ensuring conformity among adjacent sub-volumes, and
then apply the straightforward sweeping to generate the final hex-mesh.
However, up to now, the critical issue still lies in how to robustly decompose
a 3D model into sweepable sub-volumes while ensuring the necessary conformity
and feature preservation at the interfaces of different parts. Industry
resolved this issue by putting the user in the loop, relying on manual block
decomposition and automatic sweeping as a workhorse for hex mesh generation in
commercial software (Altair, 2021; ANSYS, 2021). Automatizing the process and
freeing the user from tedious critical work is an open challenge for future
methods in this family.

Figure 12. Dual methods regularize the valence of hanging nodes (black dots)
by connecting them pairwise along triangular bridges, so that the dual is a
hex-mesh. Top: the flat transition firstly introduced in (Maréchal, 2009).
Bottom: the set of atomic schemes introduced in (Livesu et al., 2021) to
handle all possible transitions in strongly and weakly balanced grids. Images
partly from (Livesu et al., 2021).

###  4.6. Grid based

A hex-mesh can be trivially created by voxelizing the interior of a closed
surface and then projecting its boundary onto the target geometry (Schneiders,
1996a). Geometric fidelity can be controlled by tuning the resolution of the
voxelization. Since the size of regular grids grows cubically, to reduce
element count a set of adaptive spatial partitioning approaches that rely on
hierarchical structures have been proposed. However, adaptive grids do not
define a conforming hex-mesh because adjacent grid elements may have different
size, generating spurious (_hanging_) nodes. Grid-based methods differ to each
other for the refinement policy they use, for the technique used to suppress
hanging nodes, or for the method used to project the mesh on the target
geometry.
Methods in this class are among the firsts that were introduced in the field.
From a mesh quality standpoint, they are typically considered inferior to
other methods because: (i) the grid is fixed in space and the result depends
on the orientation of the model; (ii) the connectivity they generate is
intricate and rich of singular edges with high valence (Livesu et al., 2021);
(iii) the meshes they generate are highly unstructured and do not endow a
coarse block decomposition (see Fig. 1 and Fig. 21 in (Livesu et al., 2020)).
Nevertheless, when compared with alternative options grid-based methods really
stand out in terms of robustness. To date, they are the only fully automatic
methods capable of successfully hex-meshing any input shape, regardless of its
geometric or topological complexity. For this reason, they are the only
automatic methods currently implemented in professional software (Distene SAS,
2020; CoreForm, 2021a; CUBIT, 2021). Despite the most prominent methods were
developed more than 10 years ago and the field remained quiet for some years,
major improvements have been proposed in recent years, also opening avenues
for further research.

##### Refinement.

Grids should satisfy both local and global criteria. At a local level, cell
size must be compatible with the local size of the input object, ensuring
geometric fidelity. At a global level, it must be possible to select a subset
of grid elements (e.g., the ones completely internal to the input shape) such
that the topology of this arrangement matches the one of the original object.
In case the grid and the input mesh are not homotopic, a bijective mapping
between them is not possible. Local criteria are easier to enforce. The most
typical split rules used in the literature are normal similarity (Ito et al.,
2009), local thickness (Livesu et al., 2021; Pitzalis et al., 2021; Maréchal,
2009), surface approximation (Gao et al., 2019) or a combination of these and
other indicators (Bawin et al., 2021). The fulfillment of global criteria is
more complex and demands to preprocess the input shape (Mitchell and Vavasis,
1992). For this reason, the vast majority of methods do not guarantee that the
output hex-mesh will have the same genus and number of connected components of
the input model (Maréchal, 2009; Livesu et al., 2021; Pitzalis et al., 2021),
or ensure this property at the cost of severe over refinement (e.g.,
iteratively splitting all grid elements until topological equivalence is
obtained (Gao et al., 2019)). Refined cells can be split in two alternative
ways: _2-refinement_ splits each edge in two, thus obtaining 8 sub-cells for
each adjacent hexahedron; _3-refinement_ splits each edge in three, thus
obtaining 27 sub-cells. In both cases, the sequence of splits is encoded in a
hierarchical tree structure, which corresponds to an octree for the
2-refinement, and to a 27-tree for the 3-refinement. Approaching this body of
literature for the first time may be confusing, because all methods generally
refer to these data structures as octrees, even though this is not always
correct. The use of 27-trees for 3-refinement is explicitly mentioned in
(Schneiders et al., 1996) and a few other articles, and is only implicitly
assumed in other articles that refer to these ones.

Figure 13. The set of 6 out of 22 schemes of 3-refinement introduced in
(Schneiders, 1996b). The 1 and the 5 are respectively the empty and the fully
refined hexahedra, while the other templates are used to manage the transition
between them. Image taken from (Schneiders, 2000).

##### Hanging nodes.

The removal of hanging nodes is obtained by substituting elements of the grid
with templated topological transitions that locally restore mesh conformity.
If adjacent grid elements differ by at most one level of refinement there
exist $2^{8}$ alternative configurations which, discarding symmetries, reduce
to 20 unique cases (Weiler et al., 1996). Existing methods can be broadly
categorized into two families: _primal_ methods aim to directly incorporate
the hanging nodes in the output hex-mesh; _dual_ methods aim to modify the
input grid such that its dual mesh contains only hexahedral cells.
Primal methods often operate on 3-refined grids and 27-trees, because it is
easier to suppress their hanging nodes (Schneiders et al., 1996). However,
handling all the possible 20 configurations is provably impossible, because
many concave transitions are bounded by an odd number of quadrilateral
elements, a condition for which it is known that a hexahedralization of the
interior does not exist (Mitchell, 1996). Transition schemes for 4 flat and
convex transitions (see Fig. 13) appeared in multiple articles (Schneiders,
2000, 1999, 1997; Tack et al., 1994) and were successfully used to compute
hexahedral meshes, prescribing additional refinement to convert unsupported
transitions into the supported ones. Over the years additional schemes were
introduced to handle concave edges (Ito et al., 2009; Zhang and Bajaj, 2006;
Elsheikh and Elsheikh, 2014), but a correct handling of concave corners
remains elusive. Several works, like (Ebeida et al., 2011; Zhang et al., 2013;
Owen et al., 2017), exploit the 2-refinement schemes introduced in (Schneiders
et al., 1996) to remove hanging nodes. Unlike from the 3-refinement
approaches, the grid needs to satisfy more strict constraint as those
described below for dual methods. Note that, as for the 3-refinement case, the
schemes in (Schneiders et al., 1996) do not allow to address all the possible
configurations, often leading to an excessive over-refinement of the grid.
Dual methods operate on 2-refined grids and octrees, and are superior to
primal methods because they can handle all possible transitions. All known
schemes operate on _balanced_ grids, that is, grids where the refinement
mismatch between adjacent elements is at most one. However, not all methods
agree on the definition of adjacent. For the majority of methods two cells are
adjacent if they share one face, edge or vertex (_strong balancing_). In
(Livesu et al., 2021) the authors relaxed this formulation, enlarging the
class of balanced grids and limiting restrictions to size mismatch only for
cells sharing a face (_weak balancing_). Weakly balanced grids permit to
greatly reduce refinement (up to $64\%$ less elements in their experiments),
but require a slightly more complex scheme set. Maréchal was the first to
observe that if all grid vertices have valence 6 and all grid edges have
valence 4, the dual of the grid is a pure hexahedral mesh (Maréchal, 2009).
Based on this observation he proposed a set of cutting schemes that,
regularizing the valence of grid elements, allow to obtain a pure hexahedral
mesh via dualization (Fig. 12).

Figure 14. 2D pipeline of the feature preservation octree-based hex-meshing.
Top row: adaptive quadtree constructed from the input, dual of the quadtree,
and quadrilateral (quad-) mesh including a scaffold mesh. Bottom row:
topological matching of feature graphs, variational padding of both the target
mesh and the scaffold, mesh deformation to fit the input, and the final pure
quad-mesh. Image from (Gao et al., 2019). Figure 15. If a quad maps two of its
four edges onto a linear feature line it becomes locally degenerate (left).
Splitting it into 5 sub-quads ensures enough degrees of freedom to produce all
valid mesh elements. Similar configurations may also occur on 3D meshes, and
can be resolved with a special padding scheme that splits a hexahedron into 6
sub elements (right). Image from (Gao et al., 2019).

Since the valence of hanging nodes is fixed pairwise, dual methods also
require that the grid is _pair_ , that is, for each cluster of grid elements
with same amount of refinement the number of hanging nodes must be even across
all grid directions. Differently from balancing, the pairing condition is non
local, hence difficult to enforce. Pairing is typically enforced directly in
the octree, fully splitting parent nodes if their siblings have been split
(Maréchal, 2009; Gao et al., 2019; Hu et al., 2013; Livesu et al., 2021). As
shown in (Pitzalis et al., 2021) all these methods operate in a restricted
space of solutions and tend to severely over refine the input grid, even if it
is already pair. The authors showed that pairing can be enforced directly in
the grid by solving a sequence of linear problems, obtaining coarser grids
that approximately halve the number of elements. Despite superior to tree-
based methods, also this method does not cover the whole space of solutions,
and may occasionally refine an already pair input grid (see Sec. 7 in
(Pitzalis et al., 2021)). Even though dual approaches exist since 2009, the
transition schemes they use were only vaguely described in the literature,
making these methods hardly reproducible. Maréchal (2009) pioneered this
technique, but his paper describes in detail only one specific transition
(Fig. 12, top). Gao and colleagues proposed three alternative schemes based on
similar ideas (Gao et al., 2019), also releasing their code, but these schemes
were recently shown to be not fully exhaustive and may fail to produce a
conforming hex-mesh starting from a balanced and paired grid (Livesu et al.,
2021). In (Livesu et al., 2021) the authors propose a comprehensive study of
dual schemes, clarifying ambiguities and implementative choices, and
ultimately deriving an exhaustive optimal set of transitions for both strongly
and weakly balanced grids (Fig. 12, bottom). CinoLib (Livesu, 2019) hosts an
open source implementation of all such schemes, as well as the code necessary
to install them in a given adaptive grid.

##### Projection.

Considering the axis-aligned nature of grid-based methods, to approximate the
input object well the boundary vertices have to be projected onto the target
geometry. To this end, maintaining the inversion-free property of a hex-mesh
poses a great challenge. While (Maréchal, 2009; Lin et al., 2015) rely on
iterative vertex smoothing to slowly move the vertices onto the boundary so
that a local smoothing can be backtracked if it causes flipped hexahedra, (Gao
et al., 2019) presents a global deformation method that can robustly align the
generated hex-mesh with the input surface (including sharp features) within a
distance bound. Fig. 14 shows the 2D pipeline of the method presented in (Gao
et al., 2019). After grid refinement and removal of hanging nodes, the grid is
partitioned into two sub-meshes: an inside target mesh that will be optimized
to be the final output, and an outside scaffold mesh that ensures the
bijectivity of the map throughout the optimization process. Geometric fidelity
is achieved by first building a topological bijectivity mapping between the
input mesh and the boundary of the target mesh, and then geometrically
deforming the target mesh towards the input surface shape using a locally
injective mapping technique (Rabinovich et al., 2017). Note that a variational
padding technique (see Sec. 5.4) is also introduced for both the target mesh
and the scaffold, so as to increase the number of degrees of freedom for
optimization. The approach can robustly produce an all-hexahedral mesh with
several guarantees: 1) the output is manifold and its boundary surface has the
same genus with the input, (2) all hexahedral elements have positive scaled
Jacobian (3) the boundary of the hex-mesh is error-bounded, i.e., within
$\epsilon$ distance from the input mesh, and (4) the boundary of the mesh has
no self-intersections thanks to the scaffold mesh. All of this is obtained by
trading robustness for efficiency, thus computational cost and memory
resources can be prohibitive for commodity hardware. On the other hand,
iterative methods such as (Maréchal, 2009; Lin et al., 2015) are quite
efficient, although may occasionally fail to preserve the shape well. Further
research is needed to devise an algorithm that optimally combines robustness,
efficiency and geometric fidelity.

Figure 16. The pipeline for polycube based hexmeshing generates a locally
injective simplicial map through volumetric deformation, and then uses it as a
medium to transfer a regular grid sampling of the polycube to the target
shape.

##### Features.

The preservation of sharp surface features is both geometrically and
topologically challenging for grid-based approaches. First of all, since the
mesh connectivity is derived by the underlying grid, surface vertices may not
have enough incident edges to reproduce high valence feature points in the
target mesh. Therefore only a subset of all possible feature networks can be
faithfully reproduced. Moreover, hexahedra that have more than one facet
exposed on the surface may easily be traversed by feature lines across more
than one edge, becoming ill-shaped or even degenerate once projected onto the
target geometry. To make sure that each element has at most one feature edge,
specific padding schemes are used (Fig. 15 and Sec. 5.4). Finally, despite the
fact that it works well in most cases, current algorithms for feature mapping
are heuristic and do not offer guarantees. The most recent methods are based
on ideas expressed in (Gao et al., 2019), and operate by iteratively
processing each feature separately, projecting its endpoints to the closest
vertices in the hex-mesh, and then finding the discrete path that connects
them with a Dijkstra search that operates on a scalar field that encodes the
euclidean distance from the input feature. Depending on the ordering of the
features and the combinatorics of the hex-mesh, there can be conflicting
configurations where a path that connects the two endpoints of a feature and
does not conflict with any previously inserted feature does not exist.
Furthermore, even if such a path exists, there may be cases in which the
previously inserted features force a path to deviate from its geometric target
significantly.

##### Assemblies and multiple materials.

While all methods described so far assume as input a single model composed of
a single material, grid-based techniques have been successfully extended to
the multi material case (Su et al., 2004; Zhang et al., 2010), and can also
handle complex non manifold CAD assemblies (Qian and Zhang, 2012). From a grid
processing perspective, these method rely on the processing techniques
described in the previous paragraphs.

###  4.7. Polycube maps

A successful line of algorithms works by volumetrically mapping a shape to an
orthogonal polyhedron (or _polycube_ (Tarini et al., 2004)) embedded in
$\mathbb{Z}^{3}$. Sampling the polycube at a dense integer lattice gives a
regular all hex connectivity, whose nodes can be positioned inside the initial
object following the inverse map (Fig. 16). Polycube methods are based on two
fundamental building blocks: the definition of the polycube structure, and the
generation of the volumetric map. These two objectives can be pursued
separately (i.e., defining a valid polycube structure first, and then
computing the map) or together, using mesh deformation to explore the space of
shapes and find the orthogonal polyhedron closest to the input object.

##### Structure.

The structure of the polycube is typically defined by assigning to each
surface element of the mesh a label that represents one of the six global axes
($\pm X,\pm Y,\pm Z$). Clusters of adjacent elements with same the label
identify the facets of the polycube. A purely local naive approach consists of
assigning to each surface element the axis closest to its normal (Gregson et
al., 2011). However, local surface orientation often does not produce a
globally consistent structure, and heuristic post-processing is necessary to
fix the labeling. Hu and Zhang (2016) use a modified CVT labeling in the space
of normals, but still require manual post-processing to sanitize topological
inconsistencies that may arise (Yu et al., 2020; Yu et al., 2021). The authors
of PolyCut (Livesu et al., 2013) listed sufficient topological conditions for
the existence of a valid polycube (following (Eppstein and Mumford, 2010)) and
– starting from naive labeling – systematically explore the space of labelings
in order to find the closest label assignment that satisfies all such
criteria. Also this method is not guaranteed to converge to a valid polycube
structure. Interactive tools for user assisted polycube construction also
exist (Li et al., 2021; Yu et al., 2020; Yu et al., 2021).

Figure 17. Top: a set of pathological shapes having sharp creases and surface
normals that push polycube deformation energies towards the generation of
globally inconsistent orthogonal polyhedra. Bottom: the associated polycubes.
Additional corners were strategically introduced to define topologically sound
global structures. Image from (Sokolov, 2016).

##### Mapping.

The volumetric map is obtained using dedicated deformation energies that
iteratively rotate surface normals until they _snap_ to the global axes. These
methods may input a pre-computed polycube segmentation (Gregson et al., 2011;
Livesu et al., 2013), freely deform the shape until the polycube structure
reveals itself (Huang et al., 2014; Fang et al., 2016), or interleave the two
operations, updating the reference labeling after each iteration (Fu et al.,
2016). Techniques that are not driven by a guiding (valid) polycube labeling
seldom converge to a topologically consistent polyhedron. As reported in
multiple articles, these pipelines include heuristic post-processing steps
that aim to remove topological artifacts (e.g., removing facets with less than
4 sides). In (2015) Sokolov and Ray observe that degenerate conditions are
both local and global, hence difficult to sanitize. In particular, since axis-
aligned features are naturally preserved during deformation, the presence of
long and slightly diagonal creases may push the deformation in the wrong
direction, resulting in globally inconsistent configurations which are
extremely hard to recover from (Fig. 17). Tiny features such as protuberances,
tunnels and handles are also critical, and may become intrinsically impossible
to resolve if the mesh lacks the sufficient resolution to define a valid
labeling.

Polycube deformation operates on a supporting tetrahedralization of the
object. Early deformation energies defined in (Gregson et al., 2011; Huang et
al., 2014) did not sufficiently penalize distorted and flipped elements,
producing non locally injective maps containing various foldovers, especially
in the vicinity of concavities. Fu and colleagues (Fu et al., 2016) introduced
a deformation energy that incorporates the AMIPS term (Fu et al., 2015), which
grows to infinity in the presence of degenerate or inverted elements. Various
similarly flip-preventing energies have been introduced in recent years (Fu
and Liu, 2016; Rabinovich et al., 2017), and could be adopted in this setting.
It is important to note, however, that the prevention of flips may reduce the
deformation space to an extent that no map that respects the boundary-
alignment constraints can be found. Alternatively to volumetric deformation,
one could in principle use a surface method to define a polycube-surface map
(e.g., with (Yang et al., 2019)) and then solve for a compatible volumetric
mapping between the two shapes. Again, however, despite the high level of
practical robustness showcased by recent approaches (Garanzha et al., 2021; Du
et al., 2020), the fully reliable automatic generation of constrained
volumetric maps without flips remains an open problem (Fu et al., 2021)
Motivated by this difficulty, an interactive polycube construction pipeline
that puts the user in the loop has been recently proposed (Li et al., 2021).
Users are allowed extensive control over each stage, such as editing the
polycube structure, positioning vertices, and exploring the trade-off among
competing quality metrics, while also providing automatic alternatives. The
robust mapping energy proposed in (Garanzha et al., 2021) is internally used
to discourage the generation of flipped elements.
The use of alternative mesh representations has also proved useful to robustly
construct volumetric mappings. In (Paillé et al., 2015) the authors represent
a tetrahedral mesh as a collection of dihedral angles, and propose a robust
spectral reconstruction method to generate an explicit mesh up to a global
similarity transformation. The use of reduced coordinates to represent and
manipulate meshes (e.g. via curvature or edge lengths) is a broad topic and
has been widely studied, especially for the surface case (Crane et al., 2011;
Campen et al., 2021). Specifically, the aforementioned paper shows that any
input polycube segmentation can be translated into a set of prescribed
dihedral angles that encode the change of normal orientation along the
surface. Using the proposed reconstruction method allows to convert such
angles into an explicit mesh, obtaining a locally injective polycube map.

##### Sampling.

The generation of the hex-mesh connectivity is based on a grid sampling in
polycube space. This is typically done by snapping polycube corners to integer
coordinates, and then sampling the space at the corresponding dense integer
lattice. Current methods assume an _ideal_ mapping that is virtually free from
distortion. In reality, hexahedra may undergo severe deformation, and
occasionally flip their orientation even if the underlying simplicial map is
injective. A valid trilinear mapping is guaranteed to exist, but it may
require a very dense polycube sampling. In practice, concrete implementations
often rely on a possibly imperfect coarse hex-mesh, and then apply untangling
techniques in case such mesh contains inverted elements (Sec. 6). No
guarantees of correctness can be provided in this case.

Regarding sampling frequency, the density of the integer grid must be fine
enough to ensure that no element of the polycube collapses. Cherchi and
colleagues (2016) showed that different samplings produce different meshes and
that a coarse sampling yields a mesh with a simpler block-structure, because
corners are more likely to snap to the same integer isolines. Finding a good
balance between these ingredients is subject to optimization (Cherchi et al.,
2016; Protais et al., 2020; Chen et al., 2019). Since snapping to coarser
grids involves bigger displacements, it becomes more complicated to preserve
the injectivity of the simplicial map. All existing methods are heuristic and
can break the map in complex cases.

Figure 18. Adaptively sampling a coarse polycube allows to restore major
features that are not explicitly encoded in parametric space (right), and that
could not be obtained with a regular sampling (left), not even with a dense
one (middle). Image from (Pitzalis et al., 2021).

While regular grids are widespread, smarter sampling schemes can be used to
control element size and anisotropy, or even to counterbalance map distortion.
In recent literature there have been a few attempts to address size control,
obtained either thickening a region of interest in polycube space prior
sampling (Xu et al., 2017) or adapting octree-based meshing to polycube space
(Hu and Zhang, 2016). More recently, Pitzalis and colleagues (2021) showed how
adaptive sampling can be unlinked from rigid octrees and extended to
generalized grids of any shape or topology. Adaptively sampled polycubes can
be used to improve geometric fidelity while keeping the mesh resolution low,
or it can also be used to restore major features that were missing in polycube
space and that could not be reproduced otherwise, not even with extremely
dense regular samplings (Fig. 18). These preliminary results suggest that a
tighter integration of adaptive sampling in polycube space could benefit the
whole pipeline. Current methods heavily rely on the ability of the polycube
finding module to catch _all_ the features of the object (at all scales) so as
to secure the proper mesh connectivity. Finding a better balance between
features that are explicitly encoded in the polycube and features that will be
reproduced with adaptive sampling, algorithms may be able to better distribute
the complexity throughout the whole pipeline, possibly increasing their
robustness. On the negative side, any non regular sampling will introduce
additional irregularities in the mesh, reducing its level of structural
regularity (Sec. 2.4) and possibly making it unusable for applications that
exploit the coarse block decomposition endowed in the mesh connectivity, such
as IGA methods (Hughes et al., 2005).

##### Features.

Desirable properties such as curvature and feature alignment depend on how the
polycube map orients these entities in $\mathbb{Z}^{3}$. In particular, since
sharp creases are preserved only if they map to integer isolines in polycube
space, there are intrinsic topological limitations to the class of feature
networks that can be correctly reproduced (e.g., a convex vertex with more
than three incoming feature lines cannot be correctly meshed). Since geometric
features are often characterized by surface normal discontinuities, labeling
methods such as (Livesu et al., 2013; Hu and Zhang, 2016) intrinsically
promote their positioning along polycube edges. Nevertheless, these methods do
not explicitly handle surface features, and may often fail to preserve them
(Guo et al., 2020). To our knowledge, the only method that explicitly promotes
feature alignment is CE-PolyCubeMaps (Guo et al., 2020). Given an input
network, the authors attempt to transform each feature into a piece-wise
linear curve that aligns with the global axes. Features that do not align (or
conflict with other features) are discarded; the others are included in the
polycube structure generation, with a feature-aware variant of PolyCut (Livesu
et al., 2013). Despite practically superior to previous approaches, also this
method does not provide strict guarantees. Furthermore, features are only
mapped to polycube edges, and the possibility to align to integer isolines
that are internal to polycube faces is not exploited.

##### Maturity.

Polycube methods have received increasing attention from the meshing community
and have now reached a discrete maturity level. The most recent algorithms
allow to blindly process datasets composed of more than a hundred shapes,
producing hex-meshes of good quality (Fu et al., 2016). In terms of mesh
structure, these methods typically produce valence semi-regular meshes, and
may occasionally produce semi-regular meshes if singularities (i.e., polycube
corners) align (Cherchi et al., 2016). The singular structure of a polycube-
based hex-mesh is fully exposed on the surface, and consists of all polycube
edges and corners. This inability to position singularities in the interior
inherently limits the map, and may occasionally be the source of unnecessary
distortion. A recent work of Guo and colleagues (2020) proposes to enhance the
singular structure with diagonal cut surfaces that penetrate the interior of
the polycube, permitting further distortion reduction. Intuitively, these cuts
can be thought of as analogous to cone singularities in surface mesh
parameterization (Soliman et al., 2018), although they are more constrained
because the two copies of each cut surface must still obey to the constrained
polycube structure. Another typical improvement consists in pushing the
external singular structure one layer inside the volume, adding a global
padding layer that avoids over-constrained hexahedra with more than one facet
exposed on the surface (see Sec. 5.4). More sophisticated padding schemes that
directly exploit the polycube map to optimally balance distortion with mesh
growth are also available (Cherchi et al., 2019a).

##### Abstract polycubes

Various techniques use skeletons or segmentation to partition a shape,
generating atlas of maps to face-adjacent cuboidal domains (Livesu et al.,
2016; Li et al., 2010, 2013; Liu et al., 2015; Usai et al., 2015). Some of
these authors explicitly refer to the structures they produce as _abstract_ or
_generalized_ polycubes, where these words mean that such structures may not
have an embedding in $\mathbb{Z}^{3}$. The major difference between embedded
and abstract polycubes is that in the former the whole polycube is considered
as a global entity that lives in a unique space, whereas in the latter each
cuboidal domain lives in its own space and has a dedicated map. For this
reason, despite the name similarity, these techniques are more similar to
domain decomposition approaches, and are discussed in Sec. 4.5.

Figure 19. Frame field based hexahedral mesh generation: Given an input
tetrahedral mesh (left), first a boundary aligned and smooth frame field is
generated (middle, yellow). The frame field serves as a proxy for the local
orientation of hexahedra and thus enables the efficient generation of an
integer-grid map (middle, red), which induces a hexahedral mesh (right). Image
from (Liu et al., 2018).

###  4.8. Frame fields

Frame fields offer a promising research direction for general hexahedral mesh
generation. They can be seen as a continuous relaxation of an integer-grid map
(cf. Sec. 2.5) and thus a key component of the corresponding hex-mesh
generation algorithms. A prototypical algorithm (see Fig. 19) consists of
three major steps: (i) optimization of a boundary-aligned frame field, (ii)
generation of an integer-grid map that resembles the frame field (Nieser et
al., 2011), and (iii) extraction of the integer level-sets which form an
explicit hexahedral mesh (Lyon et al., 2016).

##### Frame Field.

A frame field can be seen as a generalization of a vector field to a quantity
that locally describes the shape of a (linearly deformed) cube. Locally, a
frame consists of three linearly independent vectors, which represent a
parallelepiped, i.e., the orientation and shape of a linearly deformed cube.
It is important to understand that globally a frame field is significantly
more complex than three superimposed vector fields since it can contain
singularities where topologically the vector fields are nontrivially
interconnected on a branched covering (Nieser et al., 2011). As a consequence,
the connection induced by a frame field might exhibit nonzero monodromy, i.e.,
a vector does not return to itself when transported along a cycle around a
singularity.

####  4.8.1. Local frame representation.

There are various different representations to locally encode a frame. A
straightforward choice that is capable of fully describing the shape of a
linearly deformed hexahedron are three explicit vectors
$u,v,w\in\mathbb{R}^{3}$ bundled into a matrix
$F=\left(u,v,w\right)\in\mathbb{R}^{3\times 3}$. The local shape of the
hexahedron then simply corresponds to the parallelepiped formed by $u$, $v$
and $w$. In practice, often reduced representations are preferable, where $F$
is orthonormal or orthogonal, meaning that only rotations or rotations and
scaling along the principal axes are encoded. For example, unit quaternions
have been used as a compact representation for orthonormal frames (Gao et al.,
2017b; Liu et al., 2018). Another common choice for orthonormal frames are
Euler angles w.r.t. either a global coordinate system (Huang et al., 2011), or
alternatively a local coordinate system (Palmer et al., 2020; Ray et al.,
2016) to avoid gimbal locks.

##### Handling Cube Symmetries.

The explicit tangent vector representation of frames has one major
disadvantage; it is not unique. For example, the matrix $(u,v,w)$ describes an
identical parallelepiped as the matrix $(v,w,-u)$. Since there are $6$
potential permutations of the three vectors and $2^{3}$ potential choices of
sign, in total, there are $48$ different matrices, which encode a single
frame. Formally, equivalence is established by the binary octahedral group
$BO$ with 48 symmetry transformations. Since the octahedron is dual to the
cube, their symmetry transformations are identical. By fixing the orientation
(the sign of the determinant of $F$), it is possible to reduce the elements in
one equivalence class to $24$ elements with octahedral symmetry. This explains
why frame fields are often called _octahedral fields_ (Solomon et al., 2017),
while _3D cross fields_ is another common name. The non-uniqueness of the
$(u,v,w)$-representation significantly complicates the optimization of frame
fields by inducing discrete variables with values from the octahedral group
$O$. A remedy is offered by alternative representations that are equipped with
a build-in symmetry invariance and thus a unique representation of frames.

The seminal work of Huang and colleagues (2011), expressed orthonormal frames
as rotations of the polynomial $x^{4}+y^{4}+z^{4}$, which are by construction
invariant under transformations by elements of $O$. By restricting the
polynomial to a sphere, it can be expressed in the spherical harmonics basis,
lifting a single orthonormal frame to a $9$-dimensional representation vector.
Since rotations only possess three degrees of freedom, it is clear that the
spherical harmonics representation is a relaxation, i.e. not all
$9$-dimensional vectors correspond to a rotation of the polynomial
$x^{4}+y^{4}+z^{4}$. An identical representation can be derived from the
perspective of 4th-order symmetric tensors (Chemin et al., 2018; Golovaty et
al., 2019), which have been further generalized to a 15-dimensional
representation of orthogonal frames offering independent scaling of axes
(Palmer et al., 2020).

##### Differential Frame Representation.

Instead of representing a frame field by pointwise specification of frames,
one valuable alternative consists in encoding its derivative, i.e. the change
of frames. After specifying one frame in the domain, the entire frame-field
can then be re-constructed by integration. Such integration is path-
independent given that the specified derivative is integrable, i.e. all
fundamental monodromies are elements of the octahedral group when expressed in
the coordinate system of the frame itself. Corman and Crane (2019) employ such
a differential representation in a frame-field optimization setting. They
extend the theory of moving frames to frame-fields with cube symmetry.
Conceptually, the setting is analog to the 2D setting of trivial connections
(Crane et al., 2010). In both cases, all singularities and thus all
fundamental monodromies are specified as input. However, while trivial
connections in 2D are generated with a simple linear solve, the non-
commutativity of 3D rotations requires a (continuous) non-linear least-squares
optimization. Interestingly, despite the non-convex objective function,
experiments suggest that the resulting field is independent from the chosen
initial configuration. Leveraging a differential frame-field representation
for optimizing fields with unconstrained singularities has not been done so
far but is an interesting direction for future work. It would require
replacing the fixed monodromies by the feasible set of discrete choices from
the octahedral group.

####  4.8.2. Frame field optimization.

The optimization problem usually consists in finding the best frame field,
which aligns to the boundary of the domain. Best often means as-smooth-as-
possible and is measured with the Dirichlet energy
$\int_{\Omega}||\nabla\phi||^{2}dx$ where $\phi$ is a frame representation.
Some earlier methods constrain the complete boundary field to a pre-computed
solution (Li et al., 2012; Kowalski et al., 2016). However, in general, it is
preferable to only require boundary alignment and let the rest of the field
emerge freely. Typically, the domain is discretized into a tetrahedral mesh,
where frames are located at vertices (Ray et al., 2016; Palmer et al., 2020),
faces (Huang et al., 2011), or cells (Liu et al., 2018). Alternatively,
boundary element discretizations have been explored (Solomon et al., 2017),
where a triangulation of the boundary is sufficient. Algorithms that are based
on a lifted frame representation ($9$\- or $15$-dimensional) need to ensure
that they do not leave the sub-manifold of frames. This is done either with a
projection operator, or with some kind of manifold optimization. More
specifically, Huang et al. (2011) compute an initial field by optimizing a
convex relaxation of the actual problem, i.e., minimization of the Dirichlet
energy in $\mathbb{R}^{9}$, followed by a local projection onto closest
frames. The boundary alignment of the field is approximated by a single linear
constraint per boundary face. The initial field is further improved by a
nonlinear optimization restricted to the frame-manifold via Euler angles. Ray
et al. (2016) follow a very similar strategy but discretize the field on
vertices, tighten the boundary constraints and improve the performance of the
projection. Palmer et al. (2020) observed that a modified
Merriman$-$Bence$-$Osher (MBO) algorithm is beneficial because it is often
able to avoid local minima that induce global inconsistencies in the
singularity graph. The MBO algorithm alternatively diffuses the $9$\- or
$15$-dimensional coefficient space and locally projects the coefficients onto
frames. The idea of the modified MBO algorithm consists in starting with a
large diffusion kernel and iteratively shrinking it in subsequent steps. The
rationale behind this strategy is that large diffusion steps sufficiently
leave the (non-convex) manifold of frames and thus avoid local minima, while
small diffusion steps are required for the accuracy of the solution. The
projection of frames can be done approximately with gradient descent (Ray et
al., 2016), or exactly with a semidefinite relaxation (Palmer et al., 2020).
More rapid convergence than the MBO algorithm is offered by Riemannian trust-
region manifold optimization (Palmer et al., 2020), which on the downside has
a higher risk of getting trapped in local minima.

####  4.8.3. Generality

Besides providing directional guidance for hexahedral mesh elements, a key
property of a frame field is its network of singularities – which one commonly
aims to adopt for a hexahedral mesh generated based on the frame field. In
this context, frame fields are general enough that the singularity network of
any hexahedral mesh can be expressed. This means that, in contrast to many
other approaches, e.g. polycube mapping, sweeping, or grid-based approaches,
the output is not a priori restricted to a subclass of hexahedral meshes.
Therefore superior mesh quality can be achieved, specifically if complex
feature alignment is required. In particular, frame field-based methods are
able to express alignment not only to the boundary of a domain but also to
arbitrary internal structures, which is, for example, important in multi-
material applications or in the simulation of fluid-structure interaction.

The main drawback, on the other hand, is the fact that frame fields are
actually over-general for the purpose of mesh generation: Frame fields may
exhibit additional types of singularities that cannot occur in hexahedral
meshes, cf. (Viertel et al., 2016; Liu et al., 2018). Such singularity
configurations are said to be “non-meshable”. A key example are 3-5
singularities (Reberol et al., 2019), which frequently appear in smooth frame
fields but are not meshable. Existing approaches are able to automatically
repair some locally non-meshable configurations (Li et al., 2012; Jiang et
al., 2014) or involve the user to manually repair the singularity graph (Liu
et al., 2018) and then generate a frame field with a prescribed singularity
network (Liu et al., 2018; Corman and Crane, 2019). Another option is to
optimize a general frame-field in such a way that all singularities are pushed
towards the boundary in order to generate a generalized polycube as done in
(Fang et al., 2016). Additional research is required to enable complete repair
or to restrict the frame field generation and optimization to the space of
meshable configurations in the first place. Little can be learned in this
regard from the analogous 2D problem, as the gap, in terms of singularity
structure, between 2D frame fields and quad-meshes is significantly smaller.

####  4.8.4. Field-guided integer-grid map

Given a (meshable) frame field, one then aims to conceptually _integrate_ it
to obtain a parametrization, a map (in particular an integer-grid map) onto
part of $\mathbb{R}^{3}$. As the frame field typically is not integrable, a
map whose isocurves are precisely aligned with the frame field’s directions
does not exist. Approximate alignment, e.g. least-squares alignment, is thus
aimed for, as in the Poisson approach described by (Nieser et al., 2011),
generalizing the cross-field guided mapping used in the 2D case for
quadrilateral meshing (Kälberer et al., 2007; Bommes et al., 2013a). The frame
field’s singularities are adopted in this process and define the implied
hexahedral mesh’s singularity structure.

Unfortunately, this field-guided mapping approach does not guarantee a valid
resulting map without flips. While there are heuristics (Lyon et al., 2016) to
recover a valid hexahedral mesh even from some invalid integer-grid maps with
flips, no general guarantees are available. In the 2D setting, aiming at
quadrilateral mesh generation, a stream of recent work has shown ways to
reliably generate flip-free maps with prescribed singularities (for instance
implied by frame fields) and boundary alignment (Campen and Zorin, 2017;
Campen et al., 2019; Campen et al., 2021; Gillespie et al., 2021). It is based
on phrasing the problem as a constrained metric computation problem; in a
specific discrete conformal setting and formulated in per-vertex scale
variables, the problem becomes convex and can be solved reliably. Most
importantly, these methods employ on-demand mesh refinement (or modification)
to ensure feasibility, in the sense that the mesh offers sufficient degrees of
freedom to support a valid map, represented in a piecewise-linear manner.
Note, however, that field guidance is considered in these works only in the
form of adopting the singularities, not in the form of dedicatedly following
the field’s directions. Generalization of this general approach to the 3D
setting is not straightforward; the space of 3D conformal maps too restricted
to be useful. The work by (Paillé et al., 2015) may be viewed as a first step:
It describes the representation and optimization of discrete metrics on 3D
tetrahedral meshes in intrinsic variables (dihedral angles). Using linear
constraints in these variables, boundary alignment and singularities can
directly be prescribed, for instance those adopted from an optimized 3D frame
field. However, the resulting problem in this 3D case is non-convex, and the
issue of potentially required mesh refinement is unsolved. While there are
known ways to reliably generate flip-free maps in 3D (Campen et al., 2016),
they do not support the prescription of arbitrary singularities.

An alternative reliable approach proposed for the 2D setting is based on
decomposing the domain into regular pieces based on stream lines of a 2D frame
field (Myles et al., 2014). Also this does not generalize to a 3D stream
surface based approach (Kowalski et al., 2016) with similar guarantees.

###  4.9. Hex-dominant meshing

Automatic methods for all-hex meshing are only applicable to a subset of all
the possible inputs. In contrast, the grid-based approaches (Sec. 4.6) can
operate on intricate shapes and guarantee all-hex meshes; unfortunately, they
produce inferior quality results. High-quality, feature-aligned all-hex meshes
are still elusive, so industry still relies on semi-manual block
decomposition, a time-consuming process (Lu et al., 2017).

For this reason, other methods focus on the hexahedral-dominant meshing
instead of full-hex, aiming to reach the highest possible proportion of
hexahedra. Hexahedral-dominant meshing is a relaxation of the problem to
significantly improve robustness at the cost of introducing a small number of
generic polyhedra. The generation of hex-dominant meshes boosted the use of
those datasets in practical contexts such as FEM (Wicke et al., 2007).
Moreover, the recent advancement in the construction of higher-order bases
(Schneider et al., 2019a) may foster the adoption of hex-dominant meshes in
the mechanical analysis.

The first approach to produce hex-dominant meshes agglomerates neighboring
tetrahedrons to assemble hexahedral cells. The problem of finding a globally
optimal solution is NP-complete; hence the clustering process is usually
driven by local heuristics. Meshkat and colleagues (2000) have proposed the
first method following this idea. The clustering process relies on an
undirected graph representing tetrahedra and their connectivity. The graph is
enriched with particular arcs and labels on nodes to calculate the
agglomeration heuristic. Given the initial tetrahedral mesh, the algorithm
detects and replaces subgraphs with hexahedral nodes.

The method proposed by Yamakawa et al. (2002) takes as input a general 3D mesh
and distributes a set of nodes into the volume by the physical simulation of
crystal pattern formation. Then nodes are used to produce a mesh composed of
hexes, prisms, and tets, with around 50% hexahedral cells. This method allows
controlling element size and primary orientation. Since it does not require a
tetrahedral mesh as input, this method is less sensitive to the input
discretization than the approach of Meshkat and Talmor (2000).

Vyas and Shimada (2009) proposed a more sophisticated method that starts by
generating a volumetric tensor field to specify the anisotropy and
directionality of the elements. Then, such a field induces an advancing front
process where several hexahedral fronts contribute to cover the entire volume.

Lévy and Liu (2010) generalized Centroidal Voronoi Tessellation (Faber and
Gunzburger, 1999) for hex-dominant meshing, introducing Lp-Centroidal Voronoi
Tessellation. Unlike the standard Voronoi diagram, Lp-CVT favors the formation
of cubical cells by using a distance metric that takes into account a
predefined background tensor field. The resulting method excels in robustness
and controllability.

Despite the considerable advancements in performances, the methods mentioned
above cannot obtain a high hexa ratio for the general case. The approach
proposed by Sokolov et al. [2016a] produces higher hex ratios by using a
guiding frame field. A sampling process generates a point set organized as a
regular grid and locally aligned with the frame field. A constrained Delaunay
triangulation of the volumetric samples makes a tetrahedral mesh and is
finally clustered into hex elements. They obtain hex-dominant meshes with up
to 95% hexahedral cells (although in the worst case they earned less than
30%). Despite the result’s quality, this method might produce non-conforming
meshes containing configurations where a quadrilateral face is adjacent to two
separate triangular faces. This issue has been solved by (Ray et al., 2018).

The approach by Pellerin et al. (2017) explores the space of all possible
agglomerations of tets. Then a greedy process selects the configurations to
agglomerate tets into hexes. Their approach produces hex-dominant meshes with
a 60% ratio of cuboidal elements across all shown examples.

Figure 20. Some challenging examples of hex-dominant re-meshing using (Gao et
al., 2017b).

Gao et al. (2017b) directly generate conforming hybrid meshes using polyhedral
agglomeration. This method starts from tetrahedral mesh obtained by sampling a
guiding frame-field. An iterative process modifies the connectivity utilizing
a set of local operators to compose hexahedral cells. The local operators
grant the conformity of the final mesh. While this method excels in robustness
(as demonstrated by the complex example shown in Fig. 20 ), it cannot control
the class of created polyhedrons (they might have up to 40 facets in some
cases, see Table 1 in (Gao et al., 2017b)).

Figure 21. Some CAD models remeshed by (Livesu et al., 2020). In this case the
meshes are hex-only.

The recent method proposed by Livesu et al. (2020) produces _strongly hex-
dominant meshes_ , conforming meshes with less than 2% non-hexahedral cells.
In most cases (76% of the models tested), this method derives pure hex-meshes.
It mimics manual block decomposition. It extracts first a set of well-
distributed loops on the surface following a feature-aligned cross-field. Each
loop defines a cutting surface that decomposes the volume into simpler
polyhedral blocks. The cutting surfaces are added one by one until the quality
requirements of the polyhedral blocks are satisfied. These blocks are finally
converted into hex-dominant mesh via midpoint subdivision. As shown in Fig.
21, this method excels in the preservation of sharp features, which are
directly incorporated in the output connectivity.

Similarly to LoopyCuts (Livesu et al., 2020), HexDom (Yu et al., 2021) produce
a block decomposition where each block is either hex, prism, or tetrahedral
cell. This method extends the approach for polycube generation proposed in (Hu
and Zhang, 2016) based Voronoi tessellation (CVT) to include non-hex elements.
The cells are embedded in 3D using a variation of (Yu et al., 2020). The
segmentation and polycube definition process requires manual work. The methods
proposed in (Zhan et al., 1803) and the one used by the commercial package
Cubit (Meyers and Tautges, 1998) uses advancing front approaches to produce
meshes composed of hexahedral and tetrahedral elements. To improve the quality
of the final mesh Cubit use some sophisticated cleanup operation based on
connectivity editing and geometric measures.

The recent approach proposed by Bukenberger and colleagues (Bukenberger et
al., 2021) generate _At-Most-Hexa Meshes_. At-Most-Hexa Meshes are meshes
composed mainly of hexahedral elements, where no cell has more than six faces,
and no boundary face has more than four sides. Similarly to tetrahedral and
hexahedral meshes, the volume of each cell can be defined by trilinear
interpolation from its corners. At-Most-Hexa Meshes meshes are generated by
extending to the volume the 2D approaches that use Lloyd relaxation with non-
euclidean distance measures (Hausner, 2001). Using the $L_{\infty}$ norm
(instead of the simple euclidean distance), the cells emerging from the
volumetric Lloyd relaxation process become more cubical, converging to a hex-
dominant mesh. Similarly to most of the meshing methods based on Lloyd
relaxation, this method is very permissive on the required input. It works on
point clouds, triangular meshes and can be guided by an input orientation
field if available.

The hex-dominant mesh allows sufficient degrees of freedom to adapt the grid-
based methods to conform to sharp features. Trimmed hexahedral meshes are
created by intersecting a grid with a closed surface. Non-hexahedral elements
emerge along the surface where the surface is not aligned with the edges of
the grid. The technique recently proposed by Kim and colleagues (Kim and Kim,
2021) extends the trimmed hexahedral methods by creating particular vertices
where sharp features intersect with the grid. A feature simplification schema
is used when multiple features are concentrated in the same cell.

Another class of methods transforms hex-dominant meshes to increase the number
of hexahedral elements in the mesh. The approach proposed in (Yamakawa and
Shimada, 2003) increases the number of hexahedral and prism elements by
applying sequences of local operations that modify the connectivity.
Unfortunately, this method can generate non-conformal meshes. Instead, HexHoop
(Yamakawa and Shimada, 2002) converts a mesh composed of hexahedrons, prism
and tetrahedrons into a conformal pure-hex mesh. The conversion process is
based on the local application of two particular refinement schema, called
core and caps. Unfortunately, this method tends to insert a high number of
irregular vertices deteriorating the regularity of the tessellation.

Table 2. Summary of the main properties for each class of hex-meshing algorithms reported in Sec. 4. Some of the columns in this table correspond to items also listed in (Blacker, 2000). As a rough indicator for the extent of (ongoing) research activity, we list the number of overall works and recent (published within the last 5 years) works (referenced in this survey) dealing with each class. Method |  Type |  User interaction |  Shape class |  Feature preserv. |  Size contr. |  Mesh structure |  Element quality |  Robustness |  Orient. sensitive |  Total works |  Recent works |  Open problems
---|---|---|---|---|---|---|---|---|---|---|---|---
Advancing front |  Direct |  Automatic |  CAD oriented |  Surface features only |  No |  Unstructured |  Good at border, poorer inside |  Poor |  No |  7 |  0 |  Improve handling colliding fronts, complex topologies
Dual methods |  Both |  Automatic, semi-automatic |  CAD oriented |  Surface features only |  No |  Unstructured, semi-structured |  Good at border, poorer inside |  Poor |  No |  13 |  2 |  Robust handling of self- intersecting sheets
Sweeping, decomp. |  Both |  Semi-automatic |  CAD oriented |  Surface features only |  No |  Semi-structured |  Good |  Good (manual) |  No |  35 |  11 |  Automatic definition of sweepable sub-volumes
Grid based |  Indirect |  Automatic |  Any shape |  Yes, (limited valence) |  Yes |  Severely unstructured |  Poor at border, optimal inside |  Great (commercial product, demonstrated on many datasets) |  Yes |  18 |  3 |  Feature preservation, mesh size, mapping
Polycube maps |  Indirect |  Automatic, semi-automatic |  Any shape |  Yes, (limited valence) |  Yes |  Valence semi-structured |  Good (depends on map) |  Good (demonstrated on medium datasets) |  Yes |  18 |  8 |  Polycube topology, mapping, feature preservation
Frame fields |  Indirect |  Automatic, manual fixing |  Any shape |  Yes |  Yes |  Valence semi-structured |  Good (depends on map) |  Poor |  No |  17 |  12 |  Generation of hexable fields, field aligned mapping
Hex-dominant |  Both |  Automatic |  Any shape |  Yes |  Yes |  Valence semi- structured, hybrid |  Often good |  Good (demonstrated on medium datasets) |  No |  20 |  11 |  Hybrid elements (topological control, amount, quality)

##  5\. Topological operators

Scientific computing often demands to edit a given hexahedral mesh, e.g. to
improve the accuracy of a solution with a posterior refinement (Shen et al.,
2022), to generate boundary hex layers for CFD applications (Reberol et al.,
2021), to ensure mesh conformity across surface membranes (Staten et al.,
2010b) or to constrain the mesh size (Maréchal, 2009). Unlike tetrahedral
meshes, where changes of the mesh connectivity always have a local footprint,
editing the topology of a hexahedral mesh is often a global operation. This
makes hexmesh editing significantly more difficult than tetmesh editing. In
this section, we revise the most prominent operators for editing the topology
of a hexahedral mesh.

###  5.1. Sheet operators

As shown in Sec. 2.2, the dual of a hexahedral mesh is a simple arrangement of
surfaces. Each surface, i.e., a sheet, corresponds to a layer of hexes in the
primal mesh and two surfaces intersect along a chord, which is a column of
hexes in the primal mesh. Sheet operators consist in inserting or removing a
complete sheet or chord from the mesh. They are used to refine (Ko-Foa Tchon
and Camarero, 2002; Parrish et al., 2007) or coarsen meshes (Benzley et al.,
2005a; Shepherd et al., 2010), to capture analytic features (Merkley et al.,
2007), or to make conforming meshes involved in the assembly of parts. Those
parts can result from a volume decomposition during a user-assisted meshing
process (Borden et al., 2002b; Jankovich et al., 1999) or can correspond to
two adjacent models sharing contact surfaces (Staten et al., 2010b).

Using sheet operations to refine a mesh mainly involves inserting sheets,
which is quite easy to control if you avoid self-intersecting and self-