3D B-Rep meshing for real-time data-based geometric parametric analysis

This paper presents an effective framework to automatically construct 3D quadrilateral meshes of complicated geometry and arbitrary topology adapted for parametric studies. The input is a triangulation of the solid 3D model’s boundary provided from B-Rep CAD models or scanned geometry. The triangulated mesh is decomposed into a set of cuboids in two steps: pants decomposition and cuboid decomposition. This workflow includes an integration of a geometry-feature-aware pants-to-cuboids decomposition algorithm. This set of cuboids perfectly replicates the input surface topology. Using aligned global parameterization, patches are re-positioned on the surface in a way to achieve low overall distortion, and alignment to principal curvature directions and sharp features. Based on the cuboid decomposition and global parameterization, a 3D quadrilateral mesh is extracted. For different parametric instances with the same topology but different geometries, the MEG-IsoQuad method allows to have the same representation: isotopological meshes holding the same connectivity where each point on a mesh has an analogous one into all other meshes. Faithful 3D numerical charts of parametric geometries are then built using standard data-based techniques. Geometries are then evaluated in real-time. The efficiency and the robustness of the proposed approach are illustrated through a few parametric examples.


Introduction
Since few years, data-based models and reduced order modeling [10,11,47] have been particularly studied in many physics applications. Such approaches allows to determine complex solutions in real-time. By real-time we mean a computational cost drastically reduced in comparison with a classical approach. In most reduced order approaches, an offline step is first needed in order to build the database model. For parametric studies, this means constructing and running the model for given sets of parameters wisely chosen. In the online phase, executed in near real-time, this database is used to construct the solution of a new set of parameters. Common reduced order modeling techniques classically used in computational mechanics often require structured meshes to be efficient by avoiding inaccurate projection steps between all meshes. This is obviously encountered when geometric parameters are considered. One way to circumvent this issue is to build isotopological meshes of all the geometric instances. Previous work has been done to build the foundations of our approach [2,3]. Moreover, research on the same problematic has been recently published [40][41][42]. For different geometric instances with the same topological properties (e.g. the same shape), being able to construct isotopological meshes of all these instances makes it very easy to compare these meshes. Consequently databased approaches with geometric parameters are mainly guided by topological information throughout the offline step. Generating isotopological quad meshes therefore seems as the right way to solve this problem.
In numerous applications in domains like computational mechanics, aerodynamics, crash simulations and fluid-structure interactions, 3D quadrilateral meshes are recommended due to their highly regular structure. Generation of high quality quadrilateral meshes from B-Rep triangulated meshes has been recently extensively studied. All developed methods seek to construct such meshes with local and global specified properties. Most common desired properties are element quality, alignment with principal curvature directions or boundaries, global topology structure and respect of sharp features. Most recent techniques for quadrangulation use a global parameterization of the input triangulated surface [6,27,45,48]. Global parameterization is of particular interest for parameterizing the whole surface and handling discontinuities along seams. Constructing such global parameterization is usually done using a reliable field that serves as a guidance to drive the parameterization gradient. Cross fields can be used for that purpose, and in addition allow to capture principal curvature directions, boundaries and sharp features while understanding the topology. Different techniques exist for 2D or 3D quadrilateral meshes purposes [5,15,30,49]. Quadrilateral layouting (i.e., partitioning the triangulated mesh into generalized squares) is particularly efficient to build structured quadrilateral meshes [8,9].
In the present work, we seek to decompose objects defined by their triangulated boundary into geometry-aware quadrilateral layouts understanding geometric features as well as possible. This tedious task is first done using pants decomposition which was widely studied [12,21,34,53]. Secondly, we perform cuboid decomposition [32,35,36] on each pants patch that allow us to obtain a decomposition in quadrilateral patches. We present a novel decomposition method that splits automatically the surface into geometry-feature-aware patches. We use an optimization process of quadrilateral patches embedded on surfaces. An optimal global parameterization gives the final surface parameterization. Topological properties of the designed cross field is entirely provided from the layout. Our main contribution is a new way to design such quadrilateral layouts to ensure an optimized resulting quadrilateral patch decomposition of the surface suitable for isotopological meshing purposes.
With this generic setup, our method can be applied in automatic scenarios involving specific cuboid sketches per pant patch. These kind of automatic processes are really useful to compute the offline phase needed for parametric studies using reduced order modeling and data-based models.
The paper is organized as follows. We first extend the decomposition pipeline thought to understand the geometry in "Smart model decomposition" section. Then using global parameterization aligned with a cross field capturing geometric information, we extract feature-aligned quadrilateral meshes in "Direction field generation" and " Aligned global parameterization computation from geometric direction field" sections. Handling isotopological datas throughout the 3D quadrilateral meshing process, geometric snapshots are used to build a reduced order model in a full autonomous integrated algorithm as presented in the various examples of "Application of the MEG-IsoQuad method togeometric parametric analysis" section.

Problem statement
In this article, we want to compute 3D quadrilateral meshes whose properties belong to a very restricted class. The goal is to be able to compare different shapes, i.e., different geometries of the same topology class. This kind of meshes are required in most parametric applications such as statistical shape analysis and reduced order modeling [19,24]. Figure  1 shows a graphical illustration of the problem. The most desirable properties of the targeted meshes are: • Minimize the number of singularities, defined to be the nodes with a valency different from 4. • Alignment with features, curvature directions and boundaries.
• Determine high quality elements, as close to a square as possible.
• Constrain the number of elements, connectivity and features locations to enable relevant comparison of all studied meshes.

Proposed approach
We give a workflow to compute high quality 3D quadrilateral meshes with specific properties needed in parametric analysis. Our goal is to generate such meshes for objects with complex geometry and arbitrary topology provided from boundary representation or scanned techniques. We want our framework to be as robust as possible.
We present an integrated pipeline partitioning the input triangulated surface to relevant domains until having comparable meshes for parametric studies. The first contribution consists in understanding the geometry while decomposing the mesh. The second contribution consists in computing an optimized aligned global parameterization using a very coarse quadrilateral layout helping us to define the required quadrilateral mesh.
Reduced order modeling is a tedious task that requires isotopological meshes in most cases. All geometric instances of the population must satisfy such mesh properties to be compared and studied. The issue of having isotopological meshes is addressed by Fig. 2 The MEG-IsoQuad method. Building isotopological analogous feature-aligned quadrilateral meshes from triangulated meshes with the same topology. Resulting isotopological meshes have non-uniform isotropy generating the same parameterization for all members of the given population. Figure 2 provides a simple understanding of our approach: the MEG-IsoQuad method.

Building isotopological analogous quadrilateral meshes: the MEG-IsoQuad method
In this section we present a method thought to build 3D quadrilateral instances from triangulated meshes suitable for data-based models with geometric parameters. This method allows us to identify the main shape variations. Quadrilateral mesh topology needs to be the same among all snapshots. We also strive to place vertices in a geometry-aware manner, i.e., all isotopological meshes have analogous points in the whole mesh. These constraints can be satisfied using our framework based on topology and geometry aspects.

Topology and parameterization prerequisites
Following prerequisites introduce briefly basic notions of geometry and topology that are needed to present our work. The interested reader unfamiliar with those concepts can find relevant detailed information in the cited references and in the corresponding literature.
Topology context and related work Topology is the study of properties like continuity, connectedness and boundaries of a space that are preserved under continuous deformations, such as bending and stretching, but not tearing and gluing. A homeomorphism is an isomorphism that admits a continuous function between two topological equivalent spaces that has a continuous inverse function. We are interested on transformations that preserve all the topological properties of a given space. Homeomorphic spaces admit a homeomorphism between them, thus topological spaces are equivalent. See Hatcher [22] for more details.
A surface M is a 2-manifold, i.e., a topological space in which each point has a neighborhood homeomorphic to either the plane R 2 or the closed half plane R 2 + . Points with closed half-plane neighborhood are defined as the boundary ∂M of the surface M. In the following we investigate only connected, orientable, compact surfaces with boundaries. Genus of a connected and orientable surface M (i.e., a 2-manifold embedded in R 3 ) is the maximum number of non-intersecting closed curves which can be drawn on it without disconnecting the surface. With these definitions, taking into account a genus-g surface M possibly with b boundary components, we can now define the topological invariant. This invariant used to classify surfaces is called Euler characteristic: If M is a triangulated surface with vertices V , edges E and faces F , we can define the characteristic as χ(M) = dim(V ) − dim(E) + dim(F ). An example of Euler characteristic is provided in Fig. 3a. We note that surfaces with different Euler characteristic cannot be homeomorphic.
We define a homology basis for M to be any set of 2g cycles whose homology classes generate the first homology group H 1 [17]. Cutting along these cycles yields a genus-0 surface. The set of all handle and tunnel loops form a homology basis. Suppose a closed surface M ⊂ R 3 separates R 3 into a bounded space I and an unbounded space O. Handle and tunnel loops on M can be defined as follows. A loop a i is a tunnel if it spans a disk in the unbounded space O. A loop b i is a handle if it spans a disk in the bounded space I.
Pants decomposition was studied in Hatcher et al. [23], subsequent work has been done to find the optimal segmentation of a given surface into relevant pants patches [12] using the shortest homology basis [26]. Geometry-aware pants decomposition has been also investigated by Zhang and Li [53]. Li et al. [34] developed a pants decomposition framework for computing maps between surfaces with arbitrary topologies. An example of pants decomposition using a homology basis can be seen in Fig. 3b. More recently, Hajij et al. [21] tried to segment surfaces into pants using a morse function.
Let M g,b be a surface of genus g with b boundary components. A pants decomposition of M g,b is a collection of pairwise disjoint simple cycles that splits the surface into pants patches. Each pants patch is a genus-0 surface (topological sphere) with 3 boundaries. We assume that M is a surface with negative Euler characteristic, i.e., M is none of the surfaces M 0,0 (topological sphere), M 0,1 (topological disk), M 0,2 (topological cylinder) and M 1,0 (topological torus). In this case pants decompositions of M do exist, and each pants decomposition consists of 3g + b − 3 curves and divides M into 2g + b − 2 pants patches.
A direction field defined on a surface M is a tangent unit vector field: at each point of the surface, there exists a direction u such that u = 1 and u.n = 0, where n is the normal of M. A N-symmetry direction field U is a multivalued direction field: at each point of the surface M, there exists a N-symmetry direction u which is a set of N directions {u 1 , u 2 , ..., u N −1 , u N } preserved by rotations of 2π/N around the normal n of M. These fields are subject to the well known Poincaré-Hopf theorem for a 2-dimensional manifold M. If M has boundaries, the vector field must be pointing in the outward normal direction along them (for higher symmetries, see the boundary number theorem [49]). b Pants decomposition provides two pants (gray and blue meshes) by cutting along handle loops and using symmetry. Homology basis is composed by two handle loops and two tunnel loops respectively depicted in red and green

Used mapping techniques
Different approaches of mapping exist. In our case, we are interested in mappings or parameterizations which map a surface M embedded in R 3 to a canonical domain D in R 2 . Throughout the paper, triangulated disk-like surfaces M with vertices V , edges E and faces F are considered for parameterization. For instance, techniques to map a multiply connected surface already exist [52]. We use discrete harmonic mapping to solve parameterizations on disks. Harmonic mappings have attributes derived from conformal parameterization, but there is no guarantee on angles. To proceed we construct a harmonic function f : M −→ R such that f = 0. Harmonic maps minimize Dirichlet energy: The surface boundary ∂M is first mapped to the boundary of the parametric domain and then the parameterization for the interior vertices is obtained by solving the linear system: where v i , v j ∈ V , N i is the neighborhood of v i , and w ij is the scalar weight assigned to the oriented edge e ij (v i , v j ). Different parameterization methods assign different weights w ij . The first definition of weight was introduced by Tutte [51]. In the parameter space, each vertex is placed at the barycenter of its neighbors. Recently, Saboret et al. [50] have implemented a CGAL package handling some of the state-of-the-art surface parameterization methods such as least squares conformal maps, discrete conformal map, discrete authalic parameterization, Floater mean value coordinates or Tutte barycentric mapping.
In this context, we use mean value coordinates introduced by Floater [18]. According to the Rado-Kneser-Choquet theorem, if weights are positive and the boundary ∂M is mapped homeomorphically in a convex square parametric space, the mapping driven by the harmonic function f has to be bijective. This allows us to determine bijective mappings in order to avoid fold-overs in the smart model decomposition part, see e.g. "Smart model decomposition" section.

Smart model decomposition
In order to build valid quadrilateral layout to determine the surface parameterization, we proceed in two steps: geometry-aware pants decomposition ("Geometry-aware pants decomposition" section) and feature-aware cuboid decomposition ("Feature-aware cuboid decomposition" section). Pants decomposition and cuboid decomposition are constructed to respect the topology of the input triangulated surface and consider features. Cuboid decomposition approximates very roughly the geometry while faithfully replicating its topology taking into account sharp edges and vertices. Due to its patch regular structure, it can serve as parametric domain needed for pure 3D quadrilateral mesh computation.

Geometry-aware pants decomposition
Closely related literature approaches and algorithmic prerequisites We focus on the decomposition of a triangulated surface into a set of pants patches {T i }. Pants decomposition provides a canonical decomposition scheme for common surfaces. χ = −1 for a pants patch, therefore a pant is the simplest topology after the sphere, disk, cylinder and torus. Many algorithms take as input handle or tunnel loops for genus-g surfaces to segment into pants patches [34]. The handle loops required for automatic pants decomposition are computed using the technique presented by K. Dey et al. [25]. An improvement was also developed by the same group [16]. In the following, we reuse the pants decomposition algorithm of Al-Akhras et al. [4] in order to integrate a most recent method.
New extensions According to previous work presented by Al-Akhras [3], we introduce an extension to enumerate the entire space of topological pants decomposition possibilities at each step. In other words, instead of picking 2 boundaries in an arbitrary manner, we suggest to take all couples of 2 boundaries among all non-repeating and commutative combinations when determining a pants patch relative to the considered step. Besides, giving some feature points, our new algorithm is able to determine a pant by slicing the mesh along a cycle passing across these locations.
User inputs User manual inputs are reduced to feature points selection in order to guide the decomposition passing through these points of interest. For other used geometric criteria, they are suggested to the user.
Feature-aware cuboid decomposition algorithm Given the homology basis formed by handle and tunnel loops, we can take a subset H composed of g simple pairwise disjoint handle loops {h 1 , h 2 , ..., h g−1 , h g }. Slicing surface M with b boundary components along its g handle loops will lead to a genus-0 surface with 2g + b boundary components denoted as W = {w 1 , w 2 , ..., w 2g+b−1 , w 2g+b }. We iteratively pick two boundaries w i and w j from W (see algorithm 1) taking into account combinations and compute a new simple cycle w ij to bound them, i.e., w ij is homotopic to w i • w j . The three cycles w i , w j and w ij bound a pants patch T k . We remove this pants patch T k from M. The remaining patch is still genus-0 but its boundary number reduces by 1: the two cycles w i and w j are removed, and one new cycle w ij is inserted. This is iteratively performed until |W | = 3. This idea is formulated in Algorithm 1, and the operation that traces a cycle w ij homotopic to cycle w i • w j is formulated in Algorithms 2 and 3.
Choice of the geometric criterion Shortest length of cycles on weighted triangulated meshes are computed using a Dijkstra's algorithm based on weighted tree graphs. A CGAL package exists to compute these paths [29]. Different geometric criteria can be used to guide the pants decomposition. Geometric criterion can be adapted to the mesh M. In algorithm 2 and also in Algorithm 1 to sort loops in L, this can be changed for areas Compute

07:
For all couples [w i , w j ] in N c : 08: Compute a cycle w ij homotopic to w i • w j . 09: Add loop to L. 10: End For 11: Sort relevant loops in L = {l 1 , ..., l dim(N c ) } using a global geometric criterion. 12: The optimal w ij cycle is classified in L. 13: Remove w i and w j from W , and add w ij into W . 15: k ← k + 1.  Genus-2 plate pants decomposition. Triangulated surface (a) with its handle and tunnel loops in red and green, respectively (b). Pants decomposition using loops with shortest distance (c) and with symmetry (d)

Fig. 5
Geometry-aware pants decomposition with feature points. Double T (a) is decomposed into pants using 4 feature points given by the user or determined automatically. A 2-torus pants decomposition passing through feature points (b). Notice that using shortest length loops, the result will be the same for (b) of minimum curvature [31], minimum length or symmetry. For instance in Fig. 4c pants decomposition is performed using loops with minimum length, whereas in Fig. 4d decomposition is made by symmetry. In Fig. 5, user selected feature points are given to guide the pants decomposition. All of these guidings yield a geometry-aware decomposition.
Discussion on robustness A common surface admits infinitely many pants decompositions. In general, not all pants decomposition results are suitable for the next step of our algorithm. By correctly defining feature points (in general a very small set of points) to lead the segmentation, the algorithm succeeds. If the pants decomposition is guided by the different geometric criteria presented above, we obtain satisfying results for all our test cases.

Feature-aware cuboid decomposition
Closely related literature approaches Pants patches provide a very simple topology and hence each pants patch can be treated separately. Li et al. [32] presented a method that generates cuboids per pants patch with a given user data input. Al-Akhras [3] improved it in order to perform the decomposition in a way that it does not require any input from the user. Nevertheless, developed techniques do not consider sharp features and deal with a simple cuboid decomposition scheme. They both try to generate corners and polyedges on each pants patch T i and decompose it into a set of cuboids, each having 8 corners and 12 polyedges in order to construct a volumetric parameterization. In addition, these two previous approaches do not handle complicated cuboid configurations and generated arcs and nodes are not necessarily positioned onto features of the input surface mesh. Finally, prior automatic processes restrict the placement of the quadrilateral layout nodes in specific areas, e.g. using harmonic functions.
New extensions Our work presents an extension to handle more complicated geometry especially with sharp features. First, on each pants patch the number of quadrilateral patches n Q is set to respect sharp vertices and edges: we aim to perform a feature-aware cuboid decomposition. Indeed, optimal number of quadrilateral patches depends on the features embedded into the considered pants patch. Secondly, irregular vertices of the layout with a valence v = 4 are placed onto relevant areas with high concentrated Gaussian curvature and/or with significant features.
Algorithmic prerequisites The space of valid quadrilateral layouts Q is defined in equation (3). A quadrilateral layout is composed by arcs and n n nodes of valence v i . Cuboid configurations C decompose a surface into a set of quadrilateral patches which are ready for hexahedral meshing due their cubic structure. In this paper, the cuboid decomposition is dedicated only to surface decomposition. Indeed, this work only focuses on structured quadrilateral meshes.
User inputs To reach the feature-aware cuboid decomposition goal, manual user inputs are necessary. Depending on the pant geometry and features to replicate in the cuboid decomposition, 2 inputs are mandatory. Using a cuboid configuration templates library, a segmentation layout is choosen by the user; thus involving the wish to suit the desired features of the pant into the decomposition. Then, the structure of the requested template requires the manual selection of its asscociated feature points. Nature of these features is automatically suggested and the user has just to specify their location.
Feature-aware cuboid decomposition algorithm The 3 boundaries of a given pants patch will be arbitrarily denoted by B 1 , B 2 and B 3 . We process these pants patches one by one in an arbitrary order. To guarantee cuboid corner alignment, when we determine one pants patch's result, we transfer its corners on the boundaries of the adjacent pants patches if they are not processed yet. Notice that just a general algorithm with the most relevant specifications is provided using a simple decomposition scheme of 4 cuboids per pant patch.
Step 1. We generate 3 feature cutting curves W S1 , W S2 and W S3 : Then a temporary non-feature cutting curve W Ti is defined as the isoparametric curve of the function f i for the value f i . Figure 6b.
[C] We remove a long branch, i.e., by cutting along its temporary cutting curve W Tk .
After filling the cutting hole, the resulting patch is a topological cylinder with 2 boundaries B i and B j . We denote this temporary patch by P Tk . Figure 7a-c.
[D] Using previously computed temporary patches P Tk , we first set u = 0 for vertices on B i and u = 1 for vertices on B j , then solve u = 0. Using an approach close to Zeng et al. [52], among all iso-v curves along ∇u starting from B i and B j , we automatically select two relevant curves who are intersecting the filled boundary B k , one starting Notice that triangles from filling are removed during U Tk computation. Figure 8a.
[F] Points in C are mapped into each unit-square map U Tk and automatically classified using parametric coordinates u and v. Then we perform partial line inverse mapping between each classified cutting points c i related to relevant U Tk and merge computed lines to obtain valid feature cutting curves W S1 , W S2 and W S3 . Figure 8b, c.
Step 2. We generate 3 feature-aware oriented unit-square maps denoted U k : [A] We remove a long branch by cutting along its feature cutting curve W Sk . We proceed like Step 1 [C] to obtain a feature-aware cutting patch P k .
[B] Using the same workflow of Step 1 [D] feature seed points s ik and s jk are determined.
[C] Using the same workflow of Step 1 [E], we obtain 6 feature seed points s ij , s ik , s ji , s jk , s ki and s kj . Then 3 feature-aware oriented unit-square maps U k are extracted.
Step 3. We generate all arcs of Q: [A] According to Step 1 [F] feature cutting arcs W Sk have already been computed. We keep W S3 passing through all cutting feature points.  With feature points in C (see Step 3 [C]) or other feature points in D or O, a huge variety of shapes can be decomposed into cuboids and hence handling more complex geometry cases is possible. To go further, the pants-to-cuboids decomposition can be divided in two main parts: the biological and mechanical parts, see e.g. Fig. 10. For instance, specific biological quadrilateral layouts have been studied by Al-Akhras et al. [2]. Details for the mechanical workflow are explained in Fig. 11. Given a consistent pants decomposition, our algorithm is able to decompose a triangulated surface into quadrilateral patches suitable for surface parameterization. This previous task is done only by specifying relevant feature points relative to a given quadrilateral layout template.

Direction field generation
Introduced in "Topology and parameterization prerequisites" section, 4-symmetry directions fields (i.e., cross fields) are widely used to determine an aligned global parameterization [6][7][8]45]. In other words, these fields are useful to guide the parameterization to determine quality quadrilateral surfaces. Designing a smooth cross field C is done with a given set of constraints. We categorize these constraints into two groups: topological and geometric constraints. Topological constraints are imposed singularities and numbers induced by the surface topology χ(M). Geometric constraints are intrinsically embedded on the surface. Thus we seek a cross field C which is smooth, aligned with the local geometry and topologically compatible. The quadrilateral layout Q or cuboid configuration C contains topological properties whereas the surface M hold the geometric information. We consider a triangulation of the surface M, assumed to be a 2-dimensional manifold of genus-g with b boundary components. Directions will be stored at faces F to avoid definition of supplementary tangent planes. The first step of the discretization consists in finding a reference direction in each triangle. This is done by choosing a local orthonormal frame (x, y) attached on each face f . The vector x is an unitary vector along one of the oriented edges of face f , and y = n × x, where n is the normal of f . A direction u on f can be formulated in terms of polar coordinates. Due to the unit norm of such directions, it is completely parameterized by the polar angle α it forms with x. We denote α the direction angle. By unfolding adjacent triangles isometrically to a plane along their common edge, angles can be expressed in a common coordinate frame. To specify the number of N th turns the direction u A undergoes to match with u B when passing from A to B, the period jump [33] is adopted. Other works introduce an angle ω named connection angle to solve this ambiguity [14,15]. In the following, the connection angle based discretization is used for topological design and period jump based discretization is used for geometric design.

Topological cross field generation
The first discretization computes an angle c α j i from face f i expressed in adjacent face f j which is in general equal to α j j , see e.g. equation (4). The angle κ ij ∈ [−π, π] represents We seek to fix all topological degrees of freedom to restrict the cross field topologically compatible to our quadrilateral layout Q. Technically speaking, we target a smooth cross field C that is singular only at specified vertices: the position of irregular nodes of Q. The approach of Campen and Kobbelt [8] is adopted to determine the cross field topological degrees of freedom from the input quadrilateral layout. Figure 12b provides a topological cross field calculated using the minimization of energy introduced by Crane et al. [15]. The generated field respect the Gauss-Bonnet theorem by distributing the Gaussian curvature in a smart manner [14].

Geometric cross field generation
The second discretization computes an angle p α j i from face f i expressed in adjacent face f j which is in general different to α j j , see e.g. equation (5). Angles α i i and α j j are respectively expressed in their native faces. The integer p ij ∈ Z is the period jump when passing from face f i to face f j . The integer N is the symmetry of the field, i.e., 4 for cross fields.
A smooth cross field that interpolates relevant principal curvature directions, sharp features and boundaries restricted to a given topology is targeted. We search the smoothest field taking into account constrained directions given by geometric features. The period Fig. 12 Cross field generation. a Input triangulated pant with sharp features. b Topological cross field that doesn't respect the sharp features directions. c Geometric cross field with boundaries and sharp features interpolation where the field total curvature is null jumps of the field are fixed. A technique introduced by Bommes et al. [6] and later used by Al-Akhras et al. [2] is adopted to compute such field. Figure 12c shows a geometric cross field that interpolates sharp features and boundaries.

Aligned global parameterization computation from geometric direction field
We now compute an aligned global parameterization, i.e., a map from the mesh M to a disk-like surface parameter domain ∈ R 2 . We assign a couple (u, v) of parameter values on each vertex of M. The parameterization should be locally aligned with the features of the mesh, this is done using the guiding geometric cross field previously computed. Such parameterization implies that the gradients ∇u and ∇v of the discrete scalar field must follow the cross field directions on each face.

Seamless parameterization
A planar parameterization of a mesh M embedded in R 3 into a parametric domain embedded in R 2 is in general done by computing a cut graph. A cut graph G is a connected graph formed by edges of M that splits the mesh into a disk-like surface mesh M c . Seams are then defined as duplicated paths of G. Transitions across seams need to belong to a very restricted class. We search for rigid transformations with a rotation angle multiple of π 2 . Moreover across each seam edge or vertex, the corresponding transition must be integral, i.e., relative to an integer. Thus we talk about integral seamless parameterization [27]. We target the cross field first and second directions u C and v C for the gradients of the parametric coordinates ∇u and ∇v. The parameterization is then computed as the solution of a constrained minimization problem: The integer t ij ∈ Z 2 is a translation across cut edge e ij . Aligned global parameterization can also integrate sharp edges and boundary constraints to ensure that a common relevant edge is mapped to an isoparametric curve.

Node connection constraints
Once a valid quadrilateral layout and consistent geometric cross field are provided, we wish to restrict each arc in a way that two incident nodes lie on a common isoparametric curve. With isoparametric we mean that either the u or v parameter is constant along the curve when taking transitions into account. The minimization problem (6) is then constrained using node connection constraints. These constraints are derived from the quadrilateral layout Q. Typically, each arc must lie on a common isoparametric curve taking seams transitions into account [46]. We are now able to compute a suitable global parameterization of the mesh M, given the associated quadrilateral layout Q.

Quadrilateral layout embedding optimization
Depending on geometry and position of Q's nodes, a better parameterization can be found using node relocation. Due to the nature of the quadrilateral layout, fixed nodes rise to large distortions or even local non-injectivities. Mathematically speaking we are going to optimize equation (6). The nodes are re-positioned based on the gradient of the parameterization's objective functional with respect to their positions. This is done iteratively until a global embedding quality is reached. We follow the method developed by Campen and Kobbelt [8] to perform such optimization. Figure 13 shows an abdominal aorta scanned mesh with its global parameterization optimization using an arbitrary valid quadrilateral layout. Observe the final global parameterization quality and nodes relocation. The final parameterization is then used to construct feature-aligned quadrilateral meshes. Due to patch structure of the quadrilateral layout, such quadrilateral meshes can be computed patch by patch according to connectivity constraints.

Isotopological analogous quadrilateral meshes construction from aligned global parameterization
Let us reformulate the previous defined problematic in section 1.2 in a more technical manner. Given a set of input triangulated meshes we strive to find 3D quadrilateral meshes which respect the four following properties: • Pure quadrilaterals with low distortion. • Feature aligned.
• Isotopological with analogous points into other geometric instances.
Pure quadrilaterals with low distortion With the definition of a quadrilateral layout replicating perfectly the topology of the input triangulated surface mesh, it is possible to define a pure quadrilateral mesh. Low distortion is achievied by minimizing the alignment of the global parameterization with the guiding field.
Feature aligned Sharp edges, vertices and principal curvature directions are inportant features of meshes. Alignment to these features is required on all instances to fully conserve significant features. Thus the input quadrilateral layout must contain all the sharp edges and vertices to ensure right feature alignment into all input meshes.
Isotopological with analogous points into other geometric instances A sufficient condition to have isotopological meshes is to hold the same connectivity. A representative quadrilateral mesh of the population with a correct sampling is taken as reference to set the mesh connectivity. It should be noted here that it is not a sufficient condition to compare meshes because of orientation. This is why we want each point of the computed mesh to have an analogous point into all other geometric instances of the considered set. It is achieved globally by patch connectivity and locally by classifying sharp features and applying patch parameterization.
Non-uniform isotropy Mesh non-uniform isotropy provides the only degrees of freedom to have isotopological meshes from different geometric instances. Non-uniform isotropy is intrinsically set because of patch connectivity and discretization.

Reduced order models
As mentioned in the introduction, Reduced Order Models (ROMs) enable real-time analysis due to their very low computational cost during evaluation. In the context of geometric parametric studies, their construction requires accumulating a certain number of geometries depending each on a set of parameters. This set of geometric instances are called snapshots. For ROM construction we need solution vectors with the same dimension associated to all snapshots in order to build matrices required to perform a Singular Value Decomposition (SVD). Isotopological and analogous meshes are thus appropriate.

Algorithmic details
The software Rhinoceros 5 [43] is used for visualization and RhinoCommon.dll [44] for handling various Non-Uniform Rational B-Spline (NURBS) objects. Our algorithm is entirely incorporated into a Rhinoceros 5 Plug-In implemented in VB.NET. C++ processes are called from the Plug-In. Taking a B-Rep object provided from a file or constructed from the Rhinoceros 5 graphical interface, a geometric reduced order model can be built using our fully integrated workflow. We use for that purpose the ROM builder proprietary software developed by ANSYS [1]. For instance, this builder has been used to prevent excessive compression of buttock's soft tissues by bony structures [39] in real-time for paraplegic persons. All computations have been performed on an Intel Core TM i7-6820HQ CPU @ 2.70Ghz with 16Go RAM.

T-shape part
We first introduce the T-shaped part for a rapid understanding of our approach. Using this very simple shape, we define the number and the range of geometric parameters for an arbitrary study, see e.g. Table 1. These meshes come from a standard CAD software and each input triangulated mesh has an attached set of 9 geometric parameters based on the zero reference (Fig. 14a). Then a feature-aware cuboid decomposition is performed involving a pant mesh as topological input (Fig. 14b). Thereafter an optimized aligned global parameterization is computed with the quadrilateral layout and hence gives us the desired quadrilateral mesh with required geometric and connectivity properties. We proceed the workflow on all 37 input geometric instances to generate all isotopological analogous meshes (Fig. 14c).
To avoid large amount of triangulated meshes as input in order to compute an accurate ROM with 9 geometric parameters, we operate a very simple sparse grid sampling of the sets [37,38] (we took 4 snapshots per parameters, involving 37 meshes). ROM is then built using the isotopological snapshots. Once the reduced order model is built, we can generate in real-time all desired geometries handling our ROM evaluation interface (Fig.  15).

Casting part
This casting part is a classical one that can be found in mechanical systems for coupling shafts together. Starting with a feature-aware pants decomposition, we apply a specific cuboid configuration on each pant (Fig. 16c, d). Then the worflow is very similar to the T-shaped model presented in the previous paragraph. The sampling of the 6 input parameters (Fig. 16a, b) was made involving an optimal space filling technique so as to generate 65 triangulated meshes covering as best as possible the R 6 space. Parameters ranges have been defined so as to replicate a real industrial case with large and short variations, see e.g. Table 2. We generate in real-time all desired geometries handling the same ROM evaluation interface (Fig. 17) adapted to the casting part study.

Heart ventricle
Ventricle geometric reduction is dedicated to a specific biological application. Scanned data from medical imaging related to one patient is provided during one cardiac cycle uniformly sampled by 10 meshes. Meshes are called "phases" and they start from 0 to 90 (Fig. 18). Ventricule volumetric capacity during a cardiac cycle can be seen in Fig. 19. Using a specific pipeline well-suited for these particular meshes, a quadrilateral layout  is computed (Fig. 10e).
To be more precise, we modify the topology of scanned meshes from disk to cylinder by adding a little hole at the bottom of the ventricule. Afterwards an optimized aligned global parameterization is extracted using leaflet (i.e., the leaflet of the aortic valve) principal curvature direction interpolation for the geometric cross field. A reduced order model is built using the phase (i.e., time) parameter and isotopological analogous snapshots. For instance two geometries have been evaluated at phase 38.7% (systole, i.e., contraction phase) and 81% (diastole, i.e., expansion phase) in Fig. 20. Mesh linear interpolation can be simply performed if quadrilateral meshes have the required properties. An example of consistent interpolation is given in Fig. 21.

Run times and statistics
Any quadrilateral mesh included in the numerical charts range of parameters can be evaluated in seconds. Table 3 shows some statistics related to the whole workflow. B-Rep parametric geometries are turned into structured quadrilateral meshes in minutes. Then, the isotopological analogous constraints are applied to the population of meshes in order to obtain the same representation. Thanks to the same data structure and analogous properties, an accurate reduced order model is built with a prescribed number of modes. For each study, evaluation time is constant because equivalent numerical operations are performed.

Accuracy of the MEG-IsoQuad Method and geometric parametric analysis
Fidelity of the method can be measured in two steps. The first focuses on the quality of the generated quadrilateral meshes with respect to the input triangulated B-Rep. We refer the readers to Campen and Kobbelt [8] and Bommes et al. [6] to have more details on quality espects of computed quadrilaterals. In Fig. 22(a, left), we show that the quadrilateral mesh replicates the B-Rep with fidelity. Viewed contrasted shapes are relative to numerical error. Sharp features are conserved during surface conversion. Notice that for a shape with less localized curvature, the structured output mesh is still satisfactory for a correct sampling of quadrilateral elements. A geometry is evaluated in an unfavorable case and compared with a validation snapshot in Fig. 22(a, middle) which was not used for the reduced order model construction. The studied case is located in a corner of the 9-dimensional sparse grid, far from known information. Despite this difficutly, the evaluated geometry looks quite good with only 12 modes and 37 snapshots to cover this high-dimensional space.
Unwanted curvature effects are related to the lack of input geometries and dependent geometric parameters. Accuracy can be significantly improved by adding snapshots to the reduced order model and modes in the truncated basis. Standard kriging responses surfaces will be thus more faithful. Error is given in Fig. 22(a, right). Secondly, we investigate the truncated basis dimension focalized on a quarter of the casting part. For that purpose, a snapshot is taken and compared with several evaluated geometries approximated by an increased number of modes. In Fig. 22(b, left) 5 modes are used to represent the geometry until 25 in Fig. 22(b, right). Higher is the number of modes, better are the estimated shapes. While global distance error decreases, the error becomes more and more localized.

Comparison and contribution for reduced order modeling with geometric parameters
Most of current approaches to build reduced order models with geometric parameters are based on the 2 following methods: mesh morphing and immersed methods. Remark that, for elementary geometric parameters manual techniques are legitimate, see e.g. Lu et al. [37]. To our best knowledge, no recent methods using surface segmentation and global parameterization tools exist to construct accurate geometric numerical charts. Our method employs techniques from the graphics community and adapt them to the needs of reduced order modeling applications with geometric parameters. MEG-IsoQuad technique preserves input mesh features and build isotopological mesh structures directly usable in various SVD algorithms.
Firstable, mesh morphing is for instance operated to study biological shapes. Femur based Principal Component Analysis (PCA) requires mesh projection and hence morphing tools to compare instances with the same data structure, see e.g. Hraiech [24] and Grassi et al. [20]. Morphing is associated with distorted and overlapped elements and as a consequence this method is not robust for large variation of parameters.
On the other hand, researchers employ immersed methods to deal with the same mesh connectivity but different geometries. However, these approaches [13,28] also suffer from some drawbacks such as a clear definition of the integration rule for the cut finite elements. Geometric comparison between the reduced order model evaluation and a validation snapshot for the same set of parameters P acc T −shape (middle). Distance error analysis of (a, middle) mapped onto the reduced order evaluation (right). b Reduced order model truncated basis influence on the casting part. Used snapshot has parameters P acc Casting = {L = 80, Rb = 120, Rh = 35, Ri = 20, Hb = 15, Hh = 140}. Distance error is analyzed between the reduced order model evaluation and the corresponding snapshot and mapped on the reduced order model evaluation. Error with a basis constructed from the first 5 modes (left). Error using 15 modes (middle). Error using 25 modes (right)

Limitations
Surface segmentation is highly dependent on the mesh discretization and also on the intrinsic geometry. Hence, the smart model decomposition step may have inputs which are problematic. For instance, pants decomposition is known to produce highly deformed segmentations or unwanted geometries [53]. To overcome partially this issue, we propose to use both topological tools and geometric criteria. Cuboid decomposition deals with consistent pants. The quadrilateral layout computed from the cuboid decomposition step is generally coherent but it may fail when exotic layouts are requested: this is the limitation of the constrained border parameterization to determine layout arcs.
Global parameterization step, including cross fields is entirely performed by techniques from the paper of Campen and Kobbelt [8]. We refer the readers to the limitation section of their work.

Guarantees
In spite of limitations expressed above, once quadrilateral meshes are generated, the reduced order modeling operation can not fail.

Conclusion
We have presented in this work a method for designing 3D quadrilateral snapshots adapted for data-based applications. Quadrilateral layouts are computed from a geometry-featureaware pants-to-cuboids decomposition. These layouts replicate the input mesh topology and contain relevant isotopological information. Depending on the geometry, the pants cuboid configurations are chosen by the user to lie on the key features. Using aligned global parameterization, arcs and nodes of the suitable quadrilateral layout are optimized to achieve low overall patch distortion and alignment with sharp features. A key contribution of the present work is the automatic construction of a relevant quadrilateral layout per pant, setting automatically well-adapted cross field singularities in order to compute a good aligned global parameterization.
This allows us to extract a quadrilateral mesh satisfying sharp features, principal curvature directions and possibly boundary constraints. With such features, our method allows to produce isotopological meshes of all the geometric instances. They are created during the offline phase of ROM based parametric studies. Throughout design process of quadrilateral meshes, we have handled data of all geometric instances in a isotopological manner, no matter the input geometry. This helps us to define a set of isotopological analogous meshes using one mesh as reference. The previous sets were successfully used to construct accurate geometric parametric models. These data-based models are evaluated for any given parameter set in near real-time to give the specified geometry.
Abbreviations CAD: Computer aided design.