Indirect all-quadrilateral meshing based on bipartite topological labeling

Quadrilateral meshes offer certain advantages compared to triangular ones, such as reduced number of elements and alignment with problem-specific directions. We present a pipeline for the generation of quadrilateral meshes on complex geometries. It is based on two key components: robust surface meshing and efficient indirect conversion of a triangular mesh to an all-quad one. The input is a valid geometric surface mesh, i.e., a triangulation that accurately represents the geometry of the model. A right-angled triangular surface mesh is initially created by continuously modifying the input mesh while always preserving its topological validity. The main advantages of our local mesh modification-based approach are to (i) allow the generation of meshes that are globally aligned with a given direction field and (ii) to reliably handle non-manifold feature edges (in multi-volume models) and small features. The final quadrilateral mesh is obtained by merging pairs of triangles into quadrilaterals. We develop a novel bipartite labeling scheme in order to identify and correct inconsistent pairings. The procedure is based on local operations and is much more efficient than previous global strategies for tri-to-quad conversion. The whole pipeline is tested on a large number of models with diverse characteristics.


Introduction
Quadrilateral meshes are often preferable to triangular ones for numerical simulations. They have fewer elements for the same number of vertices, they are ideally capable of providing a block-structure and they can provide better alignment to geometric features, as well as to problem-specific features, providing better numerical behavior for specific physical phenomena (a typical example is the demand for structured and aligned boundary layers in Computational Fluid Dynamics). Yet, the automatic generation of quadrilateral meshes is still regarded as a challenging problem in mesh generation. Even though a lot of different approaches exist, there is not to date a conclusive method, analogous to triangular meshing which is considered highly mature and for which there exist robust algorithms based on strong mathematical foundations.
The purpose of this work is to address the problem of generating quadrilateral meshes for complex 3D models. We strive for generality in our approach; our input is simply a triangulation of the model. The triangulation can be an STL representation of the geometry, triangulation of scanned data or a mesh generated from a CAD model with standard meshing techniques (Fig. 1). The input meshes may be of bad quality and contain non-manifold feature edges. Our goal is to design a pipeline satisfying the following design goals: (i) robustness, i.e., guarantee of termination regardless of the complexity/bad quality of input data, (ii) feature preservation, i.e., the persistence of user-defined internal and boundary curves, (iii) high element quality, and (iv) efficiency, i.e., providing a quad meshing algorithm with a running time comparable to or faster than conventional trito-quad technologies. Maxence Reberol and Jean-François Remacle are contributing authors.

3
To generate a high-quality unstructured quadrilateral mesh that preserves user-defined features, we develop three independent steps: (i) we sample points on the model curves and surfaces, guided from a metric field (cross-field and an associated size-map), (ii) we use local mesh modifications to continuously remove points of the initial triangulation and add the new ones, leading to a right-angled triangular mesh and (iii) we convert the triangular mesh into an allquad mesh using a novel approach based on a topological labeling scheme.
When it comes to mesh generation, the expected reliability rate of industrial-grade algorithms is essentially 100%. The input data supplied to our meshing algorithms are often noisy. Inputs may be huge, with a wide range of scales. It is not surprising that the most complex/tricky part of our work is related to robustness. We have given special attention to ensure that this mesh generator provides results regardless of the complexity of the input data, as long as it is correct (i.e., a watertight but possibly non-manifold triangulation, in the sense that no folded elements and edge intersections may be present).
In order to ensure reproducibility of this work, the whole implementation will be available in Gmsh, the open-source mesh generator [1]. To demonstrate the robustness of our algorithm, we applied it to a large number of models found in various datasets.

Related work
Surface meshing Surface mesh generation poses various difficulties related to robustness and efficiency. Of the several methods proposed in the bibliography, we can identify two main categories [2,3]: (i) Parametric approaches, where the surface mesh is generated in the parametric space, and (ii) Nonparametric (direct) approaches, where the surface mesh is generated directly in the 3D space.
In parametric approaches, the 3D surface is mapped to a 2D parametric space [2,[4][5][6][7]. Since the CAD surfaces (typically NURBS patches) have underlying u, v representation, it can be efficient to generate a mesh in the plane with standard meshing techniques and afterward map it back to the 3D space. Generating planar meshes in the parameter space is a robust approach that is usually able to provide high-quality meshes. Yet, approaches that use the parameter plane are able to consider surfaces that are isomorphic to a punctured disk. Meshing complex models with a parameter space approach do not allow to globally align a mesh with a cross-field, since each discrete patch of a CAD model may be equipped with an independent parametrization and the feature edges that separate those patches are not necessary aligned with the cross-field. Parametrization techniques can also be used to remesh triangulations [8,9].
Nonparametric (also referred to as direct surface meshing in the literature, creating some ambiguity with the terms of direct/indirect approaches commonly used in quad meshing) approaches to surface meshing can be based on quad trees [10,11], advancing front [12,13] or Delaunay strategies [14,15]. One of the main advantages of the direct approach is its usefulness for models where an underlying parametrization is not available or when it is degraded. Local mesh modification strategies to remesh models described by STL triangulations are proposed in Refs. [16,17]. One of the main difficulties of this class of methods is that the manipulation of geometry directly on the 3D space is a challenging task that may lead to geometric or topological ambiguities.
Quadrilateral meshing Initial efforts to automatically generate quadrilateral meshes include grid-based and paving algorithms. Grid-based methods start with the generation of a background Cartesian or a quad-tree grid with the Fig. 1 Quadrilateral meshes produced by our algorithm with input a CAD model (left) and an STL triangulation (right) subsequent snapping of elements to the domain boundaries [18][19][20]. The paving algorithm, first introduced by Ref.
[21] generates quadrilateral elements in an advancing-front fashion, propagating from the boundary to the interior. Both classes of methods suffer from a degraded quadrilateral quality and high node irregularity on specific regions of the domain: the latter on the domain boundaries and the former on the front collisions on the interior. In Ref. [22], a bichromatic Delaunay quadrangulation method is presented, with our current work building upon a similar concept.
On the other hand, quad conversion or indirect methods are based on the merging of pairs of adjacent triangles of an input mesh to quadrilaterals [23,24]. In Q-Morph [25], triangles are transformed into quadrilaterals with an advancing-front algorithm. Blossom-Quad algorithm [26] computes a perfect matching to optimally pair triangles, while [27,28] produce meshes better suited for triangle pairing by generating aligned right-angled triangles.
Cross-fields Cross-fields are nowadays commonly used in the context of mesh generation, a line of research stemming from the computer graphics community and global parametrization methods [29]. For quad/hex meshing, crossfields define preferred orthogonal directions on the domain to guide creation of optimal elements [30][31][32]. Cross-fields should be as smooth as possible (except at singular points) and aligned with the boundary of the domain.

Contributions
The quad meshing pipeline that is proposed follows a modular approach, with each of the steps being an independent algorithm that can be re-used in various situations (Fig. 2). As stated before, our main concern is to provide a reliable solution, i.e., we want the quad meshing pipeline to be resilient to complex/ill-conditioned inputs.
The two main contributions of this work are: 1. Robust surface meshing. In our procedure, we follow the idea of Ref. [28] of separating the generation of points and the creation of the elements. The main novelty of this work is our direct approach. In this work, an input triangular mesh is continuously modified through robust local mesh modifications. The word continuous is chosen on purpose: each local mesh modification that is performed guarantees the topological integrity of the current triangular mesh. At the end of the remeshing process, most of the points of the initial triangular mesh, and in most cases all those points are removed from the triangulation and replaced by the ones created to accommodate the cross-field and size field characteristics. 2. Straightforward and efficient all-quad meshing. We propose a bipartite labeling scheme that propagates topology information on the vertices during point gen-eration. By using this information on the right-angled triangular mesh, we are able to optimally place new Steiner points to fix topological inconsistencies (oddbounded regions on the interior) and recover a bipartite, all-quad mesh. The approach is very efficient since it converts a global problem to a local one.

Overview
Our algorithm takes as input a watertight, possibly non-manifold triangulation T 0 . The input triangulation T 0 is classified: the word classified indicates the fact that the triangles are grouped into colors and interfaces between colors are considered as feature edges that must be conserved in the remeshing process. Internal feature edges may also exist that lie inside a group of triangles with the same color. Note that feature edges can also be detected based on the dihedral angle, given a user-defined threshold. We take into account two special categories of feature edges in T 0 : (i) non-manifold feature edges with two or more adjacent surfaces and (ii) internal (embedded) feature edges inside surfaces (e.g., a crack for solid modeling) All the necessary topological information is available in the initial mesh. Surfaces are bounded by closed feature edges and those feature edges by feature points. Consecutive feature edges form feature curves. Feature curves must be preserved during the continuous mesh modification process. An important property of our method is the ability to handle the model as a whole and thus take advantage of the global nature of the guiding cross-field. We do not follow a patchwise approach where we handle each surface independently, followed by a connection of curves to ensure conformity. A half-edge data structure [33] is used to get connectivity information, and an array of boundary edges that define feature curves is stored. We extend this data structure by storing the type of boundary edge (manifold or non-manifold), along with the triangles connected to each (one triangle for open boundary edges, two triangles for manifold edges, and a larger than two number for non-manifold ones). This feature enables us to efficiently treat boundaries during surface meshing.
The scheme provides the flexibility to preserve the topological characteristics of the input mesh, such as 'hard' edges or user-defined feature curves, without relying on extensive a priori knowledge of domain characteristics or a complicated feature recognition preprocess [34][35][36][37]. Furthermore, by utilizing an appropriate data structure on top of the halfedge one, we can handle non-manifold configurations that may occur in industrial multi-volume CAD models.
The input triangulation T 0 is the geometric model. Two other inputs are required for running the algorithm: (i) a unit cross-field and (ii) a size field h( ) that are both used for guiding the point insertion process. A cross-field is a field defined on a surface S with values in the quotient space S 1 ∕Q , where S 1 is the circle group and Q is the group of quadrilateral symmetry. It associates to each point of a surface S to be meshed a cross made of two unit vectors orthogonal to one another in the tangent plane of the surface and their opposites (Fig. 3a). Although T 0 , and h can be independent, it is beneficial to have (i) a cross-field that is aligned with the feature edges of T 0 and (ii) a size field h that takes into account both the local change of direction of the cross-field and small features of the geometry. In this work, and h are precomputed using T 0 as support with the algorithms described in Ref. [38].
The transformation of the triangulation T 0 into the final quad mesh T q is done in three sequential steps: 1. Point generation (Sect. 3) Points of the final mesh T q are generated in a frontal fashion starting from the feature curves and guided by both the cross and the size field. This set of points is embedded on the triangles of the base/initial mesh T 0 . 2. Point replacement (Sect. 4) The initial mesh T 0 is continuously transformed into another triangular mesh T 1 by connecting the newly generated points on T 0 and subsequently removing the initial mesh points, utilizing local mesh modifications. 3. From triangles to quads (Sect. 5) Mesh T 1 is transformed into a quad mesh by combining pairs of triangles. Using a binary labeling scheme for the points during frontal generation allows us to instantly extract a valid all-quad topology.
The simple model of Fig. 3 is used as an example to describe those three steps.

Point generation
In this work, we take a standard surface-to-volume point of view of mesh generation on Gmsh [1] which essentially consists of a bottom-up procedure. Model curves are discretized at first. Mesh edges on the model curves are used as boundaries of model faces and mesh triangles on model faces are used as the boundary of model volumes if they exist.
A set P of points on surfaces is generated using a frontal point propagation algorithm that is similar to Ref. [28]. The main difference with Ref. [28] is that all the operations are performed here directly on T 0 without using a parametric space. The point sampling scheme has been implemented for the special case of the sphere [39]. It is extended here for general surfaces and we reiterate the whole process for completeness. Our frontal approach is enhanced with the use of a cross-field that allows to structure the quad mesh and a size field h that allows taking into account the various feature sizes of a model as well as changes of directions of ( Fig. 3a).

Curve point generation
Each feature curve is uniquely defined from a list of nonintersecting connected edges. Given a size field h defined by a value at each point, we mesh the discrete representation of each curve by following the general guidelines of Ref. [40]. This leads to the set of points It is important to note that at this step, we can control the generation to have an even number of points for each feature curve. This provides us with a topologically necessary condition for all-quadrilateral meshing.

Surface point generation
Starting now from the N gc generated points on mesh feature curves, we want to spawn a set of points P s = { i | i = 1, … , N gs } on the surfaces in the directions provided by the cross-field ( ) and with respect to the underlying size field h( ) . The point set P s , along P c , will be used to generate a right-angled triangulation T 1 that is well suited for combining triangles into quadrilaterals and form T q . The cross-field gives N d = 2 tangent orthogonal directions and their opposites. A priority queue is initially filled with the N gc points ordered along the curves. The point i at the top of the queue then tries to insert 4 new points ij in the j = 1, … , 2N d directions defined by j i and at a distance h( i ) . In order to have points inserted 'by layers,' the priority queue that is chosen is a first-in, first-out queue. Ordering the boundary points along the domain allows smooth propagation on the interior.
Each seed point i tries to spawn ij , j = 1, … , 2N d neighbor points on T 0 . Yet, there is no guarantee that point ij is not too close to another point of the queue. Points ij are hence filtered. A rectangular exclusion zone is defined by the cross-field orientation and the size field around each vertex i in such a way that any point lying in this zone will not be inserted in the queue. Finally, seed points are removed from the queue, and accepted points are added to the end of the queue as well as in P . The procedure terminates when the queue is empty. Algorithm 1 describes the procedure. Following, we will focus in detail on the two main operations, i.e., the point insertion and the point filtering.

Intersection with triangulation
Assume a point i that lies on one of the triangles of T 0 , a direction ij = j i , i.e., a unit vector tangent to the surface and the mesh size h i = h( i ) at that point. The aim is to create an edge of size h i . Therefore, the new point ij can be computed as the intersection of the triangulated surface T 0 and a circle C i with center i and radius h i . C i lies on the plane P i that is formed by the direction vector ij and the normal to the triangulation at our origin point, i (Fig. 4).
To compute ij our goal is to find the intersection point of circle C i with the triangulation T 0 (Fig. 4). We start from the triangle of the base mesh T 0i on which i lies. First, we compute the intersection line of the plane P i and the plane of the triangle P T 0 i . Then, we find the intersection points of this line with the circle C i and choose the one that lies in direction ij . Finally, the barycentric coordinates 0 , 1 , 2 of this point with respect to the current triangle are calculated. In this way, we determine whether the intersection point lies on the triangle, and therefore if we have a successful intersection with this triangle.
In the case where the current triangle is not intersected, we move forward to another triangle. Since we have already computed the barycentric coordinates with respect to the current triangle T 0i , we know where on the plane P T 0 i the intersection point lies (Fig. 5, left). The next triangle to be searched for an intersection is thus given from the computed barycentric coordinates and the adjacency information of the input mesh. This procedure continues until a valid intersection point is retrieved.
Essentially, we perform a walk in the triangulation [41] in the desired direction until we obtain the intersection point. Our experience shows that this method is efficient since it utilizes the underlying mesh as a space searching structure. For the same order of magnitude of mesh sizes on input and desired meshes, intersection points are found after a little less than two triangle visits on average.

Filtering procedure
Each point generates ij points for j = 1, … , 2N d directions. We have to ensure that new points are not too close to already generated points. Therefore, after each point ij is generated, a filtering procedure should follow. To this end, we use RTrees as a spatial search structure [42].
For every candidate point ij , we define a large enough search box (typically 2 times the mesh size). We find then the set of points Since our objective is to create right-angled triangles, i.e., equilateral triangles in the L ∞ norm, we compute the distance between the candidate point and its surrounding ones as

Surface meshing
The objective now is to create a surface mesh with the set of optimal points P = { i | i = 1, … , N g } (where N g = N gc + N gs ) that have been generated on curves and surfaces. The idea is straightforward: connect the generated points P on the initial mesh T 0 , and subsequently remove the initial mesh points (Fig. 6). A similar idea but in a different context has been used in Ref. [43], leading to a quaddominant mesh.
Robustness is of crucial interest in our method, since modifying geometric aspects of general surfaces in 3D space is a delicate task. With our approach, the main goal is to preserve the topological integrity of the mesh through each step of the process, while not compromising the accuracy of the geometric representation of the surface.
We can identify two complementary sets of generated points i : i N gc points lying on curves (feature edges) ii N gs points lying purely on surface triangles with a corresponding unique parent element (feature edge or mesh triangle, respectively) already stored for each of these points. Correspondingly, there are N r points from the initial mesh: N rc points of the feature edges and N rs points of the 'interior' surface. This division enables us to perform a bottom-up procedure where topological entities are handled independently (first mesh curves and then surfaces). We can therefore ensure that each step will have a well-defined 'predecessor' mesh to build upon.

Local mesh modifications
The basic operations utilized to locally modify the mesh follow:

Split triangles
Given a surface point and its parent triangle, split it by replacing it by three triangles. This operation is trivial to implement since it cannot change the geometry or the topology of the mesh.

Split edges
Given a point that lies on a mesh edge, split this edge. For points on boundary edges, we already know the parent edge, while for points on triangles it is easy to compute if the point lies on a triangle edge (given a user-defined threshold value ). For non-manifold boundary edges, we split all the corresponding triangles connected to this edge (Fig. 7, left). This set of triangles is readily accessible from the extended data structure for boundary edges.

Collapse edges
Remove points from the triangulation by collapsing an edge. Collapsing an edge is not always a valid operation since it can create flipped or degenerate elements. We check the fan of n triangles connected to the vertex in question and choose an edge that can be collapsed. The resulting n − 2 triangles should not intersect with neighboring ones. Again, non-manifold edges can be collapsed if all corresponding 'half-fans' pass the validity test (Fig. 7, right).

Swap edges
Given an appropriate quality criterion, swap the edge if the topology is not violated. Obviously, a feature edge cannot be swapped. Edge swaps have a significant role during collapsing of edges, since it is not always feasible to collapse all unnecessary vertices at the first pass. Edge swaps serve at this point to create improved conditions for the next  A quality improvement with edge swaps can be performed here, though it is not a necessary condition to continue to the next step.

Outline
The procedure consists of the following (Algorithm 2). All generated points are flagged to be inserted while all initial points are flagged to be removed. Starting from the model curves, each generated point splits its corresponding parent feature edge, which can be manifold or non-manifold. Each point has a unique parent feature edge and the splitting is done based on its parameter t ∈ [0, 1] , thus defined in an unambiguous way. If a point to be inserted is 'close' (w.r.t to ) to an initial mesh point, we flag the latter to be preserved instead of the former in order to avoid small-scale geometric ambiguities that may occur. The procedure follows until all points on model curves are inserted into the mesh. Following, points on surfaces are inserted by splitting triangles. These points have a unique parent triangle, and the splitting is based on the barycentric coordinates. If a point is 'close' w.r.t to on one of the triangle's edges, we split the edge, and if a point lies 'close' to an initial point, it gets disregarded (as explained for points on lines).
At this stage, we have an 'enhanced' intermediate mesh T i containing the initial (red) and generated (green) points (Fig. 3c). This mesh is of low quality, but our interest here is that it represents as accurately as possible the underlying surface and respects the topology of the initial mesh.
We want now to remove the initial mesh points. Starting again from model curves, we iterate for all feature points to be removed and examine if one of the two connected boundary edges can be collapsed. Similarly with the procedure of splitting feature edges, we can collapse non-manifold feature edges without compromising the conformity of the mesh. Subsequently, points on surfaces are removed by collapsing interior edges. Since not all points are guaranteed to be removed in the first pass, the two discrete loops for points are encapsulated in a while loop that terminates when no point to be removed remains or breaks if the remaining 'red' points cannot be removed at all. This is a crucial part since it prevents the algorithm from producing invalid topology. When an initial mesh vertex cannot be collapsed, it can simply remain in the final mesh. In practice, we observe that no more than two iterations are necessary to remove unwanted vertices for most of the models tested. While inserting points in the triangulation by split edge/triangle operators is trivial, robustly implementing the collapsing step proved be a fairly strenuous task for the general case.

Bipartite quadrilateral labeling
At this point, we have obtained a mesh T 1 that is right-angled and ready to be transformed to an all-quad mesh by triangle pairing. A necessary condition for an all-quad mesh to exist is to have an even number of edges on the boundary of each surface of the model. In this work, we trivially impose this condition on the feature curves that bound the surfaces during the first step of the pipeline (Sect. 3.1). Dividing each curve into an even number of segments is not sufficient to ensure the existence of a perfect matching, i.e., a triangle pairing that involves all triangles of T 1 . This is a typical obstacle of triangle-pairing methods, leading to isolated triangles that cannot be locally removed. A naive way to eliminate those triangles is to perform a Catmull-Clark subdivision [44], converting a quad-dominant mesh to an all-quad at the expense of vastly increasing the number of elements as well as the number of irregular vertices. In Blossom-Quad [26] a more sophisticated solution is proposed by computing a minimum-cost perfect matching of the dual graph to offer a global solution. Still, ensuring a graph to contain at least one perfect matching remains quite complicated.
We propose here a new idea for constructing one perfect matching in a triangular mesh: all edges of the matching will be used to combine their two neighboring triangles and form a quad mesh.
Let us first recall some elements of graph theory. A graph labeling is the assignment of labels to the nodes of the graph. A graph coloring is a special case of graph labeling; it is an assignment of labels traditionally called colors to the nodes of the graph with the constraint that the graph contains no monochromatic edge, i.e., no edge connecting two nodes with the same color.
A bipartite graph G is a graph that can be 2-colored. An important property of a bipartite graph is that it does not contain any odd-cycles (an odd-cycle is a cycle of odd length).
Let us now look at the quadrilateral mesh of a domain Ω as a graph G . We assume here that the boundary Ω of Ω is divided into n separated sub-boundaries Ω i , i = 1 … n , with each of the sub-boundaries Ω i forming a boundary cycle C i in G . We assume that each C i is an even-cycle. Under these conditions, it is easy to see that G is bipartite. We consider a cycle C of G : C bounds a quadrilateral mesh and it is known [32] that every quadrilateral mesh has an even number of boundary vertices. If the domain D that is enclosed by C is simply connected, then C is its only boundary and C is an even-cycle. If D is not simply connected, its boundary contains other boundary cycles C l that are by hypothesis even-cycles. Thus, ensuring that every boundary cycle C k is even is sufficient to ensure that any other cycle C is even as well. Figure 8 illustrates the aforementioned reasoning. It should be noted that not all quad meshes are bipartite: it only holds when every boundary cycle is even. Fig. 8 Every cycle C is even so the graph is bipartite and is 2-colorable. Starting from one 2-colored boundary C k , it is possible to find the color of any vertex v ∉ C k by 2-coloring any path P between C k and v (color figure online) Fig. 9 Bipartite topological correction. a Initial labeling on the even-sided boundary. b Splitting topologically inconsistent triangle. c reposition of vertex according to opposite labeled stencil. d recovery of quadrilaterals (defined by same-labeled edges) It should be noted as well that, starting from the 2-coloring of one of the boundary cycle C k , the 2-coloring of the rest of the graph is unique: the color of any vertex v ∉ C k is found by coloring any path P between any vertex of C k and v (see Fig. 8).
We consider now a triangulation of Ω and a bipartite labeling of its vertices. The graph of a triangulation is obviously not a bipartite graph since every triangle forms an odd-cycle and therefore it cannot be 2-colored.
Assume that the labeling is done in such a way that there exists no monochromatic triangle (a triangle is monochromatic if its three vertices have the same label). Then, the set of monochromatic edges of the triangulation forms a perfect matching and removing monochromatic edges leads to the desired quadrangulation. The proof is simple: every triangle of the mesh has exactly one monochromatic edge, which means that each triangle will be chosen exactly once for creating a quad: this is by definition a perfect matching. The only possible issue would be the existence of monochromatic edges on Ω . This issue is avoidable if every boundary cycle is even: coloring should be done at first on Ω and then be propagated inside. Now the last question: is it possible to avoid monochromatic triangles? The answer is no, at least not without modifying the triangulation. We start by analyzing the simplest possible case of an even triangulation without perfect matching (see Fig. 9a). In this example, all labelings leave the internal triangle monochromatic so there exists no perfect matching in this graph. In other words, merging any of the three possible pairs of triangles would leave two non-connected triangles 'hanging.' In general, picking edges in the graph while constructing a matching may split the graph into disconnected sub-graphs with an odd number of edges on their boundaries, thus not allowing a perfect matching. This condition is a special form of a well-known graph theory theorem from Tutte [45]. By inserting an additional Steiner point of the opposite label (by splitting one edge of the triangle in question) our condition is met (Fig. 9b). A perfect matching actually exists and is composed of all monochromatic edges as demonstrated above (see Fig. 9d).

Random labeling
We start by 2-coloring the boundary cycles C k of Ω and assign a random label to all internal vertices. Figure 10a shows this initial random labeling. Random labeling is not the worst-case scenario: it typically produces 25% of monochromatic triangles. Yet, this number is too big because the number of Steiner points that should be added to remove all the monochromatic triangles is sufficiently large to damage the mesh size field.
It is possible to dramatically decrease the number of monochromatic triangles by applying the following smoothing algorithm. Each internal vertex is re-labeled if changing its label reduces the number of (its adjacent) monochromatic triangles. The smoothing usually ends with monochromatic triangles that are either isolated or form adjacent pairs (see Fig. 10b). In this specific example, no pairs are observed but they cannot be avoided in the general case.
Then, bipartite topological correction is performed. When monochromatic triangles appear in pairs, their common edge is split. When they are isolated, we choose to split their longest edge, to not degrade the accordance with the size field. This process continues until all monochromatic triangles are eliminated. Figure 10c shows the final labeling without any monochromatic triangle.

3
The triangular mesh of Fig. 10c can then be transformed into a quadrilateral mesh by removing its monochromatic edges as depicted in Fig. 10d. Mesh of Fig. 10d is not of high quality. We had to use all Gmsh's heavy optimization weaponry to obtain the final mesh of Fig. 10e.

Cross-field guided labeling
The triangular mesh of Fig. 10a is not suited to be transformed into a quadrilateral mesh; obtaining a good mesh at the end heavily relies on advanced optimization. In Sect. 3.2, we proposed an advancing-front scheme where each point tries to add other points guided by the orthogonal directions provided by a cross-field. The cross-field guided point insertion process produces an excellent labeling scheme: new vertices added by a vertex v have the opposite label of the one of v (since the edges formed between these vertices will be the cross-field aligned ones that we want to have in the final quad mesh). Figure 11a shows the initial labeling after frontal point insertion. As few as 14 monochromatic triangles are present in the triangular mesh, leading to a very small number of Steiner points. More importantly, no Steiner point has been inserted at the vicinity of the boundary, leading to several perfect layers of quads (see Fig. 11d).
By inserting the Steiner points, we have constructed the underlying topological connectivity of an all-quad mesh (Fig. 11b). Mesh vertices can be repositioned at the center of the stencil of opposite labeled neighbors (essentially a single iteration Laplacian smoothing, e.g., Fig. 11c). In the case of non-planar surfaces, the repositioned vertices are projected onto the original triangulation. We can then trivially extract an all-quad mesh, since each pair of triangles whose common edge has two same labels defines a quadrilateral (Fig. 11d). At this point, the advancing labeling scheme is converted to the true 2-coloring of the final bipartite quadrilateral mesh. Another way to see this is that the set of same-colored edges is perpendicular to the set of the edges of the perfect matching of the dual. It is interesting to note that forming the quads this way before inserting the additional vertices would lead to a quad-dominant mesh with non-quad polygonal regions that are even-sided due to the bipartite condition.
Finally, it must be noted that our algorithm is straightforward, and linear in time compared to the computation of an optimal perfect matching, which is quadratic in time and very complex. For the sake of completeness, the algorithm is presented on Algorithm 3.

Mesh optimization
The algorithm described is always able by construction to produce a topological quad mesh T q . Moreover, T q is also of high geometric quality (measured by the scaled Jacobian Q [46]) since it follows the cross-field directions.
In few cases, especially when the requested size field is much coarser than the geometric characteristics of the model, we may encounter minor defects where elements with Q ≤ 0 may occur, such as doublet quads (two quads connected only by one vertex) or flat quads on the boundary (since boundary vertices cannot be repositioned). Both cases can be locally handled, the former by merging the doublet into one quad and the latter by inserting additional vertices and thus 'pushing' the boundary irregular vertex inside the domain.
The only smoothing procedure used in this work is the Laplacian smoothing described in Sect. 5.2 for the repositioning of vertices. Further and more suited smoothing (e.g., Winslow smoothing [47]) can be performed on the final quadrilateral mesh. The next step in our developments will be to complement our approach with the more sophisticated optimization procedure proposed in [38] that takes as input the quad mesh T q , the cross-field as well as the size field h. The main features of this optimization process are: • All vertices of high valence are eliminated. • All boundary vertices have their optimal valence in such a way that we can create a boundary layer mesh without effort.
• Isolated irregular vertices corresponding to the singularities of the cross-field are preserved. • Advanced vertex relocation schemes are performed to obtain a geometry regular quad mesh.

Results
To demonstrate the capabilities of our methodology, we have tested the whole pipeline on a variety of cases with diverse characteristics (Fig. 12). One of our main interests is to be able to produce meshes on any input, regardless of its complexity (Fig. 13). Reviewing the relevant literature, we observe that the majority of the works present results on relatively simple models. We work on large model databases (that are recently becoming the common standard for validation purposes on bulk numbers of models), and we aim to have a close to 100% success rate on them. The models are obtained from the following sources: the ABC [48], MAMBO [49] and Thingi10k [50] datasets, and the supplementary test data of LoopyCuts software [51]. We present statistics for all the models from MAMBO and LoopyCuts, since all their models are directly suited for meshing (Table 1). MAMBO contains three categories of CAD models: basic, simple, and medium, of which we use the latter two since they have more complex features of interest. LoopyCuts models are in general simpler models in terms of surface complexity (smooth surfaces, absence of very small features) oriented mostly for computer graphics research.
A desired property of our mesh generation process is to keep user intervention at a minimum. Input parameters are a threshold value to suppress potential discrepancies during local mesh modifications and the angle threshold to detect hard edges. The angle threshold may be disregarded, and the algorithm can initiate from a random edge of the initial mesh instead of the model curves. We keep in mind that this may result in low-quality meshing on hard edges.
Cross-fields are precomputed and considered input to the algorithm, or they can be taken as the unit axis directions for practical purposes. Size fields can be a uniform input size or computed from the curvature and geometric characteristics of the input mesh. Both cross and size fields impact considerably the output mesh quality, and we can observe that the directionality of the elements matches the input cross-field directions. It must be noted that our point generation module can have any kind of metric input. For example, we can use an anisotropic metric or a prescribed user-defined size field that is 'incompatible' with the computed cross-field. In Fig. 14 (right), a size field of sinusoidal form is prescribed and we can observe that our algorithm can successfully handle it at the cost of lower element quality and an excess of irregular vertices to accommodate for the size transitions.
Our pipeline can succesfully produce a quad mesh in all the models that can be processed by GMSH for the initial Fig. 12 All-quad meshes on various models obtained from the following sources: ABC [48], MAMBO [49], Thingi10k [50] and LoopyCuts [51] triangulation. One of the most challenging aspects of our methodology is the removal of the initial mesh points by collapsing and testing multiple 'extreme' cases on that is an ongoing process. As desired, non-manifold features are well handled and conformity is preserved (Fig. 14, left). The final conversion to an all-quad mesh proves to be straightforward and efficient, given that the input respects the fundamental topological condition (i.e., even number of edges on each curve). The quad meshes exhibit high quality Q, computed as the scaled Jacobian (Table 1). Minimum quality over all LoopyCuts models is 0.22 with the vast majority of the models having over 0.50. On the other hand, we observe similar outputs on the MAMBO models except of a few cases where very thin sliver regions result to almost flat quad elements with Q ≃ 0 . These kind of features are common in 'real-life' models and element quality there could be improved with further optimization. The output meshes shown here consist of the 'raw' quadrilateral meshes, without heavy further optimization (an example of further optimizing a quad mesh is presented in Fig. 15).
One of our main interests is the computational efficiency of the method. The majority of operations are local and the algorithms of linear complexity. The most time-consuming parts of our pipeline is the (inevitable) use of spatial search structures (involved in filtering and also in projection of points onto the triangulation during surface meshing and all-quad conversion). The performance of the algorithm is heavily dependent on the input mesh, specifically its number of triangles and colored surfaces. More triangles lead to more spatial searches during point generation and filtering as well as more local mesh modifications during surface meshing. Spatial searches and projections are the most  important bottlenecks on the performance of our pipeline and we consider that we can speed up with more efficient procedures. In Table 1, we present the average time for each part of our pipeline, using a single-thread on a laptop with 2.6 GHz CPU.

Conclusions
We have presented a method for the creation of quadrilateral surface meshes on general complex geometries. One of the important aspects of our pipeline is its modular nature. Each step can be independently used, for example, a more efficient point sampling method could be used, or a typical CAD meshing algorithm to generate the mesh with these points. The surface meshing algorithm proposed is used to produce right-angled triangles, but it can be repurposed as a remeshing tool for triangulations, given appropriate input (i.e., an 'asterisk' field [32,39]). The proposed two-step method for surface meshing aims at balancing the trade-off between output quality and robustness. One of its strengths is the ability to handle the model as a whole and thus create cross-aligned elements globally. While producing highquality elements in this aspect, we always respect the topology and succeed at meshing small-scale and non-manifold features as shown in the results demonstrated (Sect. 6).
As stated, scaled cross-fields are beneficial for quad mesh generation. We have mainly used here size fields obtained from the scaling of cross-fields. In practice though, a user would commonly prefer prescribing a size field that is a better fit for the problem in hand, and unfortunately this is currently an open problem. One of the advantages of our approach is that our pipeline can produce meshes regardless of the topological integrity of the input metric, for example with a user-defined size field that is incompatible with the computed cross-field (Fig. 14, left). Nevertheless, adapting a direction field to an input size field would vastly improve the results. Regarding highly anisotropic metrics (for example for the generation of structured boundary layer meshes), it is straightforward to extend the point generation part of our pipeline for such a use. Nevertheless, we consider the triangulation on such a metric a much more difficult task (where the aim is having structured layers of anisotropic triangles to be combined). We believe that the best strategy would be to create structured layers of boundary quadrilaterals and then subdivide them.
A novel method to obtain an all-quad mesh from any quad-dominant mesh with minimal, if any, post-processing clean-up is introduced. By using a bipartite labeling scheme, we simplify the global treatment required for quad meshing to a set of localized operations that can be performed on a quad-dominant mesh. We are thus enabled to produce always a high-quality quad mesh, regardless of the complexity of the model. Our meshes can be further optimized with a strategy to remove the majority of irregular vertices, presented on [38].
We are currently in the process of incorporating this quad surface mesher as a component of a general high-quality hex-dominant meshing pipeline, where it will be used as the volume bounding surface. An interesting open problem to this end is if hexahedral element generation should be constrained by a quadrilateral mesh, and if yes, what would be its optimal characteristics. Besides that, we are investigating ways to take advantage of the topological labels used here for quadrilateral elements to the generation of hexahedra, possibly with a labeling scheme that can codify more topological information. Fig. 15 Left: unstructured quad mesh produced with our pipeline, where we can observe the excess of irregular vertices. Right: quasi-structured quad mesh obtained with the topological optimization of Ref. [38]