Mesh generation from CAD — mesh_cad_surface

mesh_cad_surface triangulates the faces of a CAD model into one conforming surface mesh, with edges of about one target size h, and needs no CAD kernel. The model is the B-rep a Kratos CAD JSON read keeps in mesh.info["brep"] (the C++ BrepInfo), or a STEP or IGES read through OpenCASCADE. Every B-rep edge is sampled once and shared by the faces that use it, so the faces meet exactly and a closed solid comes out watertight. Each face is a constrained Delaunay triangulation of its parameter domain, lifted onto its surface. It is the first milestone of mesh generation in meshio++ (roadmap §6), with one global size. Since v17.14.0 it also meshes a face's embedded edges, a closed edge whose uses start at different points, and loops that cross only at a coarse size. Since v17.15.0 min_angle refines every face by Delaunay refinement until its triangles' smallest angles, measured on the surface, reach the bound, and smooth_iterations smooths the inside under a guard; with neither, the mesh is the unrefined one byte for byte. Since v17.16.0 order 2 or 3 writes curved triangle6 or triangle10 elements whose edge and face nodes lie on the exact surfaces and edges, each proved valid everywhere inside or kept straight-sided; order=1, the default, is the linear mesh byte for byte. Since v17.18.0 curvature_deviation and cells_across_gaps size the mesh by one gradation-limited field over the whole model (grading bounds its growth), pole_rings grades the triangles around a pole, and a triangle that lifts folded is refined away rather than written; with the sizing options and pole_rings at 0, a model that does not fold is meshed byte for byte as before. Since v17.19.0 sizes can be given per face (face_sizes), per edge (edge_sizes) or by a background mesh (background), joining the same field, and the report names the edges that hold a face finer than its given size; virtual_edges skips the slivers of edges a model splits its rims into; and metric="varying" meshes each face in its surface's own metric rather than one average map. With all three off, the mesh is as before byte for byte. Since v17.20.0 reparam_condition_number meshes a face whose parameters are unusable (a side collapsed onto a point, a pole, a corner its parameters stretch thirtyfold) in a flattening of its own coarse mesh, virtual_faces meshes a sliver face together with the neighbour it meets tangentially, and model_tolerances accepts the gaps a STEP or IGES model declares edge by edge; with all three off, the mesh is as before byte for byte.

import meshioplusplus as mp

model = mp.read("part.cad.json")  # or part.step / part.iges with OpenCASCADE
surface, report = mp.mesh_cad_surface(model, 0.1, return_report=True)
mp.write("part.vtu", surface)
report["num_failed"], report["num_open_edges"], report["min_angle_deg"]

meshioplusplus mesh-cad-surface part.cad.json part.vtu --size 0.1 --deflection 0.001 --json
meshioplusplus mesh-cad-surface part.cad.json part.vtu --size 0.1 --min-angle 25 --smooth 2
meshioplusplus mesh-cad-surface part.cad.json part.msh --size 0.2 --order 3
meshioplusplus mesh-cad-surface part.step part.vtu --size 2 --curvature-deviation 0.01 --cells-across-gaps 2 --min-angle 25
meshioplusplus mesh-cad-surface part.step part.vtu --size 2 --face-sizes 3=0.5,7=0.5 --background sizes.vtu --metric varying --min-angle 25
meshioplusplus mesh-cad-surface part.step part.vtu --size 2 --metric-background hessian.vtu --metric-array f:hessian --metric-error 0.001

# v17.19.0: face 3 at 0.5, edge 12 at 0.2, a background mesh's point data "size" elsewhere.
surface, report = mp.mesh_cad_surface(
    model, 2.0, min_angle=25, face_sizes={3: 0.5}, edge_sizes={12: 0.2},
    background=mp.read("sizes.vtu"), metric="varying", return_report=True,
)
[f["limited_by_edges"] for f in report["faces"]]  # the edges holding each face finer

# v17.20.0: faces whose parameters stretch more than fourfold meshed in a flattening,
# sliver faces meshed with their neighbours, the model's own tolerances accepted.
surface, report = mp.mesh_cad_surface(
    model, 2.0, min_angle=25, reparam_condition_number=4, virtual_faces=True,
    model_tolerances=True, return_report=True,
)
report["reparametrised"], report["merged_faces"], report["tolerance_uses"]

# Acute corners lopped (Pav and Walkington 2005), each with the bound they prove.
surface, report = mp.mesh_cad_surface(model, 2.0, min_angle=20, lop_corners=True,
                                      return_report=True)
[(c["stated_bound_deg"], c["guaranteed"]) for c in report["lopped_corners"]]

# An anisotropic metric from a field's Hessian on a coarse mesh: triangles stretch along
# the field's layers (roadmap 6.10.2).
coarse = mp.mesh_cad_surface(model, 2.0, edges=False)
coarse.point_data["f"] = field_at(coarse.points)  # any scalar field
coarse = mp.hessian(coarse, "f", location="point")
surface = mp.mesh_cad_surface(model, 2.0, metric_background=coarse,
                              metric_array="f:hessian", metric_error=1e-3)

The planar triangulator underneath is public too: triangulate_polygon_2d(outer, holes=(), target_size=0, min_angle=0, max_points=1000000) triangulates a polygon with holes in the plane z = 0, refined by Ruppert's algorithm to min_angle when it is set. Its coordinates must be zero or of magnitude within [2^-128, 2^128], where the exact predicates stay exact (a product inside them would otherwise overflow or underflow); anything else is refused by name.

How it works

1. Topology. Every edge's uses (BrepEdge.topology: a face id or name, a trim index and a relative direction) are resolved to the face trims they name. A trim no edge names becomes its own edge; with match_unreferenced_trims (the default), unreferenced trims whose 3-D images coincide (both ends and the middle within tolerance) are paired, so models without edge topology, such as the untrimmed patches of the Surfaces dialect, still conform; with model_tolerances they pair within the model's own tolerance (item 21). An untrimmed face is its patch's parameter rectangle. A loop's trims are walked head to tail by their ends in the face's parameters, whatever order they are listed in.
2. Corners from topology. A trim's end and the next trim's start are one corner, and so are the ends of every use of one edge, joined by a union-find rather than by comparing coordinates. Two faces whose trims disagree by the model's tolerance still share their corner points exactly; the disagreement is measured and reported as max_edge_gap. A closed edge whose uses start at different points (a cap's circle listed from another angle than the cylinder's rim) is anchored on the use that cannot move, one that runs between other trims of its loop, and every use that is a loop of its own and closed in its face's parameters is rotated to start at the point of it closest to the anchor's start (a scan by arc length refined by Newton), so the uses meet at one corner again.
3. Edges, sampled once. Each edge is sampled on its first use, by arc length (Gauss–Legendre on the trim's knot spans), into ceil(L / h) equal pieces between the trim's corners (a polyline trim keeps its vertices), at least min_points_per_edge points inside. With max_chordal_deviation, a piece is halved while its sagitta, the distance of the edge from its chord, exceeds the cap (OpenCASCADE's "linear deflection"). Every other use, the neighbouring face or the same face across a seam, projects the same 3-D points onto its own trim, and the points' ids are shared: this is what makes the faces conform; a rotated use runs past its end and on from its start. An edge shorter than the tolerance, a pole, keeps its two ends, which are one corner. With a size field (item 13) the pieces are equal in ∫ ds / h(s) instead, the integral by adaptive Simpson along the arc length, so an edge is finer where the field is.
4. Crossings at a coarse size. The sampled loops of a face can cross in its parameters where the curves themselves do not: two close loops, a thin wall, a hole nearer its rim than h. Each face's loops are searched for segments that meet anywhere but at a shared sample (a crossing, a touch, an overlap) with exact orient2d tests over a bucket grid of the plane the triangulation is given. Every piece involved is halved at its arc-length midpoint in every face that uses its edge, so the neighbours take the same new points, and the faces those edges touch are searched again until nothing crosses. The requests of all faces are merged in (edge, piece) order, so the result does not depend on the thread count. A crossing whose pieces are already within the tolerance is where the curves themselves meet: that face stays self_intersecting, as does one still crossing after 64 rounds, and one whose splits would pass max_points is budget_exceeded.
5. Faces, in their parameters. Each face's loops become closed polylines in its parameters, mapped by one affine transformation (the Cholesky factor of the face's mean first fundamental form) so that a Delaunay triangulation there follows the surface's own metric on average rather than raw parameters, whose scales can differ by orders of magnitude (the Scordelis–Lo roof's differ fortyfold). The constrained Delaunay triangulation inserts the boundary points with Lawson's flips, recovers every boundary segment by re-triangulating the cavity it crosses (Anglada 1997), and keeps the triangles inside by the even-odd rule, so holes need no orientation. Interior points come from an equilateral lattice whose step is at most h on the surface, kept only farther than h / 2 from the boundary points and the points already kept and 0.4 h from the boundary segments. The points are lifted onto the surface, triangles a pole collapses are dropped, and the winding follows dS/du × dS/dv, reversed for a face whose normal is swapped.
6. Embedded edges. A face's embedded edges (a coupling or load curve inside it) are edges like any other: named by a BrepEdge or paired by their images, sampled once and sharing their points with every face that uses the same edge (a fin standing on a plate along one). An embedded edge's end within the tolerance of a corner of its face's loops is that corner, and a closed one is one corner. After the loops are classified, its samples are inserted and its segments recovered by the same cavity re-triangulation as the loops (Anglada 1997); since the even-odd rule has already counted the loops alone, an embedded segment never changes which triangles are inside. The lattice keeps its spacing from it as from the boundary.
7. Refinement (min_angle). Each face is refined in its plane by Ruppert's algorithm with Chew's circumcentres (Ruppert 1995; Chew 1989): a triangle whose circumradius-to-shortest-edge ratio exceeds B = 1 / (2 sin(min_angle)) gets a new point, unless that point lies strictly inside the diametral circle of a boundary piece (the exact diametral predicate) or beyond one, which is split instead; a piece a vertex encroaches is split first. Bad triangles are taken worst first from a priority queue whose ties are broken by the triangle's sorted vertex ids, never by a slot or a counter (Shewchuk 2002), so a face refines the same way on every run. The triangle is judged on the surface: its ratio is that of the lifted triangle, and the point is the lifted triangle's own circumcentre carried to the plane by the triangle's affine map, where the plane's circumcentre would not round a triangle the face's metric stretches (Chew 1993; Boissonnat and Oudot 2005). Where the plane's Delaunay connectivity is not the surface's, an edge whose flip raises its two triangles' smallest angle on the surface is flipped (when its quadrilateral is convex, by exact orient2d) and refinement runs again. A point closer to its cavity's vertices than half Ruppert's insertion radius (B l for a circumcentre, B l / sqrt 2 for a split, l the triangle's shortest edge), which a plane that is the surface never produces, is declined and its triangle reported, so points do not pile up where one affine map cannot follow the surface. Where the surface stretches the plane more than twofold (the condition number of J L^-1 above 2, J the patch's Jacobian and L the face's map), a point placed for the surface can leave smaller triangles as bad on the surface as the one it replaced, and refinement need not stop: there no point lands closer than h / 4 to its cavity and no triangle shorter than h / 4 is refined, so finitely many points fit, and a triangle left under the bound is warned. With metric="varying" (item 18) that region is refined like any other.
8. Splits shared across faces. A boundary piece belongs to an edge, shared with the neighbouring face, so a face asks for its splits rather than making them. Every piece asked for is split by arc length in every use of its edge, the crossings are searched again, and every face is meshed again, until no face asks; the requests are merged in (edge, piece) order, so the result does not depend on the thread count. Splitting by arc length rather than by parameter keeps the halves close to straight, which Ruppert's encroachment argument assumes; the chordal sagitta of the cylinder's rims stays under 3% of the piece at min_angle = 30. After 16 rounds (the fixtures need at most 7), or when the splits would pass max_points, refinement stops with a warning and the faces report what is left. Each face refines with at most its share of max_points and twice its unrefined points summed over the rounds (1024 at least; in a sized run, sixteen times the points its size field asks for if that is more, and four times them in any one meshing, since a face sized down to a narrow layer re-spends its points each round and needs ten times them in all; a face that converges adds at most 1.46 times its points over all rounds where it has 256 or more, 1.07 times on the fixtures, and at most 521 points on any face, 3.07 times on the coarse capped cylinder's 136 at 30 degrees): a face that spends them, which on a surface past 20.7 degrees can happen without end, stops the rounds and is meshed unrefined on its current boundary, for its size as for its angles, with a warning, rather than every round spending them again.
9. Small input angles and poles. Ruppert's algorithm need not terminate where two boundary pieces meet at less than 60 degrees (the exact lens_angle predicate decides). There a piece at the corner splits on concentric shells, a power-of-two distance from the corner, so the pieces of the cluster stay equal; a skinny triangle that holds the angle, or whose shortest edge spans it, is left as the angle forces it; and a piece of the cluster is not split when its halves would be shorter than the skinny triangle's shortest edge (Shewchuk 2002; Miller, Pav and Walkington 2003). Each such corner is reported with its angle in the face's plane and the smallest angle, on the surface, of the triangles it forces. A pole (a degenerate trim, the apex of a fan) is reported the same way, with its reach: the triangles within h of it, or farther where the face's metric departs more than twofold from its plane's, are left to it. With pole_rings = m (at most 8) the pole is graded instead (item 14).
10. Guarded smoothing (smooth_iterations). Each pass moves every interior point toward the centroid of its neighbours (laplacian) or its star's area-weighted circumcentre (odt, the optimal Delaunay move of Chen and Xu 2004), in the plane, lifted onto the surface. A move is kept only when no triangle of the star inverts (exact orient2d) and the star's smallest angle on the surface does not drop, since unguarded smoothing can invert or worsen triangles (Field 1988); Lawson's flips follow each pass, judged on the surface. Boundary points never move, so shared edges stay shared. A face that asks for boundary splits is meshed again in the next round, so it is smoothed only in the meshing that is kept; when the rounds stop with its splits not made, that meshing is made again, smoothed, from the same boundary and refinement budget.
11. Curved elements (order). The linear mesh is made first, refinement and smoothing included, and its triangles keep their corners; order 2 adds a node on each edge, order 3 two on each edge and one inside. A node on an edge inside a face is the surface at the parameter fraction (1/2, or 1/3 and 2/3) of the edge's two corners, computed once per face and shared by its two triangles; a node on a B-rep edge is the master trace at the arc-length fraction of its piece, computed once per piece and shared by every face that uses the edge, so neighbours conform as the linear mesh does. At a pole, whose parameters a degenerate trim spans, the pole's corner is taken at the point of that segment toward the edge's other end. The face node is the surface at the parameter centroid. Each element is then checked against its face's normal S_u × S_v at that centroid: its area element J = (x_ξ × x_η) · N, a polynomial of degree 2 or 4, is expanded in Bernstein form; all coefficients positive proves J > 0 everywhere, a value not positive proves it folds, and otherwise the triangle is split in four and each part judged, to depth 8, undecided counting as folded (Johnen, Remacle and Geuzaine 2013). A positive J at the nodes is necessary, not sufficient: a curved element can fold between them. An element that folds is made straight-sided, its edge nodes on its chords and its face node at its corners' centroid; since the nodes are shared its neighbours are checked again, to a fixed point, and every one made straight is counted in num_curved_straightened and warned. A straightened boundary node leaves the curve for every face that uses it alike, so the mesh stays conforming.
12. Exact predicates. Every orientation and in-circle decision is Shewchuk's adaptive exact predicate, re-implemented from the paper, as are refinement's encroachment (diametral) and small-angle (lens_angle) tests, built the same way. No coordinate is compared to another with a tolerance in the triangulation; the quality measures (ratios, angles) are the only floating-point judgements, and they choose where to refine, never whether the triangulation is valid.
13. Size fields (curvature_deviation, cells_across_gaps, grading). A first, uniform meshing at h gives the samples: every point of every face, with the surface's unit normal and its largest absolute principal curvature k there, from the first and second fundamental forms (I⁻¹ II, the NURBS second derivatives by the quotient rule of Piegl and Tiller's A4.4). Each asks for the curvature size 2 sqrt(2 d / k - d²), the chord of a circle of radius 1 / k whose sagitta is d = curvature_deviation, and the gap size t / n, where t is the distance a ray along the normal travels, inward to the far wall of a thin part or outward across a gap, to a face that is neither the point's own nor a neighbour (on the first mesh's triangles, through a bounding-volume hierarchy) and n = cells_across_gaps; both are divided by a slack of 1.5, the most an edge may exceed the local size. The field is their gradation-limited envelope h(x) = min(h, min_i h_i + (g - 1) max(0, |x - x_i| - r)), Persson's (2006) continuous form of gradation, exact for point sources, with g = grading (1.5 when 0) and r = 3h/4 the reach of a sample's own size; it is never above h nor below h / 64, and looked up through a uniform grid. One function of space for the whole model, so the two faces of an edge and the edge itself agree on it by construction. The second meshing, the result, samples the edges by it (item 3), spaces the lattice's admission radii by h(x), and refines (item 7, with no angle sought when min_angle is 0) every triangle whose longest edge on the surface exceeds 1.5 h at its centroid, worst first; a point inserted for size alone keeps a third of the local size from its cavity's vertices, and asks a boundary piece to split only when the piece is itself half again too long. With both sources at 0 the field is not built; when no sample asks for less than h (a flat model with no thin part) the mesh is the uniform one byte for byte.
14. Graded poles (pole_rings). At a pole the parameters collapse onto a segment and the lattice leaves a fan of thin triangles. With pole_rings = m, each pole gets m rings of points at 3-D distances h, 2h, … from it: along the segment's inward perpendicular, a bisection on |S(u, v) - P| = r finds each ring's parameters, and the ring gets round(L / h) points equally spaced by its own traced length L, skipping any near an existing point. The lattice is kept off the band they cover, out to the last ring plus h / 2; the pole's exemption shrinks to its first ring, so refinement and the surface flips work everywhere else. On the unit sphere at min_angle = 25 two rings lift the smallest angle from 10.16 to 26.14 degrees at h = 0.2 and from 6.33 to 25.12 at h = 0.1, every pole with valence 6.
15. Folds. On a face smaller than h, or strongly curved within one triangle (a blade, a fillet, a band swept past half a turn with one boundary piece per arc), a triangle valid in the plane can lift onto the surface folded, its normal against the surface's. A triangle is folded when its normal turns more than 120 degrees from S_u × S_v at its parameter centroid, or at all three of its corners (the margin keeps a sliver across a narrow strip from counting on rounding). A face that folds is refined even when no angle or size is sought: a folded triangle is bad first of all, split at its centroid in the plane, and one whose centroid encroaches a boundary piece asks that piece to split, in every face that uses its edge, as refinement's splits are (item 8). A fold left after the last round, past the face's points or out of the splits' reach, fails the face as degenerate rather than being written. None of the fixtures folds; on the real STEP corpus several hundred faces did.
16. Given sizes (face_sizes, edge_sizes, background). Sizes given by the caller join the field of item 13 as samples on the first mesh: a face's size at every first-mesh point of that face, an edge's at every point of its line cells, and, when curvature_deviation is on, the edge's own bend there, 2 sqrt(2 d / k - d²) from the curvature k of the circle through three consecutive points of it. A background mesh (a surface, or a volume, whose skin is taken) is evaluated directly wherever the field is: the closest point on it (through a bounding-volume hierarchy), its one-component point data background_array (default "size") interpolated barycentrically there. Sizes only refine: the ceiling stays target_size, so a region is coarsened by raising target_size and giving finer sizes elsewhere. The gradation of item 13 still holds, so a face given a size grades its neighbours away from it, and a face beside a finer one is held finer along the edges they share. After meshing, a face given a size whose shared edges carry a median field size smaller than its size over the slack 1.5 lists those edges in its limited_by_edges, and each is warned ("face f: size limited by edge e"). An index out of range, a size that is not finite and positive, and a background without its array or with one of more than one component are refused before meshing.
17. Virtual edges (virtual_edges). A model often splits a rim into a long edge and a sliver beside it, and the sliver's corner forces a point there, so a triangle on each side has an edge of the sliver's length. With virtual_edges, a corner where exactly two edges meet is suppressed when both are used by the same faces, the rim turns there by under 30 degrees, one of the two edges is shorter than h / 2, and the chord across the corner, between its two neighbouring samples, stays within 1.5 h (the local size, with a field) and a tenth of its length of the corner. A corner beside one already suppressed across a single piece is kept, so no chord skips two. A suppressed corner is in no face's loop and no line cell: the two line cells beside it become one, straight across it, and the faces stay conforming since every face of both edges skips it alike. The chains of edges joined through suppressed corners are reported (virtual_edges, each a list of BrepEdge indices in walking order); the model is not changed, and the corner stays in the points as an orphan, which clean(remove_orphans=True) drops (as subdivide's orphans are left). Curved elements follow each edge's own curve, so virtual_edges is refused with order 2 or 3. It is off by default; on the real STEP corpus it forms 67 chains, meshes one more face and closes the Ventilator's 73 open edges, but leaves 24 more faces under the angle bound (Sheffer et al. 2000 build the same virtual topology, for faces as well).
18. A varying metric (metric="varying"). One affine map per face follows the surface's metric only on average. With metric="varying" the face is still triangulated in its plane, but every Lawson flip after an insertion is decided by the in-circumellipse test of Borouchaki et al. (1997) in the surface's own metric, the first fundamental form in the plane's coordinates M = (J L⁻¹)ᵀ (J L⁻¹), evaluated at the shared edge's midpoint: the four points are mapped by its Cholesky factor and the exact incircle is taken on the mapped coordinates. A flip still needs a strictly convex quadrilateral (exact orient2d), and only the new vertex's star is tested, each flip adding an edge to it that is never flipped away, so the flips end. Refinement measures encroachment in the same metric (the diametral circle at the boundary piece's midpoint), so a piece short on the surface but long in the plane is not split by the points its own size puts beside it; the stretched-region exemption of item 7 is off. A face still folded (item 15) after three meshings in a row, or when no split is left, is meshed again in its mean metric and warned, since the metric at an edge's midpoint is the surface's to first order only. On a cone frustum whose radius varies tenfold, a torus band whose parallels vary eightfold and the stretched bilinear patch, every angle meets 25 degrees with fewer points than the mean metric, which leaves 39 to 153 triangles under it; unrefined, the smallest angle rises from 9.8 to 29.8 degrees on the cone and from 4.8 to 25.8 on the torus, and the sphere's poles meet 25 degrees without rings.
19. Discrete reparametrisation (reparam_condition_number). A face's own parameters can be unusable for the plane of item 5: a patch whose corner control points coincide collapses a side onto a point, a hemisphere's fan meets at its pole, a trim ends short of a collapsed side where the parameters stretch thirtyfold. With reparam_condition_number = K (0 is off, otherwise above 1), each face's unrefined triangulation, the lattice with its corners and boundary samples, is kept as its coarse mesh (consecutive samples a pole collapses merged into one vertex, a seam's two sides kept apart), and its condition number measured: the largest singular value of J L⁻¹ at any coarse triangle's centroid over the smallest at any. A face above K is flattened. The coarse mesh is first mapped into the plane by least-squares conformal maps (Lévy et al. 2002), kept when no triangle flips, the orientation holds and the loops' images stay simple (the exact crossing search of item 4); otherwise by mean-value coordinates, the outer loop on a circle by chord length and each hole closed by a virtual fan (Floater 2003), only when every boundary corner is 60 degrees or more, since the circle keeps the samples' spacing but flattens their corners. A band closed without a seam (a cone's frustum, a cylinder trimmed by two full circles) has a coarse mesh with two loops and Euler characteristic 0, an annulus rather than a disc, and a fan would crush its second rim into the middle of the circle: it is laid between two concentric circles instead (method 3, "annulus"), the longer loop on the outer one, the other on a circle exp(-2π w / L) times as large, the conformal modulus of a cylinder of circumference L (the loops' mean length) and height w (the area over L), each by arc length and the inner one starting at its vertex nearest the outer's first so the map does not twist; the interior is the same mean-value solve, and the same exact checks decide. The planar triangulation then meshes the inner loop as a hole. Inside a coarse triangle the parameters are the barycentric blend of its corners', and in a triangle with a corner at a pole they are collapsed coordinates (Duffy 1982), the pole's corner taken at the point of its segment toward the opposite side, so the parameters along an edge out of the pole agree from both sides. The face is then meshed in the flattening as it would be in its parameters, with the flattening's mean metric in place of L, every point evaluated on its surface at the parameters its position maps to, located through a bucket grid; the boundary samples keep their exact parameters, so the face conforms with its neighbours unchanged. The flattening is taken only when its own condition number is below the face's; otherwise, or when it cannot be built, the face is meshed in its parameters and the reason is reported and warned. The decision is made at the face's first meshing and kept; the coarse mesh and its flattening are rebuilt every meshing, since splits move the boundary. In a flattening a pole is a point of the boundary, not a fan's apex, so pole_rings and the pole's exemption skip the face, and it is meshed in its mean metric whatever metric says. It is Remacle et al.'s (2010) route for faces whose parameters are unusable, with the coarse mesh as the discrete surface; order 1 only.
20. Virtual faces (virtual_faces). A sliver face, a strip between two long faces, forces triangles as narrow as it is, and its neighbours are refined against it. With virtual_faces, after the first meshing a face with no interior point whose width, twice its area over its perimeter, is under h (the local size, with a field) is a sliver, and is merged into a host: of its neighbours across manifold edges, the one sharing the most edge length (ties to the lower index) that is not itself a sliver or a host, meets it tangentially (the normals turn under 30 degrees at every shared sample, the virtual edges' threshold) and shares one run of edges with it in both loops; the next neighbour is tried when one is not. A sliver left over, whose neighbours are all slivers or hosts (a fillet split into strips, a strip between two merged ones), then joins a group it borders, in rounds until none joins: the group of the host it shares the most edge length with, the same checks made against the group's merged loop, so a chain of strips grows from the face it starts on. A group holds at most eight faces. The coarse meshes (item 19) are glued along the shared samples, in the order the slivers joined, into one disc, flattened as in item 19 with each triangle's own surface, and meshed as one face whose loop leaves the shared edges out, as virtual_edges leaves out a suppressed corner. Each triangle belongs to the face whose coarse triangle holds its centroid, so cad:face and face:<k> keep each model face's triangles and each face's record counts its own. The shared edges' line cells are dropped and their samples stay as orphan points, which clean(remove_orphans=True) drops. A group still folded after three meshings in a row, or when no split is left, is split back into its faces, each meshed on its own from then on, so a merge does not fail its own faces; a neighbour can still fold on the boundary the pair leaves it (two faces of the Ventilator model on the corpus). Each sliver found is reported (merged_faces: sliver, host (its group's), width, the BrepEdge indices it shares with the group, and why not when it is not merged). Order 1 only.
21. Model tolerances (model_tolerances). A STEP or IGES read through OpenCASCADE writes each edge's own tolerance and its two vertices' into the edge's extra JSON ({"tolerance": t, "end_tolerances": [t0, t1]}, the approximation tolerance added where an edge was approximated), which CAD JSON carries through a write and a read; any reader may write the same keys. A BrepVertex declares its own ({"tolerance": t}, which an OpenCASCADE read writes for every vertex with the trims it starts and its parameters on each face), and a Kratos CAD JSON file its model's at top level ("tolerances": {"model_tolerance": t}, kept in BrepInfo::mExtra). With model_tolerances, the model's tolerance is the floor of every declared one, an edge's tolerance is the larger of that and its own, and a corner's the larger of that, its edges' ends' and its vertices' (each vertex at the end of the trims it names nearer its parameters, the start without them). The trims no edge names (item 1) pair within the model's tolerance. They accept only what the global tolerance refuses: a use of a closed edge that cannot be rotated onto its anchor's start (item 2) is taken to start there when it lies within the edge's tolerance, instead of failing its face; a corner's points may lie apart, before its face is an open_loop, and an edge's uses, before it is projection_failed, by h / 4 or ten times the tolerance, whichever is larger. A use that can be rotated is rotated as without the option, since keeping a start a declared gap away was measured to fold triangles. Whether an edge is degenerate and the sampling stay on the global tolerance, so a declared tolerance never makes a real edge a pole. Each place a declared tolerance accepted what the global one refused is reported (tolerance_uses: edge, declared tolerance, gap, and closed_start, corner, projection or pairing, a pair having no edge, -1), and a declared tolerance that is not finite and non-negative is ignored with a warning naming its edge, vertex or the model.
22. Corner lopping (lop_corners). Concentric shells keep refinement finite at a corner under 60 degrees, but only measure what the corner forces. Corner lopping (S. E. Pav and N. J. Walkington, Delaunay refinement by corner lopping, IMR 2005, Algorithm 2) replaces them and proves a bound. The input points are the ends of the face's edges and any kink inside one (a pole's aside, which keeps its own handling), and every one at an acute corner (an inside wedge under 90 degrees in the face's plane) gets a ball of radius d = (sqrt 2 - 1) lfs0, where lfs0 is its distance in 3-D to the face's nearest other corner. Each of its edges is split where it leaves the ball. Then, in priority: a vertex inside a ball shrinks it to half (SPLITBALL(p, 1/2)) and splits the edges there again; an encroached piece is split, a piece of a ball's own by shrinking that ball; and a bad triangle is refined unless its point (the surface circumcentre of item 8) lies in a ball, where it is left, to be tried again when the ball shrinks. Concentric shells and their terminator are off. Balls are measured in 3-D, carried between rounds and faces by corner position (the smallest radius any face or split left), and asked of the edges like every other split: a piece is cut where its distance from the corner, bisected along the trace, is the radius, which then becomes that split's own distance. The paper assumes no curve turns more than pi / 3, so a piece that turns more (its tangents summed at nine points) is halved first. Corollary 3 proves termination for min_angle under arcsin(1 / (2 sqrt 2 mu)) and arcsin(1 / (2 zeta eta)), with eta = (1 + tan(dt/2)) / (1 - tan(dt/2)), mu = 1 / (1 - tan(dt/2)) and zeta = sqrt 2 sin(dt/2) / (sqrt(1 + sin dt) - 1) for the most any piece turns, dt. That is 20.7 degrees for straight edges, falling as the edges turn. Lemma 7 says every triangle left under it has its vertices within 2 d of a corner. Each corner is reported with its angle, lfs0, its final radius, the face's turning, the bound (on a planar face only, since the theorem is the plane's), the triangles left in its ball and their smallest angle, whether Lemma 7 held (checked at the end, not assumed), and guaranteed, or a reason: a curved face, min_angle not under the bound, refinement that did not finish, or a triangle past 2 d. A corner still under 60 degrees is listed among the face's corners as before. The lattice of item 7 is no part of the paper's input, and its points inside a ball shrink it to about the local spacing: on a planar lens of two circular trims meeting at 5 degrees, refined to 20 degrees, both corners end at radius 0.026 at h = 0.05 (272 points) and 0.013 at h = 0.02 (942 points), each keeping only its own corner triangle, under a stated bound of 20.68 and 20.69 degrees. Off by default; the mesh is then the shells' byte for byte.
23. A metric source (metric_background, metric_array, metric_error). A size field is isotropic: a layer that varies across but not along needs triangles fine across it only. With a metric_background mesh (a surface, or a volume whose skin is taken) whose 9-component point data metric_array is a field's Hessian (what hessian writes as <array>:hessian) and a target interpolation error metric_error, every background point carries the metric M = R |Λ| Rᵀ / (8 ε) of its symmetrised Hessian R Λ Rᵀ, under which linear interpolation errs by about ε on a unit edge (A. Loseille and F. Alauzet 2011, the bound h ≤ sqrt(8 ε / |λ|) along each eigenvector). At a point of a face the metric is read at the nearest point of the background's surface, interpolated barycentrically (a convex combination stays positive semi-definite), projected on that surface's tangent plane, and its sizes clamped between a 64th of target_size and the local size of item 13's field: like every size source it only refines. It composes with the surface as (J L⁻¹)ᵀ M (J L⁻¹), which item 18's flip and encroachment tests take unchanged, so metric_background meshes in the varying route whatever metric says; and refinement judges every triangle in it: its lengths against the unit by the slack 1.5, its circumradius-to-shortest-edge ratio and smallest angle on its corners mapped by M^½ at its centroid, its circumcentre in that space carried to the plane by barycentrics, the distances that guard an insertion, the boundary pieces it splits, and the point budget, area · sqrt(det M) over an equilateral unit triangle. Lawson's flips after refinement and guarded smoothing raise the smallest angle in the metric too. The lattice and the edges are sampled at the isotropic size, never finer than the metric's coarsest direction, and refinement adds the fine one. On the 6.11 layer (tanh((z - 0.3) / 0.05) on the unit sphere, its Hessian from a background at h = 0.1, ε = 8e-3), the mesh has 1,008 points with an L2 interpolation error of 1.53e-2 and a maximum of 5.5e-2, its triangles in the layer stretched 5.4 to one along it, against the uniform mesh at h = 0.05 (11,281 points, 1.58e-2, 9.7e-2) and the isotropic size from the same Hessian (3,848 points, 1.39e-2, 0.14). The reported angles stay the surface's, so a stretched mesh reports small ones. A missing or wrong-sized array, an error that is not finite and positive, and order 2 or 3 are refused by name; without a metric background the mesh is byte-identical.

Notebook 42 renders these meshes and the Scordelis–Lo roof, checks their convergence, and measures refinement, smoothing, curved elements and the interpolation error.

Options

Python	C++ (MeshCadSurfaceOptions)	Meaning
target_size	mTargetSize	The target edge length h; required, > 0
max_chordal_deviation=0	mMaxChordalDeviation	The sagitta cap of the boundary segments; 0 is off
min_points_per_edge=0	mMinPointsPerEdge	At least this many points inside every edge
tolerance=0	mGeometricTolerance	Edges shorter than this are degenerate, and unreferenced trims closer than this are paired; 0 is 1e-7 times the diagonal of the control nets' box
record_face_ids=True	mRecordFaceIds	cad:face cell data and a face:<k> cell region per face (failed faces keep an empty one)
edges=True	mEdges	A line block of the discretised edges with cad:edge, the BrepEdge index (-1 for an edge no BrepEdge names)
match_unreferenced_trims=True	mMatchUnreferencedTrims	Pair unreferenced trims by their images
strict=False	mStrict	Raise (RuntimeError; C++ std::runtime_error) when any face fails, instead of reporting it
max_points=50000000	mMaxPoints	The most points in the result; a face or edge that would exceed it fails, and refinement shares what the boundary leaves equally among the faces, each at most twice its unrefined points
min_angle=0	mMinAngleDeg	Refine until every triangle's smallest angle on the surface is at least this (degrees, under 60); 0 keeps the unrefined mesh. Guaranteed up to 20.7 degrees away from corners under 60 degrees and poles; above, a target whose shortfall is reported (it usually reaches about 30 here)
smooth_iterations=0	mSmoothIterations	Passes of guarded smoothing
smooth_method="laplacian"	mSmoothMethod	laplacian or odt
order=1	mOrder	1 linear triangle and line; 2 curved triangle6 and line3; 3 curved triangle10 and line4
curvature_deviation=0	mCurvatureDeviation	Size curved regions so a triangle sits at most this far from the surface (its chordal error); 0 is off
cells_across_gaps=0	mCellsAcrossGaps	At least this many elements across a thin wall or a gap; 0 is off
grading=0	mGrading	How fast the size may grow, at most grading - 1 per unit length (a ratio above 1); 0 is 1.5, used only when a size source is on
pole_rings=0	mPoleRings	Rings of points about h apart around each pole (0 to 8); 0 is the pole's fan
metric="mean"	mMetric	mean, one affine map per face, or varying, the surface's own metric wherever a flip or an encroachment is decided (item 18)
face_sizes=None	MeshCadSizeSources::mFaceSizes	Sizes given per face: a dict {face index: size} or (index, size) pairs; they only refine (item 16)
edge_sizes=None	MeshCadSizeSources::mEdgeSizes	Sizes given per edge (BrepEdge index), the same way
background=None	MeshCadSizeSources::pBackground	A mesh whose point data is the size at its points, interpolated at the closest point
background_array="size"	MeshCadSizeSources::mBackgroundArray	The background's one-component point data array
virtual_edges=False	MeshCadSizeSources::mVirtualEdges	Suppress the corners between slivers of edges and their neighbours (item 17); order 1 only
reparam_condition_number=0	mReparamConditionNumber	Mesh a face whose condition number exceeds this in a flattening of its coarse mesh (item 19); 0 is off, otherwise above 1; order 1 only
virtual_faces=False	MeshCadSurfaceExtras::mVirtualFaces	Mesh a sliver face with the neighbour it meets tangentially (item 20); order 1 only
model_tolerances=False	MeshCadSurfaceExtras::mModelTolerances	Accept the gaps each edge's own declared tolerance allows (item 21)
lop_corners=False	MeshCadSurfaceExtras::mLopCorners	Lop acute corners in place of concentric shells and state the bound Pav and Walkington prove (item 22); only with refinement
metric_background=None	MeshCadSizeSources::pMetricBackground	A mesh whose point data is a field's Hessian: an anisotropic metric from it sizes the triangles (item 23); order 1 only
metric_array=""	MeshCadSizeSources::mMetricArray	The background's 9-component Hessian array
metric_error=0	MeshCadSizeSources::mMetricError	The target interpolation error, positive with a metric background
—	mInsertion	Laid out for off-centres (6.2.9); refused unless left at its default. Off-centres were measured and parked: on lattice-seeded faces they change nothing (roadmap 6.2.9)

In C++ the sources travel beside the options, mesh_cad_surface(const BrepInfo&, const MeshCadSurfaceOptions&, const MeshCadSizeSources&, MeshCadSizeReport* = nullptr), since MeshCadSurfaceOptions and MeshCadSurfaceResult keep their ABI 34 layouts (MeshCadSizeSources itself gained the metric source in ABI 35); the two-argument overload is the new one with empty sources, byte for byte. MeshCadSizeReport holds mLimitedByEdges per face and the mVirtualEdges chains. Since v17.20.0 mesh_cad_surface(const BrepInfo&, const MeshCadSurfaceOptions&, const MeshCadSizeSources&, const MeshCadSurfaceExtras&, MeshCadSizeReport* = nullptr, MeshCadSurfaceExtrasReport* = nullptr) takes the inputs that have no room left in those structs; the four-argument overload is it with empty extras, byte for byte. MeshCadSurfaceExtrasReport holds mReparametrised (MeshCadReparamFace: the face, its condition number in its parameters and in the flattening, whether it was taken, the method and the reason), mMergedFaces (MeshCadMergedFace) and mToleranceUses (MeshCadToleranceUse); mLoppedCorners is laid out and stays empty.

Output and report

The mesh has the points and one triangle block (and the line block of the edges, embedded edges included). It writes to VTU, Gmsh, STL and every other surface format unchanged. Gmsh gets each face as a named physical surface face:<k>: 2.2 tags every triangle, and 4.1, whose groups live on entities, splits the triangle block into one block per face (since v17.17.0), so a box opens in Gmsh with six physical surfaces; the faces also travel as the cad:face data. With order 2 the blocks are triangle6 and line3, which VTU, VTK, Gmsh, Medit, XDMF, Exodus, MED, CGNS, Abaqus, Nastran, UNV, MDPA, Elmer, FEBio and XPLT carry; with order 3 they are triangle10 and line4, which Gmsh, Medit (TrianglesP3), Elmer and XPLT carry, and MDPA under their meshio names, which Kratos itself does not define (VTU has no triangle10). The points are the linear mesh's, with their ids, followed by the edge nodes and then each face's own nodes; the regions and cad:face / cad:edge index the elements as before.

The report (return_report=True; MeshCadSurfaceResult in C++) has the counts (num_points, num_triangles, num_edge_segments, num_faces, num_failed, num_edges), the conformity of the result (num_open_edges, triangle edges used once: zero on a closed solid, the model's own boundary on an open one; num_non_manifold_edges; num_misoriented_edges, used twice the same way round), num_degenerate_dropped, the extreme angles and edge lengths, tolerance, max_edge_gap, one record per face (status, reason, counts, extremes) and the warnings. Refinement adds num_refinement_points (inside the faces and on their edges), num_faces_below_bound, num_smoothed and angle_histogram, every angle of every triangle on the surface in 18 bins of 10 degrees; each face record adds num_refinement_points, num_smoothed, below_bound (a triangle under min_angle), bound_deg (the smallest angle of its triangles that no corner forces and no stretched region leaves: at least min_angle where refinement reached its bound), min_input_angle_deg and corners, each corner under 60 degrees or pole with its point, parameters, position, angle_deg, min_forced_deg, reach and pole. num_curved_straightened counts the curved elements kept straight-sided; the angles, edge lengths and histogram are those of the straight triangles through the elements' corners. Since v17.19.0 each face record adds limited_by_edges, the edges that hold the face finer than its given size (item 16, empty without given sizes), and the report adds virtual_edges, the chains of item 17. Since v17.20.0 the report adds reparametrised (one entry per face over reparam_condition_number: face, condition_number, condition_number_disc, taken, method (conformal, mean_value or none) and reason), merged_faces (sliver, host, width, edges, merged, reason), tolerance_uses (edge, declared, gap, where) and lopped_corners (face, point, angle_deg, lfs0, lop_radius, turning_deg, stated_bound_deg, num_forced, min_forced_deg, within_two_radii, guaranteed, reason; item 22).

A face that cannot be meshed is left out and reported, never meshed wrong. Its status names why: unsupported (an invalid surface or trim, a trim that leaves its surface's knot range by more than a thousandth of it, where the patch would be evaluated clamped and points would lie off the face, or parameters that, mapped to the face's plane, fall outside the range [2^-128, 2^128] the exact predicates hold in), open_loop (its loops do not close, or a corner the topology joins it to lies elsewhere), self_intersecting (its boundary meets itself in its parameters where the curves themselves meet, such as a hole tangent to the rim, or still crosses after 64 rounds of splits), degenerate, projection_failed (an edge's uses lie farther apart than h / 4, or a closed edge's uses start at different points and none of them can be started elsewhere), budget_exceeded, empty. The face's edges stay in the line block, and num_open_edges shows the hole.

An embedded edge that cannot be inserted is left out with a warning naming the face and the trim, and the face is meshed without it: one that leaves the face's domain or ends on its boundary away from a corner, and one that crosses or touches the boundary or another embedded edge at this size. Those are not split, since refining the loops around a curve that really crosses them would only crowd boundaries the neighbouring faces share. Its line cells stay in the line block.

Volume meshes

A closed result is the skin a tetrahedral mesher needs, and octree_mesh fills it natively: octree_mesh(mesh_cad_surface(model, h, min_angle=25), skin_grading="uniform", surface_size=h) gives conforming tetrahedra of the inside, the line block of the edges skipped. The tetrahedra fit the skin by isosurface stuffing (cut points on the skin, lattice vertices warped onto them, Labelle and Shewchuk 2007), not by keeping its triangles: the boundary is the skin's within a measured distance, and in exchange every dihedral angle is bounded away from 0 and 180 degrees by construction, 10.7 to 164.8 degrees where the report's guarantee_applies holds (uniform skin grading, no warp refused, no node moved onto a feature). On the unit sphere at h = 0.2, 0.1 and 0.05 the enclosed volume is the skin's within 0.30%, 0.075% and 0.019%, second order, so the sphere's own 4/3 π is met to the skin's chordal error; the boundary is within 0.036 h, 0.028 h and 0.03 h of the skin (hausdorff_to_skin); and the dihedral angles lie between 25.6 and 132.4 degrees, then 21.6 and 144.2. A box is filled exactly, by BCC tetrahedra only (35.3 to 125.3 degrees). Uniform grading matters: the default graded skin cuts transition shapes along the skin, which on the same sphere left angles of 0.6 to 179 degrees. The curved elements of order 2 and 3 are not a skin octree_mesh takes (it refuses triangle6 by name); mesh the volume from the linear result. tests/cpp/test_mesh_cad.cpp (MeshCad.OctreeMeshFillsTheSkinWithTetrahedra) and tests/python/test_mesh_cad.py pin these. A tetrahedral mesher that keeps the triangles exactly is an optional backend still to come (roadmap 6.7).

Limitations

This first milestone deliberately leaves out:

• Quality guarantees. Unrefined, the lattice gives near-equilateral triangles inside and edges within about [0.5 h, 1.5 h] on planar faces, with small angles near sharp corners and poles. min_angle bounds the rest, guaranteed to 20.7 degrees (Ruppert's B = sqrt 2) on a plane away from corners under 60 degrees, and measured past it; on a curved face the bound is the surface's, where termination is no longer the plane's theorem: max_points, the insertion-radius declines and the report are the guard. The smallest angle a corner under 60 degrees forces is measured and reported; with lop_corners the corner is lopped instead and the bound Pav and Walkington (2005) prove is stated on a planar face, checked against the mesh (item 22). On a curved face it is not stated: their theorem is the plane's.
• Sizing. The field's sources are the model's curvature and its gaps, measured on the first mesh, and the sizes given per face, per edge or by a background mesh; given sizes only refine below target_size, and a size from a solution's error estimate is a background the caller builds, or a metric from its Hessian (item 23). The size field is isotropic; the metric source's anisotropy is judged by refinement, which only adds points, so a stretched triangle comes from splitting the lattice's, not from a mesh built in the metric from the start, and the lattice and edges stay at the isotropic size. Sampled at the first mesh's points, a feature smaller than h that no point lands on (a hole narrower than h between two samples) is sized only as far as the curvature at the points around it reaches. A sized run meshes the model twice.
• Metrics. By default one affine map per face: a face whose parameterisation stretches very differently from place to place is meshed by its average. Where it stretches the plane more than twofold, refinement keeps to h / 4 and warns of the triangles it leaves under min_angle; a planar bilinear patch whose corner stretches 2.4-fold still meets 25 degrees at h from 0.2 to 0.05. metric="varying" follows the surface's metric there, but it is isotropic (the surface's own lengths, scaled by the size field), not an anisotropic metric a solution asks for, and it is opt-in: on the corpus it leaves 59 fewer faces under the bound, yet 6 fewer faces mesh, their folds outlasting the 16 rounds (on the RC buggy's suspension). A fold the varying metric keeps is meshed in the mean one.
• Virtual topology. virtual_edges merges slivers of edges and virtual_faces merges sliver faces into a neighbour, a chain of slivers growing from the face it starts on, at most eight faces a group; a chain no face that is not a sliver borders, or a sliver whose neighbours all meet it at a crease, is meshed face by face as before; on the corpus 37 slivers merged before chains did, and two faces beside merged ones fold on the boundary the pairs leave them. Suppressed corners and shared edges stay as orphan points, and both are order 1 only.
• Reparametrisation. reparam_condition_number flattens one face at a time, or a sliver with its host; it is order 1 only, so a curved element's nodes still sit at its face's own parameters. A face whose coarse mesh has a pinched boundary, or more than one hole that the conformal map folds, or whose flattening folds at a corner under 60 degrees, is declined and meshed in its parameters, with the reason reported; a band closed without a seam is a disc in its coarse mesh (the seam's copies stay apart), and a coarse mesh with one hole is laid between two circles when the conformal map folds it. On the CAD corpus at a threshold of 4 it flattens 10 faces and declines 15, and one face fails where a split its flattened neighbour asked for finds its trim leaving its patch, so it stays opt-in.
• Poles. Without pole_rings a pole is a fan apex, ungraded; refinement leaves its reach to it and reports it, so on the sphere the triangles near the poles stay below min_angle. The reach is where the face's latitude–longitude parameters stretch more than twofold, so it does not shrink with h: 0.375 on the unit sphere, 7.5% to 10% of its area at h from 0.2 to 0.05, where 40% to 47% of the triangles stay under 25 degrees at min_angle = 25. With rings the angles are met, but the curved elements of the pole's own fan keep their nodes at the degenerate parameters' placement, so their maximum interpolation error converges slower than elsewhere.
• Folds. A fold that the boundary rounds cannot reach before their limit of 16, or within the face's share of points, fails its face (one of the Ventilator model's 305 faces at min_angle = 25, none at 0).
• Embedded edges that meet. Embedded edges that cross each other (stiffeners in a grid) or end on the boundary between its corners are left out with a warning, not split where they meet.
• Touching curves. Loops whose curves touch, such as a hole tangent to the rim, stay self_intersecting.
• Curved elements. The cubic elements' area converges at order 4, not 6. A node inside a face sits at a parameter fraction, and placing it at the surface point nearest the chord instead was measured and made the cubic elements worse, not better (the parameter placement keeps every node the image of one lattice under one smooth map; roadmap 6.6.2), so the shortfall is not the placement. An element that folds is made straight-sided, not untangled by moving its nodes (6.6.3).
• Other elements. No quadrilaterals and no boundary layers; volume meshes come from octree_mesh, fitted to the skin rather than keeping its triangles.

Roadmap §6 lists what comes next. Everything is serial within a face; faces and edges run in parallel and the output is byte-identical across backends and thread counts.

Bindings

Surface	Entry points
Python	meshioplusplus.mesh_cad_surface(mesh_or_brep, target_size, …, min_angle=0, smooth_iterations=0, smooth_method="laplacian", order=1, curvature_deviation=0, cells_across_gaps=0, grading=0, pole_rings=0, face_sizes=None, edge_sizes=None, background=None, background_array="size", virtual_edges=False, metric="mean", reparam_condition_number=0, model_tolerances=False, virtual_faces=False, lop_corners=False, metric_background=None, metric_array="", metric_error=0, return_report=False), triangulate_polygon_2d(outer, holes=(), target_size=0, min_angle=0, max_points=1000000) (core only: without the compiled core both raise NotImplementedError)
C++	operations/mesh_cad.hpp: mesh_cad_surface(const BrepInfo&, const MeshCadSurfaceOptions&) and, since v17.19.0, mesh_cad_surface(…, const MeshCadSizeSources&, MeshCadSizeReport*) and, since v17.20.0, mesh_cad_surface(…, const MeshCadSizeSources&, const MeshCadSurfaceExtras&, MeshCadSizeReport*, MeshCadSurfaceExtrasReport*), triangulate_polygon_2d (and its TriangulatePolygon2dOptions overload)
C	mio_mesh_cad_surface on a mio_format_info, mio_mesh_cad_surface_opts / _report (with min_angle_deg, smooth_iterations, smooth_method, order, pole_rings, grading, curvature_deviation, cells_across_gaps, num_curved_straightened), per-face mio_mesh_cad_surface_face_status / _reason, mio_mesh_cad_surface_warnings, mio_mesh_cad_surface_angle_histogram, mio_mesh_cad_surface_num_corners / _corner; since v17.19.0 mio_mesh_cad_surface2(info, opts, opts_ex, report) with mio_mesh_cad_surface_opts_ex (metric, virtual_edges, the borrowed face_size_faces / face_sizes / num_face_sizes and edge arrays, background, background_array, since v17.20.0 reparam_condition_number, model_tolerances, virtual_faces, lop_corners and, from ABI 35, metric_background, metric_array and metric_error, a reserved tail), mio_mesh_cad_surface_face_size_limits, mio_mesh_cad_surface_num_virtual_edges / _virtual_edge; since v17.20.0 mio_mesh_cad_surface_num_reparametrised / _reparametrised / _reparam_reason, _num_merged_faces / _merged_face / _merged_face_edges / _merged_face_reason, _num_tolerance_uses / _tolerance_use (with mio_mesh_cad_tolerance_where) and _num_lopped_corners / _lopped_corner / _lopped_corner_reason. The options struct has no room left, so the new inputs are a second struct and entry point, as mio_fit_cad2 is; mio_mesh_cad_surface is mio_mesh_cad_surface2 with a NULL opts_ex
Fortran	info%mesh_cad_surface(target_size, …, min_angle, smooth_iterations, smooth_method, order, grading, curvature_deviation, cells_across_gaps, pole_rings, face_size_faces, face_sizes, edge_size_edges, edge_sizes, background, background_array, metric, virtual_edges, reparam_condition_number, model_tolerances, virtual_faces, lop_corners, metric_background, metric_array, metric_error) (indices 1-based), mio_mesh_cad_face_size_limits, mio_mesh_cad_num_reparametrised / mio_mesh_cad_get_reparametrised, mio_mesh_cad_num_merged_faces / mio_mesh_cad_get_merged_face, mio_mesh_cad_num_tolerance_uses / mio_mesh_cad_get_tolerance_use (MIO_MESH_CAD_TOLERANCE_*), mio_mesh_cad_num_lopped_corners / mio_mesh_cad_get_lopped_corner / mio_mesh_cad_lopped_corner_reason, mio_mesh_cad_surface_report, mio_mesh_cad_face_status / _reason, mio_mesh_cad_angle_histogram, mio_mesh_cad_num_corners / mio_mesh_cad_get_corner
Julia	mesh_cad_surface(info; target_size, …, min_angle, smooth_iterations, smooth_method, order, grading, curvature_deviation, cells_across_gaps, pole_rings, face_sizes, edge_sizes, background, background_array, metric, virtual_edges, reparam_condition_number, model_tolerances, virtual_faces, lop_corners, metric_background, metric_array, metric_error) (1-based indices; limited_by_edges per face, virtual_edges, reparametrised, merged_faces, tolerance_uses and lopped_corners in the result)
R	mio_mesh_cad_surface(info, target_size, …, min_angle, smooth_iterations, smooth_method, order, grading, curvature_deviation, cells_across_gaps, pole_rings, face_sizes, edge_sizes, background, background_array, metric, virtual_edges, reparam_condition_number, model_tolerances, virtual_faces, lop_corners, metric_background, metric_array, metric_error) (sizes as a named vector or a two-column matrix of 1-based indices; face_size_limits, virtual_edges and the data frames reparametrised, merged_faces, tolerance_uses and lopped_corners in the result)
WebAssembly	meshCadSurface(info, {targetSize, …, minAngle, smoothIterations, smoothMethod, order, curvatureDeviation, cellsAcrossGaps, grading, poleRings, faceSizes, edgeSizes, background, backgroundArray, metric, virtualEdges, reparamConditionNumber, modelTolerances, virtualFaces, lopCorners, metricBackground, metricArray, metricError}) (sizes as an object {index: size}; limitedByEdges per face, virtualEdges, reparametrised, mergedFaces, toleranceUses and loppedCorners in the result) on readMeshSelective(path, {info: true}).info, in both artifacts
CLIs	mesh-cad-surface IN OUT --size H [--deflection D] [--min-points-per-edge N] [--tolerance T] [--no-face-ids] [--no-edges] [--no-match] [--strict] [--min-angle DEG] [--smooth N] [--smooth-method laplacian|odt] [--order 1|2|3] [--curvature-deviation D] [--cells-across-gaps N] [--grading G] [--pole-rings M] [--face-sizes K=H,…] [--edge-sizes K=H,…] [--background PATH [--background-array NAME]] [--metric mean|varying] [--virtual-edges] [--reparam-condition-number K] [--model-tolerances] [--virtual-faces] [--lop-corners] [--metric-background PATH --metric-array NAME --metric-error E] [--json] in both CLIs
MCP	mesh_cad_surface (face_sizes and edge_sizes as JSON objects, background_path, background_array, metric, virtual_edges, reparam_condition_number, model_tolerances, virtual_faces, lop_corners, metric_background_path, metric_array, metric_error)

The pipeline has a MeshCadSurface step (v17.17.0): the B-rep a CAD JSON, STEP or IGES Input reads travels to it, through IgaRefine or from a FitCad step's fitted solid, and is dropped by the first step that changes the mesh, so {"Input": {"Path": "part.step"}, "Operations": [{"Op": "MeshCadSurface", "TargetSize": 2, "MinAngle": 25}, {"Op": "Quality"}], "Output": {"Path": "part.vtu"}} meshes a STEP solid in a run. Since v17.18.0 the step also takes CurvatureDeviation, CellsAcrossGaps, Grading and PoleRings, and since v17.19.0 FaceSizes and EdgeSizes (flat [index, size, …] lists), Background (a path) and BackgroundArray, Metric and VirtualEdges, and since v17.20.0 ReparamConditionNumber, ModelTolerances, VirtualFaces and LopCorners (the step's report adds NumReparametrised, NumMergedFaces and NumLoppedCorners), in both engines.

Verification

The predicates are checked against 128-bit integer determinants, including collinear, cocircular, coplanar and cospherical configurations and Kettner et al.'s classroom grid; the triangulation on squares, L-shapes, holes, collinear and near-duplicate points, slivers and random polygons for exact area, segment recovery, Euler characteristic and the Delaunay property of every unconstrained edge, and segments recovered after classification for an unchanged inside; the crossing search case by case, on one-ulp near misses and against a test of every pair on random polygons; the edge sampler on rational arcs for arc length to 1e-12, on reparametrised, reversed uses and on closed uses started a quarter, a half and three quarters of a turn later. The meshes of boxes, a capped cylinder through a seam, a sphere through its poles, an untrimmed rational patch, the trimmed plate and the CAD JSON fixtures are checked for watertightness, area and enclosed volume converging at second order, every point on the exact surface, the Hausdorff distance to a dense reference, determinism over thread counts and VTU and Gmsh round trips; a capped cylinder whose cap starts elsewhere, a plate with a hole 0.02 from its rim at h = 1.1 (and the same plate as a closed solid), a hole tangent to its rim, the two-patch beam's embedded curve, a fin sharing a plate's embedded edge and a closed embedded circle; with OpenCASCADE, a STEP solid against BRepMesh's own triangulation (tests/cpp/test_predicates.cpp, test_cdt2d.cpp, test_edge_sampler.cpp, test_mesh_cad.cpp, tests/python/test_mesh_cad.py). Refinement is checked both sides of its preconditions (test_cdt2d.cpp, test_mesh_cad_quality.cpp): on planar domains whose input angles are 60 degrees or more, every triangle meets 20.7, 25 and 30 degrees; adversarial domains (slivers down to a width of 1e-6, a needle-shaped hole, random star polygons) end under max_points with every triangle left below the bound attributed to a corner under 60 degrees; diametral and lens_angle agree with 128-bit integer evaluations, on their boundaries too. On the trimmed plate, a capped cylinder and a sphere at h = 0.5, 0.2 and 0.1 (and 0.05 measured), every triangle meets 25 degrees except those a reported pole forces (the sphere's smallest angle away from its poles is 36.8, 26.1, 26.0 and 25.2 degrees); area and enclosed volume still converge at second order; every boundary point a face's refinement adds is shared by every face of its edge; the mesh is the same over thread counts and a shuffled face order; smoothing never moves a boundary point nor lowers a star's smallest angle; and min_angle = 0 with no smoothing reproduces the unrefined mesh byte for byte, which the operation golden pins. The fuzz targets, run with refinement on, found five inputs refinement never finished or spent max_points on, each fixed and replayed in ctest (the fifth, a cylinder with a control point pulled out, by its input alone): a sliver whose split the triangulation refuses, a planar bilinear patch whose stretched corner ran refinement away, a reweighted roof that runs away at 30 degrees, and a face whose trim leaves its surface's knot range (Cdt2dRefine.ASplitTheTriangulationRefusesDeclinesItsTriangle, MeshCadQuality.AStretchedCornerDoesNotRunAway, MeshCadQuality.AFaceWhoseRefinementRunsAwayIsMeshedUnrefined, MeshCadQuality.ATrimLeavingItsSurfaceIsReported).

Refinement to 25 degrees costs 1.0 to 5.4 times the unrefined meshing on the CAD JSON fixtures (about 4 times on most, the sphere at h = 0.05 164 ms against 30) and adds under 8% more points (22% on the near-hole plate at h = 0.3); two smoothing passes about double that. The meshioplusplus_bench_ops row mesh_cad_surface_refined is hashed by the determinism job.

Interpolation error, the measure a mesh is for (Shewchuk 2002), is checked too (tests/cpp/test_mesh_cad_interpolation.cpp): the piecewise-linear interpolant of a smooth field is compared with the field at the points of the 7-point degree-5 rule of Dunavant (1985) on each triangle, the field evaluated in space there, as h halves. On the sphere (sin 2x cos 1.5y + exp(z/2)) the L2 error's observed orders are 2.01 and 2.03 unrefined and 2.06 and 2.03 at 25 degrees, and the maximum error's, away from the poles, 2.12 and 1.97; on the Scordelis–Lo roof the L2 orders are 2.05 and 2.01 and the maximum's 2.08 and 1.98; a tanh layer across the sphere converges at 1.90: the size field follows the geometry, not a solution, so a field's layer is left to a background size. That background is built from the field itself in MeshCadInterpolation.AHessianSizedMeshReachesTheLayerWithFewerPoints: on a mesh of the sphere at h = 0.1, hessian gives the layer's second derivatives, projected on the tangent plane (P H P, P = I - n nᵀ), and the size at a point is sqrt(8 ε / |λ|max), the length at which linear interpolation's error is about ε there (Loseille and Alauzet 2011), between h / 4 and 0.2, with ε = 0.008. Meshed from target_size = 0.2 with that background, both meshes at min_angle = 25 with two pole rings, the layer's L2 error is 1.39e-2 with 3,848 points, where the uniform mesh at h = 0.05 has 1.58e-2 with 11,281: the same error with a third of the points. The test pins a smaller error with under 0.4 times the points. The test pins the L2 orders to 2 ± 0.2 and the maximum's to 2 ± 0.35, beside the angle histogram that the report gives for the same meshes. The curved elements are measured the same way, their own Lagrange interpolant on the curved element integrated by an 8 by 8 collapsed Gauss–Legendre rule: on the roof the L2 orders are 3.09 and 3.05 for quadratic and 4.09 and 3.93 for cubic elements; on the sphere, away from the poles, 3.00 and 3.04, and 3.73 and 3.76, the cubic short of 4 as its area is short of 6 (placing the interior nodes by the surface rather than its parameters was measured and parked: it lowered the cubic orders, roadmap 6.6.2); over the whole sphere the poles' fans hold them to 2.76 and 2.53, and 3.73 and 3.43. With pole_rings = 2 they are 2.74 and 3.11, and 3.05 and 4.41, and the linear maximum's over the whole sphere 1.93 and 1.87 (1.34 at the last halving without rings); the curved maximum stays at the pole's own fan (MeshCadInterpolation.PoleRingsLiftTheWholeSphere). The test pins these with margin, and every order below the one beneath it at every size. Two fuzz targets feed mesh_cad_surface mutated CAD JSON models and triangulate_polygon_2d near-duplicate, near-collinear and touching polygons, and expect a result or a reported face, never a crash or a sanitizer report (fuzzing).

Curved elements are checked both sides of the validity test (tests/cpp/test_curved_triangle.cpp): the shape functions interpolate their nodes and the identity element has J = 1; the Bernstein coefficients reproduce known polynomials of degree 2 and 4; a quadratic and a cubic element whose J is positive at every node but negative inside (sampled on a lattice of 4186 points) is proved invalid, and so is one whose mid-edge node is pushed across the opposite edge; an element bent onto a cylinder is proved valid without reaching the depth limit. A quadratic element with a single mid-edge node moved was found valid wherever its nodes pass, over a 121 by 121 grid of positions, so the probe moves three. The meshes (tests/cpp/test_mesh_cad_curved.cpp, tests/python/test_mesh_cad.py) put every node of the sphere's quadratic and cubic elements on the sphere to 1e-12, share every edge node between exactly the two elements or the element and the line that use it, keep the linear mesh's points as their first points, and are the same over thread counts and a shuffled face order. Area and enclosed volume, integrated over the curved elements by an 8 by 8 collapsed Gauss–Legendre rule, converge on the capped cylinder and the sphere (h = 0.4, 0.2, 0.1) at observed orders of 2.0 to 2.1 for linear elements, 3.9 to 4.3 for quadratic and 4.0 to 4.3 for cubic, the cubic errors 6 to 20 times smaller than the quadratic; the test pins 1.8 and 3.5. The coarse capped cylinder at h = 1.75 with cubic elements folds some, which are made straight-sided, warned and counted, every element written then valid; at h = 0.3 none folds.

Sizing, poles and folds have their own probes (tests/cpp/test_mesh_cad_size.cpp, tests/cpp/test_edge_sampler.cpp, tests/python/test_mesh_cad.py). The NURBS second derivatives match central differences of the first on every face of a sphere, a capped cylinder (both rational) and a box, off their knots, and the principal curvatures are exact on the analytic surfaces: 1 / R on the sphere, 1 / R and 0 on a cylinder, 0 on a plane. On the capped cylinder of radius 1 at h = 0.2, curvature_deviation = 0.002 keeps every triangle within 0.00197 of the surface (0.000498 at 0.0005), with 5,226 triangles where a uniform mesh of the same chordal error needs 5,618 (19,518 against 21,360); a plate 0.05 thick at h = 0.25 gets at most 0.025 across its rim with cells_across_gaps = 2 and 0.0166 with 3; the field equals a brute-force evaluation of its envelope and is Lipschitz with constant grading - 1; gradings of 2, 1.5 and 1.2 spend 5,138, 5,222 and 5,370 triangles; a sized mesh is the same over thread counts and a shuffled face order. On the sphere at min_angle = 25 two pole rings give 26.14 and 25.12 degrees at h = 0.2 and 0.1 with valence 6 at both poles. Narrow torus bands do not fold; a band swept 270 to 340 degrees, meshed coarser than itself, does, and is unfolded with every triangle facing out at its corners (MeshCadFolds.ABandSweptPastHalfATurnIsUnfolded). With every new option at 0 the operation goldens do not move. The mesh_cad_surface_sized bench row (a slab two lattice steps thick, two cells across its walls) is hashed by the determinism job. The CAD corpus meshes real STEP models in a scheduled job: its first runs found the fold above on several hundred faces, an arc-length table that halved towards 2^30 pieces on short traces far from the origin, and a STEP read that leaked OpenCASCADE's exception, each fixed with its own test. Its v17.18.0 failures were dominated by two classes, both fixed in v17.19.0 with synthetic repros: 95 faces whose parameters lay a hair off zero, outside the exact predicates' range, now flushed to zero when (and only when) the face would otherwise be refused (MeshCad.ParametersAHairOffZeroAreMeshed), and 118 surfaces of revolution smaller than h whose every triangle collapsed onto the seam or a pole, now meshed again with a small interior grid (test_a_surface_of_revolution_smaller_than_h_is_meshed).

Given sizes, virtual edges and the varying metric have their own probes (tests/cpp/test_mesh_cad_size.cpp, test_mesh_cad_quality.cpp, test_cdt2d.cpp, tests/python/test_mesh_cad.py). A face of the unit cube given 0.1 at target_size = 0.4 meshes with a median edge of 0.1 within 15% while its neighbours stay coarser; beside a face given 0.4, exactly the shared edge is reported as limiting it; an edge given 0.05 is sampled at most that apart; a background plane mesh whose size grows linearly in x sizes a box face within 25% per strip; empty sources are the two-argument overload byte for byte, and bad ones are refused. Two plates folded along a rim split by an edge of h / 50 keep a triangle of 1.0 degree unrefined and need 152 triangles at 25 degrees; with virtual_edges the smallest unrefined angle is 23.8 degrees and 25 degrees takes 88 triangles, the plates conform across the rim and the chain is reported. Cdt2d with the identity metric is the plane's triangulation exactly, and with a stretched one is valid and rounder in that metric. Under metric="varying" the cone frustum, the torus band and the stretched bilinear patch meet 25 degrees unwarned with at most the mean metric's points, the same at every thread count, and the box, the capped cylinder and the sphere stay watertight; a band swept past half a turn is unfolded in either metric.

The reparametrisation, virtual faces and model tolerances have their own probes (tests/cpp/test_mesh_cad_reparam.cpp, test_mesh_cad_quality.cpp, tests/python/test_mesh_cad.py). A planar patch whose side v = 1 collapses onto a point (its corner 26.6 degrees) and a dome over the same collapse leave triangles under 25 degrees at min_angle = 25 in both the mean and the varying metric at h = 0.2 and 0.1; with reparam_condition_number = 4 each is flattened conformally and meets 25 degrees unwarned (26.6 and 28.0 degrees at h = 0.2), the plane's area exact to 1e-12. The same patch trimmed short of its collapse, a corner its parameters stretch thirtyfold, keeps a triangle under 5 degrees and a warning in the mean metric; flattened it meets 25 degrees with fewer points than the varying metric, every point on the patch. Both hemispheres of the unit sphere are flattened through their poles and meet 25 degrees at h = 0.2, 0.1 and 0.05 (25.3, 25.4 and 25.1 degrees, where two pole rings leave 7.4 at h = 0.05), watertight, every point on the sphere to 1e-12, the area within 2 h² of 4π, and the same over thread counts and a shuffled face order. A threshold no face passes is the option off byte for byte; a box is never flattened, and a capped cylinder, whose rational circles vary their speed to a condition number of about 1.4, not at 2. The collapsed coordinates are continuous across an edge out of a pole. Three coplanar plates whose middle one is h / 20 wide need 796 triangles at h = 0.2 and 25 degrees face by face; with virtual_faces the sliver is merged into a neighbour, 170 triangles meet 25 degrees (25.35), every triangle lies in the face its cad:face names, the strip is open only along its outside, the area is exact and the mesh is the same over thread counts. With its far neighbour turned up 90 degrees, the sliver merges with the flat one only. Two stacked bands whose shared rim's uses start 1e-5 apart, more than the global tolerance, fail a face; with the rim declaring 1e-4 and model_tolerances, both mesh, and the use is reported as closed_start; a gap past the declaration still fails, a negative declaration is warned, and a box with nothing declared is unchanged. With OpenCASCADE, every edge of a STEP read declares its tolerance. The meshioplusplus_bench_ops rows mesh_cad_surface_reparam and mesh_cad_surface_virtual_faces are hashed by the determinism job, and the CAD JSON fuzz target turns the three options on in a share of its inputs.

References

• J. R. Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, Discrete & Computational Geometry 18 (1997).
• L. P. Chew, Guaranteed-quality mesh generation for curved surfaces, 9th ACM Symposium on Computational Geometry (1993).
• G. L. Miller, S. E. Pav and N. J. Walkington, When and why Ruppert's algorithm works, 12th International Meshing Roundtable (2003).
• S. E. Pav and N. J. Walkington, Delaunay refinement by corner lopping, 14th International Meshing Roundtable (2005).
• J.-D. Boissonnat and S. Oudot, Provably good sampling and meshing of surfaces, Graphical Models 67 (2005).
• D. A. Field, Laplacian smoothing and Delaunay triangulations, Communications in Applied Numerical Methods 4 (1988).
• L. Chen and J. Xu, Optimal Delaunay triangulations, Journal of Computational Mathematics 22 (2004).
• J. R. Shewchuk, What is a good linear element? Interpolation, conditioning, and quality measures, 11th International Meshing Roundtable (2002).
• D. A. Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the triangle, IJNME 21 (1985).
• L. Kettner, K. Mehlhorn, S. Pion, S. Schirra and C. Yap, Classroom examples of robustness problems in geometric computations, Computational Geometry 40 (2008).
• J. R. Shewchuk, Delaunay Refinement Algorithms for Triangular Mesh Generation, Computational Geometry 22 (2002).
• C. L. Lawson, Software for C1 surface interpolation, Mathematical Software III (1977).
• M. V. Anglada, An improved incremental algorithm for constructing restricted Delaunay triangulations, Computers & Graphics 21 (1997).
• J. Ruppert, A Delaunay refinement algorithm for quality 2-dimensional mesh generation, Journal of Algorithms 18 (1995).
• L. P. Chew, Guaranteed-quality triangular meshes, Cornell TR 89-983 (1989).
• C. Geuzaine and J.-F. Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, IJNME 79 (2009).
• J. Schöberl, NETGEN: an advancing front 2D/3D-mesh generator based on abstract rules, Computing and Visualization in Science 1 (1997).
• H. Si, TetGen, a Delaunay-based quality tetrahedral mesh generator, ACM TOMS 41 (2015).
• S.-W. Cheng, T. K. Dey and J. R. Shewchuk, Delaunay Mesh Generation, CRC Press (2012).
• A. Johnen, J.-F. Remacle and C. Geuzaine, Geometrical validity of curvilinear finite elements, Journal of Computational Physics 233 (2013).
• S. J. Owen, A survey of unstructured mesh generation technology, 7th International Meshing Roundtable (1998).
• Open CASCADE Technology, Mesh user's guide (BRepMesh, linear and angular deflection).
• P.-O. Persson, Mesh size functions for implicit geometries and PDE-based gradient limiting, Engineering with Computers 22 (2006).
• H. Borouchaki, F. Hecht and P. J. Frey, Mesh gradation control, IJNME 43 (1998).
• L. Piegl and W. Tiller, The NURBS Book, 2nd ed., Springer (1997): A2.3 (basis function derivatives) and A4.4 (derivatives of a rational surface).
• H. Borouchaki, P. L. George, F. Hecht, P. Laug and E. Saltel, Delaunay mesh generation governed by metric specifications. Part I. Algorithms, Finite Elements in Analysis and Design 25 (1997).
• A. Sheffer, T. D. Blacker and M. Bercovier, Virtual topology operators for meshing, International Journal of Computational Geometry & Applications 10 (2000).
• J.-F. Remacle, C. Geuzaine, G. Compère and E. Marchandise, High-quality surface remeshing using harmonic maps, IJNME 83 (2010).
• B. Lévy, S. Petitjean, N. Ray and J. Maillot, Least squares conformal maps for automatic texture atlas generation, ACM Transactions on Graphics 21 (SIGGRAPH 2002).
• M. S. Floater, Mean value coordinates, Computer Aided Geometric Design 20 (2003).
• M. G. Duffy, Quadrature over a pyramid or cube of integrands with a singularity at a vertex, SIAM Journal on Numerical Analysis 19 (1982).
• A. Loseille and F. Alauzet, Continuous mesh framework part I: well-posed continuous interpolation error, SIAM Journal on Numerical Analysis 49 (2011).

No code from Triangle, TetGen, Gmsh, CGAL or the CDT library is read or vendored; see THIRD_PARTY_NOTICES.md.