Computer-Aided Design 109 (2019) 49–59

Contents lists available at ScienceDirect

Computer-Aided Design
journal homepage: www.elsevier.com/locate/cad

Automatic and high-quality surface mesh generation for CAD models✩

∗
Jianwei Guo a , Fan Ding a , Xiaohong Jia b,c , , Dong-Ming Yan a
a
National Laboratory of Pattern Recognition, Institute of Automation, Chinese Academy of Sciences, Beijing 100190, China
b
KLMM, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing 100190, China
c
School of Mathematical Science, University of Chinese Academy of Sciences, Beijing 100049, China

article info a b s t r a c t

Article history: In this paper, we present a fully automatic framework that tessellates industrial computer-aided design
Received 16 May 2018 (CAD) models into high-quality triangular meshes. In contrast to previous approaches that are purely
Accepted 19 December 2018 parametric or performed directly in 3D space, our method is based on a remeshing algorithm that
can achieve accuracy and high-quality simultaneously. Given an input standard CAD model, which is
Keywords:
CAD tessellation represented by B-rep format, we first rebuild the parametric domain for each surface patch according to
High-quality mesh an initial triangulation, and then the boundaries of the parametric domain are retriangulated by exploiting
Remeshing the Constrained Delaunay Triangulation (CDT). In the second stage, the 2D triangulation is projected back
to the 3D space, and a modified global isotropic remeshing process is applied, which further improves
the regularity and angle quality of the tessellated meshes. Experiments demonstrate the validity of the
proposed approach and its ability to generate high-quality meshes. Moreover, we evaluate our technique
and compare it with state-of-the-art CAD model tessellation approaches.
© 2018 Elsevier Ltd. All rights reserved.

1. Introduction algorithms (e.g., MeshFix [8] and PolyMender [9]) experience dif-
ficulties in dealing with these problems. Moreover, even some
Tessellating Computer-Aided Design (CAD) models into discrete CAD models have correct topology and geometry, and also contain
representations is important in various downstream applications many small-scale features (e.g., tiny surfaces or thin blending
[1,2] ranging from engineering to scientific research, such as strips, as shown in Fig. 13) that do not contribute to downstream
numerical simulations, rapid prototyping, computational fluid dy- applications. Preserving such small features creates dense triangu-
namics, and visualization, to name a few. The tessellation pro- lations locally.
cess aims to generate a discrete mesh that approximates a CAD From a practical point of view, the ability to generate high-
model with simple discrete elements, e.g., triangles/quads for a quality meshes is the key to the success of numerical analyses
surface mesh and tetrahedra/pyramids/hexahedra for a volumetric [10]. A good mesh usually achieves balance between quality and
mesh. Among all types of tessellations, triangular surface mesh computation time, whereas low-quality (extremely thin or badly
generation attracts the most attention because of its simplicity distorted) elements often hinder the solution convergence and
and flexibility, which is also a prerequisite for many volumetric increase analysis errors. Furthermore, surface meshes are used as
tetrahedral meshing procedures [3]. input to volume mesh generators, thus considerably influencing
Although numerous meshing methods are available, automated the generation of quality volumetric meshes and furthering the
generation of high-quality meshes remains a challenge [4,5]. Even numerical results.
with modern commercial software (e.g., Ansys and Hypermesh) In this paper, we focus on generating isotropic surface meshes
and open-source packages (e.g., NetGen [6] and Gmsh [7]), gen- from CAD models that enhance existing approaches in several as-
erating correct and satisfying meshes is still a time consuming pects. The goal is to obtain well-shaped tessellations of CAD models
process that involves an excessive amount of human effort. As the with high accuracy, while considering the simplicity of patch lay-
complexity of the constructed CAD products increases, existing out and feature preservation. On the one hand, the system robustly
methods cannot always achieve accurate output due to incorrect, handles complex CAD models even with the (near)-degenerate
degenerate or ambiguous geometric designs in the CAD system, cases (Section 4.1) that commonly occur in the conceptual-to-
as illustrated in Fig. 1. Even post-processing with mesh repair early design of industrial products. We address this issue by first
detecting such configurations in an initial triangulation and then
✩ This paper has been recommended for acceptance by S.J. Owen removing non-manifold edges or inserting ‘‘anchor points’’ in the
∗ Corresponding author. parametric domain. Then the use of constrained Delaunay triangu-
E-mail addresses: jianwei.guo@nlpr.ia.ac.cn (J. Guo), xhjia@amss.ac.cn (X. Jia). lation (CDT) can generate correct results. Moreover, we produce

https://doi.org/10.1016/j.cad.2018.12.005
0010-4485/© 2018 Elsevier Ltd. All rights reserved.
50 J. Guo, F. Ding, X. Jia et al. / Computer-Aided Design 109 (2019) 49–59

triangulated by using any 2D mesh generation procedure, and then

the mesh is mapped back onto the original surface. However, this
approach introduces size stretching and shape distortion of the
elements in degenerated or badly parameterized surfaces when
the mapping is not affine.
(3) Hybrid approaches. Considering the advantages and disad-
vantages of the previous approaches, recent algorithms utilize
3D geometric and 2D parametric information to produce high-
quality meshes. Starting from an approximated boundary repre-
sentation (using STL file format), Béchet et al. [34] propose an
adaptive surface mesh generating method. The principle used to
generate the mesh is based on the Delaunay method, which is
associated with refinement and smoothing operations. Wang et al.
[35] propose a remeshing-based algorithm, called EQSM, to gen-
Fig. 1. (a) An input CAD model; (b) Tessellation result by using an open-source erate unstructured triangular meshes. An initial triangulation of
toolkit called Open Cascade Technology (OCCT) [13]; (c) Tessellation result of NetGen the points defining the intersection curves of each surface region
[6]; (d) Mesh repair result [9] taking (b) as input. is generated in the parametric domain. These meshes are then
mapped onto 3D space and refined by basic remeshing operations.
However, the input of this method requires several restrictions,
high-quality meshes for complex geometric shapes with smooth e.g., the mapping of each patch should be bijective everywhere,
sizing gradation control. In this step, we first perform patch clus- and the boundary should be closed. Marchandise et al. [36,37]
tering and constraint simplification to suppress undesired internal present several different parameterization techniques for quality
features that lead to low-quality elements. Then, we exploit an surface remeshing, where the implementation of those mappings,
improved version of a real-time isotropic remeshing technique as well as the boundary conditions, have been presented in a
[11,12] that applies a series of local operators for mesh optimiza- comprehensive and unified manner. Aubry et al. [38] present a
tion. In summary, the main contributions of this work include the novel parametric surface meshing technique, which includes semi-
following: structured boundary-layer surface mesh generation and advanc-
ing front point creation approaches. Thereafter, [39] extends this
• an automatic and robust surface tessellation framework that approach to obtain a robust anisotropic tessellation technique.
could generate correct surface meshes from complex indus- However, this approach only handles NURBS-based geometry. Our
trial CAD models, especially by handling two common cases approach also falls into this category. Compared with the afore-
of degeneracies encountered in many industrial CAD geome- mentioned methods, our approach can deal with more complex
tries during mesh generation. composite parametric surfaces, as well as achieve higher quality.
• improvement of a state-of-the-art remeshing approach to
generate high-quality surface meshes, particularly for angles Model repair. Although various mesh generation methods exist,
and regularity. the tessellated CAD models may still exhibit topological or geo-
• introduction of simple yet effective heuristics to merge the metrical defects. Thus, mesh repair, which converts an imperfect
small patches and remove the internal feature boundaries, CAD model into a clean and manifold closed triangular mesh, is a
which are basically irrelevant to an overall simulation of the highly important task. This topic has been extensively studied [40–
CAD model. Furthermore, this approach allows generating a 42] and a website featuring freely obtainable implementations of
smooth geometry for improved meshing. several repairing methods is available at http://www.meshrepair.
org. According to the approach employed, repairing algorithms can
be distinguished by local correction or global remeshing. Typical
2. Related work local approaches [8,43–46] are performed directly on the input
tessellation and only modify the mesh in a small region around
CAD mesh generation. Over the past decades, numerous tessel- individual defects and flaws. By contrast, global methods [47–49]
lation methods have been available for meshing analytical 3D are typically based on a complete remeshing of the input by using
surfaces by different geometric primitives, e.g., triangles, quadri- some intermediate data structures (such as volume representation
laterals or general n-sided polygons. In this section, we restrict the [9,50,51]). However, no available approach can address all the
discussion to the most related works, which focus on triangular defects because the type of defects depends on the upstream and
mesh generation. For more comprehensive discussions, we refer downstream applications in a given scenario. Most repairing algo-
the reader to the survey [14] and textbooks [2,15,16]. rithms focus on certain defect types and ignore or even introduce
A recent survey by Shimada [17] discusses current trends and other flaws. Furthermore, they do not consider improving the
issues in meshing and geometric processing for computational triangle quality of the output mesh, and usually a post-process is
engineering analyses. In general, these triangulation algorithms required in practical applications.
can be categorized into three types:
(1) Direct approaches, which work directly on the 3D surface 3. Preliminaries and overview
in the physical space. The commonly used approaches include
Delaunay triangulation [18–20], Advancing Front Technique (AFT) In this section, we present an overview of the proposed ap-
[21–23], and octree-based approaches [24,25]. However, direct proach. To precisely explain our techniques, we first provide the
methods cannot guarantee the validity of the mesh, e.g., Delaunay- related terminology used in the paper.
based methods have difficulty in recovering surface boundaries
[26], whereas methods using AFT can encounter problems associ- 3.1. Definitions
ated with colliding fronts [27].
(2) Purely parametric approaches (also called indirect approaches) In solid modeling systems, a CAD model is typically represented
[28–33]. This type of method utilizes bijective mapping between by a collection of trimmed parametric surface patches {Pi }ni=1 .
the surface and a parametric domain. The parametric space is Each patch Pi is represented by an analytically defined mapping,
 J. Guo, F. Ding, X. Jia et al. / Computer-Aided Design 109 (2019) 49–59 51

Fig. 2. Our algorithm converts an input standard CAD model (881 patches) into a high-quality triangular mesh. At the bottom, we show the histograms of the angle (blue)
and triangle quality (green) distributions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. The pipeline of our framework based on a two-stage remeshing approach.

denoted as S(u, v ) = (x(u, v ), y(u, v ), z(u, v )), from a bounded 3. Parametric remeshing: Based on initial triangulation, we re-
domain of R2 , called parametric space, into a 3D coordinate system build the parametric domain for each patch and reapply
R3 , called object space. Furthermore, each patch Pi is uniquely CDT to preserve the boundaries for each parametric domain.
identified by its patch ID i and equipped with a surface type Ci , such Then, the 2D triangulation is projected back to 3D space.
as plane, cylinder, sphere, and B-spline surface. 4. Global remeshing: To improve the mesh quality, a global
To explore the topologic and geometric entities in the CAD isotropic remeshing process is applied. The user-specified
model, several third-party CAD kernels are available, e.g., OCCT meshing size is also satisfied in this stage.
[13], Computational Analysis Programming Interface (CAPRI) [52]. In
our implementation, we use OCCT to access the geometry of shapes Each step involves a precise task while gathering enough in-
represented in B-rep format. It also provides access to the topol- formation for the next step. Details of each step are given in the
ogy of the model, including the connectivity information between following sections.
different patches. However, to operate our mesher, the user is not
required to have a deep knowledge of OCCT. The only parameter for 4. Methodology
OCCT that can be specified by the user is the tessellation precision,
and this parameter is set as 100 (the default value) in all our 4.1. Initial triangulation
experiments.
The core idea of our algorithm is that we only perform surface
3.2. Overview remeshing on a discrete model to achieve high efficiency because
the edges and facets that depict the discrete surface mesh have lin-
The input to our framework is a standard CAD model, G = ear equations. Therefore, a suitable initial triangulation is required,
{Pi }ni=1 , which consists of n trimmed parametric surface patches but it should be generated quickly, and the approximation error
that can meet non-smoothly on boundaries. In addition, the user is should be small.
required to specify the target edge length, ltarget , as the desired size In our approach, we exploit OCCT to generate the initial trian-
of the approximation mesh. The output is a 2-manifold watertight gulation. Although this toolkit is designed to generate triangles for
triangular mesh M = {ti }m i=1 , consisting of as much as possible visualization purposes only (leading to non-conforming meshes at
of quasi-equilateral triangles ti . As shown in Fig. 3, our algorithm face boundaries) and has no gradation in element size, it can con-
undergoes four main stages:
trol the distance between the initial tessellation and the original
1. Initial triangulation: Initially, the CAD model G is converted analytical surfaces. Thus, geometric fidelity is preserved, i.e., the
into a simple triangulation by using OCCT, where the geo- approximation error is small.
metric fidelity (measured as approximation error) is guar- The initial triangulation is not closed yet because each patch Pi
anteed. is tessellated into a patch mesh Mi independently. This condition is
2. (Near-)degeneracy handling (optional): We detect the incor- favorable because our parametric remeshing stage is still operated
rect triangulation caused by poorly designed geometry and in a patch-wise manner. Besides, we collect the original vertices of
repair them by removing non-manifold edges or inserting the input geometry from the B-rep model. These vertices are then
‘‘anchor points’’, and then compute CDT in the parametric identified in the initial triangulation and marked with a ‘‘locked’’
domain. flag. Thus, if a vertex is marked as ‘‘locked’’, it remains fixed in later
52 J. Guo, F. Ding, X. Jia et al. / Computer-Aided Design 109 (2019) 49–59

Fig. 4. 2D illustration of degeneracy caused by tangency.

Fig. 6. An example of a highly distorted initial parametric space.

a threshold (we set it to be 0.04ltarget where ltarget is the target edge

length), we insert an ‘‘anchor point’’, D, on curve ⌢ BC, which
is the closest point to point A (Fig. 5(c)). With the help of this
‘‘anchor point’’, an edge AD prohibits the generation of intersecting
Fig. 5. Solution to problem caused by curve proximity. triangles.

4.2. Parametric remeshing

steps to preserve the geometric feature unless its flag is changed
at a certain stage. We explain this phenomenon and the moment a As each patch Pi has a 2D parametric representation as de-
feature vertex is unlocked in Section 4.3. signed, therefore a straightforward way to tessellate is to mesh the
(Near-)degeneracy handling. As the precision of industrial prod- 2D domain in this parametric space and then map it back to the 3D
ucts increases, a complex CAD model usually contains a large object space. However, in many cases, such surface mesh genera-
number of surface patches; for example, even the shell of a cell- tion is easily corrupted. The reason is twofold: (1) the patches in
phone can contain hundreds of patches (see Fig. 12). During the a CAD model are often not topologically connected [53], i.e., small
construction of the CAD model some operations may generate gaps or overlaps between the neighboring patches exist, and the
small artifacts that may not be noticed by the designers. Moreover, parameterization may be different along one common boundary.
the incorrect CAD data arising from translation errors between These conditions will result in non-conformal triangles; and (2) the
different CAD products and formats, or even from limitations inside parametric space of one patch can be highly/extremely distorted
CAD kernels, can cause meshing failures. In this paper, we address (see Fig. 6). Therefore, a valid mesh in the parametric space may be
two cases of degeneracies encountered in a wide range of complex invalid or at low quality when mapped in a 3D space [54].
CAD models in practical engineering applications. To address these problems, we apply a planar parameterization
(1) Tangency case: Non-manifold vertices and edges usually method to rebuild the parametric domain based on the initial
appear due to patch tangency, in which an intersection between triangulation. Thereafter, we recover the boundary edges corre-
the inner and outer boundaries occur. Fig. 4(a) shows such tan- sponding to the original curves and perform edge sampling. Finally,
gency degeneracy, which has barely attracted attention in previous the CDT is used to finish the 2D triangulation. Our parametric
methods. In this case, a planar patch P1 is tangent to a cylindrical remeshing is detailed as follows.
patch P2 . At the position of P1 where the inner and outer bound-
Patch reparameterization. We first rebuild a one-to-one mapping
aries are tangent, two vertices that coincide or are very close exist,
from each patch mesh Mi to a simple 2D domain, using the Scalable
as shown in segment BD in Figs. 4(b) and (c). As a result, this
Locally Injective Mappings (SLIM) parameterization [55]. The result
intersection forms non-manifold triangles because OCCT cannot
is a pair of parameter coordinates (u, v ) for each vertex of the initial
distinguish the inner and outer regions at this position. To avoid
triangulation. SLIM is robust and fast by efficiently minimizing a
such non-manifold triangles while keeping the shape boundaries,
nonlinear, conformal or isometric distortion energy. It is guaran-
we first detect the planar patch P1 where the tangency happens by
teed to produce optimized maps without any flipped triangles. Fur-
sorting out the original vertices that are linked to multiple (more
thermore, it corresponds to a locally injective mapping with a free
than two) boundary edges. Then segment BD is removed from the
border. Therefore, it is suitable for our method because a complex
constraints when building CDT for this patch, so that only one
patch mesh with arbitrarily shaped borders can be parameterized.
outer boundary A − B − C − D − E − F − G − A exists. In other
adjacent patches P2 and P3 , segment BD is retained to preserve Boundary curve recovery and sampling. From the initial triangu-
the boundary. Therefore, they can form a closed manifold mesh in lation, we recover the 3D boundary edges that correspond to the
the later stage. original curves of the input CAD model. This step aims to preserve
(2) Curve proximity case: Another common case is that one point the geometry of boundaries in the final mesh. For each patch mesh
of the patch is too close to one curve, as shown in Fig. 5(a). Although Mi , we use a tracing scheme for recovery. Here we suppose that
this case is not a type of degeneracy, it results in intersecting the mesh is represented by the half-edge data structure. As shown
triangles when the boundary sampling rate is sparse (Fig. 5(b)). in Fig. 7, starting from a half-edge h1 whose vertex is marked as
Thus, it requires dense but unnecessary samples to overcome this ‘‘locked’’, we recursively visit its previous border half-edge, until
problem. In our approach, we calculate the Euclidean distance another half-edge whose vertex has a flag of ‘‘locked’’ is met. Then
between a concave point A on one border and its nearest border these piecewise linear segments form a boundary curve ci =
edge to avoid such dense sampling. If the distance is shorter than {h1 , . . . , hk }. After processing all half-edges, we collect a set of
 J. Guo, F. Ding, X. Jia et al. / Computer-Aided Design 109 (2019) 49–59 53

Fig. 8. (a) Initial triangulation using OCCT causes gaps or overlaps between neigh-
boring patch meshes. (b) Our boundary curve recovery and sampling could generate
consistent vertices along the shared curves.

coarse triangulation by a parametric remeshing approach. In the

2D space, we organize the sampled points on boundary curves into
a bounded region, which is defined by a Planar Straight Line Graph
Fig. 7. Our tracing scheme for recovering boundary edges. (a) An input parametric
(PSLG). The segments of the PSLG are considered as constraints and
curve formed by two points; (b) Half-edge data structure for this curve (we show
only one side); (c) Half-edge h1 whose vertex is marked as ‘‘locked’’ is first identified,
are fed to our algorithm to build a CDT. Then, an improved mesh
shown in red; (d) Other half-edges belonging to this curve are recursively traced is obtained by invoking the Delaunay refinement method with
until h7 is met. (For interpretation of the references to color in this figure legend, default shape and size criteria1 (default shape criterion B = 0.125,
the reader is referred to the web version of this article.) size criterion s = 0.5). We do not need to tune the criteria be-
cause we will satisfy the user-specified meshing size in the global
Algorithm 1: Boundary curve recovery remeshing stage. Finally, the 2D mesh of each patch is mapped
back to the 3D space, and an improved tessellation with correct
1 Input: a coarse patch mesh Mi ; conformity is obtained, as shown in Fig. 8(b).
2 Output: recovered boundary curves set, C = {ci }ni=1 ;
3 foreach halfedge h do 4.3. Global remeshing
4 if h is a border edge and its vertex is marked as "locked"
then Although 2D Delaunay refinement is performed, the quality
5 construct an empty boundary curve ci = ∅; of parametric remeshing is still unsatisfied. To improve the reg-
6 halfedge h′ ←− h; ularity and angle quality of the final mesh, another high-quality
7 do remeshing is required. In our approach, we leverage a real-time
8 add h′ into ci ; isotropic remeshing technique inspired by [11] to refine the mesh.
9 h′ ←− previous halfedge of h′ ; Besides, an industrial CAD product contains many small surfaces
10 while the vertex of h′ is not marked as "locked"; and internal ‘‘feature’’ boundaries because it is usually modeled
11 add ci into C ; via various Solid Boolean Operations, translations or corrections.
Thus, before remeshing, we first remove these small geometric
features, in which most unintended discontinuities reside. This
global patch-independent remeshing can achieve better quality
elements by ignoring the internal CAD patch limitations.
boundary curves, C = {ci }ri=1 . The pseudo-code of this step is
illustrated in Algorithm 1. Patch mesh clustering. As illustrated in Fig. 9, a complex CAD
Currently, each curve c of the input geometry has been identi- model often contains many small surface patches that is not related
fied in its two incident patch meshes (the neighboring relations to the geometry of the model, thus, the arrangement of these
between different patches have been stored in advance in the patches cannot reflect the intention of the designer. The proximity
initial triangulation), corresponding to ci and ci′ , respectively. How- features among these patches lead to low quality elements around
ever, one curve can be discretized differently in different patches them because they contribute to constraining the resulting mesh.
because OCCT tessellates each patch independently (see Fig. 8). As To suppress these features, we cluster such patches into more
a result, this discretization causes gaps or overlaps when merging desirable subregions using a hierarchical agglomerative clustering
adjacent patch meshes. To achieve consistent discretizations for algorithm. First, for each patch mesh Mi consisting∑ of a list of∑tri-
each shared curve c, we produce the same number of uniform angles {tj }, we compute its weighted centroid qi = |tj |btj / |tj |
samples on ci and ci′ , respectively. We denote the length of curve ntj / |tj |, by averaging the facet
∑ ∑
and weighted normal n ⃗i = |tj |⃗
c as lc , and the curve sampling interval as δ = α ltarget (α is barycenter btj and normal n ⃗ tj of each triangle with area |tj |. Then,
set to 0.6 −⌈ 0.⌉8). Thus, the number of points to be sampled is we group neighboring patch meshes within a close distance and
N = max(2, lδc ), where the max(a, b) function returns the greater with similar normals. Specifically, two patch meshes, Mi and Mj ,
value of a and b. Finally, we generate N samples uniformly on the are grouped together iff:
piecewise linear segments of ci and ci′ . The parameter coordinates qi − qj  < εD ,
 
(u, v ) of each sample is interpolated using the values defined on
the vertices of ci and ci′ . ⃗ j > εA ,
⃗i · n
n (1)
Triangulation. Based on the rebuilt parametric space and recov-
ered boundary curves, we now refine each patch mesh of the 1 http://doc.cgal.org/latest/Mesh_2/.
54 J. Guo, F. Ding, X. Jia et al. / Computer-Aided Design 109 (2019) 49–59

Fig. 10. Results before (left) and after (right) edge flipping by angle. Green circles
indicate the entangled triangles caused by splitting long edges.

Fig. 9. Illustration of how clustering approach provides us a mesh that is simpler

than merely mimicking the CAD topology: (a) patch meshes before grouping, and
of the vertices on feature edges. Thus, we propose valence
(b) patch meshes after grouping. optimization for feature vertices. We observe that a low-
quality triangle with one obtuse angle is usually situated
around a feature vertex with valence smaller than 6 (actu-
ally, the valence is 4 or 5, because valence 3 is impossible
where · is the dot product between two normal vectors. We set
on feature edges). Besides, a feature vertex with a valence
εD = 2 and εA = 0.8 by default.
larger than 6 results in small angles. To improve the valence,
After this step, all patch meshes are assembled by merging
we propagate the valence optimization of feature vertices
common vertices to yield a single mesh. In this merging process,
to inner vertices by edge splitting and collapsing. Fig. 11
each triangle edge shared by two groups is marked by a ‘‘feature’’
illustrates our solution. Vertex A is a feature vertex with non-
flag.
6 valence and its two incident half-edges h1 and h2 are on
Constraint simplification. The above clustering approach only the same side of A. We denote θ as the angle between h1
works well for small patches. Many undesired internal constraints and h2 , nt as the number of triangles between h1 and h2 . The
remain between large patches, which lead to stretched triangles. desired number of triangles is computed as nt ′ = ⌊ 60θ ◦ + 0.5⌋.
To remove them, we consider all the ‘‘feature’’ edges along one Then two cases are considered: (1) If nt < nt ′ , as shown
common boundary between two groups, and compute the dihedral in Fig. 11(a) where θ = 180◦ , nt = 2, nt ′ = 3, we split
angle between the two incident triangles of each ‘‘feature’’ edge. If the opposite edge BD of A. The reason for selecting BD is
all of the dihedral angles are lower than a threshold (default value that the squared difference between the valence of vertex C
is 5◦ ), then the ‘‘feature’’ edges on this boundary, as well as their and the optimal value 6 is smaller than that of E. However,
‘‘locked’’ vertices, are released. splitting introduces another vertex G with valence 4, thus,
a local edge flip is performed to improve the valence of G
Isotropic remeshing. The general pipeline of our global optimiza- (see Fig. 11(b)). (2) If nt > nt ′ , as shown in Fig. 11(c) where
tion is similar to that in [11,12], but we further improve the quality nt = 4, nt ′ = 3, we collapse the opposite edge CD of the
of this algorithm, especially for CAD models. To make the paper smallest angle iteratively until nt = nt ′ . But note that the
self-contained, each local operation is outlined sequentially as case of nt > nt ′ is rare because above edge collapse and
follows (⋆ indicates the newly added or improved operations): valence optimization for inner vertices can largely prevent
this case from happening. Finally, after processing each side
• Edge split. Given a target edge length ltarget , we split all edges of the feature vertices, we can significantly improve the mesh
at their midpoint that are longer than 43 ltarget . quality near feature edges.
• Angle optimization⋆. We flip an edge if the sum of its two • Tangential Laplacian smoothing. To optimize vertex distribu-
opposite angles are larger than 180◦ . This operation is per- tion, we relocate the positions by computing the Optimal
formed because the initial tessellation forms many edges that Delaunay Triangulation (ODT). Each vertex is moved to the
are considerably longer than ltarget . As a result, the above edge average pi of the barycenter btj of its incident triangles tj ∈ T ,
split operation generates entangled triangles around these weighted by the triangle area:
long edges (see Fig. 10), which can lead the algorithm to be ∑
t ∈T |tj |btj
stuck in a loop of edge collapsing and splitting. By contrast, pi = ∑j , (2)
our angle optimization could generate better shaped triangles tj ∈T |tj |
to avoid this problem and could further accelerate the entire
remeshing process. This operation is only performed in the Then the new position pi can be optionally projected back
first iteration, because the inputs to later iterations are al- onto the original surface patch to decrease approximation
ready well-shaped. Therefore, this operation will not conflict errors. Since projecting any 3D point back to the underlying
surface patch is very slow in OpenCASCADE, to achieve a
with the following valence optimization in later iterations.
good balance between the meshing speed and approximation
• Edge collapse. We collapse all edges shorter than 45 ltarget .
error, we only apply the projections in the first iteration.
• Valence optimization for inner vertices⋆. A non-feature edge is
flipped if this operation meets two conditions: (1) decreases Typically, 4 − 8 iterations are enough to achieve a good tradeoff
the squared difference of the valences of the four vertices of between the mesh quality and the performance of the algorithm.
the two incident triangles to their optimal value 6 (assuming Finally, after all iterations, we again perform edge flipping by
that no border edges for closed CAD models exist); (2) does angles. We found that this operation improves the angle quality,
not increase the projection distance from the midpoint of this without violating the valence regularization.
edge to the input CAD model. The former improves valence
regularization whereas the latter maintains the approxima- 4.4. Adaptive mesh generation
tion error.
• Valence optimization for feature vertices⋆. Since the feature In engineering and scientific computation, an adaptive tes-
edges cannot be flipped, it is difficult to improve the valence sellation is preferred because the mesh size is smaller near the
 J. Guo, F. Ding, X. Jia et al. / Computer-Aided Design 109 (2019) 49–59 55

are inherently non-smooth. Thus, we use the curvature controlled

metric to generate adaptive surface mesh.
In our pipeline, we only need to modify the isotropic remesh-
ing stage by replacing the constant target edge length with the
adaptive sizing field ρ (x) that is sensitive to local curvatures. To
compute the local target edge length (for edge split and collapse)
and weighted areas (for tangential smoothing), we apply a warp
method for mapping. First, we compute the average value E(ρ ) of
the sizing field. Then, for any edge e = (xa , xb ), the desired local
target edge length is defined as:
2E(ρ )
l′target = ltarget . (3)
ρ (xa ) + ρ (xb )
For the tangential relaxation, the barycenter computation in Eq. (4)
Fig. 11. Illustration of the valence optimization for feature vertex. is replaced by:
∑
t ∈T |tj |w(btj )btj
pi = ∑j , (4)
tj ∈T |tj |w(btj )
small features to achieve high geometrical accuracy and is larger
ρ (x )+ρ (x )+ρ (x )
elsewhere to avoid increasing mesh counts unnecessarily. Our al- where w (btj ) = a b
3E(ρ )
c
is the sizing field at the barycenter
gorithm can be modified to generate adaptive meshes by applying btj of triangle tj = (xa , xb , xc ).
a sizing function, ρ (x), that defines a target element size at every
point on the model. The sizing function can be specified by the 5. Experimental results
user input, or derived from geometrical features. In the remeshing
In this section, we first present a number of experimental re-
literature, the preference is local feature size (the distance to the sults that demonstrate the effectiveness of the proposed algo-
medial axis of the model) or curvature. However, the local feature rithm. Then we provide a complete comparison with state-of-the-
size is suitable for smooth shapes, whereas most CAD models art approaches in terms of the meshing validity and quality. Our

Fig. 12. A gallery of examples generated by our framework. For each model, we show the input CAD geometry, tessellation result, zoom-in views, and histograms of the
angle (blue) and triangle quality (green) distributions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this
article.)
56 J. Guo, F. Ding, X. Jia et al. / Computer-Aided Design 109 (2019) 49–59

Fig. 13. Tessellation results without (middle) and with (right) the operations of patch mesh clustering and constraint simplification.

Fig. 15. Comparison on ‘‘Receiver’’. Each row shows the surface mesh, cutaway
view of tetrahedral meshing, and slivers with dihedral angles smaller than 10◦ .

Fig. 14. Comparison with NetGen [6] and Gmsh [7] on the ‘‘Airplane bracket’’
dataset. From the output of Gmsh, we failed to generate a tetrahedral mesh. For
NetGen and our method, we show the surface mesh, cutaway view of tetrahedral
meshing, and livers with dihedral angles smaller than 10◦ .

algorithm is implemented in C++ using the CGAL library [56] for

CDT in 2D domain. All results presented in this paper are obtained
on an Intel i7-3770, 3.40 GHz desktop with 16 GB memory.
Fig. 16. Remeshing comparison with RAR [12] on ‘‘Iron man mask’’. The top row
5.1. Tessellation results (from left to right) is the input CAD, our remeshing result, and RAR’s result. The red
color indicates obtuse triangles and the green color shows triangles with θmin < 30◦ .
The bottom row shows histograms of the angle (left) and triangle quality (right)
We have evaluated our method on a wide range of CAD models distributions. (For interpretation of the references to color in this figure legend, the
of different complexities, ranging from CAD components to me- reader is referred to the web version of this article.)
chanical assemblies. Figs. 2, 3, 12, and 13 illustrate some selected
tessellation results, as well as the histograms of the angle and
triangle quality distributions. The quality of a triangle t is measured minimal angles of all triangles; θ<30◦ % is the percentage of triangles
by Q (t) = √6 pSht , where St is the area of t, pt is the half-perimeter with angles smaller than 30◦ ; θ>90◦ % and θ>120◦ % are the percent-
3 t t
of t and ht is the longest edge length of t [57]. More numerical ages of triangles
∑n with angles larger than 90◦ and 120◦ , respectively;
statistics of the meshing quality are presented in Table 1. Here Qav g ξ = n i=1 |δi − 6| is the measure of mesh regularity, where δi
1

is the average triangle quality in the triangulation; θmin and θmax represents the valence of the ith vertex. The value of ξ is positive,
are the minimal and maximal angles, and θ̄min is the average of the and a smaller value indicates a more regular triangulation. dRMS is
 J. Guo, F. Ding, X. Jia et al. / Computer-Aided Design 109 (2019) 49–59 57

Table 1
Our tessellation qualities. |P | is the number of surface patches in the input; |M| is the number of triangles in the output. The other measurements are explained in Section 5.
We use ⋆ after the model name to indicate that the timing also includes the time cost of degeneracy repair.
Type Model |P | |M| Qav g θmin θ̄min θmax θ<30◦ % θ>90◦ % θ>120◦ % ξ dRMS dH (×10−2 ) Time (s)
Gear (Fig. 12) 229 77k 0.87 2.42 48.81 157.25 0.43 1.68 0.03 0.28 0.017 0.074 8.5
Uniform Helicopter Cockpit (Fig. 12) 72 18k 0.88 7.05 50.23 153.63 1.45 0.68 0.21 0.27 0.006 0.236 5.2
Component⋆ (Fig. 12) 183 42k 0.87 2.55 49.39 168.08 0.29 1.54 0.04 0.34 0.005 0.649 5.9
Car (Fig. 2) 881 93k 0.81 0.15 45.32 178.67 1.28 3.05 0.07 0.41 0.074 0.575 15.7
Cover (Fig. 3) 130 48k 0.86 0.78 48.22 130.82 0.32 1.85 0.01 0.31 0.037 0.146 5.6
Suspension Component (Fig. 3) 249 37k 0.85 3.08 47.60 156.53 0.53 2.75 0.02 0.30 0.102 0.512 6.8
Adaptive Seat⋆ (Fig. 12) 257 23k 0.84 0.01 46.92 178.88 0.44 2.76 0.04 0.33 0.029 0.217 10.2
Mini Servo (Fig. 12) 168 35k 0.86 0.88 48.13 176.04 0.27 2.18 0.03 0.32 0.036 0.232 4.5
Cellphone⋆ (Fig. 12) 985 58k 0.83 0.63 45.42 177.81 1.58 3.91 1.04 0.43 0.031 0.177 20.1
Heart (Fig. 13) 93 26k 0.87 0.01 49.05 177.42 1.62 1.19 0.05 0.32 0.006 0.227 8.2

the root mean square distance, and dH is the Hausdorff distance generates intersecting triangles near very narrow sharp features,
between the output mesh and input CAD model, measured with as illustrated in Fig. 14. By contrast, our algorithm is more robust
the Metro tool [58]. Here the input CAD model is approximated in obtaining correct tessellation results by carefully handling tan-
with a dense triangular mesh. gency degeneracy and performing boundary curve recovery before
From these results, we observe that our approach is robust global remeshing.
to tessellate CAD models and to generate high-quality triangular Next, we conduct a comparison in terms of tessellation qual-
meshes. Furthermore, we can remove the internal feature bound- ity. Besides using surface mesh quality metrics, we also evaluate
aries to achieve better quality, as shown in Fig. 13. Notice that the
these three methods by generating tetrahedral meshes. Then we
output meshes may still have a few low-quality triangles. The rea-
compute two criteria for evaluating the quality of a tetrahedral
son is that a complex CAD model usually contains various features
mesh: dihedral angle ϑ , and radius ratio κ . The radius ratio, κ ,
that affect surface mesh quality, i.e., short curves, curve proximities
is defined as κ = r in , where rin and rcir are the inscribed and
3r
and small input angles. An ill-shaped and jagged boundary will re- cir

sult in poor quality mesh elements [59]. This result is unavoidable circumscribed sphere radius of a tetrahedron, respectively. In our
because the features formed by input boundaries cannot be fully approach, we use TetGen [61] to generate tetrahedral meshes from
suppressed by our clustering approach. Cleaning such features on each surface mesh. Figs. 14 and 15 present the comparison for
the geometry level is another topic about CAD defeaturing [60], tetrahedral meshing. The quality is listed in Table 2, where NetGen
which is beyond the scope of this paper. However, histograms of runs by default in parallel mode. Comparison results verify that our
the angle and triangle quality distributions demonstrate that we tessellation technique yields higher quality surface and tetrahedral
can generate high-quality meshes generally. meshes. Furthermore, NetGen and Gmsh often lead to high-density
triangles (marked with green circles) around tiny features, which
Performance evaluation. Our algorithm can achieve greater ef-
ficiency than approaches that directly work on smooth surfaces, are undesirable in practical simulation.
because most of our geometry computations are performed on Finally, to demonstrate the effectivity of our improved remesh-
the discrete model. Specifically, the processing time mainly com- ing approach, we compare it with original real-time adaptive
prises three components: initial triangulation, parametric remesh- remeshing (RAR) method [12]. Fig. 16 shows the remeshing results,
ing, and global remeshing. Thanks to OCCT, the first step can be in which we use the same initial triangulation and 5 iterations
performed very efficiently. Parametric remeshing is the most time for our method and RAR. We achieve a higher quality mesh,
consuming step in our framework, because if we process each sur- especially near the feature boundaries, due to our operations of
face patch one-by-one, the total processing time is linear with the angle optimization and valence optimization for feature vertices.
number of input surface patches. Fortunately, we found that our Besides, the values of ξ of our method and RAR are 0.35 and 0.48,
parametric remeshing works on each patch mesh independently. respectively.
Thus this step can be parallelized. We use the OpenMP library for
threaded parametric remeshing. Comparison with commercial meshers. Many commercial soft-
In the last step, although we have to carefully deal with the ware, such as Altair HyperMesh,2 Pro/ENGINEER,3 MeshGems,4
sharp features that are common in CAD models, our algorithm are available in the market. These software packages are fully
is still considerably faster by utilizing real-time remeshing algo- optimized and provide various state-of-the-art surface meshing
rithms [11,12]. Table 1 lists the timings for tessellating each model. algorithms to generate high-quality meshes. However, they fail
In addition, performance can be further improved by disabling the when handling complex models with degeneracies, as shown in
back-to-surface projection without considerable loss of quality. Fig. 17. In this example, for the tangency case (labeled by (1) in the
figure), MeshGems and HyperMesh generate degenerate triangles.
5.2. Comparisons Although Pro/ENGINEER can correctly handle this tangency case,
it fails to tessellate the model where the curve proximity is met,
Comparison with open-source methods. We now compare our as labeled by (2) in the figure. Compared with these packages,
approach with two popular open-source software, NetGen [6] and our tessellation framework is much simpler, and it outperforms
Gmsh [7], which are both readily available and widely used in them by successfully tessellating this complex model. Therefore,
the meshing and simulation community. For fair comparison, we our algorithm can be used as a drop-in replacement for the task
tune the parameters of each algorithm to generate surface meshes involving CAD tessellation.
with a similar number of vertices and triangles. In addition, the
sizing function used in all methods is defined based on the cur-
vature. First, we compare the robustness of these three methods. 2 http://www.altairhyperworks.com.au/product/HyperMesh.
As shown in Fig. 1, NetGen cannot handle the degeneracy caused 3 https://www.ptc.com/en/products/cad/pro-engineer.
by tangency. While Gmsh is able to deal with such a case, it easily 4 http://www.meshgems.com/.
58 J. Guo, F. Ding, X. Jia et al. / Computer-Aided Design 109 (2019) 49–59

Table 2
Comparison of graded tetrahedral meshing qualities. The best result of each measurement is marked in bold font. |P | is the number of input patches; |M| is the number
of surface triangles; |τ | is the number of generated tetrahedra.
Model Method |M| |τ | Time (s) Triangle Mesh Tetrahedral Mesh
Qav g θ̄min θ<30◦ % ξ dRMS dH (×10−2 ) ϑmin ϑmax ϑ̄min κmin κ̄ ϑ<10◦ % ϑ<20◦ %
NetGen 26k 146k 16.5 0.824 45.52 3.36 0.35 0.084 0.339 2.17 173.27 69.58 0.005 0.75 0.11 1.38
Receiver (Fig. 14), |P |= 112 Gmsh 20k 106k 34.8 0.829 46.82 2.15 0.61 0.017 0.045 0.67 177.91 69.65 0 0.73 0.14 1.51
Ours 26k 142k 7.6 0.827 46.28 1.01 0.31 0.013 0.161 1.98 172.64 70.04 0.018 0.77 0.09 1.24
NetGen 55k 136k 18.3 0.799 44.14 3.1 0.32 0.165 0.796 0.07 179.28 69.68 0 0.72 0.35 1.71
Airplane bracket (Fig. 15), |P |= 228 Gmsh 63k × 89.5 0.631 38.14 16.28 0.76 0.128 0.473 × × × × × × ×
Ours 49k 125k 9.8 0.817 45.78 2.61 0.34 0.031 0.136 0.27 178.89 69.54 0 0.76 0.29 1.41

Fig. 17. Comparison with some commercial meshers. The green lines represent border edges, while red points represent the vertices of degenerate triangles. (For
interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

5.3. Limitations (e.g., tight radius fillets, turbulators, and flow diversion devices)
and preserve them properly during mesh generation.
Although the tessellation results of our algorithm satisfy the
requirements for most simulation applications, we have not tack- Acknowledgments
led the exact approximation errors. We hope to address this issue
by integrating another parameter (approximation tolerance) into This work is partially funded by the National Natural Science
the remeshing stage in the future. Next, the SLIM parameterization Foundation of China (61802406, 61872354, 61772523 and
method we used only grants local injectivity, thus a global overlap 61620106003), the Beijing Natural Science Foundation, China
could still happen depending on the shape of the patch. Although (4184102), the Open Funding Project of State Key Laboratory of Vir-
we did not meet this case in our experiments, it remains prob- tual Reality Technology and Systems of Beihang University, China
lematic. We would like to explore other more advanced parame- (Grant No. BUAA-VR-17KF-06), and the Open Projects Program of
terization techniques. In addition, since our framework does not National Laboratory of Pattern Recognition (NLPR), China (Grant
involve any geometry repair mechanism, we cannot handle geom- No. 201700020).
etry errors existing in original CAD models, e.g., self-intersection.
Thus, designers are required to modify the geometry manually or References
automatically by using other repairing algorithms.
[1] Klingner BM, Feldman BE, Chentanez N, O’Brien JF. Fluid animation with
6. Conclusion and future work dynamic meshes. ACM Trans on Graphics (Proc SIGGRAPH) 2006;25(3):820–
5.
[2] Frey PJ, George P-L. Mesh generation: application to finite elements. ISTE;
This paper has proposed a novel method for automated mesh 2007.
generation for CAD models. Our algorithm relies on two classical [3] George P-L, Borouchaki H. Delaunay triangulation and meshing.
remeshing approaches to improve the regularity and angle quality [4] Baker TJ. Mesh generation: art or science? Prog Aerosp Sci 2005;41(1):29–63.
of the mesh. Robustness has been emphasized through the entire [5] Ito Y. Challenges in unstructured mesh generation for practical and efficient
meshing process, and by taking care of the degeneracies, unin- computational fluid dynamics simulations. Comput & Fluids 2013;85:47–52.
[6] Schöberl J. Netgen an advancing front 2d/3d-mesh generator based on ab-
tended construction artifacts, and shape boundary preservation. stract rules. Comput Vis Sci 1997;1(1):41–52.
The experiments on a variety of complex models and the compar- [7] Geuzaine C, Remacle J-F. Gmsh: a three-dimensional finite element mesh
ison with the state-of-the-art tessellation packages, demonstrate generator with built-in pre- and post-processing facilities. Internat J Numer
that our approach is simple yet efficient in generating high-quality Methods Engrg 2009;79(11):1309–31,.
surface meshes. [8] Attene M. A lightweight approach to repairing digitized polygon meshes. Vis
Comput 2010;26(11):1393–406.
In our future work, we would like to extend our approach to [9] Ju T. Robust repair of polygonal models. ACM Trans on Graphics (Proc SIG-
a mesh generator that can control anisotropy and directionality GRAPH) 2004;23(3):888–95.
via a metric tensor field. We are also interested in exploring more [10] Du Q, Wang D, Zhu L. On mesh geometry and stiffness matrix conditioning for
accurate sizing functions rather than the simple curvature to de- general finite element spaces. SIAM J Numer Anal 2009;47(2):1421–44.
fine adaptive mesh with smoothed gradient of element scales. [11] Botsch M, Kobbelt L. A remeshing approach to multiresolution modeling. In:
Eurographics symposium on geometry processing. 2004, p. 189–96.
Lastly, although we considered filtering small features, we indeed [12] Dunyach M, Vanderhaeghe D, Barthe L, Botsch M. Adaptive remeshing for
did not distinguish between important and unimportant features. real-time mesh deformation.. In: Eurographics (Short Papers). 2013, p. 29–32.
One good research topic is to recognize such important features [13] Open Cascade Technology, http://www.opencascade.com.
 J. Guo, F. Ding, X. Jia et al. / Computer-Aided Design 109 (2019) 49–59 59

[14] Owen SJ. A survey of unstructured mesh generation technology. In: Interna- [39] Aubry R, Dey S, Mestreau EL, Karamete BK, Gayman D. A robust conforming
tional meshing roundtable. 1998, p. 239–67. nurbs tessellation for industrial applications based on a mesh generation
[15] Bern M, Plassmann P. Mesh generation. Handbook of computational geome- approach. Comput Aided Des 2015;63:26–38.
try. Elsevier Science; 2000. [40] Ju T. Fixing geometric errors on polygonal models: a survey. J Comput Sci Tech
[16] Lo DS. Finite element mesh generation. CRC Press; 2014. 2009;24(1):19–29.
[17] Shimada K. Current issues and trends in meshing and geometric process- [41] Campen M, Attene M, Kobbelt L. A practical guide to polygon mesh repairing..
ing for computational engineering analyses. ASME J Comput Inf Sci Eng In: Eurographics (Tutorials). 2012.
2011;11(2):021008. [42] Attene M, Campen M, Kobbelt L. Polygon mesh repairing: an application
[18] Caendish JC, Field DA, Frey WH. An apporach to automatic three- perspective. ACM Comput Surv 2013;45(2):15.
dimensional finite element mesh generation. Internat J Numer Methods Engrg
[43] Turk G, Levoy M. Zippered polygon meshes from range images. In: Proceed-
1985;21(2):329–47.
ings of the 21st annual conference on computer graphics and interactive
[19] Lévy B, Liu Y. Lp centroidal Voronoi tesselation and its applications. ACM Trans
techniques. SIGGRAPH ’94, ACM; 1994, p. 311–8.
on Graphics (Proc SIGGRAPH) 2010;29(4):119:1–11.
[44] Borodin P, Novotni M, Klein R. Progressive gap closing for meshrepairing. In:
[20] Lévy B, Bonneel N. Variational anisotropic surface meshing with voronoi Par-
Advances in modelling, animation and rendering. Springer; 2002, p. 201–13.
allel Linear Enumeration. In: Proceedings of the 21st international meshing
roundtable. 2012, p. 349–66. [45] Liepa P. Filling holes in meshes. In: Proceedings of the 2003 Eurograph-
[21] Nakahashi K, Sharov D. Direct surface triangulation using the advancing front ics/ACM SIGGRAPH symposium on geometry processing. Eurographics Asso-
method. In: 12th computational fluid dynamics conference. 1995, p. 1686. ciation; 2003, p. 200–5.
[22] Lan T, Lo S. Finite element mesh generation over analytical curved surfaces. [46] Campen M, Kobbelt L. Exact and robust (self-) intersections for polygonal
Comput Struct 1996;59(2):301–9. meshes. In: Computer graphics forum, vol. 29. Wiley Online Library; 2010,
[23] Tristano JR, Owen SJ, Canann SA. Advancing front surface mesh generation p. 397–406.
in parametric space using a riemannian surface definition.. In: International [47] Andújar C, Brunet P, Ayala D. Topology-reducing surface simplification using
meshing roundtable. 1998, p. 429–45. a discrete solid representation. ACM Trans Graph. 2002;21(2):88–105.
[24] Yerry MA, Shephard MS. Automatic three-dimensional mesh genera- [48] Bischoff S, Pavic D, Kobbelt L. Automatic restoration of polygon models. ACM
tion by the modified-octree technique. Internat J Numer Methods Engrg Trans Graph. 2005;24(4):1332–52.
1984;20(11):1965–90. [49] Zhou Q-Y, Ju T, Hu S-M. Topology repair of solid models using skeletons. IEEE
[25] Shostko AA, Löhner R, Sandberg WC. Surface triangulation over intersecting Trans Vis Comp Graph 2007;13(4).
geometries. Internat J Numer Methods Engrg 1999;44(9):1359–76. [50] Nooruddin FS, Turk G. Simplification and repair of polygonal models using
[26] Liu Y, Lo S, Guan Z-Q, Zhang H-W. Boundary recovery for 3d delaunay trian- volumetric techniques. IEEE Trans Vis Comp Graph 2003;9(2):191–205.
gulation. Finite Elem Anal Des 2014;84:32–43. [51] Bischoff S, Kobbelt L. Structure preserving cad model repair. In: Computer
[27] Shewchuk JR. Delaunay refinement mesh generation. In: Tech. rep.. DTIC graphics forum (Proc. EUROGRAPHICS), vol. 24. Wiley Online Library; 2005,
Document; 1997.
p. 527–36.
[28] Sheng X, Hirsch BE. Triangulation of trimmed surfaces in parametric space.
[52] Haimes R. Capri: computational analysis programming interface. a solid mod-
Comput Aided Des 1992;24(8):437–44.
eling based infra-structure for engineering analysis and design. revision 2.0.
[29] Cuillière J-C. An adaptive method for the automatic triangulation of 3d para-
Massachusetts Inst. of Technology, Cambridge, MA; 2004.
metric surfaces. Comput Aided Des 1998;30(2):139–49.
[53] Dassi F, Mola A, Si H. Curvature-adapted remeshing of cad surfaces. Eng
[30] Borouchaki H, Laug P, George P-L. Parametric surface meshing using a com-
bined advancing-front generalized delaunay approach. Internat J Numer Comput 2017.
Methods Engrg 2000;49(1–2):233–59. [54] Loseille A. Unstructured mesh generation and adaptation. Handb Numer Anal
[31] Zheng Y, Weatherill NP, Hassan O. Topology abstraction of surface models for 2017;18:263–302.
three-dimensional grid generation. Eng Comput 2001;17(1):28–38. [55] Rabinovich M, Poranne R, Panozzo D, Sorkine -Hornung O. Scalable locally
[32] Cripps RJ, Parwana S. A robust efficient tracing scheme for triangulating injective mappings. ACM Trans Graph 2017;36(2):16.
trimmed parametric surfaces. Comput Aided Des 2011;43(1):12–20. [56] CGAL, Computational Geometry Algorithms Library, http://www.cgal.org.
[33] Laug P. Some aspects of parametric surface meshing. Finite Elem Anal Des [57] Frey P, Borouchaki H. Surface mesh evaluation. In: 6th Intl. Meshing
2010;46(1):216–26. Roundtable. 1997, p. 363–74.
[34] Béchet E, Cuilliere J-C, Trochu F. Generation of a finite element mesh from [58] Cignoni P, Rocchini C, Scopigno R. Metro: measuring error on simplified
stereolithography (stl) files. Comput Aided Des 2002;34(1):1–17. surfaces. Comput Graph Forum 1998;17(2):167–74.
[35] Wang D, Hassan O, Morgan K, Weatherill N. Eqsm: an efficient high quality [59] Shewchuk JR. Mesh generation for domains with small angles. In: Proceedings
surface grid generation method based on remeshing. Comput Methods Appl of the sixteenth annual symposium on computational geometry. ACM; 2000,
Mech Engrg 2006;195(41):5621–33. p. 1–10.
[36] Marchandise E, Remacle J-F, Geuzaine C. Quality surface meshing using dis- [60] Inoue K, Itoh T, Yamada A, Furuhata T, Shimada K. Face clustering of
crete parametrizations. In: Proceedings of the 20th international meshing a large-scale cad model for surface mesh generation. Comput Aided Des
roundtable. Springer; 2011, p. 21–39. 2001;33(3):251–61.
[37] Marchandise E, Remacle J, Geuzaine C. Optimal parametrizations for surface [61] Si H. TetGen, a Delaunay-based quality tetrahedral mesh generator. ACM
remeshing. Eng Comput 2014;30(3):383–402.
Trans Math Software 2015;41(2):11:1–36.
[38] Aubry R, Karamete BK, Mestreau EL, Dey S. A three-dimensional para-
metric mesher with surface boundary-layer capability. J Comput Phys
2014;270:161–81.