From cd3bd4d90e467b156cf249ad492e7b5d31d5939e Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 6 Jul 2026 09:42:26 +0930 Subject: [PATCH 01/12] feat: preserve edge constraints through cuts, add corefine and implicit-cut snapping - Track fixed-edge constraints by geometry (endpoint coords) instead of by Edge_index so they survive clip()/collect_garbage() and the exact-kernel round trip, via snapshot_constraint_coords()/rebuild_fixed_edges_from_coords(). Re-resolve constraints after cutWithSurface and clipWithPlane. - Add TriMesh::corefine: insert the shared intersection polyline into both meshes so patches from each side stitch together watertight. - cut_with_implicit_function: add snap_tol to reuse an existing endpoint when a crossing lands within that fraction of an edge, avoiding sliver triangles; nudge exactly-on-value properties to the positive side so every vertex is unambiguously classified. - export_mesh: guard non-positive duplicate_vertex_threshold (1/0 quantisation was scaling coords to non-finite and silently dropping the whole mesh). - Add regression tests for user constraints surviving cut + protected remesh. Co-Authored-By: Claude Opus 4.8 (1M context) --- loop_cgal/bindings.cpp | 12 +- src/mesh.cpp | 182 ++++++++++++++++++++++++++- src/mesh.h | 29 ++++- src/meshutils.cpp | 17 ++- tests/test_fixed_edge_constraints.py | 153 ++++++++++++++++++++++ 5 files changed, 386 insertions(+), 7 deletions(-) create mode 100644 tests/test_fixed_edge_constraints.py diff --git a/loop_cgal/bindings.cpp b/loop_cgal/bindings.cpp index 4eee326..e9bbabd 100644 --- a/loop_cgal/bindings.cpp +++ b/loop_cgal/bindings.cpp @@ -31,6 +31,13 @@ PYBIND11_MODULE(_loop_cgal, m) "Clip the mesh with the halfspace ax+by+cz+d <= 0. " "Uses PMP::clip(mesh, Plane_3) directly — no corefinement, no skirt construction. " "Returns the number of faces removed (0 = no-op).") + .def("corefine", &TriMesh::corefine, py::arg("other"), + py::arg("use_exact_kernel") = true, + "Corefine this mesh with another, inserting the shared intersection " + "polyline into BOTH meshes in place. After the call the two meshes " + "carry coincident vertices/edges along the intersection curve, so " + "patches taken from each side stitch together watertight. " + "Returns the number of vertices added to this mesh.") .def("cut_with_surface", &TriMesh::cutWithSurface, py::arg("surface"), py::arg("preserve_intersection") = false, py::arg("preserve_intersection_clipper") = false, @@ -51,7 +58,10 @@ PYBIND11_MODULE(_loop_cgal, m) "Vertex index pairs defining edges to be fixed in mesh when remeshing.") .def("cut_with_implicit_function", &TriMesh::cut_with_implicit_function, py::arg("property"), py::arg("value"),py::arg("cutmode") = ImplicitCutMode::KEEP_POSITIVE_SIDE, - "Cut the mesh with an implicit function defined by vertex properties.") + py::arg("snap_tol") = 1e-4, + "Cut the mesh with an implicit function defined by vertex properties. " + "snap_tol reuses an existing vertex when a crossing lands within that " + "fraction of an edge from it, preventing sliver triangles (0 = disable).") .def("_cgal_area", &TriMesh::area, "Surface area computed directly from CGAL mesh.") .def("_cgal_n_faces", &TriMesh::n_faces, "Number of faces in the CGAL mesh.") .def("_cgal_n_vertices", &TriMesh::n_vertices, "Number of vertices in the CGAL mesh.") diff --git a/src/mesh.cpp b/src/mesh.cpp index 7940391..7cf673c 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -121,6 +121,58 @@ void TriMesh::init() _edge_is_constrained_map = CGAL::make_boolean_property_map(_fixedEdges); } +std::vector> TriMesh::snapshot_constraint_coords() const +{ + std::vector> coords; + coords.reserve(_fixedEdges.size()); + for (auto e : _fixedEdges) + { + // _fixedEdges may hold edges that clip() has since tombstoned (removed but + // not yet garbage-collected); skip those so we never dereference a dead + // halfedge or capture a stale position. + if (!_mesh.has_valid_index(e) || _mesh.is_removed(e)) + continue; + auto h = _mesh.halfedge(e); + coords.emplace_back(_mesh.point(_mesh.source(h)), _mesh.point(_mesh.target(h))); + } + return coords; +} + +void TriMesh::rebuild_fixed_edges_from_coords( + const std::vector> &saved) +{ + // Border edges of the (possibly newly built) mesh are always constraints — + // this re-derives the outer boundary plus any fresh cut boundary left by the + // clip, exactly as init() would. + _fixedEdges = collect_border_edges(_mesh); + + // Coordinate → vertex lookup for the current mesh. Positions are compared by + // exact double value: surviving vertices keep their coordinates bit-for-bit + // across clip()+collect_garbage() and across the exact round trip (an exact + // number built from a double converts back to that same double), so no + // tolerance is needed to re-find them. + std::map, TriangleMesh::Vertex_index> lut; + for (auto v : _mesh.vertices()) + { + const auto &p = _mesh.point(v); + lut.emplace(std::make_tuple(p.x(), p.y(), p.z()), v); + } + + for (const auto &pr : saved) + { + auto it0 = lut.find(std::make_tuple(pr.first.x(), pr.first.y(), pr.first.z())); + auto it1 = lut.find(std::make_tuple(pr.second.x(), pr.second.y(), pr.second.z())); + if (it0 == lut.end() || it1 == lut.end()) + continue; // an endpoint no longer exists — constraint did not survive + auto h = _mesh.halfedge(it0->second, it1->second); + if (h == TriangleMesh::null_halfedge()) + continue; // endpoints exist but are no longer directly connected + _fixedEdges.insert(_mesh.edge(h)); + } + + _edge_is_constrained_map = CGAL::make_boolean_property_map(_fixedEdges); +} + // Move constructor: _edge_is_constrained_map stores a raw pointer into // _fixedEdges, so after relocating _fixedEdges we must re-bind the map. TriMesh::TriMesh(TriMesh&& other) noexcept @@ -479,13 +531,38 @@ int TriMesh::cutWithSurface(TriMesh &clipper, try { if (use_exact_kernel){ + // The exact round trip rebuilds _mesh from scratch (convert_to_double_mesh), + // so every Edge_index in _fixedEdges is invalidated and the constrained + // map cannot be threaded through the exact clip. Capture constraints by + // geometry beforehand and re-resolve them against the rebuilt mesh after. + auto saved_constraints = snapshot_constraint_coords(); Exact_Mesh exact_clipper = convert_to_exact(scaled_clipper); Exact_Mesh exact_mesh = convert_to_exact(*this); flag = PMP::clip(exact_mesh, exact_clipper, CGAL::parameters::clip_volume(false)); set_mesh(convert_to_double_mesh(exact_mesh)); + if (flag) + rebuild_fixed_edges_from_coords(saved_constraints); } else{ - flag = PMP::clip(_mesh, scaled_clipper._mesh, CGAL::parameters::clip_volume(false)); + // Pass the constrained-edge map for both meshes: CGAL reads existing + // constraints on input and, on output, marks the intersection edges and + // any surviving/split constraint edges. After the clip, _fixedEdges holds + // valid (pre-garbage-collection) indices, so snapshot it by geometry, + // collect garbage, then re-resolve — Edge_index values do not survive + // collect_garbage(). + flag = PMP::clip( + _mesh, scaled_clipper._mesh, + CGAL::parameters::edge_is_constrained_map(_edge_is_constrained_map) + .clip_volume(false), + CGAL::parameters::edge_is_constrained_map( + scaled_clipper._edge_is_constrained_map)); + if (flag) + { + auto saved_constraints = snapshot_constraint_coords(); + if (_mesh.has_garbage()) + _mesh.collect_garbage(); + rebuild_fixed_edges_from_coords(saved_constraints); + } } } catch (const std::exception &e) @@ -530,6 +607,62 @@ int TriMesh::cutWithSurface(TriMesh &clipper, return faces_before - faces_after; } +int TriMesh::corefine(TriMesh &other, bool use_exact_kernel) +{ + if (_mesh.number_of_faces() == 0 || other._mesh.number_of_faces() == 0) + { + std::cerr << "Error: corefine called on an empty mesh." << std::endl; + return 0; + } + + // Merge collocated border vertices first — corefinement needs valid, + // duplicate-free borders to compute a clean intersection polyline. + PMP::stitch_borders(_mesh); + PMP::remove_isolated_vertices(_mesh); + PMP::stitch_borders(other._mesh); + PMP::remove_isolated_vertices(other._mesh); + + const std::size_t v_before = _mesh.number_of_vertices(); + + try + { + if (use_exact_kernel) + { + Exact_Mesh em1 = convert_to_exact(*this); + Exact_Mesh em2 = convert_to_exact(other); + PMP::corefine(em1, em2); + set_mesh(convert_to_double_mesh(em1)); + other.set_mesh(convert_to_double_mesh(em2)); + } + else + { + PMP::corefine(_mesh, other._mesh); + } + } + catch (const std::exception &e) + { + std::cerr << "corefine failed: " << e.what() << std::endl; + return 0; + } + + if (_mesh.has_garbage()) + _mesh.collect_garbage(); + if (other._mesh.has_garbage()) + other._mesh.collect_garbage(); + + // Re-derive border edges (the intersection split changed topology on both). + init(); + other.init(); + + const std::size_t v_after = _mesh.number_of_vertices(); + if (LoopCGAL::verbose) + { + std::cout << "corefine: this mesh " << v_before << " -> " << v_after + << " vertices." << std::endl; + } + return static_cast(v_after) - static_cast(v_before); +} + int TriMesh::clipWithPlane(double a, double b, double c, double d, bool use_exact_kernel) { if (_mesh.number_of_vertices() == 0 || _mesh.number_of_faces() == 0) @@ -540,6 +673,12 @@ int TriMesh::clipWithPlane(double a, double b, double c, double d, bool use_exac const int faces_before = static_cast(_mesh.number_of_faces()); + // Capture constraints by geometry before the clip. The plane overload of + // PMP::clip does not accept an edge_is_constrained_map (it uses an internal + // one), and both the exact round trip and collect_garbage() below invalidate + // every Edge_index, so _fixedEdges must be re-resolved by coordinate after. + auto saved_constraints = snapshot_constraint_coords(); + try { if (use_exact_kernel) @@ -572,6 +711,13 @@ int TriMesh::clipWithPlane(double a, double b, double c, double d, bool use_exac // can accumulate enough tombstones to make get_points() take seconds. if (_mesh.has_garbage()) _mesh.collect_garbage(); + + // Re-derive the constraint set: border edges of the clipped mesh (including + // the new cut boundary) plus any user-added interior constraints that + // survived, re-resolved by coordinate. Without this, stale Edge_index values + // in _fixedEdges would mark unrelated (recycled) edges as protected during a + // later remesh(). + rebuild_fixed_edges_from_coords(saved_constraints); return faces_before - faces_after; } @@ -581,7 +727,7 @@ NumpyMesh TriMesh::save(double area_threshold, return export_mesh(_mesh, area_threshold, duplicate_vertex_threshold); } -void TriMesh::cut_with_implicit_function(const std::vector &property, double value, ImplicitCutMode cutmode) +void TriMesh::cut_with_implicit_function(const std::vector &property, double value, ImplicitCutMode cutmode, double snap_tol) { std::cout << "Cutting mesh with implicit function at value " << value << std::endl; std::cout << "Mesh has " << _mesh.number_of_vertices() << " vertices and " @@ -600,6 +746,17 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do { vertex_properties[vim[v]] = property[vim[v]]; } + + // The sign-based case table below uses strict >/< comparisons, so a vertex + // whose property equals `value` exactly matches neither side: its triangle is + // left un-split and the isocontour is torn where it should pass through that + // vertex. Nudge any exactly-on-value property to the next representable + // double (consistently to the positive side) so every vertex is unambiguously + // classified; NaN (off-extent) compares false and is left untouched. The + // shift is sub-ULP, so it does not move the seam for any non-degenerate input. + for (auto &vp : vertex_properties) + if (vp == value) + vp = std::nextafter(value, HUGE_VAL); auto property_map = boost::make_iterator_property_map( vertex_properties.begin(), vim); @@ -713,6 +870,27 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do ratio = 0.5; else ratio = (value - f0) / denom; + + // Snap the crossing to an existing endpoint when it lands within snap_tol + // (measured as a fraction of the edge) of that endpoint. Otherwise the + // crossing point sits an infinitesimal distance from a real vertex, and + // the sub-triangles built from it in the case table below are slivers that + // a downstream pass would have to collapse. Snapping reuses the existing + // vertex instead: any sub-triangle that collapses to a repeated index is + // dropped by the degeneracy check when the new mesh is assembled, so no + // sliver is ever created. newverts starts as a copy of verts, so the + // original vertex indices ends.first/ends.second are valid indices into it. + if (ratio <= snap_tol) + { + new_point_on_edge[eid] = ends.first; + continue; + } + if (ratio >= 1.0 - snap_tol) + { + new_point_on_edge[eid] = ends.second; + continue; + } + Point p0 = verts[ends.first]; Point p1 = verts[ends.second]; Point np = Point(p0.x() + ratio * (p1.x() - p0.x()), p0.y() + ratio * (p1.y() - p0.y()), p0.z() + ratio * (p1.z() - p0.z())); diff --git a/src/mesh.h b/src/mesh.h index e153e33..b7e3c32 100644 --- a/src/mesh.h +++ b/src/mesh.h @@ -48,12 +48,25 @@ class TriMesh int clipWithPlane(double a, double b, double c, double d, bool use_exact_kernel = true); + // Corefine this mesh with another, inserting the SHARED intersection + // polyline into BOTH meshes. After the call the two meshes carry + // coincident vertices and edges along the intersection curve, so patches + // taken from each (e.g. a contact cut by a fault and the matching fault + // patch) stitch together watertight. Both meshes are mutated in place. + // Returns the number of intersection polyline segments inserted. + int corefine(TriMesh &other, bool use_exact_kernel = true); + // Method to remesh the triangle mesh void remesh(bool split_long_edges, double target_edge_length, int number_of_iterations, bool protect_constraints, bool relax_constraints); void init(); - void cut_with_implicit_function(const std::vector& property, double value, ImplicitCutMode cutmode = ImplicitCutMode::KEEP_POSITIVE_SIDE); + // Cut the mesh along the isocontour property==value. + // snap_tol: when an edge crossing lands within this fraction of an edge + // from an existing endpoint, reuse that endpoint instead of inserting a + // near-coincident vertex. This prevents the sliver triangles that a raw + // insertion produces (0 = never snap, exact old behaviour). + void cut_with_implicit_function(const std::vector& property, double value, ImplicitCutMode cutmode = ImplicitCutMode::KEEP_POSITIVE_SIDE, double snap_tol = 1e-4); // Getters for mesh properties void reverseFaceOrientation(); NumpyMesh save(double area_threshold, double duplicate_vertex_threshold); @@ -72,6 +85,20 @@ class TriMesh // Internal constructor used by cutWithSurface to wrap a scaled copy. explicit TriMesh(TriangleMesh m) : _mesh(std::move(m)) { init(); } + // Capture the endpoint coordinates of every edge currently in + // _fixedEdges that is still live in _mesh. Edge_index handles are not + // stable across clip()/collect_garbage() or the exact-kernel round + // trip, so constraints must be tracked by geometry, not by index. + std::vector> snapshot_constraint_coords() const; + + // Rebuild _fixedEdges (and re-bind the constrained-edge map) from a + // coordinate snapshot: border edges of the current _mesh are always + // re-derived, and each saved endpoint pair is re-resolved to an edge + // descriptor in the current _mesh when both endpoints still exist and + // are directly connected. Use after any op that invalidates indices. + void rebuild_fixed_edges_from_coords( + const std::vector> &saved); + std::set _fixedEdges; TriangleMesh _mesh; // The underlying CGAL surface mesh CGAL::Boolean_property_map> diff --git a/src/meshutils.cpp b/src/meshutils.cpp index 7bffc50..6b39e58 100644 --- a/src/meshutils.cpp +++ b/src/meshutils.cpp @@ -84,20 +84,31 @@ NumpyMesh export_mesh(const TriangleMesh &tm, double area_threshold, } }; - const double inv = 1.0 / duplicate_vertex_threshold; // quantisation + // A non-positive threshold means "do not merge duplicates". Quantising with + // inv = 1/threshold would then be +inf, scaling every coordinate to a + // non-finite value that the isfinite guard below rejects — silently dropping + // the entire mesh. Guard it explicitly and emit every vertex as unique. + const bool dedup = duplicate_vertex_threshold > 0.0; + const double inv = dedup ? 1.0 / duplicate_vertex_threshold : 0.0; // quantisation std::unordered_map qmap; // grid → index int next_idx = 0; for (VIndex v : tm.vertices()) { const auto &p = tm.point(v); - + // Validate vertex coordinates if (!std::isfinite(p.x()) || !std::isfinite(p.y()) || !std::isfinite(p.z())) { if (LoopCGAL::verbose) std::cout << "Warning: Non-finite vertex coordinates, skipping vertex\n"; continue; } - + + if (!dedup) { // merging disabled → keep every vertex distinct + vertices.push_back({p.x(), p.y(), p.z()}); + vertex_index_map[v] = next_idx++; + continue; + } + // Compute quantized key with overflow protection double x_scaled = p.x() * inv; double y_scaled = p.y() * inv; diff --git a/tests/test_fixed_edge_constraints.py b/tests/test_fixed_edge_constraints.py new file mode 100644 index 0000000..9a20c53 --- /dev/null +++ b/tests/test_fixed_edge_constraints.py @@ -0,0 +1,153 @@ +"""Regression tests: user-added edge constraints must survive a cut and be +honoured by a subsequent protected remesh. + +Background +---------- +``add_fixed_edges`` records interior (non-border) edges that must be preserved +when the mesh is remeshed. ``cut_with_surface`` / ``clip_with_plane`` mutate the +mesh via CGAL's clip, which tombstones/garbage-collects elements and (for the +exact-kernel path) rebuilds the mesh from scratch — both of which invalidate the +``Edge_index`` handles stored for the constraints. Before the fix the constraint +set was never refreshed, so after a cut it either silently pointed at recycled +edges or was thrown away entirely, and the subsequent ``remesh`` freely erased +the constrained polyline. + +Discriminator +------------- +A constrained edge may legitimately be *split* into collinear sub-edges by +remesh (CGAL keeps the sub-edges constrained), but the polyline itself must stay +exactly on its original line with both endpoints intact. When the constraint is +lost the region is remeshed freely and no mesh edges remain on the line. We +therefore count mesh edges lying exactly on the constrained segment: preserved → +a connected chain covering the full span; lost → zero. +""" +from __future__ import annotations +import numpy as np +import pytest +import pyvista as pv +import loop_cgal + + +NX, NY = 11, 4 # grid columns x=0..10, rows y=0..3 (rows 1,2 are interior) +# Fully interior vertical edge at x=8 spanning y=1..2 (both endpoints interior). +EDGE_X = 8.0 +EDGE_Y0, EDGE_Y1 = 1.0, 2.0 +A_COORD = np.array([EDGE_X, EDGE_Y0, 0.0]) +B_COORD = np.array([EDGE_X, EDGE_Y1, 0.0]) + + +def _build_grid(): + verts = np.array( + [[float(i), float(j), 0.0] for j in range(NY) for i in range(NX)], + dtype=np.float64, + ) + + def vid(i, j): + return j * NX + i + + tris = [] + for j in range(NY - 1): + for i in range(NX - 1): + a, b = vid(i, j), vid(i + 1, j) + c, d = vid(i + 1, j + 1), vid(i, j + 1) + tris += [[a, b, c], [a, c, d]] + return verts, np.array(tris, dtype=np.int32), vid + + +def _on_line(p): + """True if point p lies exactly on the constrained segment (x=8, z=0, 1<=y<=2).""" + return ( + abs(p[0] - EDGE_X) < 1e-9 + and abs(p[2]) < 1e-9 + and (EDGE_Y0 - 1e-9) <= p[1] <= (EDGE_Y1 + 1e-9) + ) + + +def _line_edges(verts, tris): + """Return the set of mesh edges (as index pairs) lying on the constrained segment.""" + edges = set() + for t in tris: + for a, b in ((t[0], t[1]), (t[1], t[2]), (t[2], t[0])): + if _on_line(verts[a]) and _on_line(verts[b]): + edges.add((min(int(a), int(b)), max(int(a), int(b)))) + return edges + + +def _has_coord(verts, coord): + return bool(np.any(np.all(np.isclose(verts, coord), axis=1))) + + +def _cut_and_remesh(add_constraint, method, use_exact_kernel): + verts, tris, vid = _build_grid() + tm = loop_cgal.TriMesh.from_vertices_and_triangles(verts, tris) + if add_constraint: + tm.add_fixed_edges(np.array([[vid(8, 1), vid(8, 2)]], dtype=np.int32)) + + if method == "surface": + # Vertical clipper at x=5, normal -x, so the side x>5 (containing the + # x=8 constrained edge) is kept. The cut boundary is near the edge but + # does not cross it. + clipper = pv.Plane( + center=(5.0, 1.5, 0.0), direction=(-1.0, 0.0, 0.0), + i_size=8.0, j_size=8.0, i_resolution=4, j_resolution=4, + ).triangulate() + removed = tm.cut_with_surface( + loop_cgal.TriMesh(clipper), use_exact_kernel=use_exact_kernel + ) + assert removed > 0, "cut removed nothing — clipper did not intersect" + else: # plane: keep -x+5 < 0 <=> x > 5 + removed = tm.clip_with_plane(-1.0, 0.0, 0.0, 5.0, use_exact_kernel=use_exact_kernel) + assert removed > 0, "plane clip removed nothing" + + # target_edge_length (0.3) is shorter than the 1.0-long constrained edge, so + # a lost constraint would let remesh freely re-triangulate the region. + tm.remesh(True, 0.3, 3, True, False) + return tm.get_vertices_and_triangles() + + +@pytest.mark.parametrize("use_exact_kernel", [True, False]) +@pytest.mark.parametrize("method", ["surface", "plane"]) +def test_user_constraint_survives_cut_and_remesh(method, use_exact_kernel): + """A user-added interior constraint must survive the cut and be honoured by remesh.""" + verts, tris = _cut_and_remesh(True, method, use_exact_kernel) + + # 1. Both original endpoints are still present, exactly (constrained edge + # endpoints are never moved by smoothing). + assert _has_coord(verts, A_COORD), f"[{method}/{use_exact_kernel}] endpoint (8,1,0) lost" + assert _has_coord(verts, B_COORD), f"[{method}/{use_exact_kernel}] endpoint (8,2,0) lost" + + # 2. The constrained polyline is tiled by mesh edges lying exactly on the + # line — a lost constraint leaves zero. + edges = _line_edges(verts, tris) + assert len(edges) >= 1, ( + f"[{method}/{use_exact_kernel}] no mesh edges on the constrained line — " + f"constraint was not preserved through the cut" + ) + + # 3. The chain covers the full span y=1..2 (endpoints connected through + # collinear sub-vertices, all exactly on the line). + on_line_ys = sorted( + p[1] for p in verts if _on_line(p) + ) + assert on_line_ys[0] == pytest.approx(EDGE_Y0) + assert on_line_ys[-1] == pytest.approx(EDGE_Y1) + + +@pytest.mark.parametrize("use_exact_kernel", [True, False]) +@pytest.mark.parametrize("method", ["surface", "plane"]) +def test_constraint_makes_a_difference(method, use_exact_kernel): + """Control: the same cut+remesh without the constraint erases the polyline. + + This pins the discriminating power of the test above — the only difference + between the two runs is the ``add_fixed_edges`` call. + """ + _, tris_with = _cut_and_remesh(True, method, use_exact_kernel) + verts_with, tris_with = _cut_and_remesh(True, method, use_exact_kernel) + verts_without, tris_without = _cut_and_remesh(False, method, use_exact_kernel) + + with_edges = len(_line_edges(verts_with, tris_with)) + without_edges = len(_line_edges(verts_without, tris_without)) + assert with_edges > without_edges, ( + f"[{method}/{use_exact_kernel}] constraint made no difference " + f"(with={with_edges}, without={without_edges}) — it is not being honoured" + ) From 6c9e17585aa459c715356aee691b571bca62800a Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 6 Jul 2026 09:43:00 +0930 Subject: [PATCH 02/12] fix: remove duplicate face-insertion loop in array constructor MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The numpy array constructor looped over the triangle array twice: the first loop validated indices and skipped degenerate/out-of-range faces, the second re-added every face with no validation. The second loop was dead work for valid meshes (re-adding an existing face returns null_face in CGAL) and, worse, re-introduced the out-of-bounds vector access the first loop was written to guard against — a triangle with an out-of-range index skipped by loop 1 would hit undefined behavior in loop 2. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/mesh.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 7cf673c..8c8635b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -95,11 +95,6 @@ TriMesh::TriMesh(const pybind11::array_t &vertices, _mesh.add_face(vertex_indices[v0], vertex_indices[v1], vertex_indices[v2]); } - for (ssize_t i = 0; i < tris.shape(0); ++i) - { - _mesh.add_face(vertex_indices[tris(i, 0)], vertex_indices[tris(i, 1)], - vertex_indices[tris(i, 2)]); - } if (LoopCGAL::verbose) { std::cout << "Loaded mesh with " << _mesh.number_of_vertices() From 3378aa9b46332ec0af37e603c333d170a374e4d2 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 6 Jul 2026 09:43:26 +0930 Subject: [PATCH 03/12] fix: compact mesh before cut_with_implicit_function to avoid OOB vertex_properties and the caller's property array are indexed by vertex.idx(). When the mesh carries tombstoned vertices from a prior clip, idx() values can exceed number_of_vertices(), so indexing either array reads/writes out of bounds. collect_garbage() at entry renumbers live vertices 0..n-1, keeping all idx() values in range and making the property-size check meaningful. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/mesh.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/mesh.cpp b/src/mesh.cpp index 8c8635b..bcf859b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -728,6 +728,13 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do std::cout << "Mesh has " << _mesh.number_of_vertices() << " vertices and " << _mesh.number_of_faces() << " faces." << std::endl; std::cout << "Property size: " << property.size() << std::endl; + // Compact first: vertex_properties/property are indexed by vertex.idx(), and a + // mesh carrying tombstoned vertices (garbage from a prior clip) has idx() values + // that exceed number_of_vertices() — indexing the property arrays with them + // would read/write out of bounds. After collect_garbage() live vertices are + // renumbered 0..n-1 contiguously, so idx() is always in range. + if (_mesh.has_garbage()) + _mesh.collect_garbage(); if (property.size() != _mesh.number_of_vertices()) { std::cerr << "Error: Property size does not match number of vertices." << std::endl; From 40cfc8264d6d67c4815fdcaf62b1528219e3ae6f Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 6 Jul 2026 09:43:48 +0930 Subject: [PATCH 04/12] fix: gate cut_with_implicit_function debug prints on verbose The three status prints at the top of the function wrote to stdout unconditionally, unlike every other diagnostic in the file. Wrap them in the LoopCGAL::verbose guard. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/mesh.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index bcf859b..219c88a 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -724,10 +724,13 @@ NumpyMesh TriMesh::save(double area_threshold, void TriMesh::cut_with_implicit_function(const std::vector &property, double value, ImplicitCutMode cutmode, double snap_tol) { - std::cout << "Cutting mesh with implicit function at value " << value << std::endl; - std::cout << "Mesh has " << _mesh.number_of_vertices() << " vertices and " - << _mesh.number_of_faces() << " faces." << std::endl; - std::cout << "Property size: " << property.size() << std::endl; + if (LoopCGAL::verbose) + { + std::cout << "Cutting mesh with implicit function at value " << value << std::endl; + std::cout << "Mesh has " << _mesh.number_of_vertices() << " vertices and " + << _mesh.number_of_faces() << " faces." << std::endl; + std::cout << "Property size: " << property.size() << std::endl; + } // Compact first: vertex_properties/property are indexed by vertex.idx(), and a // mesh carrying tombstoned vertices (garbage from a prior clip) has idx() values // that exceed number_of_vertices() — indexing the property arrays with them From 221be62a1ec1de6d009346bfa29e32de2103e74f Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 6 Jul 2026 09:44:15 +0930 Subject: [PATCH 05/12] fix: validate stream state and face indices in read_from_file MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Face indices read from disk were used as verts[i0] with no bounds check, and stream read failures (truncated file) went unchecked — a corrupt or truncated LCMESH file caused out-of-bounds access instead of a clean error. Check the stream after each read and validate every face index against the vertex count, throwing runtime_error on failure. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/mesh.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 219c88a..6518201 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1178,12 +1178,14 @@ TriMesh TriMesh::read_from_file(const std::string& path) char magic[6]; in.read(magic, 6); - if (std::string(magic, 6) != "LCMESH") + if (!in || std::string(magic, 6) != "LCMESH") throw std::runtime_error("TriMesh::read_from_file — bad magic in: " + path); uint32_t nv, nf; in.read(reinterpret_cast(&nv), 4); in.read(reinterpret_cast(&nf), 4); + if (!in) + throw std::runtime_error("TriMesh::read_from_file — truncated header in: " + path); TriangleMesh mesh; std::vector verts(nv); @@ -1192,6 +1194,8 @@ TriMesh TriMesh::read_from_file(const std::string& path) in.read(reinterpret_cast(&x), 8); in.read(reinterpret_cast(&y), 8); in.read(reinterpret_cast(&z), 8); + if (!in) + throw std::runtime_error("TriMesh::read_from_file — truncated vertex data in: " + path); verts[i] = mesh.add_vertex(Point(x, y, z)); } @@ -1200,6 +1204,10 @@ TriMesh TriMesh::read_from_file(const std::string& path) in.read(reinterpret_cast(&i0), 4); in.read(reinterpret_cast(&i1), 4); in.read(reinterpret_cast(&i2), 4); + if (!in) + throw std::runtime_error("TriMesh::read_from_file — truncated face data in: " + path); + if (i0 >= nv || i1 >= nv || i2 >= nv) + throw std::runtime_error("TriMesh::read_from_file — face index out of range in: " + path); mesh.add_face(verts[i0], verts[i1], verts[i2]); } From f5f933deff1e5b01238adc3030234b62db1edcbe Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 6 Jul 2026 09:44:58 +0930 Subject: [PATCH 06/12] docs: correct corefine return-value comment and document clip-split limitation - corefine's header comment claimed it returns "intersection polyline segments inserted", but it returns the number of vertices added to this mesh (matching the pybind docstring). Fix the header to match the implementation. - Note in clipWithPlane that a constraint the plane cuts through is not recovered, since the plane overload takes no edge_is_constrained_map and the split breaks endpoint adjacency. - Fix the include-guard end comment (MESH_HANDLER_H -> MESH_H). Co-Authored-By: Claude Opus 4.8 (1M context) --- src/mesh.cpp | 6 ++++++ src/mesh.h | 4 ++-- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 6518201..fbb7769 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -712,6 +712,12 @@ int TriMesh::clipWithPlane(double a, double b, double c, double d, bool use_exac // survived, re-resolved by coordinate. Without this, stale Edge_index values // in _fixedEdges would mark unrelated (recycled) edges as protected during a // later remesh(). + // + // Limitation: a constraint that the plane cuts *through* is not recovered. + // rebuild_fixed_edges_from_coords only re-resolves a constraint when its two + // original endpoints are still directly connected; the plane overload of + // PMP::clip takes no edge_is_constrained_map, so a split introduces a midpoint + // that breaks that adjacency and the sub-edges are left unconstrained. rebuild_fixed_edges_from_coords(saved_constraints); return faces_before - faces_after; } diff --git a/src/mesh.h b/src/mesh.h index b7e3c32..5c11084 100644 --- a/src/mesh.h +++ b/src/mesh.h @@ -53,7 +53,7 @@ class TriMesh // coincident vertices and edges along the intersection curve, so patches // taken from each (e.g. a contact cut by a fault and the matching fault // patch) stitch together watertight. Both meshes are mutated in place. - // Returns the number of intersection polyline segments inserted. + // Returns the number of vertices added to this mesh by the corefinement. int corefine(TriMesh &other, bool use_exact_kernel = true); // Method to remesh the triangle mesh @@ -105,4 +105,4 @@ class TriMesh _edge_is_constrained_map; }; -#endif // MESH_HANDLER_H \ No newline at end of file +#endif // MESH_H \ No newline at end of file From 85dd1832af3be291d84faafd909e8b7484f5ab83 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 6 Jul 2026 09:45:53 +0930 Subject: [PATCH 07/12] refactor: collapse redundant nested try/catch in cutWithSurface The clip was wrapped in two nested try blocks both catching std::exception; the outer only guarded non-throwing flag checks and logging, so its catch was unreachable. Flatten to a single try around the clip, with the success/failure reporting after it. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/mesh.cpp | 92 +++++++++++++++++++++++----------------------------- 1 file changed, 41 insertions(+), 51 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index fbb7769..df51646 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -520,61 +520,41 @@ int TriMesh::cutWithSurface(TriMesh &clipper, std::cout << "Clipping tm with clipper." << std::endl; } + bool flag = false; try { - bool flag = false; - try - { - if (use_exact_kernel){ - // The exact round trip rebuilds _mesh from scratch (convert_to_double_mesh), - // so every Edge_index in _fixedEdges is invalidated and the constrained - // map cannot be threaded through the exact clip. Capture constraints by - // geometry beforehand and re-resolve them against the rebuilt mesh after. - auto saved_constraints = snapshot_constraint_coords(); - Exact_Mesh exact_clipper = convert_to_exact(scaled_clipper); - Exact_Mesh exact_mesh = convert_to_exact(*this); - flag = PMP::clip(exact_mesh, exact_clipper, CGAL::parameters::clip_volume(false)); - set_mesh(convert_to_double_mesh(exact_mesh)); - if (flag) - rebuild_fixed_edges_from_coords(saved_constraints); - } - else{ - // Pass the constrained-edge map for both meshes: CGAL reads existing - // constraints on input and, on output, marks the intersection edges and - // any surviving/split constraint edges. After the clip, _fixedEdges holds - // valid (pre-garbage-collection) indices, so snapshot it by geometry, - // collect garbage, then re-resolve — Edge_index values do not survive - // collect_garbage(). - flag = PMP::clip( - _mesh, scaled_clipper._mesh, - CGAL::parameters::edge_is_constrained_map(_edge_is_constrained_map) - .clip_volume(false), - CGAL::parameters::edge_is_constrained_map( - scaled_clipper._edge_is_constrained_map)); - if (flag) - { - auto saved_constraints = snapshot_constraint_coords(); - if (_mesh.has_garbage()) - _mesh.collect_garbage(); - rebuild_fixed_edges_from_coords(saved_constraints); - } - } - } - catch (const std::exception &e) - { - std::cerr << "Corefinement failed: " << e.what() << std::endl; - } - if (!flag) - { - std::cerr << "Warning: Clipping operation failed." << std::endl; + if (use_exact_kernel){ + // The exact round trip rebuilds _mesh from scratch (convert_to_double_mesh), + // so every Edge_index in _fixedEdges is invalidated and the constrained + // map cannot be threaded through the exact clip. Capture constraints by + // geometry beforehand and re-resolve them against the rebuilt mesh after. + auto saved_constraints = snapshot_constraint_coords(); + Exact_Mesh exact_clipper = convert_to_exact(scaled_clipper); + Exact_Mesh exact_mesh = convert_to_exact(*this); + flag = PMP::clip(exact_mesh, exact_clipper, CGAL::parameters::clip_volume(false)); + set_mesh(convert_to_double_mesh(exact_mesh)); + if (flag) + rebuild_fixed_edges_from_coords(saved_constraints); } - else - { - if (LoopCGAL::verbose) + else{ + // Pass the constrained-edge map for both meshes: CGAL reads existing + // constraints on input and, on output, marks the intersection edges and + // any surviving/split constraint edges. After the clip, _fixedEdges holds + // valid (pre-garbage-collection) indices, so snapshot it by geometry, + // collect garbage, then re-resolve — Edge_index values do not survive + // collect_garbage(). + flag = PMP::clip( + _mesh, scaled_clipper._mesh, + CGAL::parameters::edge_is_constrained_map(_edge_is_constrained_map) + .clip_volume(false), + CGAL::parameters::edge_is_constrained_map( + scaled_clipper._edge_is_constrained_map)); + if (flag) { - std::cout << "Clipping successful. Result has " - << _mesh.number_of_vertices() << " vertices and " - << _mesh.number_of_faces() << " faces." << std::endl; + auto saved_constraints = snapshot_constraint_coords(); + if (_mesh.has_garbage()) + _mesh.collect_garbage(); + rebuild_fixed_edges_from_coords(saved_constraints); } } } @@ -582,6 +562,16 @@ int TriMesh::cutWithSurface(TriMesh &clipper, { std::cerr << "Error during clipping: " << e.what() << std::endl; } + if (!flag) + { + std::cerr << "Warning: Clipping operation failed." << std::endl; + } + else if (LoopCGAL::verbose) + { + std::cout << "Clipping successful. Result has " + << _mesh.number_of_vertices() << " vertices and " + << _mesh.number_of_faces() << " faces." << std::endl; + } } else { From 078148694ae7ed61d06bd4217b575b9c7e1b8f15 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 6 Jul 2026 09:46:18 +0930 Subject: [PATCH 08/12] cleanup: consistent vertex naming in add_fixed_edges v0/v1 were bound to columns 1/0 for the validity check while the halfedge lookup used columns 0/1 directly, so the error messages printed mismatched vertex pairs. Bind v0=col0, v1=col1 once and reuse for the lookup; the resolved halfedge (columns 0->1) is unchanged. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/mesh.cpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index df51646..355d5bd 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -215,20 +215,18 @@ void TriMesh::add_fixed_edges(const pybind11::array_t &pairs) for (ssize_t i = 0; i < pairs_buf.shape(0); ++i) { - TriangleMesh::Vertex_index v0 = TriangleMesh::Vertex_index(pairs_buf(i, 1)); - TriangleMesh::Vertex_index v1 = TriangleMesh::Vertex_index(pairs_buf(i, 0)); + TriangleMesh::Vertex_index v0 = TriangleMesh::Vertex_index(pairs_buf(i, 0)); + TriangleMesh::Vertex_index v1 = TriangleMesh::Vertex_index(pairs_buf(i, 1)); if (!_mesh.is_valid(v0) || !_mesh.is_valid(v1)) { std::cerr << "Invalid vertex indices: (" << v0 << ", " << v1 << ")" << std::endl; continue; // Skip invalid vertex pairs } - TriangleMesh::Halfedge_index edge = - _mesh.halfedge(TriangleMesh::Vertex_index(pairs_buf(i, 0)), - TriangleMesh::Vertex_index(pairs_buf(i, 1))); + TriangleMesh::Halfedge_index edge = _mesh.halfedge(v0, v1); if (edge == TriangleMesh::null_halfedge()) { - std::cerr << "Half-edge is null for vertices (" << v1 << ", " << v0 << ")" + std::cerr << "Half-edge is null for vertices (" << v0 << ", " << v1 << ")" << std::endl; continue; } From b627380be2297ceac4deb86cab968f06aac15fdf Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 6 Jul 2026 10:47:08 +0930 Subject: [PATCH 09/12] test: harden implicit-cut coverage, add corefine tests, cover recent fixes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 1 of the improvement plan — raise the assertion floor and cover the newest features. - Replace the smoke test for cut_with_implicit_function (which only asserted n_verts/n_faces >= 0 — true for any mesh) with tests/test_implicit_cut.py that pin real geometry: seam inserted exactly on the isovalue, no degenerate faces, PRESERVE_INTERSECTION keeps full area, NaN triangles kept untouched, snap_tol reuses vertices. Three tests assert the *tight* cut and are marked xfail (strict) to document a discovered bug: KEEP_POSITIVE/NEGATIVE_SIDE retain a one-triangle fringe of wrong-side sub-triangles at the seam. Fixed in Phase 3. - Add tests/test_corefine.py: the exact-kernel path inserts a coincident shared polyline into both meshes (watertight) and preserves area. The use_exact_kernel=False path hard-crashes (SIGBUS) on open surfaces, so it is documented and skipped rather than executed (a native crash cannot be xfailed). - Add a truncated-file regression test for read_from_file and a raw-constructor test proving out-of-range/degenerate triangles are skipped without crashing and no face is double-inserted. Co-Authored-By: Claude Opus 4.8 (1M context) --- tests/test_corefine.py | 113 +++++++++++++++++++++++ tests/test_degenerate.py | 25 ++++- tests/test_file_io.py | 24 +++++ tests/test_implicit_cut.py | 169 ++++++++++++++++++++++++++++++++++ tests/test_mesh_operations.py | 16 +--- 5 files changed, 333 insertions(+), 14 deletions(-) create mode 100644 tests/test_corefine.py create mode 100644 tests/test_implicit_cut.py diff --git a/tests/test_corefine.py b/tests/test_corefine.py new file mode 100644 index 0000000..001da35 --- /dev/null +++ b/tests/test_corefine.py @@ -0,0 +1,113 @@ +"""Tests for ``TriMesh.corefine``. + +Corefinement inserts the shared intersection polyline into *both* meshes so that +patches taken from each side (e.g. a contact cut by a fault and the matching +fault patch) stitch together watertight. The contract we pin here is exactly +that: after corefining two surfaces that cross along a line, both meshes carry +the *same* set of vertices along that line. + +Two crossing grids are used: ``A`` lies in the z=0 plane, ``B`` in the y=0 +plane, so their intersection is the segment y=0, z=0 parametrised by x. + +Note: only the default exact-kernel path is exercised. ``use_exact_kernel=False`` +runs CGAL corefinement on the inexact kernel, which hard-crashes (SIGBUS) on +these open surfaces — see ``test_corefine_inexact_kernel_is_unsafe`` below. A +native crash cannot be trapped by ``xfail`` (it takes down the whole pytest +process), so that path is documented and skipped rather than executed. +""" +from __future__ import annotations + +import numpy as np +import pytest + +import loop_cgal + +EXT = 1.0 +N = 6 + + +def _grid_plane(fixed_axis: str): + """A triangulated grid in [-EXT, EXT]^2, embedded in the plane fixed_axis==0.""" + a = np.linspace(-EXT, EXT, N) + b = np.linspace(-EXT, EXT, N) + if fixed_axis == "z": + verts = np.array([[x, y, 0.0] for y in b for x in a], dtype=np.float64) + elif fixed_axis == "y": + verts = np.array([[x, 0.0, z] for z in b for x in a], dtype=np.float64) + else: + raise ValueError(fixed_axis) + + def vid(i, j): + return j * N + i + + tris = [] + for j in range(N - 1): + for i in range(N - 1): + tris += [ + [vid(i, j), vid(i + 1, j), vid(i + 1, j + 1)], + [vid(i, j), vid(i + 1, j + 1), vid(i, j + 1)], + ] + return verts, np.array(tris, dtype=np.int32) + + +def _seam_x(points): + """x-coordinates of vertices lying on the intersection line y=0, z=0.""" + on = (np.abs(points[:, 1]) < 1e-9) & (np.abs(points[:, 2]) < 1e-9) + return set(np.round(points[on, 0], 9)) + + +def _crossing_pair(): + A = loop_cgal.TriMesh.from_vertices_and_triangles(*_grid_plane("z")) + B = loop_cgal.TriMesh.from_vertices_and_triangles(*_grid_plane("y")) + return A, B + + +def test_corefine_adds_vertices_and_returns_the_count(): + """corefine returns the number of vertices it added to this mesh.""" + A, B = _crossing_pair() + before = A.n_points + added = A.corefine(B) + assert added == A.n_points - before + assert added > 0, "crossing surfaces must gain seam vertices" + + +def test_corefine_inserts_coincident_polyline_into_both_meshes(): + """The shared intersection polyline must be identical on both meshes.""" + A, B = _crossing_pair() + A.corefine(B) + pa, _ = A.get_vertices_and_triangles() + pb, _ = B.get_vertices_and_triangles() + seam_a = _seam_x(pa) + seam_b = _seam_x(pb) + assert len(seam_a) > 2, "expected an inserted polyline, not just the original corners" + assert seam_a == seam_b, ( + f"seam vertices differ between meshes — not watertight:\n" + f" A: {sorted(seam_a)}\n B: {sorted(seam_b)}" + ) + + +def test_corefine_mutates_both_meshes(): + """Both meshes are corefined in place (the second argument is not read-only).""" + A, B = _crossing_pair() + b_before = B.n_points + A.corefine(B) + assert B.n_points > b_before, "the other mesh must also receive the shared polyline" + + +def test_corefine_preserves_area(): + """Corefinement only subdivides — it must not change either surface's area.""" + A, B = _crossing_pair() + area_a, area_b = A.area, B.area + A.corefine(B) + assert A.area == pytest.approx(area_a, rel=1e-6) + assert B.area == pytest.approx(area_b, rel=1e-6) + + +@pytest.mark.skip( + reason="corefine(use_exact_kernel=False) hard-crashes (SIGBUS) on open " + "intersecting surfaces; a native crash would take down the pytest process. " + "Tracked separately — the inexact path needs input validation or removal." +) +def test_corefine_inexact_kernel_is_unsafe(): + A, B = _crossing_pair() + A.corefine(B, use_exact_kernel=False) diff --git a/tests/test_degenerate.py b/tests/test_degenerate.py index f263c5d..c6fc199 100644 --- a/tests/test_degenerate.py +++ b/tests/test_degenerate.py @@ -27,4 +27,27 @@ def test_invalid_index_raises_value_error(): def test_degenerate_triangle_raises_value_error(): surface = make_polydata_degenerate_triangle() with pytest.raises(ValueError, match="degenerate triangles"): - _ = loop_cgal.TriMesh(surface) \ No newline at end of file + _ = loop_cgal.TriMesh(surface) + + +def test_raw_constructor_skips_bad_triangles_without_crashing(): + """The C++ array constructor (bypassing Python validation) must skip + out-of-range and degenerate triangles gracefully. + + Regression: the constructor looped over the triangle array twice — the + second loop re-added every face with no validation, dereferencing the + vertex-index vector out of bounds for any triangle the first loop skipped. + """ + from loop_cgal._loop_cgal import TriMesh as RawTriMesh + + verts = np.array( + [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]], dtype=np.float64 + ) + # one valid triangle, one with an out-of-range index, one degenerate + tris = np.array([[0, 1, 2], [0, 2, 99], [1, 1, 3]], dtype=np.int32) + + m = RawTriMesh(verts, tris) # must not crash + assert m._cgal_n_vertices() == 4 + # only the single valid triangle survives; the bad ones are skipped and, + # crucially, no face is inserted twice. + assert m._cgal_n_faces() == 1 \ No newline at end of file diff --git a/tests/test_file_io.py b/tests/test_file_io.py index 945ee00..9f17de8 100644 --- a/tests/test_file_io.py +++ b/tests/test_file_io.py @@ -118,6 +118,30 @@ def test_read_missing_file(tmp_path): TriMesh.read_from_file(str(tmp_path / "nonexistent.lcm")) +def test_read_truncated_raises(tmp_path, unit_mesh): + """A file with a valid header but truncated body must raise, not read OOB. + + Regression: face indices and vertex data were consumed with no stream-state + or bounds check, so a truncated LCMESH file caused out-of-bounds access + instead of a clean error. + """ + good = tmp_path / "good.lcm" + unit_mesh.write_to_file(str(good)) + raw = good.read_bytes() + + # Keep magic + the two 4-byte counts (14 bytes) but drop all vertex/face data. + truncated = tmp_path / "trunc.lcm" + truncated.write_bytes(raw[:14]) + with pytest.raises(Exception, match="truncated"): + TriMesh.read_from_file(str(truncated)) + + # Truncated partway through the face block must also be caught. + mid = tmp_path / "mid.lcm" + mid.write_bytes(raw[:-4]) + with pytest.raises(Exception, match="truncated"): + TriMesh.read_from_file(str(mid)) + + # --------------------------------------------------------------------------- # Loaded mesh is fully usable (no dangling state) # --------------------------------------------------------------------------- diff --git a/tests/test_implicit_cut.py b/tests/test_implicit_cut.py new file mode 100644 index 0000000..5b1dc2a --- /dev/null +++ b/tests/test_implicit_cut.py @@ -0,0 +1,169 @@ +"""Behavioural tests for ``cut_with_implicit_function``. + +These replace an earlier smoke test whose only assertions were ``n_verts >= 0`` +and ``n_faces >= 0`` — conditions that are true for *any* mesh and so could +never fail. The tests below pin the actual geometry of the isocontour cut: + +* the seam is inserted exactly on the isovalue, +* ``KEEP_POSITIVE_SIDE`` / ``KEEP_NEGATIVE_SIDE`` retain only their side + (no fringe of wrong-side triangles clinging to the seam), +* ``PRESERVE_INTERSECTION`` keeps the whole surface but still inserts the seam, +* NaN (off-extent) property values leave their triangles untouched, +* ``snap_tol`` reuses an existing vertex instead of inserting a near-duplicate. + +The property is chosen to equal the vertex x-coordinate, so the isocontour +``property == value`` is the plane ``x == value`` and every geometric claim can +be checked directly against vertex coordinates. +""" +from __future__ import annotations + +import numpy as np +import pytest + +import loop_cgal +from loop_cgal._loop_cgal import ImplicitCutMode + +SIZE = 10.0 +EPS = 1e-6 + + +def _grid(nx: int = 11, ny: int = 11, size: float = SIZE): + """A flat (z=0) triangulated grid spanning [0, size] in x and y.""" + verts = np.array( + [[size * i / (nx - 1), size * j / (ny - 1), 0.0] + for j in range(ny) for i in range(nx)], + dtype=np.float64, + ) + + def vid(i, j): + return j * nx + i + + tris = [] + for j in range(ny - 1): + for i in range(nx - 1): + tris += [ + [vid(i, j), vid(i + 1, j), vid(i + 1, j + 1)], + [vid(i, j), vid(i + 1, j + 1), vid(i, j + 1)], + ] + return verts, np.array(tris, dtype=np.int32) + + +def _cut(prop, value, mode, snap_tol=1e-4): + """Run the cut on a fresh grid whose property equals its x-coordinate.""" + verts, tris = _grid() + tm = loop_cgal.TriMesh.from_vertices_and_triangles(verts, tris) + if prop == "x": + prop = verts[:, 0].tolist() + tm.cut_with_implicit_function(prop, value, mode, snap_tol) + v, f = tm.get_vertices_and_triangles() + return tm, v, f + + +def _used_x(v, f): + """x-coordinates of only the vertices actually referenced by a face. + + ``get_vertices_and_triangles`` also returns vertices left isolated by the + cut, so filtering to referenced vertices is what tells us where the surface + actually reaches. + """ + return v[np.unique(f), 0] + + +def test_seam_inserted_exactly_on_isovalue(): + """A cut between grid lines inserts seam vertices lying exactly on the plane x==value.""" + _, v, f = _cut("x", 4.5, ImplicitCutMode.KEEP_POSITIVE_SIDE) + xs = _used_x(v, f) + assert np.any(np.isclose(xs, 4.5, atol=EPS)), "no seam vertices were inserted on the isovalue" + + +def test_no_degenerate_faces(): + """The cut must never emit a face with a repeated vertex index.""" + _, _, f = _cut("x", 4.5, ImplicitCutMode.KEEP_POSITIVE_SIDE) + assert len(f) > 0 + for tri in f: + assert tri[0] != tri[1] and tri[1] != tri[2] and tri[0] != tri[2], ( + f"degenerate face emitted: {tri}" + ) + + +def test_preserve_intersection_keeps_full_area_and_inserts_seam(): + """PRESERVE_INTERSECTION splits along the seam but keeps the whole surface.""" + tm, v, f = _cut("x", 4.5, ImplicitCutMode.PRESERVE_INTERSECTION) + assert tm.area == pytest.approx(SIZE * SIZE, rel=1e-6), "PRESERVE_INTERSECTION must keep the full area" + xs = _used_x(v, f) + assert np.any(np.isclose(xs, 4.5, atol=EPS)), "seam must still be inserted" + assert xs.min() < 4.5 < xs.max(), "both sides must be retained" + + +def test_nan_triangles_are_kept_untouched(): + """Triangles with a NaN (off-extent) property survive regardless of side/value. + + The far-negative region (x < 3, marked NaN) must remain even under + KEEP_POSITIVE_SIDE with a value well above it. + """ + verts, tris = _grid() + prop = verts[:, 0].copy() + prop[verts[:, 0] < 3.0] = np.nan + tm = loop_cgal.TriMesh.from_vertices_and_triangles(verts, tris) + tm.cut_with_implicit_function(prop.tolist(), 4.5, ImplicitCutMode.KEEP_POSITIVE_SIDE) + v, f = tm.get_vertices_and_triangles() + xs = _used_x(v, f) + assert xs.min() < EPS, "NaN-property region (x≈0) was dropped but should be kept untouched" + + +def test_snap_tol_reuses_nearby_vertex(): + """A value within snap_tol of a grid line inserts no near-duplicate vertices.""" + verts, tris = _grid() + tm = loop_cgal.TriMesh.from_vertices_and_triangles(verts, tris) + n_before = tm.n_points + # 4.0 is a grid line; 4.0 + 1e-6 is far inside snap_tol (1e-4 of a unit edge). + tm.cut_with_implicit_function( + verts[:, 0].tolist(), 4.0 + 1e-6, ImplicitCutMode.KEEP_POSITIVE_SIDE, 1e-4 + ) + assert tm.n_points == n_before, ( + f"snap_tol should reuse existing vertices, but point count grew " + f"{n_before} -> {tm.n_points}" + ) + + +@pytest.mark.xfail( + reason="known bug: wrong-side seam sub-triangles are retained; fixed by the " + "data-driven case-table refactor", + strict=True, +) +def test_keep_positive_side_is_tight(): + """KEEP_POSITIVE_SIDE must retain exactly the region x >= value.""" + tm, v, f = _cut("x", 4.5, ImplicitCutMode.KEEP_POSITIVE_SIDE) + xs = _used_x(v, f) + assert xs.min() >= 4.5 - EPS, ( + f"kept a wrong-side fringe: min referenced x = {xs.min():.4f} (< 4.5)" + ) + assert tm.area == pytest.approx((SIZE - 4.5) * SIZE, rel=1e-3) + + +@pytest.mark.xfail( + reason="known bug: wrong-side seam sub-triangles are retained; fixed by the " + "data-driven case-table refactor", + strict=True, +) +def test_keep_negative_side_is_tight(): + """KEEP_NEGATIVE_SIDE must retain exactly the region x <= value.""" + tm, v, f = _cut("x", 4.5, ImplicitCutMode.KEEP_NEGATIVE_SIDE) + xs = _used_x(v, f) + assert xs.max() <= 4.5 + EPS, ( + f"kept a wrong-side fringe: max referenced x = {xs.max():.4f} (> 4.5)" + ) + assert tm.area == pytest.approx(4.5 * SIZE, rel=1e-3) + + +@pytest.mark.xfail( + reason="known bug: wrong-side seam sub-triangles are retained; fixed by the " + "data-driven case-table refactor", + strict=True, +) +def test_sides_partition_the_original(): + """Positive and negative sides must partition the original area (they only + share the zero-area seam), so their areas sum to the whole.""" + pos, _, _ = _cut("x", 4.5, ImplicitCutMode.KEEP_POSITIVE_SIDE) + neg, _, _ = _cut("x", 4.5, ImplicitCutMode.KEEP_NEGATIVE_SIDE) + assert pos.area + neg.area == pytest.approx(SIZE * SIZE, rel=1e-3) diff --git a/tests/test_mesh_operations.py b/tests/test_mesh_operations.py index 922bf6c..bc485dd 100644 --- a/tests/test_mesh_operations.py +++ b/tests/test_mesh_operations.py @@ -406,16 +406,6 @@ def _add_collocated_seam(pv_mesh: pv.PolyData) -> pv.PolyData: ) -def test_cut_with_implicit_function(square_surface): - tm = loop_cgal.TriMesh(square_surface) - # create a scalar property that varies across vertices - saved = tm.save() - nverts = np.array(saved.vertices).shape[0] - prop = [float(i) / max(1, (nverts - 1)) for i in range(nverts)] - # cut at 0.5 keeping positive side - tm.cut_with_implicit_function(prop, 0.5, ImplicitCutMode.KEEP_POSITIVE_SIDE) - res = tm.save() - v = np.array(res.vertices).shape[0] - f = np.array(res.triangles).shape[0] - assert v >= 0 - assert f >= 0 +# cut_with_implicit_function is covered thoroughly in tests/test_implicit_cut.py +# (the previous smoke test here only asserted n_verts/n_faces >= 0, which is true +# for any mesh and so could never fail). From 4ccdb15d9bebca9d9313a9876b612583bd972a73 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 6 Jul 2026 10:55:11 +0930 Subject: [PATCH 10/12] refactor: data-driven marching-triangles table + fix wrong-side fringe MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 3 of the improvement plan. Refactor: the six hand-written split cases in cut_with_implicit_function were sign-paired duplicates — CASE 1≡5, 2≡6, 3≡7 produce identical oriented sub-triangles because the split topology depends only on *which* vertex is the "lone" one (the vertex whose side differs), not on its sign. Collapse them into one table keyed on the lone vertex (~130 lines -> ~30), and drop the dead all-positive skip (already covered by the `active` mask) and per-case verbose prints. Fix: KEEP_POSITIVE_SIDE / KEEP_NEGATIVE_SIDE kept a one-triangle fringe past the isovalue. When a straddling triangle is split, the wrong-side sub-triangle has two crossing vertices whose value is exactly `value`; the assembly filter's strict < / > comparison spared it. Use inclusive <= / >= so it is dropped and the cut stops exactly at the seam. Verified: for property=x on a 10x10 grid, POS+NEG areas now sum to the original and each matches the analytic half-area. Removes the strict-xfail markers from the three tightness tests added in Phase 1, which now pass. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/mesh.cpp | 176 +++++++++---------------------------- tests/test_implicit_cut.py | 15 ---- 2 files changed, 42 insertions(+), 149 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 355d5bd..33dd4d3 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -846,9 +846,6 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do continue; } auto tri = tri_array[t]; - // if all > value skip (hanging_wall in python) - if (vertex_properties[tri[0]] > value && vertex_properties[tri[1]] > value && vertex_properties[tri[2]] > value) - continue; // for each edge of tri, check if edge crosses for (auto eid : tri2edge[t]) { @@ -898,138 +895,43 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do new_point_on_edge[eid] = newverts.size() - 1; } - double v1 = vertex_properties[tri[0]]; - double v2 = vertex_properties[tri[1]]; - double v3 = vertex_properties[tri[2]]; - // replicate python cases - // convert tri to vector of 3 original indices and 2 new points - std::array extended = {tri[0], tri[1], tri[2], 0, 0}; - // retrieve relevant edges indices + // Edge ids of the three triangle sides, and the crossing vertex inserted on + // each (SIZE_MAX where that side does not cross the isovalue). std::size_t e01 = edge_index_map[std::make_pair(std::min(tri[0], tri[1]), std::max(tri[0], tri[1]))]; std::size_t e12 = edge_index_map[std::make_pair(std::min(tri[1], tri[2]), std::max(tri[1], tri[2]))]; std::size_t e20 = edge_index_map[std::make_pair(std::min(tri[2], tri[0]), std::max(tri[2], tri[0]))]; - // Get new points where available - std::size_t np_e01 = new_point_on_edge.count(e01) ? new_point_on_edge[e01] : SIZE_MAX; - std::size_t np_e12 = new_point_on_edge.count(e12) ? new_point_on_edge[e12] : SIZE_MAX; - std::size_t np_e20 = new_point_on_edge.count(e20) ? new_point_on_edge[e20] : SIZE_MAX; - // Helper to append triangle - auto append_tri = [&](std::array tarr) - { newtris.push_back(tarr); }; - - // CASE 1: v1 > value and v2 > value and v3 value && v2 > value && v3 < value) - { - std::size_t p1 = np_e12; - std::size_t p2 = np_e20; - extended[3] = p1; - extended[4] = p2; - std::array m1 = {extended[0], extended[1], extended[3]}; - std::array m2 = {extended[0], extended[3], extended[4]}; - std::array m3 = {extended[4], extended[3], extended[2]}; - newtris[t] = m1; - append_tri(m2); - append_tri(m3); - if (LoopCGAL::verbose) - { - std::cout << "CASE 1 executed" << std::endl; - } - continue; - } - // CASE 2 - if (v1 > value && v2 < value && v3 > value) - { - std::size_t p1 = np_e01; - std::size_t p2 = np_e12; - extended[3] = p1; - extended[4] = p2; - std::array m1 = {extended[0], extended[3], extended[2]}; - std::array m2 = {extended[3], extended[4], extended[2]}; - std::array m3 = {extended[3], extended[1], extended[4]}; - newtris[t] = m1; - append_tri(m2); - append_tri(m3); - if (LoopCGAL::verbose) - { - std::cout << "CASE 2 executed" << std::endl; - } - continue; - } - // CASE 3 - if (v1 < value && v2 > value && v3 > value) - { - std::size_t p1 = np_e01; - std::size_t p2 = np_e20; - extended[3] = p1; - extended[4] = p2; - std::array m1 = {extended[0], extended[3], extended[4]}; - std::array m2 = {extended[3], extended[1], extended[2]}; - std::array m3 = {extended[4], extended[3], extended[2]}; - newtris[t] = m1; - append_tri(m2); - append_tri(m3); - if (LoopCGAL::verbose) - { - std::cout << "CASE 3 executed" << std::endl; - } - continue; - } - // CASE 5 - if (v1 < value && v2 < value && v3 > value) - { - std::size_t p1 = np_e12; - std::size_t p2 = np_e20; - extended[3] = p1; - extended[4] = p2; - std::array m1 = {extended[0], extended[1], extended[3]}; - std::array m2 = {extended[0], extended[3], extended[4]}; - std::array m3 = {extended[4], extended[3], extended[2]}; - newtris[t] = m1; - append_tri(m2); - append_tri(m3); - if (LoopCGAL::verbose) - { - std::cout << "CASE 5 executed" << std::endl; - } - continue; - } - // CASE 6 - if (v1 < value && v2 > value && v3 < value) - { - std::size_t p1 = np_e01; - std::size_t p2 = np_e12; - extended[3] = p1; - extended[4] = p2; - std::array m1 = {extended[0], extended[3], extended[2]}; - std::array m2 = {extended[3], extended[4], extended[2]}; - std::array m3 = {extended[3], extended[1], extended[4]}; - newtris[t] = m1; - append_tri(m2); - append_tri(m3); - if (LoopCGAL::verbose) - { - std::cout << "CASE 6 executed" << std::endl; - } - continue; - } - // CASE 7 - if (v1 > value && v2 < value && v3 < value) - { - std::size_t p1 = np_e01; - std::size_t p2 = np_e20; - extended[3] = p1; - extended[4] = p2; - std::array m1 = {extended[0], extended[3], extended[4]}; - std::array m2 = {extended[3], extended[2], extended[4]}; - std::array m3 = {extended[3], extended[1], extended[2]}; - newtris[t] = m1; - append_tri(m2); - append_tri(m3); - if (LoopCGAL::verbose) - { - std::cout << "CASE 7 executed" << std::endl; - } - continue; - } + const std::size_t c01 = new_point_on_edge.count(e01) ? new_point_on_edge[e01] : SIZE_MAX; + const std::size_t c12 = new_point_on_edge.count(e12) ? new_point_on_edge[e12] : SIZE_MAX; + const std::size_t c20 = new_point_on_edge.count(e20) ? new_point_on_edge[e20] : SIZE_MAX; + + // Split the straddling triangle along the seam. Exactly one vertex (the + // "lone" vertex) lies on the opposite side of the isovalue from the other + // two — the nudge above guarantees no vertex is exactly on it — so the two + // triangle sides incident to the lone vertex are the ones that cross. The + // three sub-triangles produced depend only on *which* vertex is lone, not on + // its sign; the sign only decides, in the assembly filter below, which + // sub-triangles survive. This collapses the six near-identical hand-written + // cases (whose sign-paired variants produced identical oriented triangles) + // into one table keyed on the lone vertex. + const bool p0 = vertex_properties[tri[0]] > value; + const bool p1 = vertex_properties[tri[1]] > value; + const bool p2 = vertex_properties[tri[2]] > value; + int lone; + if (p1 == p2) lone = 0; // tri[0] differs from tri[1] == tri[2] + else if (p0 == p2) lone = 1; // tri[1] differs + else lone = 2; // tri[2] differs + + std::array, 3> sub; + if (lone == 0) + sub = {{{tri[0], c01, c20}, {c01, tri[1], tri[2]}, {c20, c01, tri[2]}}}; + else if (lone == 1) + sub = {{{tri[0], c01, tri[2]}, {c01, c12, tri[2]}, {c01, tri[1], c12}}}; + else + sub = {{{tri[0], tri[1], c12}, {tri[0], c12, c20}, {c20, c12, tri[2]}}}; + + newtris[t] = sub[0]; + newtris.push_back(sub[1]); + newtris.push_back(sub[2]); } // Build new CGAL mesh from newverts and newtris @@ -1043,12 +945,18 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do // skip degenerate if (tri[0] == tri[1] || tri[1] == tri[2] || tri[0] == tri[2]) continue; + // Drop sub-triangles on the unwanted side. Comparisons are inclusive of the + // seam: a wrong-side seam sub-triangle has two crossing vertices whose value + // is exactly `value` plus one vertex strictly on the wrong side, so a strict + // comparison would spare it and leave a one-triangle fringe past the + // isovalue. Using <=/>= drops it, so the cut stops exactly at the seam. + // (NaN values compare false either way, so off-extent triangles are kept.) if (ImplicitCutMode::KEEP_NEGATIVE_SIDE == cutmode) { double v0 = newvals[tri[0]]; double v1 = newvals[tri[1]]; double v2 = newvals[tri[2]]; - if (v0 > value && v1 > value && v2 > value) + if (v0 >= value && v1 >= value && v2 >= value) { continue; } @@ -1058,7 +966,7 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do double v0 = newvals[tri[0]]; double v1 = newvals[tri[1]]; double v2 = newvals[tri[2]]; - if (v0 < value && v1 < value && v2 < value) + if (v0 <= value && v1 <= value && v2 <= value) { continue; } diff --git a/tests/test_implicit_cut.py b/tests/test_implicit_cut.py index 5b1dc2a..57d2bfe 100644 --- a/tests/test_implicit_cut.py +++ b/tests/test_implicit_cut.py @@ -126,11 +126,6 @@ def test_snap_tol_reuses_nearby_vertex(): ) -@pytest.mark.xfail( - reason="known bug: wrong-side seam sub-triangles are retained; fixed by the " - "data-driven case-table refactor", - strict=True, -) def test_keep_positive_side_is_tight(): """KEEP_POSITIVE_SIDE must retain exactly the region x >= value.""" tm, v, f = _cut("x", 4.5, ImplicitCutMode.KEEP_POSITIVE_SIDE) @@ -141,11 +136,6 @@ def test_keep_positive_side_is_tight(): assert tm.area == pytest.approx((SIZE - 4.5) * SIZE, rel=1e-3) -@pytest.mark.xfail( - reason="known bug: wrong-side seam sub-triangles are retained; fixed by the " - "data-driven case-table refactor", - strict=True, -) def test_keep_negative_side_is_tight(): """KEEP_NEGATIVE_SIDE must retain exactly the region x <= value.""" tm, v, f = _cut("x", 4.5, ImplicitCutMode.KEEP_NEGATIVE_SIDE) @@ -156,11 +146,6 @@ def test_keep_negative_side_is_tight(): assert tm.area == pytest.approx(4.5 * SIZE, rel=1e-3) -@pytest.mark.xfail( - reason="known bug: wrong-side seam sub-triangles are retained; fixed by the " - "data-driven case-table refactor", - strict=True, -) def test_sides_partition_the_original(): """Positive and negative sides must partition the original area (they only share the zero-area seam), so their areas sum to the whole.""" From 5d2eea3416935c3423a833bdd277034e6b63f2de Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 6 Jul 2026 11:08:03 +0930 Subject: [PATCH 11/12] feat!: raise on genuine failure instead of returning 0 / crashing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase 2 of the improvement plan — adopt a clear error-handling contract. A return of 0 is now reserved for the one legitimate no-op (a cut whose meshes do not intersect); every genuine failure raises instead of printing to stderr and returning a sentinel. pybind11 maps invalid_argument -> ValueError and runtime_error -> RuntimeError. - cutWithSurface: raise ValueError on empty/invalid source or clipper; raise RuntimeError if the clip of intersecting meshes throws or reports failure. Non-intersecting meshes still return 0. - clipWithPlane: raise ValueError on empty mesh; raise RuntimeError if the clip throws. - cut_with_implicit_function: raise ValueError on property/vertex size mismatch. - corefine: refuse use_exact_kernel=False up front with ValueError. That path ran CGAL corefinement on the inexact predicate kernel (Simple_cartesian), which hard-crashes (SIGBUS) even on trivial valid inputs — uncatchable, so it must be prevented rather than trapped. corefine failures now raise RuntimeError. Tests: add tests/test_error_contract.py; convert the previously-skipped corefine inexact test into one asserting it raises ValueError. BREAKING CHANGE: operations that previously returned 0 on invalid/empty input or logged-and-continued on failure now raise ValueError/RuntimeError. Co-Authored-By: Claude Opus 4.8 (1M context) --- loop_cgal/bindings.cpp | 4 +- src/mesh.cpp | 101 +++++++++++++++++------------------ src/mesh.h | 3 ++ tests/test_corefine.py | 25 ++++----- tests/test_error_contract.py | 71 ++++++++++++++++++++++++ 5 files changed, 138 insertions(+), 66 deletions(-) create mode 100644 tests/test_error_contract.py diff --git a/loop_cgal/bindings.cpp b/loop_cgal/bindings.cpp index e9bbabd..4e6aa27 100644 --- a/loop_cgal/bindings.cpp +++ b/loop_cgal/bindings.cpp @@ -37,7 +37,9 @@ PYBIND11_MODULE(_loop_cgal, m) "polyline into BOTH meshes in place. After the call the two meshes " "carry coincident vertices/edges along the intersection curve, so " "patches taken from each side stitch together watertight. " - "Returns the number of vertices added to this mesh.") + "Returns the number of vertices added to this mesh. " + "Only the exact kernel is supported; use_exact_kernel=False raises " + "ValueError (inexact-kernel corefinement crashes).") .def("cut_with_surface", &TriMesh::cutWithSurface, py::arg("surface"), py::arg("preserve_intersection") = false, py::arg("preserve_intersection_clipper") = false, diff --git a/src/mesh.cpp b/src/mesh.cpp index 33dd4d3..3e25cf6 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -387,31 +387,22 @@ int TriMesh::cutWithSurface(TriMesh &clipper, std::cout << "Cutting mesh with surface." << std::endl; } - // Validate input meshes - if (!CGAL::is_valid_polygon_mesh(_mesh, LoopCGAL::verbose)) - { - std::cerr << "Error: Source mesh is invalid!" << std::endl; - return 0; - } - - if (!CGAL::is_valid_polygon_mesh(clipper._mesh, LoopCGAL::verbose)) - { - std::cerr << "Error: Clipper mesh is invalid!" << std::endl; - return 0; - } - + // Validate input meshes. A bad or empty input is a programming error, not a + // legitimate no-op, so raise rather than return 0 (which is reserved for the + // meshes-don't-intersect case below). pybind11 maps invalid_argument -> + // Python ValueError. if (_mesh.number_of_vertices() == 0 || _mesh.number_of_faces() == 0) - { - std::cerr << "Error: Source mesh is empty!" << std::endl; - return 0; - } + throw std::invalid_argument("cutWithSurface: source mesh is empty"); if (clipper._mesh.number_of_vertices() == 0 || clipper._mesh.number_of_faces() == 0) - { - std::cerr << "Error: Clipper mesh is empty!" << std::endl; - return 0; - } + throw std::invalid_argument("cutWithSurface: clipper mesh is empty"); + + if (!CGAL::is_valid_polygon_mesh(_mesh, LoopCGAL::verbose)) + throw std::invalid_argument("cutWithSurface: source mesh is not a valid polygon mesh"); + + if (!CGAL::is_valid_polygon_mesh(clipper._mesh, LoopCGAL::verbose)) + throw std::invalid_argument("cutWithSurface: clipper mesh is not a valid polygon mesh"); // Merge any collocated border vertices on the target before clipping. // Collocated vertices produce zero-area faces whose degenerate bounding boxes @@ -558,13 +549,16 @@ int TriMesh::cutWithSurface(TriMesh &clipper, } catch (const std::exception &e) { - std::cerr << "Error during clipping: " << e.what() << std::endl; + // The meshes intersect, so the clip was expected to succeed — a thrown + // exception here is a real failure, not a no-op. Rethrow with context + // (pybind11 maps runtime_error -> Python RuntimeError). + throw std::runtime_error(std::string("cutWithSurface: clip failed: ") + e.what()); } + // The meshes intersect but PMP::clip reported failure — surface it rather + // than silently returning 0 (which callers read as "no intersection"). if (!flag) - { - std::cerr << "Warning: Clipping operation failed." << std::endl; - } - else if (LoopCGAL::verbose) + throw std::runtime_error("cutWithSurface: clip of intersecting meshes failed"); + if (LoopCGAL::verbose) { std::cout << "Clipping successful. Result has " << _mesh.number_of_vertices() << " vertices and " @@ -573,6 +567,7 @@ int TriMesh::cutWithSurface(TriMesh &clipper, } else { + // Legitimate no-op: nothing to cut. Return 0 (see faces_before-faces_after). if (LoopCGAL::verbose) { std::cout << "Meshes do not intersect. No clipping performed." @@ -592,11 +587,21 @@ int TriMesh::cutWithSurface(TriMesh &clipper, int TriMesh::corefine(TriMesh &other, bool use_exact_kernel) { + // The inexact kernel path is unsafe by construction: TriangleMesh uses + // Simple_cartesian, whose predicates are inexact, and CGAL + // corefinement requires exact predicates to make consistent branching + // decisions. On Simple_cartesian it does not merely give a wrong result — it + // hard-crashes (SIGBUS) even on trivial valid inputs, which no try/catch can + // trap. Refuse it up front with a catchable error instead. The exact path + // converts to the EPECK Exact_Mesh, which is safe. + if (!use_exact_kernel) + throw std::invalid_argument( + "corefine requires the exact kernel; use_exact_kernel=False is unsafe " + "with the inexact predicate kernel and crashes. Call with " + "use_exact_kernel=True (the default)."); + if (_mesh.number_of_faces() == 0 || other._mesh.number_of_faces() == 0) - { - std::cerr << "Error: corefine called on an empty mesh." << std::endl; - return 0; - } + throw std::invalid_argument("corefine: called on an empty mesh"); // Merge collocated border vertices first — corefinement needs valid, // duplicate-free borders to compute a clean intersection polyline. @@ -609,23 +614,15 @@ int TriMesh::corefine(TriMesh &other, bool use_exact_kernel) try { - if (use_exact_kernel) - { - Exact_Mesh em1 = convert_to_exact(*this); - Exact_Mesh em2 = convert_to_exact(other); - PMP::corefine(em1, em2); - set_mesh(convert_to_double_mesh(em1)); - other.set_mesh(convert_to_double_mesh(em2)); - } - else - { - PMP::corefine(_mesh, other._mesh); - } + Exact_Mesh em1 = convert_to_exact(*this); + Exact_Mesh em2 = convert_to_exact(other); + PMP::corefine(em1, em2); + set_mesh(convert_to_double_mesh(em1)); + other.set_mesh(convert_to_double_mesh(em2)); } catch (const std::exception &e) { - std::cerr << "corefine failed: " << e.what() << std::endl; - return 0; + throw std::runtime_error(std::string("corefine failed: ") + e.what()); } if (_mesh.has_garbage()) @@ -649,10 +646,7 @@ int TriMesh::corefine(TriMesh &other, bool use_exact_kernel) int TriMesh::clipWithPlane(double a, double b, double c, double d, bool use_exact_kernel) { if (_mesh.number_of_vertices() == 0 || _mesh.number_of_faces() == 0) - { - std::cerr << "Error: Source mesh is empty!" << std::endl; - return 0; - } + throw std::invalid_argument("clipWithPlane: source mesh is empty"); const int faces_before = static_cast(_mesh.number_of_faces()); @@ -679,7 +673,8 @@ int TriMesh::clipWithPlane(double a, double b, double c, double d, bool use_exac } catch (const std::exception &e) { - std::cerr << "Plane clip failed: " << e.what() << std::endl; + // pybind11 maps runtime_error -> Python RuntimeError. + throw std::runtime_error(std::string("clipWithPlane: clip failed: ") + e.what()); } const int faces_after = static_cast(_mesh.number_of_faces()); @@ -733,10 +728,10 @@ void TriMesh::cut_with_implicit_function(const std::vector &property, do if (_mesh.has_garbage()) _mesh.collect_garbage(); if (property.size() != _mesh.number_of_vertices()) - { - std::cerr << "Error: Property size does not match number of vertices." << std::endl; - return; - } + throw std::invalid_argument( + "cut_with_implicit_function: property size (" + + std::to_string(property.size()) + ") does not match the number of " + "mesh vertices (" + std::to_string(_mesh.number_of_vertices()) + ")"); // Create a property map for vertex properties typedef boost::property_map::type VertexIndexMap; VertexIndexMap vim = get(boost::vertex_index, _mesh); diff --git a/src/mesh.h b/src/mesh.h index 5c11084..f67362c 100644 --- a/src/mesh.h +++ b/src/mesh.h @@ -54,6 +54,9 @@ class TriMesh // taken from each (e.g. a contact cut by a fault and the matching fault // patch) stitch together watertight. Both meshes are mutated in place. // Returns the number of vertices added to this mesh by the corefinement. + // Only the exact kernel is supported: use_exact_kernel=false throws + // std::invalid_argument, because corefinement on the inexact predicate + // kernel (Simple_cartesian) crashes rather than merely misbehaving. int corefine(TriMesh &other, bool use_exact_kernel = true); // Method to remesh the triangle mesh diff --git a/tests/test_corefine.py b/tests/test_corefine.py index 001da35..0afb806 100644 --- a/tests/test_corefine.py +++ b/tests/test_corefine.py @@ -9,11 +9,10 @@ Two crossing grids are used: ``A`` lies in the z=0 plane, ``B`` in the y=0 plane, so their intersection is the segment y=0, z=0 parametrised by x. -Note: only the default exact-kernel path is exercised. ``use_exact_kernel=False`` -runs CGAL corefinement on the inexact kernel, which hard-crashes (SIGBUS) on -these open surfaces — see ``test_corefine_inexact_kernel_is_unsafe`` below. A -native crash cannot be trapped by ``xfail`` (it takes down the whole pytest -process), so that path is documented and skipped rather than executed. +Note: only the default exact-kernel path performs a corefinement. Corefinement +on the inexact kernel (Simple_cartesian) hard-crashes (SIGBUS), so +``use_exact_kernel=False`` is refused up front with a ValueError — see +``test_corefine_inexact_kernel_raises``. """ from __future__ import annotations @@ -103,11 +102,13 @@ def test_corefine_preserves_area(): assert B.area == pytest.approx(area_b, rel=1e-6) -@pytest.mark.skip( - reason="corefine(use_exact_kernel=False) hard-crashes (SIGBUS) on open " - "intersecting surfaces; a native crash would take down the pytest process. " - "Tracked separately — the inexact path needs input validation or removal." -) -def test_corefine_inexact_kernel_is_unsafe(): +def test_corefine_inexact_kernel_raises(): + """corefine(use_exact_kernel=False) must raise, not crash. + + Corefinement on the inexact predicate kernel (Simple_cartesian) hard-crashes + (SIGBUS) even on trivial valid inputs, so the inexact path is refused up + front with a catchable ValueError rather than run. + """ A, B = _crossing_pair() - A.corefine(B, use_exact_kernel=False) + with pytest.raises(ValueError, match="exact kernel"): + A.corefine(B, use_exact_kernel=False) diff --git a/tests/test_error_contract.py b/tests/test_error_contract.py new file mode 100644 index 0000000..dafe401 --- /dev/null +++ b/tests/test_error_contract.py @@ -0,0 +1,71 @@ +"""Error-handling contract for the mutating operations. + +The cut/clip/corefine operations now *raise* on a genuine failure (invalid or +empty input, or an operation that was expected to succeed but did not) instead +of printing to stderr and returning 0. A return of 0 is reserved for the one +legitimate no-op: a cut whose meshes do not intersect. + +pybind11 maps C++ std::invalid_argument -> ValueError and std::runtime_error / +other std::exception -> RuntimeError. +""" +from __future__ import annotations + +import numpy as np +import pytest + +import loop_cgal +from loop_cgal._loop_cgal import ImplicitCutMode +from loop_cgal._loop_cgal import TriMesh as RawTriMesh + + +def _unit_square(): + verts = np.array( + [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]], dtype=np.float64 + ) + tris = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int32) + return loop_cgal.TriMesh.from_vertices_and_triangles(verts, tris) + + +def _empty_mesh(): + """A TriMesh with no faces (built directly via the C++ ctor, bypassing + the Python validation that would otherwise reject it).""" + return RawTriMesh(np.empty((0, 3), dtype=np.float64), np.empty((0, 3), dtype=np.int32)) + + +def test_cut_with_surface_empty_clipper_raises(): + tm = _unit_square() + with pytest.raises(ValueError, match="empty"): + tm.cut_with_surface(_empty_mesh()) + + +def test_cut_with_surface_empty_source_raises(): + empty = _empty_mesh() + with pytest.raises(ValueError, match="empty"): + empty.cut_with_surface(_unit_square()) + + +def test_clip_with_plane_empty_mesh_raises(): + empty = _empty_mesh() + with pytest.raises(ValueError, match="empty"): + empty.clip_with_plane(1.0, 0.0, 0.0, 0.0) + + +def test_cut_with_implicit_function_property_size_mismatch_raises(): + tm = _unit_square() + wrong = [0.0, 1.0] # mesh has 4 vertices + with pytest.raises(ValueError, match="does not match"): + tm.cut_with_implicit_function(wrong, 0.5, ImplicitCutMode.KEEP_POSITIVE_SIDE) + + +def test_non_intersecting_cut_is_a_no_op_not_an_error(): + """Disjoint meshes are the one legitimate no-op: return 0, do not raise.""" + tm = _unit_square() + far = np.array( + [[10, 10, 0], [11, 10, 0], [11, 11, 0], [10, 11, 0]], dtype=np.float64 + ) + clipper = loop_cgal.TriMesh.from_vertices_and_triangles( + far, np.array([[0, 1, 2], [0, 2, 3]], dtype=np.int32) + ) + removed = tm.cut_with_surface(clipper) + assert removed == 0 + assert tm.n_cells == 2 # mesh untouched From b159d2705f89678e899e46aa320a65bd9a6d1fb5 Mon Sep 17 00:00:00 2001 From: lachlangrose Date: Mon, 28 Sep 2026 21:39:46 +0930 Subject: [PATCH 12/12] ci: build against CGAL 6 on ubuntu, and state the minimum The ubuntu-latest runner is still 24.04, whose libcgal-dev is 5.6. Passing an edge_is_constrained_map to the mesh-mesh PMP::clip does not compile there: 5.6's put() for Boolean_property_map takes a non-const reference while CGAL's own corefinement internals pass it as const, so the build dies inside Corefinement/face_graph_utils.h. Fixed in CGAL 6. Verified by compiling src/mesh.cpp against each: 5.6 (ubuntu 24.04) fails, 6.0.1 (debian trixie), 6.1.1 (ubuntu 26.04) and 6.2 (homebrew) all succeed. That is why only the ubuntu jobs were red while every macOS job passed. CGAL is header-only, so fetch the 6.0.1 release rather than wait for the runner image to catch up, and pass it through CMAKE_ARGS the same way the Windows job passes its toolchain. Pinned, not floating, and it matches the version the project Dockerfile already builds. Also state the floor in CMakeLists: find_package(CGAL 6.0 REQUIRED) turns a 200-line template error into one clear line for anyone still on 5.6. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/ci.yml | 14 ++++++++++++-- CMakeLists.txt | 5 ++++- 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 34b47d8..d990e46 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -15,7 +15,8 @@ jobs: os: [ubuntu-latest, macos-latest, macos-14] python-version: ${{ fromJSON(vars.PYTHON_VERSIONS)}} env: - MACOSX_DEPLOYMENT_TARGET: '15.0' + MACOSX_DEPLOYMENT_TARGET: '15.0' + CGAL_VERSION: '6.0.1' steps: - uses: actions/checkout@v4 - uses: actions/setup-python@v5 @@ -25,8 +26,17 @@ jobs: run: | python -m pip install --upgrade pip - if: matrix.os == 'ubuntu-latest' + # ubuntu-latest is still 24.04, whose libcgal-dev is 5.6. Passing an + # edge_is_constrained_map to the mesh-mesh PMP::clip does not compile + # there: 5.6's put() for Boolean_property_map takes a non-const + # reference and the corefinement internals pass it as const. Fixed in + # CGAL 6, so fetch the headers (CGAL is header-only) rather than wait + # for the runner image to catch up. run: | - sudo apt-get install -y libgl1-mesa-dev libeigen3-dev libcgal-dev + sudo apt-get install -y libgl1-mesa-dev libeigen3-dev libboost-dev libgmp-dev libmpfr-dev + wget -q https://github.com/CGAL/cgal/releases/download/v${CGAL_VERSION}/CGAL-${CGAL_VERSION}.tar.xz + tar -xf CGAL-${CGAL_VERSION}.tar.xz + echo "CMAKE_ARGS=-DCGAL_DIR=$PWD/CGAL-${CGAL_VERSION}/lib/cmake/CGAL" >> $GITHUB_ENV - if: ${{ startsWith(matrix.os, 'macos') }} run: | brew install cgal diff --git a/CMakeLists.txt b/CMakeLists.txt index fe9c2a8..a79f145 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,7 +10,10 @@ endif() set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) -find_package(CGAL REQUIRED) +# 6.0 minimum: the mesh-mesh PMP::clip with an edge_is_constrained_map does not +# compile against 5.6 (its put() for Boolean_property_map takes a non-const ref). +# Stating the floor turns a 200-line template error into one clear message. +find_package(CGAL 6.0 REQUIRED) find_package(pybind11 REQUIRED) include_directories(${CMAKE_SOURCE_DIR}/src)