diff --git a/src/polymesh/impl/impl_low_level_api_mutable.hh b/src/polymesh/impl/impl_low_level_api_mutable.hh
index db531f3477cedd761f6a491a37fe7050b2704544..721d3177244ddd70b60bd6ed244077cadf6a7120 100644
--- a/src/polymesh/impl/impl_low_level_api_mutable.hh
+++ b/src/polymesh/impl/impl_low_level_api_mutable.hh
@@ -487,6 +487,65 @@ inline void low_level_api_mutable::face_split(face_index f, vertex_index v) cons
     } while (h != h_begin);
 }
 
+inline halfedge_index low_level_api_mutable::face_cut(face_index f, halfedge_index h0, halfedge_index h1) const
+{
+    // must be non-adjacent halfedges
+    POLYMESH_ASSERT(h0 != h1);
+    POLYMESH_ASSERT(next_halfedge_of(h0) != h1);
+    POLYMESH_ASSERT(prev_halfedge_of(h0) != h1);
+
+    // add face and edge
+    auto nf = alloc_face();
+    auto ne = alloc_edge();
+    auto nh0 = halfedge_of(ne, 0);
+    auto nh1 = halfedge_of(ne, 1);
+
+    halfedge_of(f) = nh0;
+    halfedge_of(nf) = nh1;
+
+    auto h0_next = next_halfedge_of(h0);
+    auto h1_next = next_halfedge_of(h1);
+
+    // rewire faces
+    {
+        auto h = h0_next;
+        do
+        {
+            if (is_boundary(opposite(h)))
+                halfedge_of(nf) = h;
+            face_of(h) = nf;
+        } while (h != h1);
+
+        face_of(nh0) = f;
+        face_of(nh1) = nf;
+    }
+
+    // fix face halfedge of f (nf is already fixed)
+    {
+        auto h = h1_next;
+        do
+        {
+            if (is_boundary(opposite(h)))
+            {
+                halfedge_of(f) = h;
+                break;
+            }
+        } while (h != h0);
+    }
+
+    // vertices
+    to_vertex_of(nh0) = to_vertex_of(h1);
+    to_vertex_of(nh1) = to_vertex_of(h0);
+
+    // connect halfedges
+    connect_prev_next(h1, nh1);
+    connect_prev_next(nh1, h0_next);
+    connect_prev_next(h0, nh0);
+    connect_prev_next(nh0, h1_next);
+
+    return nh0;
+}
+
 inline vertex_index low_level_api_mutable::edge_split_and_triangulate(edge_index e) const
 {
     auto v = add_vertex();
diff --git a/src/polymesh/impl/impl_ranges.hh b/src/polymesh/impl/impl_ranges.hh
index 5ff49d9ad0969ea831b52231b7a17ce11c12465d..e834348c13a56e48558b8128f9577bbc3289b2d2 100644
--- a/src/polymesh/impl/impl_ranges.hh
+++ b/src/polymesh/impl/impl_ranges.hh
@@ -495,6 +495,7 @@ void smart_collection<mesh_ptr, tag, iterator>::reserve(int capacity) const
 template <class mesh_ptr, class tag, class iterator>
 typename smart_collection<mesh_ptr, tag, iterator>::handle smart_collection<mesh_ptr, tag, iterator>::operator[](int idx) const
 {
+    POLYMESH_ASSERT(idx < iterator::primitive_size(*m));
     return (*m)[index(idx)];
 }
 
@@ -988,6 +989,12 @@ vertex_handle face_collection<iterator>::split(face_handle f) const
     return this->m->handle_of(low_level_api(this->m).face_split(f.idx));
 }
 
+template <class iterator>
+halfedge_handle face_collection<iterator>::cut(face_handle f, halfedge_handle h0, halfedge_handle h1) const
+{
+    return this->m->handle_of(low_level_api(this->m).face_cut(f.idx, h0.idx, h1.idx));
+}
+
 template <class iterator>
 void face_collection<iterator>::split(face_handle f, vertex_handle v) const
 {
diff --git a/src/polymesh/low_level_api.hh b/src/polymesh/low_level_api.hh
index 02b2de231f696de01bc153d1cbb2f0439ee4b9d8..d3c626a2684adbdf08cf3295331f6fbae1e57bed 100644
--- a/src/polymesh/low_level_api.hh
+++ b/src/polymesh/low_level_api.hh
@@ -318,6 +318,9 @@ public:
     /// splits a half-edge at a given ISOLATED vertex
     void halfedge_split(halfedge_index h, vertex_index v) const;
 
+    /// cuts a face along two given vertices (indicated by halfedges)
+    halfedge_index face_cut(face_index f, halfedge_index h0, halfedge_index h1) const;
+
     /// fills a face
     face_index face_fill(halfedge_index h) const;
 
diff --git a/src/polymesh/ranges.hh b/src/polymesh/ranges.hh
index 4351986612c65107de4d754604f73d95b67fb89f..f44526203dc70f6aa7946ea946cd9423e69259c6 100644
--- a/src/polymesh/ranges.hh
+++ b/src/polymesh/ranges.hh
@@ -332,6 +332,13 @@ struct face_collection : smart_collection<Mesh*, face_tag, iterator>
     /// same as before but splits at a given ISOLATED vertex
     void split(face_handle f, vertex_handle v) const;
 
+    /// Cuts a face topologically into two parts
+    /// Inserts an edge between h0.vertex_to() and h1.vertex_to()
+    /// Creates a new face as well
+    /// The returned halfedge is of the new edge and goes from h0.vertex_to() to h1.vertex_to()
+    /// (and thus points towards f)
+    halfedge_handle cut(face_handle f, halfedge_handle h0, halfedge_handle h1) const;
+
     /// Fills the half-edge ring of a given boundary half-edge
     /// Returns the new face
     face_handle fill(halfedge_handle h) const;