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)
diff --git a/loop_cgal/bindings.cpp b/loop_cgal/bindings.cpp
index 4eee326..4e6aa27 100644
--- a/loop_cgal/bindings.cpp
+++ b/loop_cgal/bindings.cpp
@@ -31,6 +31,15 @@ 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. "
+              "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,
@@ -51,7 +60,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..3e25cf6 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -95,11 +95,6 @@ TriMesh::TriMesh(const pybind11::array_t<double> &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()
@@ -121,6 +116,58 @@ void TriMesh::init()
   _edge_is_constrained_map = CGAL::make_boolean_property_map(_fixedEdges);
 }
 
+std::vector<std::pair<Point, Point>> TriMesh::snapshot_constraint_coords() const
+{
+  std::vector<std::pair<Point, Point>> 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<std::pair<Point, Point>> &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<std::tuple<double, double, double>, 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
@@ -168,20 +215,18 @@ void TriMesh::add_fixed_edges(const pybind11::array_t<int> &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;
     }
@@ -342,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
@@ -473,46 +509,65 @@ 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){
-          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));
-        }
-        else{
-          flag = PMP::clip(_mesh, scaled_clipper._mesh, CGAL::parameters::clip_volume(false));
-        }
-      }
-      catch (const std::exception &e)
-      {
-        std::cerr << "Corefinement failed: " << e.what() << 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);
       }
-      if (!flag)
-      {
-        std::cerr << "Warning: Clipping operation failed." << std::endl;
-      }
-      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);
         }
       }
     }
     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)
+      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 "
+                << _mesh.number_of_faces() << " faces." << std::endl;
     }
   }
   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."
@@ -530,16 +585,77 @@ int TriMesh::cutWithSurface(TriMesh &clipper,
   return faces_before - faces_after;
 }
 
-int TriMesh::clipWithPlane(double a, double b, double c, double d, bool use_exact_kernel)
+int TriMesh::corefine(TriMesh &other, bool use_exact_kernel)
 {
-  if (_mesh.number_of_vertices() == 0 || _mesh.number_of_faces() == 0)
+  // The inexact kernel path is unsafe by construction: TriangleMesh uses
+  // Simple_cartesian<double>, 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)
+    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.
+  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
   {
-    std::cerr << "Error: Source mesh is empty!" << std::endl;
-    return 0;
+    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)
+  {
+    throw std::runtime_error(std::string("corefine failed: ") + e.what());
+  }
+
+  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<int>(v_after) - static_cast<int>(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)
+    throw std::invalid_argument("clipWithPlane: source mesh is empty");
 
   const int faces_before = static_cast<int>(_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)
@@ -557,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<int>(_mesh.number_of_faces());
@@ -572,6 +689,19 @@ 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().
+  //
+  // 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;
 }
 
@@ -581,17 +711,27 @@ NumpyMesh TriMesh::save(double area_threshold,
   return export_mesh(_mesh, area_threshold, duplicate_vertex_threshold);
 }
 
-void TriMesh::cut_with_implicit_function(const std::vector<double> &property, double value, ImplicitCutMode cutmode)
+void TriMesh::cut_with_implicit_function(const std::vector<double> &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 (property.size() != _mesh.number_of_vertices())
+  if (LoopCGAL::verbose)
   {
-    std::cerr << "Error: Property size does not match number of vertices." << std::endl;
-    return;
+    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
+  // 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())
+    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<TriangleMesh, boost::vertex_index_t>::type VertexIndexMap;
   VertexIndexMap vim = get(boost::vertex_index, _mesh);
@@ -600,6 +740,17 @@ void TriMesh::cut_with_implicit_function(const std::vector<double> &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);
 
@@ -690,9 +841,6 @@ void TriMesh::cut_with_implicit_function(const std::vector<double> &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])
     {
@@ -713,6 +861,27 @@ void TriMesh::cut_with_implicit_function(const std::vector<double> &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()));
@@ -721,138 +890,43 @@ void TriMesh::cut_with_implicit_function(const std::vector<double> &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<std::size_t, 5> 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<std::size_t, 3> tarr)
-    { newtris.push_back(tarr); };
-
-    // CASE 1: v1 > value and v2 > value and v3<value
-    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<std::size_t, 3> m1 = {extended[0], extended[1], extended[3]};
-      std::array<std::size_t, 3> m2 = {extended[0], extended[3], extended[4]};
-      std::array<std::size_t, 3> 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<std::size_t, 3> m1 = {extended[0], extended[3], extended[2]};
-      std::array<std::size_t, 3> m2 = {extended[3], extended[4], extended[2]};
-      std::array<std::size_t, 3> 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<std::size_t, 3> m1 = {extended[0], extended[3], extended[4]};
-      std::array<std::size_t, 3> m2 = {extended[3], extended[1], extended[2]};
-      std::array<std::size_t, 3> 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<std::size_t, 3> m1 = {extended[0], extended[1], extended[3]};
-      std::array<std::size_t, 3> m2 = {extended[0], extended[3], extended[4]};
-      std::array<std::size_t, 3> 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<std::size_t, 3> m1 = {extended[0], extended[3], extended[2]};
-      std::array<std::size_t, 3> m2 = {extended[3], extended[4], extended[2]};
-      std::array<std::size_t, 3> 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<std::size_t, 3> m1 = {extended[0], extended[3], extended[4]};
-      std::array<std::size_t, 3> m2 = {extended[3], extended[2], extended[4]};
-      std::array<std::size_t, 3> 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<std::array<std::size_t, 3>, 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
@@ -866,12 +940,18 @@ void TriMesh::cut_with_implicit_function(const std::vector<double> &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;
       }
@@ -881,7 +961,7 @@ void TriMesh::cut_with_implicit_function(const std::vector<double> &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;
       }
@@ -995,12 +1075,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<char*>(&nv), 4);
     in.read(reinterpret_cast<char*>(&nf), 4);
+    if (!in)
+        throw std::runtime_error("TriMesh::read_from_file — truncated header in: " + path);
 
     TriangleMesh mesh;
     std::vector<TriangleMesh::Vertex_index> verts(nv);
@@ -1009,6 +1091,8 @@ TriMesh TriMesh::read_from_file(const std::string& path)
         in.read(reinterpret_cast<char*>(&x), 8);
         in.read(reinterpret_cast<char*>(&y), 8);
         in.read(reinterpret_cast<char*>(&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));
     }
 
@@ -1017,6 +1101,10 @@ TriMesh TriMesh::read_from_file(const std::string& path)
         in.read(reinterpret_cast<char*>(&i0), 4);
         in.read(reinterpret_cast<char*>(&i1), 4);
         in.read(reinterpret_cast<char*>(&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]);
     }
 
diff --git a/src/mesh.h b/src/mesh.h
index e153e33..f67362c 100644
--- a/src/mesh.h
+++ b/src/mesh.h
@@ -48,12 +48,28 @@ 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 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
         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<double>& 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<double>& 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,10 +88,24 @@ 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<std::pair<Point, Point>> 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<std::pair<Point, Point>> &saved);
+
         std::set<TriangleMesh::Edge_index> _fixedEdges;
         TriangleMesh _mesh; // The underlying CGAL surface mesh
         CGAL::Boolean_property_map<std::set<TriangleMesh::Edge_index>>
             _edge_is_constrained_map;
 };
 
-#endif // MESH_HANDLER_H
\ No newline at end of file
+#endif // MESH_H
\ No newline at end of file
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<QKey, int, QHash> 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_corefine.py b/tests/test_corefine.py
new file mode 100644
index 0000000..0afb806
--- /dev/null
+++ b/tests/test_corefine.py
@@ -0,0 +1,114 @@
+"""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 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
+
+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)
+
+
+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()
+    with pytest.raises(ValueError, match="exact kernel"):
+        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_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
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_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"
+    )
diff --git a/tests/test_implicit_cut.py b/tests/test_implicit_cut.py
new file mode 100644
index 0000000..57d2bfe
--- /dev/null
+++ b/tests/test_implicit_cut.py
@@ -0,0 +1,154 @@
+"""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}"
+    )
+
+
+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)
+
+
+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)
+
+
+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).