From dd7d2d03555a33c4f4c78c83d63959c2180a06d9 Mon Sep 17 00:00:00 2001
From: Philip Trettner <Philip.Trettner@rwth-aachen.de>
Date: Tue, 24 Jul 2018 12:43:47 +0200
Subject: [PATCH] low level API

---
 Readme.md                                     |   5 -
 src/polymesh/Mesh.hh                          | 307 +-----
 src/polymesh/impl/impl_cursors.hh             |  74 +-
 src/polymesh/impl/impl_low_level_api_base.hh  | 294 ++++++
 .../impl/impl_low_level_api_mutable.hh        | 825 +++++++++++++++
 src/polymesh/impl/impl_mesh.hh                | 986 +-----------------
 src/polymesh/impl/impl_primitive.hh           |  24 +-
 src/polymesh/impl/impl_ranges.hh              |  82 +-
 src/polymesh/iterators.hh                     |   4 +-
 src/polymesh/low_level_api.hh                 | 275 ++++-
 src/polymesh/tmp.hh                           |  25 +
 11 files changed, 1571 insertions(+), 1330 deletions(-)
 create mode 100644 src/polymesh/impl/impl_low_level_api_base.hh
 create mode 100644 src/polymesh/impl/impl_low_level_api_mutable.hh

diff --git a/Readme.md b/Readme.md
index ac0232e..bc48ad3 100644
--- a/Readme.md
+++ b/Readme.md
@@ -13,7 +13,6 @@ Best used with glm and glow.
 * attribute transformations (also between different types)
 * Debug: store compactify generation in handles to check for invalidation
 * Debug: insert is_removed assertions into handle access
-* Switch primitives and valid_primitives, check if compact flag is inlined
 * Test self-adjacent faces
 * smart ranges: filter, map
 * vector, set, map -> range
@@ -30,7 +29,3 @@ Best used with glm and glow.
 * structure of arrays instead of AOS
 * lowlevel API that allows direct half-edge manipulation and does not fix boundaries (but also mirrors high level one)
 * primitive collection sort and sort_by functions
-
-Low-level TODO:
-
-* face_of(opposite(h)) -> opposite_face_of
diff --git a/src/polymesh/Mesh.hh b/src/polymesh/Mesh.hh
index 7ee7118..635e0a5 100644
--- a/src/polymesh/Mesh.hh
+++ b/src/polymesh/Mesh.hh
@@ -10,6 +10,7 @@
 
 // often used helper
 #include "attribute_collection.hh"
+#include "low_level_api.hh"
 
 namespace polymesh
 {
@@ -113,14 +114,18 @@ public:
     /// Creates a new mesh and calls copy_from(*this);
     SharedMesh copy() const;
 
-    // internal helper
+    // internal primitives
 private:
-    // reserves a certain number of primitives
-    void reserve_faces(int capacity);
-    void reserve_vertices(int capacity);
-    void reserve_edges(int capacity);
-    void reserve_halfedges(int capacity);
+    std::vector<halfedge_index> mFaceToHalfedge;
+    std::vector<halfedge_index> mVertexToOutgoingHalfedge;
 
+    std::vector<vertex_index> mHalfedgeToVertex;
+    std::vector<face_index> mHalfedgeToFace;
+    std::vector<halfedge_index> mHalfedgeToNextHalfedge;
+    std::vector<halfedge_index> mHalfedgeToPrevHalfedge;
+
+    // primitive size
+private:
     int size_all_faces() const { return (int)mFaceToHalfedge.size(); }
     int size_all_vertices() const { return (int)mVertexToOutgoingHalfedge.size(); }
     int size_all_edges() const { return (int)mHalfedgeToNextHalfedge.size() >> 1; }
@@ -131,121 +136,8 @@ private:
     int size_valid_edges() const { return ((int)mHalfedgeToNextHalfedge.size() - mRemovedHalfedges) >> 1; }
     int size_valid_halfedges() const { return (int)mHalfedgeToNextHalfedge.size() - mRemovedHalfedges; }
 
-    // returns the next valid idx (returns the given one if valid)
-    // NOTE: the result can be invalid if no valid one was found
-    vertex_index next_valid_idx_from(vertex_index idx) const;
-    edge_index next_valid_idx_from(edge_index idx) const;
-    face_index next_valid_idx_from(face_index idx) const;
-    halfedge_index next_valid_idx_from(halfedge_index idx) const;
-    // returns the next valid idx (returns the given one if valid) counting DOWNWARDS
-    vertex_index prev_valid_idx_from(vertex_index idx) const;
-    edge_index prev_valid_idx_from(edge_index idx) const;
-    face_index prev_valid_idx_from(face_index idx) const;
-    halfedge_index prev_valid_idx_from(halfedge_index idx) const;
-
-    /// Adds a single non-connected vertex
-    /// Does NOT invalidate iterators!
-    vertex_index add_vertex();
-
-    /// Allocates a new vertex
-    vertex_index alloc_vertex();
-    /// Allocates a new face
-    face_index alloc_face();
-    /// Allocates a new edge
-    edge_index alloc_edge();
-
-    /// Allocates a given amount of vertices, faces, and halfedges
-    /// NOTE: leaves ALL of them in an unspecified state
-    void alloc_primitives(int vertices, int faces, int halfedges);
-
-    /// Adds a face consisting of N vertices
-    /// The vertices must already be sorted in CCW order
-    /// (note: trying to add already existing halfedges triggers assertions)
-    face_index add_face(vertex_handle const *v_handles, int vcnt);
-    face_index add_face(vertex_index const *v_indices, int vcnt);
-    face_index add_face(halfedge_handle const *half_loop, int vcnt);
-    face_index add_face(halfedge_index const *half_loop, int vcnt);
-
-    /// Adds an edge between two existing, distinct vertices
-    /// if edge already exists, returns it
-    edge_index add_or_get_edge(vertex_index v_from, vertex_index v_to);
-    /// Adds an edge between to existing, distinct vertices that are pointed to by two given halfedges
-    /// if edge already exists, returns it
-    edge_index add_or_get_edge(halfedge_index h_from, halfedge_index h_to);
-
-    /// same as add_or_get_edge but returns the apattrriate half-edge
-    /// Assures:
-    ///     return_value.from_vertex == v_from
-    ///     return_value.to_vertex == v_to
-    halfedge_index add_or_get_halfedge(vertex_index v_from, vertex_index v_to);
-    /// same as add_or_get_edge but returns the apattrriate half-edge
-    /// Assures:
-    ///     return_value.from_vertex == h_from.to_vertex
-    ///     return_value.to_vertex == h_to.to_vertex
-    halfedge_index add_or_get_halfedge(halfedge_index h_from, halfedge_index h_to);
-
-    /// connect two half-edges in a prev-next relation
-    void connect_prev_next(halfedge_index prev, halfedge_index next);
-
-    /// splits a face
-    vertex_index face_split(face_index f);
-    /// splits an edge
-    vertex_index edge_split(edge_index e);
-    /// splits a half-edge
-    vertex_index halfedge_split(halfedge_index h);
-
-    /// fills a face
-    face_index face_fill(halfedge_index h);
-
-    /// attaches a given vertex to the to-vertex of a given half-edge
-    void halfedge_attach(halfedge_index h, vertex_index v);
-
-    /// collapse a vertex
-    void vertex_collapse(vertex_index v);
-    /// collapse a half-edge
-    void halfedge_collapse(halfedge_index h);
-
-    /// rotates an edge to next
-    void edge_rotate_next(edge_index e);
-    /// rotates an edge to prev
-    void edge_rotate_prev(edge_index e);
-    /// rotates a half-edge to next
-    void halfedge_rotate_next(halfedge_index h);
-    /// rotates a half-edge to prev
-    void halfedge_rotate_prev(halfedge_index h);
-
-    /// removes a face (actually sets the removed status)
-    /// modifies all adjacent vertices so that they correctly report is_boundary true
-    void remove_face(face_index f_idx);
-    /// removes both adjacent faces, then removes both half edges
-    void remove_edge(edge_index e_idx);
-    /// removes all adjacent edges, then the vertex
-    void remove_vertex(vertex_index v_idx);
-
-    /// choses a new outgoing half-edge for a given vertex, prefers boundary ones
-    /// assumes non-isolated vertex
-    void fix_boundary_state_of(vertex_index v_idx);
-    /// choses a new half-edge for a given face, prefers boundary ones
-    void fix_boundary_state_of(face_index f_idx);
-    /// choses a new half-edge for all vertices of a face, prefers boundary ones
-    void fix_boundary_state_of_vertices(face_index f_idx);
-
-    // attributes
-    bool is_free(halfedge_index idx) const;
-
-    bool is_boundary(vertex_index idx) const;
-    bool is_boundary(halfedge_index idx) const;
-    bool is_boundary(edge_index idx) const;
-    bool is_boundary(face_index idx) const;
-
-    bool is_removed(vertex_index idx) const;
-    bool is_removed(face_index idx) const;
-    bool is_removed(edge_index idx) const;
-    bool is_removed(halfedge_index idx) const;
-
-    bool is_isolated(vertex_index idx) const;
-    bool is_isolated(edge_index idx) const;
-
+    // primitive access
+private:
     vertex_index &to_vertex_of(halfedge_index idx);
     face_index &face_of(halfedge_index idx);
     halfedge_index &next_halfedge_of(halfedge_index idx);
@@ -260,45 +152,26 @@ private:
     halfedge_index halfedge_of(face_index idx) const;
     halfedge_index outgoing_halfedge_of(vertex_index idx) const;
 
-    void set_removed(vertex_index idx);
-    void set_removed(face_index idx);
-    void set_removed(edge_index idx);
-
-    /// Returns the opposite of a given valid half-edge
-    halfedge_index opposite(halfedge_index he) const;
-    face_index opposite_face_of(halfedge_index he) const;
-
-    /// Makes two half-edges adjacent
-    /// Ensures:
-    ///     * he_in->next == he_out
-    ///     * he_out->prev == he_in
-    /// Requires:
-    ///     * he_in->is_free()
-    ///     * he_out->is_free()
-    /// Only works if a free incident half-edge is available
-    void make_adjacent(halfedge_index he_in, halfedge_index he_out);
-
-    /// finds the next free incoming half-edge around a certain vertex
-    /// starting from in_begin, EXCLUDING in_end (if in_begin == in_end, the whole vertex is searched)
-    /// returns invalid index if no edge is found
-    halfedge_index find_free_incident(halfedge_index in_begin, halfedge_index in_end) const;
-    /// finds a free incident incoming half-edge around a given vertex
-    halfedge_index find_free_incident(vertex_index v) const;
-
-    /// returns half-edge going from `from`, point to `to`
-    /// returns invalid index if not exists
-    halfedge_index find_halfedge(vertex_index from, vertex_index to) const;
-
-    /// returns edge index belonging to a half-edge
-    edge_index edge_of(halfedge_index idx) const { return edge_index(idx.value >> 1); }
-    /// returns a half-edge belonging to an edge
-    halfedge_index halfedge_of(edge_index idx, int i) const { return halfedge_index((idx.value << 1) + i); }
-    vertex_index to_vertex_of(edge_index idx, int i) const { return to_vertex_of(halfedge_of(idx, i)); }
-    face_index face_of(edge_index idx, int i) const { return face_of(halfedge_of(idx, i)); }
-
-    vertex_index &from_vertex_of(halfedge_index idx);
-    vertex_index from_vertex_of(halfedge_index idx) const;
+    // primitive allocation
+private:
+    /// Allocates a new vertex
+    vertex_index alloc_vertex();
+    /// Allocates a new face
+    face_index alloc_face();
+    /// Allocates a new edge
+    edge_index alloc_edge();
+    /// Allocates a given amount of vertices, faces, and halfedges
+    /// NOTE: leaves ALL of them in an unspecified state
+    void alloc_primitives(int vertices, int faces, int halfedges);
 
+    // reserves a certain number of primitives
+    void reserve_faces(int capacity);
+    void reserve_vertices(int capacity);
+    void reserve_edges(int capacity);
+    void reserve_halfedges(int capacity);
+
+    // primitive reordering
+private:
     /// applies an index remapping to all face indices (p[curr_idx] = new_idx)
     void permute_faces(std::vector<int> const &p);
     /// applies an index remapping to all edge (and half-edge) indices (p[curr_idx] = new_idx)
@@ -306,106 +179,6 @@ private:
     /// applies an index remapping to all vertices indices (p[curr_idx] = new_idx)
     void permute_vertices(std::vector<int> const &p);
 
-    // internal datastructures
-private:
-    // 4 byte per face
-    struct face_info
-    {
-        halfedge_index halfedge; ///< one half-edge bounding this face
-
-        bool is_valid() const { return halfedge.is_valid(); }
-        void set_removed() { halfedge = halfedge_index::invalid(); }
-        // is_boundary: check if half-edge is boundary
-    };
-
-    // 4 byte per vertex
-    struct vertex_info
-    {
-        halfedge_index outgoing_halfedge;
-
-        /// a vertex can be valid even without outgoing halfedge
-        bool is_valid() const { return outgoing_halfedge.value >= -1; }
-        bool is_isolated() const { return !outgoing_halfedge.is_valid(); }
-        void set_removed() { outgoing_halfedge = halfedge_index(-2); }
-        // is_boundary: check if outgoing_halfedge is boundary
-    };
-
-    // 16 byte per edge
-    struct halfedge_info
-    {
-        vertex_index to_vertex;       ///< half-edge points towards this vertex
-        face_index face;              ///< might be invalid if boundary
-        halfedge_index next_halfedge; ///< CCW
-        halfedge_index prev_halfedge; ///< CW
-        // opposite half-edge idx is "idx ^ 1"
-        // edge idx is "idx >> 1"
-
-        bool is_valid() const { return to_vertex.is_valid(); }
-
-        /// a half-edge is free if it is a boundary, aka has no associated face
-        bool is_free() const { return !face.is_valid(); }
-
-        // CAUTION: remove both HE belonging to an edge
-        void set_removed() { to_vertex = vertex_index::invalid(); }
-    };
-
-    // internal primitives
-private:
-    // std::vector<face_info> mFaces;
-    // std::vector<vertex_info> mVertices;
-    // std::vector<halfedge_info> mHalfedges;
-
-    std::vector<halfedge_index> mFaceToHalfedge;
-    std::vector<halfedge_index> mVertexToOutgoingHalfedge;
-
-    std::vector<vertex_index> mHalfedgeToVertex;
-    std::vector<face_index> mHalfedgeToFace;
-    std::vector<halfedge_index> mHalfedgeToNextHalfedge;
-    std::vector<halfedge_index> mHalfedgeToPrevHalfedge;
-
-    /*
-    struct face_info &face(face_index i)
-    {
-        assert(i.is_valid() && i.value < size_all_faces());
-        return mFaces[i.value];
-    }
-    struct face_info const &face(face_index i) const
-    {
-        assert(i.is_valid() && i.value < size_all_faces());
-        return mFaces[i.value];
-    }
-    struct vertex_info &vertex(vertex_index i)
-    {
-        assert(i.is_valid() && i.value < size_all_vertices());
-        return mVertices[i.value];
-    }
-    struct vertex_info const &vertex(vertex_index i) const
-    {
-        assert(i.is_valid() && i.value < size_all_vertices());
-        return mVertices[i.value];
-    }
-    struct halfedge_info &halfedge(halfedge_index i)
-    {
-        assert(i.is_valid() && i.value < size_all_halfedges());
-        return mHalfedges[i.value];
-    }
-    struct halfedge_info const &halfedge(halfedge_index i) const
-    {
-        assert(i.is_valid() && i.value < size_all_halfedges());
-        return mHalfedges[i.value];
-    }
-    struct halfedge_info &halfedge(edge_index i, int h)
-    {
-        assert(i.is_valid() && i.value < size_all_edges());
-        return mHalfedges[(i.value << 1) + h];
-    }
-    struct halfedge_info const &halfedge(edge_index i, int h) const
-    {
-        assert(i.is_valid() && i.value < size_all_edges());
-        return mHalfedges[(i.value << 1) + h];
-    }
-    */
-
     // internal state
 private:
     bool mCompact = true;
@@ -434,7 +207,7 @@ private:
 
     // friends
 private:
-    friend struct vertex_handle;
+    /*friend struct vertex_handle;
     friend struct all_vertex_iterator;
     friend struct valid_vertex_iterator;
     friend struct vertex_attribute_base;
@@ -479,9 +252,17 @@ private:
     template <class iterator>
     friend struct edge_collection;
     template <class iterator>
-    friend struct halfedge_collection;
+    friend struct halfedge_collection;*/
+
+    // for attribute registration
+    template <class tag>
+    friend struct primitive_attribute_base;
 
-    friend struct low_level_api;
+    // for low level access
+    template <class MeshT>
+    friend struct low_level_api_base;
+    friend struct low_level_api_const;
+    friend struct low_level_api_mutable;
 };
 }
 
@@ -490,6 +271,8 @@ private:
 #include "impl/impl_attributes.hh"
 #include "impl/impl_cursors.hh"
 #include "impl/impl_iterators.hh"
+#include "impl/impl_low_level_api_base.hh"
+#include "impl/impl_low_level_api_mutable.hh"
 #include "impl/impl_mesh.hh"
 #include "impl/impl_primitive.hh"
 #include "impl/impl_ranges.hh"
diff --git a/src/polymesh/impl/impl_cursors.hh b/src/polymesh/impl/impl_cursors.hh
index c68761b..ec9dca8 100644
--- a/src/polymesh/impl/impl_cursors.hh
+++ b/src/polymesh/impl/impl_cursors.hh
@@ -18,39 +18,39 @@ auto primitive_handle<tag>::operator[](FuncT&& f) const -> tmp::result_type_of<F
     return f(*static_cast<typename primitive<tag>::handle const*>(this));
 }
 
-inline bool vertex_handle::is_removed() const { return idx.is_valid() && mesh->is_removed(idx); }
-inline bool face_handle::is_removed() const { return idx.is_valid() && mesh->is_removed(idx); }
-inline bool edge_handle::is_removed() const { return idx.is_valid() && mesh->is_removed(idx); }
-inline bool halfedge_handle::is_removed() const { return idx.is_valid() && mesh->is_removed(idx); }
-
-inline bool vertex_handle::is_isolated() const { return mesh->is_isolated(idx); }
-inline bool edge_handle::is_isolated() const { return mesh->is_isolated(idx); }
-
-inline bool vertex_handle::is_boundary() const { return mesh->is_boundary(idx); }
-inline bool face_handle::is_boundary() const { return mesh->is_boundary(idx); }
-inline bool edge_handle::is_boundary() const { return mesh->is_boundary(idx); }
-inline bool halfedge_handle::is_boundary() const { return mesh->is_boundary(idx); }
-
-inline vertex_handle halfedge_handle::vertex_to() const { return mesh->handle_of(mesh->to_vertex_of(idx)); }
-inline vertex_handle halfedge_handle::vertex_from() const { return mesh->handle_of(mesh->from_vertex_of(idx)); }
-inline halfedge_handle halfedge_handle::next() const { return mesh->handle_of(mesh->next_halfedge_of(idx)); }
-inline halfedge_handle halfedge_handle::prev() const { return mesh->handle_of(mesh->prev_halfedge_of(idx)); }
-inline halfedge_handle halfedge_handle::opposite() const { return mesh->handle_of(mesh->opposite(idx)); }
-inline edge_handle halfedge_handle::edge() const { return mesh->handle_of(mesh->edge_of(idx)); }
-inline face_handle halfedge_handle::face() const { return mesh->handle_of(mesh->face_of(idx)); }
-inline face_handle halfedge_handle::opposite_face() const { return mesh->handle_of(mesh->opposite_face_of(idx)); }
-
-inline halfedge_handle edge_handle::halfedgeA() const { return mesh->handle_of(mesh->halfedge_of(idx, 0)); }
-inline halfedge_handle edge_handle::halfedgeB() const { return mesh->handle_of(mesh->halfedge_of(idx, 1)); }
-inline vertex_handle edge_handle::vertexA() const { return mesh->handle_of(mesh->to_vertex_of(idx, 0)); }
-inline vertex_handle edge_handle::vertexB() const { return mesh->handle_of(mesh->to_vertex_of(idx, 1)); }
-inline face_handle edge_handle::faceA() const { return mesh->handle_of(mesh->face_of(idx, 0)); }
-inline face_handle edge_handle::faceB() const { return mesh->handle_of(mesh->face_of(idx, 1)); }
+inline bool vertex_handle::is_removed() const { return idx.is_valid() && low_level_api(mesh).is_removed(idx); }
+inline bool face_handle::is_removed() const { return idx.is_valid() && low_level_api(mesh).is_removed(idx); }
+inline bool edge_handle::is_removed() const { return idx.is_valid() && low_level_api(mesh).is_removed(idx); }
+inline bool halfedge_handle::is_removed() const { return idx.is_valid() && low_level_api(mesh).is_removed(idx); }
+
+inline bool vertex_handle::is_isolated() const { return low_level_api(mesh).is_isolated(idx); }
+inline bool edge_handle::is_isolated() const { return low_level_api(mesh).is_isolated(idx); }
+
+inline bool vertex_handle::is_boundary() const { return low_level_api(mesh).is_boundary(idx); }
+inline bool face_handle::is_boundary() const { return low_level_api(mesh).is_boundary(idx); }
+inline bool edge_handle::is_boundary() const { return low_level_api(mesh).is_boundary(idx); }
+inline bool halfedge_handle::is_boundary() const { return low_level_api(mesh).is_boundary(idx); }
+
+inline vertex_handle halfedge_handle::vertex_to() const { return mesh->handle_of(low_level_api(mesh).to_vertex_of(idx)); }
+inline vertex_handle halfedge_handle::vertex_from() const { return mesh->handle_of(low_level_api(mesh).from_vertex_of(idx)); }
+inline halfedge_handle halfedge_handle::next() const { return mesh->handle_of(low_level_api(mesh).next_halfedge_of(idx)); }
+inline halfedge_handle halfedge_handle::prev() const { return mesh->handle_of(low_level_api(mesh).prev_halfedge_of(idx)); }
+inline halfedge_handle halfedge_handle::opposite() const { return mesh->handle_of(low_level_api(mesh).opposite(idx)); }
+inline edge_handle halfedge_handle::edge() const { return mesh->handle_of(low_level_api(mesh).edge_of(idx)); }
+inline face_handle halfedge_handle::face() const { return mesh->handle_of(low_level_api(mesh).face_of(idx)); }
+inline face_handle halfedge_handle::opposite_face() const { return mesh->handle_of(low_level_api(mesh).opposite_face_of(idx)); }
+
+inline halfedge_handle edge_handle::halfedgeA() const { return mesh->handle_of(low_level_api(mesh).halfedge_of(idx, 0)); }
+inline halfedge_handle edge_handle::halfedgeB() const { return mesh->handle_of(low_level_api(mesh).halfedge_of(idx, 1)); }
+inline vertex_handle edge_handle::vertexA() const { return mesh->handle_of(low_level_api(mesh).to_vertex_of(idx, 0)); }
+inline vertex_handle edge_handle::vertexB() const { return mesh->handle_of(low_level_api(mesh).to_vertex_of(idx, 1)); }
+inline face_handle edge_handle::faceA() const { return mesh->handle_of(low_level_api(mesh).face_of(idx, 0)); }
+inline face_handle edge_handle::faceB() const { return mesh->handle_of(low_level_api(mesh).face_of(idx, 1)); }
 
 inline face_handle vertex_handle::any_face() const
 {
-    auto h = mesh->outgoing_halfedge_of(idx);
-    return mesh->handle_of(h.is_valid() ? mesh->face_of(h) : face_index::invalid());
+    auto h = low_level_api(mesh).outgoing_halfedge_of(idx);
+    return mesh->handle_of(h.is_valid() ? low_level_api(mesh).face_of(h) : face_index::invalid());
 }
 
 inline face_handle vertex_handle::any_valid_face() const
@@ -61,23 +61,23 @@ inline face_handle vertex_handle::any_valid_face() const
     return mesh->handle_of(face_index::invalid());
 }
 
-inline halfedge_handle vertex_handle::any_outgoing_halfedge() const { return mesh->handle_of(mesh->outgoing_halfedge_of(idx)); }
+inline halfedge_handle vertex_handle::any_outgoing_halfedge() const { return mesh->handle_of(low_level_api(mesh).outgoing_halfedge_of(idx)); }
 
 inline halfedge_handle vertex_handle::any_incoming_halfedge() const
 {
-    auto h = mesh->outgoing_halfedge_of(idx);
-    return mesh->handle_of(h.is_valid() ? mesh->opposite(h) : halfedge_index::invalid());
+    auto h = low_level_api(mesh).outgoing_halfedge_of(idx);
+    return mesh->handle_of(h.is_valid() ? low_level_api(mesh).opposite(h) : halfedge_index::invalid());
 }
 
 inline edge_handle vertex_handle::any_edge() const
 {
-    auto h = mesh->outgoing_halfedge_of(idx);
-    return mesh->handle_of(h.is_valid() ? mesh->edge_of(h) : edge_index::invalid());
+    auto h = low_level_api(mesh).outgoing_halfedge_of(idx);
+    return mesh->handle_of(h.is_valid() ? low_level_api(mesh).edge_of(h) : edge_index::invalid());
 }
 
-inline vertex_handle face_handle::any_vertex() const { return mesh->handle_of(mesh->to_vertex_of(mesh->halfedge_of(idx))); }
+inline vertex_handle face_handle::any_vertex() const { return mesh->handle_of(low_level_api(mesh).to_vertex_of(low_level_api(mesh).halfedge_of(idx))); }
 
-inline halfedge_handle face_handle::any_halfedge() const { return mesh->handle_of(mesh->halfedge_of(idx)); }
+inline halfedge_handle face_handle::any_halfedge() const { return mesh->handle_of(low_level_api(mesh).halfedge_of(idx)); }
 
 inline face_vertex_ring face_handle::vertices() const { return {*this}; }
 
diff --git a/src/polymesh/impl/impl_low_level_api_base.hh b/src/polymesh/impl/impl_low_level_api_base.hh
new file mode 100644
index 0000000..856d53d
--- /dev/null
+++ b/src/polymesh/impl/impl_low_level_api_base.hh
@@ -0,0 +1,294 @@
+#pragma once
+
+#include "../Mesh.hh"
+
+namespace polymesh
+{
+template <class MeshT>
+tmp::ref_if_mut<vertex_index, MeshT> low_level_api_base<MeshT>::to_vertex_of(halfedge_index idx) const
+{
+    return m.to_vertex_of(idx);
+}
+
+template <class MeshT>
+tmp::ref_if_mut<face_index, MeshT> low_level_api_base<MeshT>::face_of(halfedge_index idx) const
+{
+    return m.face_of(idx);
+}
+
+template <class MeshT>
+tmp::ref_if_mut<halfedge_index, MeshT> low_level_api_base<MeshT>::next_halfedge_of(halfedge_index idx) const
+{
+    return m.next_halfedge_of(idx);
+}
+
+template <class MeshT>
+tmp::ref_if_mut<halfedge_index, MeshT> low_level_api_base<MeshT>::prev_halfedge_of(halfedge_index idx) const
+{
+    return m.prev_halfedge_of(idx);
+}
+
+template <class MeshT>
+tmp::ref_if_mut<halfedge_index, MeshT> low_level_api_base<MeshT>::halfedge_of(face_index idx) const
+{
+    return m.halfedge_of(idx);
+}
+
+template <class MeshT>
+tmp::ref_if_mut<halfedge_index, MeshT> low_level_api_base<MeshT>::outgoing_halfedge_of(vertex_index idx) const
+{
+    return m.outgoing_halfedge_of(idx);
+}
+
+
+template<class MeshT>
+int low_level_api_base<MeshT>::size_all_faces() const
+{
+    return m.size_all_faces();
+}
+
+template<class MeshT>
+int low_level_api_base<MeshT>::size_all_vertices() const
+{
+    return m.size_all_vertices();
+}
+
+template<class MeshT>
+int low_level_api_base<MeshT>::size_all_edges() const
+{
+    return m.size_all_edges();
+}
+
+template<class MeshT>
+int low_level_api_base<MeshT>::size_all_halfedges() const
+{
+    return m.size_all_halfedges();
+}
+
+template<class MeshT>
+int low_level_api_base<MeshT>::size_valid_faces() const
+{
+    return m.size_valid_faces();
+}
+
+template<class MeshT>
+int low_level_api_base<MeshT>::size_valid_vertices() const
+{
+    return m.size_valid_vertices();
+}
+
+template<class MeshT>
+int low_level_api_base<MeshT>::size_valid_edges() const
+{
+    return m.size_valid_edges();
+}
+
+template<class MeshT>
+int low_level_api_base<MeshT>::size_valid_halfedges() const
+{
+    return m.size_valid_halfedges();
+}
+
+template <class MeshT>
+halfedge_index low_level_api_base<MeshT>::find_free_incident(halfedge_index in_begin, halfedge_index in_end) const
+{
+    assert(to_vertex_of(in_begin) == to_vertex_of(in_end));
+
+    auto he = in_begin;
+    do
+    {
+        assert(to_vertex_of(he) == to_vertex_of(in_end));
+
+        // free? found one!
+        if (is_free(he))
+            return he;
+
+        // next half-edge of vertex
+        he = opposite(next_halfedge_of(he));
+    } while (he != in_end);
+
+    return halfedge_index::invalid();
+}
+
+template <class MeshT>
+halfedge_index low_level_api_base<MeshT>::find_free_incident(vertex_index v) const
+{
+    auto in_begin = opposite(outgoing_halfedge_of(v));
+    return find_free_incident(in_begin, in_begin);
+}
+
+template <class MeshT>
+halfedge_index low_level_api_base<MeshT>::find_halfedge(vertex_index from, vertex_index to) const
+{
+    auto he_begin = outgoing_halfedge_of(from);
+    if (!he_begin.is_valid())
+        return halfedge_index::invalid(); // isolated vertex
+
+    auto he = he_begin;
+    do
+    {
+        // found?
+        if (to_vertex_of(he) == to)
+            return he;
+
+        // advance
+        he = next_halfedge_of(opposite(he));
+
+    } while (he != he_begin);
+
+    return halfedge_index::invalid(); // not found
+}
+
+template <class MeshT>
+bool low_level_api_base<MeshT>::is_boundary(vertex_index idx) const
+{
+    auto oh = outgoing_halfedge_of(idx);
+    return !oh.is_valid() || is_boundary(oh);
+}
+
+template <class MeshT>
+bool low_level_api_base<MeshT>::is_free(halfedge_index idx) const
+{
+    return face_of(idx).is_invalid();
+}
+template <class MeshT>
+bool low_level_api_base<MeshT>::is_boundary(halfedge_index idx) const
+{
+    return is_free(idx);
+}
+template <class MeshT>
+bool low_level_api_base<MeshT>::is_boundary(face_index idx) const
+{
+    return is_free(opposite(halfedge_of(idx)));
+}
+template <class MeshT>
+bool low_level_api_base<MeshT>::is_boundary(edge_index idx) const
+{
+    return is_free(halfedge_of(idx, 0)) || is_free(halfedge_of(idx, 1));
+}
+
+template <class MeshT>
+bool low_level_api_base<MeshT>::is_isolated(vertex_index idx) const
+{
+    return outgoing_halfedge_of(idx).is_invalid();
+}
+template <class MeshT>
+bool low_level_api_base<MeshT>::is_isolated(edge_index idx) const
+{
+    return is_free(halfedge_of(idx, 0)) && is_free(halfedge_of(idx, 1));
+}
+
+template <class MeshT>
+bool low_level_api_base<MeshT>::is_removed(vertex_index idx) const
+{
+    return outgoing_halfedge_of(idx).value == -2;
+}
+template <class MeshT>
+bool low_level_api_base<MeshT>::is_removed(face_index idx) const
+{
+    return halfedge_of(idx).is_invalid();
+}
+template <class MeshT>
+bool low_level_api_base<MeshT>::is_removed(edge_index idx) const
+{
+    return to_vertex_of(halfedge_of(idx, 0)).is_invalid();
+}
+template <class MeshT>
+bool low_level_api_base<MeshT>::is_removed(halfedge_index idx) const
+{
+    return to_vertex_of(idx).is_invalid();
+}
+
+template <class MeshT>
+halfedge_index low_level_api_base<MeshT>::opposite(halfedge_index he) const
+{
+    return halfedge_index(he.value ^ 1);
+}
+template <class MeshT>
+face_index low_level_api_base<MeshT>::opposite_face_of(halfedge_index he) const
+{
+    return face_of(opposite(he));
+}
+
+template <class MeshT>
+vertex_index low_level_api_base<MeshT>::from_vertex_of(halfedge_index idx) const
+{
+    return to_vertex_of(opposite(idx));
+}
+
+template <class MeshT>
+vertex_index low_level_api_base<MeshT>::next_valid_idx_from(vertex_index idx) const
+{
+    auto s = size_all_vertices();
+    auto i = idx;
+    while (i.value < s && is_removed(i))
+        i.value++;
+    return i;
+}
+
+template <class MeshT>
+vertex_index low_level_api_base<MeshT>::prev_valid_idx_from(vertex_index idx) const
+{
+    auto i = idx;
+    while (i.value >= 0 && is_removed(i))
+        i.value--;
+    return i;
+}
+
+template <class MeshT>
+edge_index low_level_api_base<MeshT>::next_valid_idx_from(edge_index idx) const
+{
+    auto s = size_all_edges();
+    auto i = idx;
+    while (i.value < s && is_removed(i))
+        i.value++;
+    return i;
+}
+
+template <class MeshT>
+edge_index low_level_api_base<MeshT>::prev_valid_idx_from(edge_index idx) const
+{
+    auto i = idx;
+    while (i.value >= 0 && is_removed(i))
+        i.value--;
+    return i;
+}
+
+template <class MeshT>
+face_index low_level_api_base<MeshT>::next_valid_idx_from(face_index idx) const
+{
+    auto s = size_all_faces();
+    auto i = idx;
+    while (i.value < s && is_removed(i))
+        i.value++;
+    return i;
+}
+
+template <class MeshT>
+face_index low_level_api_base<MeshT>::prev_valid_idx_from(face_index idx) const
+{
+    auto i = idx;
+    while (i.value >= 0 && is_removed(i))
+        i.value--;
+    return i;
+}
+
+template <class MeshT>
+halfedge_index low_level_api_base<MeshT>::next_valid_idx_from(halfedge_index idx) const
+{
+    auto s = size_all_halfedges();
+    auto i = idx;
+    while (i.value < s && is_removed(i))
+        i.value++;
+    return i;
+}
+
+template <class MeshT>
+halfedge_index low_level_api_base<MeshT>::prev_valid_idx_from(halfedge_index idx) const
+{
+    auto i = idx;
+    while (i.value >= 0 && is_removed(i))
+        i.value--;
+    return i;
+}
+}
diff --git a/src/polymesh/impl/impl_low_level_api_mutable.hh b/src/polymesh/impl/impl_low_level_api_mutable.hh
new file mode 100644
index 0000000..01b0fdb
--- /dev/null
+++ b/src/polymesh/impl/impl_low_level_api_mutable.hh
@@ -0,0 +1,825 @@
+#pragma once
+
+#include "../Mesh.hh"
+
+namespace polymesh
+{
+inline vertex_index low_level_api_mutable::add_vertex() const { return alloc_vertex(); }
+
+inline vertex_index low_level_api_mutable::alloc_vertex() const { return m.alloc_vertex(); }
+inline face_index low_level_api_mutable::alloc_face() const { return m.alloc_face(); }
+inline edge_index low_level_api_mutable::alloc_edge() const { return m.alloc_edge(); }
+inline void low_level_api_mutable::alloc_primitives(int vertices, int faces, int halfedges) const { m.alloc_primitives(vertices, faces, halfedges); }
+
+inline void low_level_api_mutable::permute_faces(const std::vector<int> &p) const { m.permute_faces(p); }
+inline void low_level_api_mutable::permute_edges(const std::vector<int> &p) const { m.permute_edges(p); }
+inline void low_level_api_mutable::permute_vertices(const std::vector<int> &p) const { m.permute_vertices(p); }
+
+inline face_index low_level_api_mutable::add_face(const vertex_handle *v_handles, int vcnt) const
+{
+    m.mFaceInsertCache.resize(vcnt);
+    for (auto i = 0; i < vcnt; ++i)
+        m.mFaceInsertCache[i] = add_or_get_halfedge(v_handles[i].idx, v_handles[(i + 1) % vcnt].idx);
+    return add_face(m.mFaceInsertCache.data(), vcnt);
+}
+
+inline face_index low_level_api_mutable::add_face(const vertex_index *v_indices, int vcnt) const
+{
+    m.mFaceInsertCache.resize(vcnt);
+    for (auto i = 0; i < vcnt; ++i)
+        m.mFaceInsertCache[i] = add_or_get_halfedge(v_indices[i], v_indices[(i + 1) % vcnt]);
+    return add_face(m.mFaceInsertCache.data(), vcnt);
+}
+
+inline face_index low_level_api_mutable::add_face(const halfedge_handle *half_loop, int vcnt) const
+{
+    m.mFaceInsertCache.resize(vcnt);
+    for (auto i = 0; i < vcnt; ++i)
+        m.mFaceInsertCache[i] = half_loop[i].idx;
+    return add_face(m.mFaceInsertCache.data(), vcnt);
+}
+
+inline face_index low_level_api_mutable::add_face(const halfedge_index *half_loop, int vcnt) const
+{
+    assert(vcnt >= 3 && "no support for less-than-triangular faces");
+
+    auto fidx = alloc_face();
+
+    // ensure that half-edges are adjacent at each vertex
+    for (auto i = 0; i < vcnt; ++i)
+    {
+        auto h0 = half_loop[i];
+        auto h1 = half_loop[(i + 1) % vcnt];
+
+        // half-edge must form a chain
+        assert(to_vertex_of(h0) == from_vertex_of(h1) && "half-edges do not form a chain");
+        // half-edge must be free, i.e. allow a new polygon
+        assert(is_free(h0) && "half-edge already contains a face");
+
+        // make them adjacent
+        make_adjacent(h0, h1);
+
+        // link face
+        face_of(h0) = fidx;
+    }
+
+    // fix boundary states
+    for (auto i = 0; i < vcnt; ++i)
+    {
+        auto h = half_loop[i];
+        auto v = to_vertex_of(h);
+        auto f = opposite_face_of(h);
+
+        // fix vertex
+        fix_boundary_state_of(v);
+
+        // fix face
+        if (f.is_valid())
+            fix_boundary_state_of(f);
+    }
+
+    // set up face data
+    halfedge_of(fidx) = half_loop[0];
+
+    // fix new face
+    fix_boundary_state_of(fidx);
+
+    return fidx;
+}
+
+inline edge_index low_level_api_mutable::add_or_get_edge(vertex_index v_from, vertex_index v_to) const
+{
+    assert(v_from != v_to);
+
+    // already exists?
+    auto he = find_halfedge(v_from, v_to);
+    if (he.is_valid())
+        return edge_of(he);
+
+    // allocate new
+    auto e = alloc_edge();
+    auto h_from_to = halfedge_of(e, 0);
+    auto h_to_from = halfedge_of(e, 1);
+
+    // setup data (self-connected edge)
+    to_vertex_of(h_from_to) = v_to;
+    to_vertex_of(h_to_from) = v_from;
+    next_halfedge_of(h_from_to) = h_to_from;
+    next_halfedge_of(h_to_from) = h_from_to;
+
+    // link from vertex
+    if (is_isolated(v_from))
+        outgoing_halfedge_of(v_from) = h_from_to;
+    else
+    {
+        auto from_in = find_free_incident(v_from);
+        assert(from_in.is_valid() && "vertex is already fully connected");
+
+        auto from_out = next_halfedge_of(from_in);
+
+        connect_prev_next(from_in, h_from_to);
+        connect_prev_next(h_to_from, from_out);
+    }
+
+    // link to vertex
+    if (is_isolated(v_to))
+        outgoing_halfedge_of(v_to) = h_to_from;
+    else
+    {
+        auto to_in = find_free_incident(v_to);
+        assert(to_in.is_valid() && "vertex is already fully connected");
+
+        auto to_out = next_halfedge_of(to_in);
+
+        connect_prev_next(to_in, h_to_from);
+        connect_prev_next(h_from_to, to_out);
+    }
+
+    return e;
+}
+
+inline halfedge_index low_level_api_mutable::add_or_get_halfedge(vertex_index v_from, vertex_index v_to) const
+{
+    auto e = add_or_get_edge(v_from, v_to);
+    auto h0 = halfedge_of(e, 0);
+    auto h1 = halfedge_of(e, 1);
+    return to_vertex_of(h0) == v_to ? h0 : h1;
+}
+
+inline edge_index low_level_api_mutable::add_or_get_edge(halfedge_index h_from, halfedge_index h_to) const
+{
+    assert(h_from != h_to);
+
+    auto v_from = to_vertex_of(h_from);
+    auto v_to = to_vertex_of(h_to);
+
+    auto ex_he = find_halfedge(v_from, v_to);
+    if (ex_he.is_valid())
+    {
+        assert(prev_halfedge_of(ex_he) == h_from && prev_halfedge_of(opposite(ex_he)) == h_to);
+
+        // TODO: Maybe try rewriting an existing halfedge that does NOT yet have the right connection.
+        return edge_of(ex_he);
+    }
+
+    assert(is_free(h_from) && is_free(h_to) && "Cannot insert into a face");
+
+    // allocate new
+    auto e = alloc_edge();
+    auto h_from_to = halfedge_of(e, 0);
+    auto h_to_from = halfedge_of(e, 1);
+
+    // setup data (self-connected edge)
+    to_vertex_of(h_from_to) = v_to;
+    to_vertex_of(h_to_from) = v_from;
+
+    // Link from side
+    auto h_from_next = next_halfedge_of(h_from);
+
+    connect_prev_next(h_from, h_from_to);
+    connect_prev_next(h_from_next, h_to_from);
+
+    // Link to side
+    auto h_to_next = next_halfedge_of(h_to);
+
+    connect_prev_next(h_to, h_to_from);
+    connect_prev_next(h_to_next, h_from_to);
+
+    return e;
+}
+
+inline halfedge_index low_level_api_mutable::add_or_get_halfedge(halfedge_index h_from, halfedge_index h_to) const
+{
+    auto e = add_or_get_edge(h_from, h_to);
+    auto h0 = halfedge_of(e, 0);
+    auto h1 = halfedge_of(e, 1);
+    return next_halfedge_of(h_from) == h0 ? h0 : h1;
+}
+
+inline void low_level_api_mutable::make_adjacent(halfedge_index he_in, halfedge_index he_out) const
+{
+    // see http://kaba.hilvi.org/homepage/blog/halfedge/halfedge.htm ::makeAdjacent
+
+    auto he_b = next_halfedge_of(he_in);
+    auto he_d = prev_halfedge_of(he_out);
+
+    // already correct
+    if (he_b == he_out)
+        return;
+
+    // find free half-edge after `out` but before `in`
+    auto he_g = find_free_incident(opposite(he_out), he_in);
+    assert(he_g.is_valid()); // unable to make adjacent
+
+    auto he_h = next_halfedge_of(he_g);
+
+    // properly rewire
+    connect_prev_next(he_in, he_out);
+    connect_prev_next(he_g, he_b);
+    connect_prev_next(he_d, he_h);
+}
+
+inline void low_level_api_mutable::remove_face(face_index f_idx) const
+{
+    assert(!is_removed(f_idx));
+
+    auto he_begin = halfedge_of(f_idx);
+    auto he = he_begin;
+    do
+    {
+        assert(face_of(he) == f_idx);
+
+        // set half-edge face to invalid
+        face_of(he) = face_index::invalid();
+
+        // fix outgoing vertex half-edge
+        // (vertex correctly reports is_boundary)
+        outgoing_halfedge_of(from_vertex_of(he)) = he;
+
+        // fix opposite face half-edge
+        auto ohe = opposite(he);
+        auto of = face_of(ohe);
+        if (of.is_valid())
+            halfedge_of(of) = ohe;
+
+        // advance
+        he = next_halfedge_of(he);
+    } while (he != he_begin);
+
+    // mark removed
+    // (at the end!)
+    set_removed(f_idx);
+
+    // bookkeeping
+    m.mRemovedFaces++;
+    m.mCompact = false;
+}
+
+inline void low_level_api_mutable::remove_edge(edge_index e_idx) const
+{
+    auto h_in = halfedge_of(e_idx, 0);
+    auto h_out = halfedge_of(e_idx, 1);
+
+    assert(!is_removed(h_in));
+    assert(!is_removed(h_out));
+
+    auto f0 = face_of(h_in);
+    auto f1 = face_of(h_out);
+
+    // remove adjacent faces
+    if (f0.is_valid())
+        remove_face(f0);
+    if (f1.is_valid())
+        remove_face(f1);
+
+    // rewire vertices
+    auto v_in_to = to_vertex_of(h_in);
+    auto v_out_to = to_vertex_of(h_out);
+
+    auto hi_out_prev = prev_halfedge_of(h_out);
+    auto hi_out_next = next_halfedge_of(h_out);
+
+    auto hi_in_prev = prev_halfedge_of(h_in);
+    auto hi_in_next = next_halfedge_of(h_in);
+
+    // modify vertex if outgoing half-edge is going to be removed
+    auto &v_in_to_out = outgoing_halfedge_of(v_in_to);
+    if (v_in_to_out == h_out)
+    {
+        if (hi_in_next == h_out) // v_in_to becomes isolated
+            v_in_to_out = halfedge_index::invalid();
+        else
+            v_in_to_out = hi_in_next;
+    }
+
+    auto &v_out_to_out = outgoing_halfedge_of(v_out_to);
+    if (v_out_to_out == h_in)
+    {
+        if (hi_out_next == h_in) // v_out_to becomes isolated
+            v_out_to_out = halfedge_index::invalid();
+        else
+            v_out_to_out = hi_out_next;
+    }
+
+    // reqire half-edges
+    connect_prev_next(hi_out_prev, hi_in_next);
+    connect_prev_next(hi_in_prev, hi_out_next);
+
+    // remove half-edges
+    set_removed(e_idx);
+
+    // bookkeeping
+    m.mRemovedHalfedges++;
+    m.mRemovedHalfedges++;
+    m.mCompact = false;
+}
+
+inline void low_level_api_mutable::remove_vertex(vertex_index v_idx) const
+{
+    assert(!is_removed(v_idx));
+
+    // remove all outgoing edges
+    while (!is_isolated(v_idx))
+        remove_edge(edge_of(outgoing_halfedge_of(v_idx)));
+
+    // mark removed
+    set_removed(v_idx);
+
+    // bookkeeping
+    m.mRemovedVertices++;
+    m.mCompact = false;
+}
+
+inline void low_level_api_mutable::fix_boundary_state_of(vertex_index v_idx) const
+{
+    assert(!is_isolated(v_idx));
+
+    auto he_begin = outgoing_halfedge_of(v_idx);
+    auto he = he_begin;
+    do
+    {
+        // if half-edge is boundary, set it
+        if (is_free(he))
+        {
+            outgoing_halfedge_of(v_idx) = he;
+            return;
+        }
+
+        // advance
+        he = next_halfedge_of(opposite(he));
+    } while (he != he_begin);
+}
+
+inline void low_level_api_mutable::fix_boundary_state_of(face_index f_idx) const
+{
+    auto he_begin = halfedge_of(f_idx);
+    auto he = he_begin;
+    do
+    {
+        // if half-edge is boundary, set it
+        if (is_free(opposite(he)))
+        {
+            halfedge_of(f_idx) = he;
+            return;
+        }
+
+        // advance
+        he = next_halfedge_of(he);
+    } while (he != he_begin);
+}
+
+inline void low_level_api_mutable::fix_boundary_state_of_vertices(face_index f_idx) const
+{
+    auto he_begin = halfedge_of(f_idx);
+    auto he = he_begin;
+    do
+    {
+        // fix vertex
+        fix_boundary_state_of(to_vertex_of(he));
+
+        // advance
+        he = next_halfedge_of(he);
+    } while (he != he_begin);
+}
+inline void low_level_api_mutable::set_removed(vertex_index idx) const { outgoing_halfedge_of(idx) = halfedge_index::invalid(); }
+inline void low_level_api_mutable::set_removed(face_index idx) const { halfedge_of(idx) = halfedge_index::invalid(); }
+inline void low_level_api_mutable::set_removed(edge_index idx) const
+{
+    to_vertex_of(halfedge_of(idx, 0)) = vertex_index::invalid();
+    to_vertex_of(halfedge_of(idx, 1)) = vertex_index::invalid();
+}
+
+inline void low_level_api_mutable::connect_prev_next(halfedge_index prev, halfedge_index next) const
+{
+    next_halfedge_of(prev) = next;
+    prev_halfedge_of(next) = prev;
+}
+
+inline vertex_index low_level_api_mutable::face_split(face_index f) const
+{
+    // TODO: can be made more performant
+
+    auto h_begin = halfedge_of(f);
+
+    // remove face
+    remove_face(f);
+
+    // add vertex
+    vertex_index vs[3];
+    vs[0] = add_vertex();
+
+    // add triangles
+    auto h = h_begin;
+    do
+    {
+        vs[1] = from_vertex_of(h);
+        vs[2] = to_vertex_of(h);
+
+        add_face(vs, 3);
+
+        // NOTE: add_face inserted a new halfedge
+        h = next_halfedge_of(opposite(next_halfedge_of(h)));
+    } while (h != h_begin);
+
+    return vs[0];
+}
+
+inline vertex_index low_level_api_mutable::edge_split(edge_index e) const
+{
+    auto h0 = halfedge_of(e, 0);
+    auto h1 = halfedge_of(e, 1);
+
+    auto v0 = to_vertex_of(h0);
+    auto v1 = to_vertex_of(h1);
+    auto f0 = face_of(h0);
+    auto f1 = face_of(h1);
+
+    // add vertex
+    auto v = add_vertex();
+
+    // add two new edges
+    auto e1 = alloc_edge();
+    auto e2 = alloc_edge();
+    auto e1h0 = halfedge_of(e1, 0);
+    auto e1h1 = halfedge_of(e1, 1);
+    auto e2h0 = halfedge_of(e2, 0);
+    auto e2h1 = halfedge_of(e2, 1);
+
+    // rewire edges
+    auto h0_prev = prev_halfedge_of(h0);
+    auto h0_next = next_halfedge_of(h0);
+    auto h1_prev = prev_halfedge_of(h1);
+    auto h1_next = next_halfedge_of(h1);
+
+    face_of(e1h0) = f0;
+    face_of(e2h0) = f0;
+    face_of(e1h1) = f1;
+    face_of(e2h1) = f1;
+
+    to_vertex_of(e1h0) = v0;
+    to_vertex_of(e2h0) = v;
+    to_vertex_of(e1h1) = v;
+    to_vertex_of(e2h1) = v1;
+
+    connect_prev_next(h0_prev, e2h0);
+    connect_prev_next(e2h0, e1h0);
+    connect_prev_next(e1h0, h0_next);
+
+    connect_prev_next(h1_prev, e1h1);
+    connect_prev_next(e1h1, e2h1);
+    connect_prev_next(e2h1, h1_next);
+
+    // rewire vertices
+    auto &v0_out = outgoing_halfedge_of(v0);
+    auto &v1_out = outgoing_halfedge_of(v1);
+    if (v0_out == h1)
+        v0_out = e1h1;
+    if (v1_out == h0)
+        v1_out = e2h0;
+
+    // rewire faces
+    if (f0.is_valid())
+    {
+        auto &f0_h = halfedge_of(f0);
+        if (f0_h == h0)
+            f0_h = e1h0;
+    }
+    if (f1.is_valid())
+    {
+        auto &f1_h = halfedge_of(f1);
+        if (f1_h == h1)
+            f1_h = e2h1;
+    }
+
+    // remove edge
+    set_removed(e);
+    m.mRemovedHalfedges += 2;
+    m.mCompact = false;
+
+    return v;
+}
+
+inline vertex_index low_level_api_mutable::halfedge_split(halfedge_index h) const
+{
+    // add vertex and edge
+    auto v = add_vertex();
+    auto e = alloc_edge();
+
+    auto h0 = h;
+    auto h1 = opposite(h);
+    auto h2 = halfedge_of(e, 0);
+    auto h3 = halfedge_of(e, 1);
+
+    auto v0 = to_vertex_of(h0);
+    auto v1 = to_vertex_of(h1);
+
+    // rewire edges
+    auto h0_next = next_halfedge_of(h0);
+    auto h1_prev = prev_halfedge_of(h1);
+
+    auto f0 = face_of(h0);
+    auto f1 = face_of(h1);
+
+    face_of(h2) = f0;
+    face_of(h3) = f1;
+
+    to_vertex_of(h0) = v;
+    to_vertex_of(h1) = v1; //< already there
+    to_vertex_of(h2) = v0;
+    to_vertex_of(h3) = v;
+
+    connect_prev_next(h0, h2);
+    connect_prev_next(h2, h0_next);
+
+    connect_prev_next(h1_prev, h3);
+    connect_prev_next(h3, h1);
+
+    // rewire vertices
+    auto &v0_out = outgoing_halfedge_of(v0);
+    if (v0_out == h1)
+        v0_out = h3;
+
+    outgoing_halfedge_of(v) = is_free(h1) ? h1 : h2; // boundary
+
+    // rewire faces
+    // -> already ok
+
+    return v;
+}
+
+inline face_index low_level_api_mutable::face_fill(halfedge_index h) const
+{
+    assert(is_boundary(h));
+
+    auto f = alloc_face();
+
+    halfedge_of(f) = h;
+
+    auto h_begin = h;
+    do
+    {
+        // set face
+        face_of(h) = f;
+
+        // fix face boundary
+        if (is_boundary(opposite(h)))
+            halfedge_of(f) = h;
+
+        // fix adj face boundary
+        auto adj_face = opposite_face_of(h);
+        if (adj_face.is_valid())
+            fix_boundary_state_of(adj_face);
+
+        // advance
+        h = next_halfedge_of(h);
+    } while (h != h_begin);
+
+    // fix vertex boundaries
+    fix_boundary_state_of_vertices(f);
+
+    return f;
+}
+
+inline void low_level_api_mutable::halfedge_attach(halfedge_index h, vertex_index v) const
+{
+    assert(is_isolated(v));
+
+    auto h_next = next_halfedge_of(h);
+    auto v_to = to_vertex_of(h);
+
+    auto f = face_of(h);
+
+    auto e = alloc_edge();
+    auto h0 = halfedge_of(e, 0);
+    auto h1 = halfedge_of(e, 1);
+
+    face_of(h0) = f;
+    to_vertex_of(h0) = v;
+
+    face_of(h1) = f;
+    to_vertex_of(h1) = v_to;
+
+    outgoing_halfedge_of(v) = h1;
+
+    connect_prev_next(h, h0);
+    connect_prev_next(h0, h1);
+    connect_prev_next(h1, h_next);
+}
+
+inline void low_level_api_mutable::vertex_collapse(vertex_index v) const
+{
+    // isolated vertices are just removed
+    if (is_isolated(v))
+    {
+        remove_vertex(v);
+    }
+    // boundary vertices are special
+    else if (is_boundary(v))
+    {
+        assert(0 && "not implemented");
+    }
+    else // interior vertex
+    {
+        auto h_begin = next_halfedge_of(outgoing_halfedge_of(v));
+
+        remove_vertex(v);
+
+        assert(is_boundary(h_begin));
+
+        // TODO: optimize
+        std::vector<halfedge_index> hs;
+        auto h = h_begin;
+        do
+        {
+            // add half-edge ring
+            hs.push_back(h);
+
+            // advance
+            h = next_halfedge_of(h);
+        } while (h != h_begin);
+
+        // add face
+        add_face(hs.data(), (int)hs.size());
+    }
+}
+
+inline void low_level_api_mutable::halfedge_collapse(halfedge_index h) const
+{
+    // TODO: collapse half-edge
+    // preserve adjacent non-triangles
+
+    assert(0 && "not implemented");
+}
+
+inline void low_level_api_mutable::edge_rotate_next(edge_index e) const
+{
+    assert(!is_boundary(e) && "does not work on boundaries");
+    assert(m.handle_of(e).vertexA().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
+    assert(m.handle_of(e).vertexB().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
+
+    auto h0 = halfedge_of(e, 0);
+    auto h1 = halfedge_of(e, 1);
+
+    auto h0_next = next_halfedge_of(h0);
+    auto h0_prev = prev_halfedge_of(h0);
+    auto h1_next = next_halfedge_of(h1);
+    auto h1_prev = prev_halfedge_of(h1);
+
+    auto h0_next_next = next_halfedge_of(h0_next);
+    auto h1_next_next = next_halfedge_of(h1_next);
+
+    // fix vertices
+    auto &v0_out = outgoing_halfedge_of(to_vertex_of(h0));
+    if (v0_out == h1)
+        v0_out = h0_next;
+    auto &v1_out = outgoing_halfedge_of(to_vertex_of(h1));
+    if (v1_out == h0)
+        v1_out = h1_next;
+
+    // fix faces
+    halfedge_of(face_of(h0)) = h0;
+    halfedge_of(face_of(h1)) = h1;
+
+    // fix half-edges
+    to_vertex_of(h0) = to_vertex_of(h0_next);
+    to_vertex_of(h1) = to_vertex_of(h1_next);
+    face_of(h0_next) = face_of(h1);
+    face_of(h1_next) = face_of(h0);
+
+    // move to next
+    connect_prev_next(h1_prev, h0_next);
+    connect_prev_next(h0_prev, h1_next);
+
+    connect_prev_next(h0_next, h1);
+    connect_prev_next(h1_next, h0);
+
+    connect_prev_next(h0, h0_next_next);
+    connect_prev_next(h1, h1_next_next);
+
+    // fix boundary state
+    fix_boundary_state_of(face_of(h0));
+    fix_boundary_state_of(face_of(h1));
+}
+
+inline void low_level_api_mutable::edge_rotate_prev(edge_index e) const
+{
+    assert(!is_boundary(e) && "does not work on boundaries");
+    assert(m.handle_of(e).vertexA().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
+    assert(m.handle_of(e).vertexB().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
+
+    auto h0 = halfedge_of(e, 0);
+    auto h1 = halfedge_of(e, 1);
+
+    auto h0_next = next_halfedge_of(h0);
+    auto h0_prev = prev_halfedge_of(h0);
+    auto h1_next = next_halfedge_of(h1);
+    auto h1_prev = prev_halfedge_of(h1);
+
+    auto h0_prev_prev = prev_halfedge_of(h0_prev);
+    auto h1_prev_prev = prev_halfedge_of(h1_prev);
+
+    // fix vertex
+    auto &v0_out = outgoing_halfedge_of(to_vertex_of(h0));
+    if (v0_out == h1)
+        v0_out = h0_next;
+    auto &v1_out = outgoing_halfedge_of(to_vertex_of(h1));
+    if (v1_out == h0)
+        v1_out = h1_next;
+
+    // fix faces
+    halfedge_of(face_of(h0)) = h0;
+    halfedge_of(face_of(h1)) = h1;
+
+    // fix half-edge
+    to_vertex_of(h1) = to_vertex_of(h0_prev_prev);
+    to_vertex_of(h0) = to_vertex_of(h1_prev_prev);
+    face_of(h0_prev) = face_of(h1);
+    face_of(h1_prev) = face_of(h0);
+
+    // move to next
+    connect_prev_next(h0_prev, h1_next);
+    connect_prev_next(h1_prev, h0_next);
+
+    connect_prev_next(h1, h0_prev);
+    connect_prev_next(h0, h1_prev);
+
+    connect_prev_next(h0_prev_prev, h0);
+    connect_prev_next(h1_prev_prev, h1);
+
+    // fix boundary state
+    fix_boundary_state_of(face_of(h0));
+    fix_boundary_state_of(face_of(h1));
+}
+
+inline void low_level_api_mutable::halfedge_rotate_next(halfedge_index h) const
+{
+    assert(m.handle_of(h).next().next().next() != h && "does not work for triangles");
+    assert(!m.handle_of(h).edge().is_boundary() && "does not work on boundaries");
+    assert(m.handle_of(h).vertex_to().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
+
+    auto h0 = h;
+    auto h1 = opposite(h);
+
+    auto h0_next = next_halfedge_of(h0);
+    auto h1_prev = prev_halfedge_of(h1);
+    auto h0_next_next = next_halfedge_of(h0_next);
+
+    // fix vertex
+    auto &v_out = outgoing_halfedge_of(to_vertex_of(h0));
+    if (v_out == h1)
+        v_out = h0_next;
+
+    // fix faces
+    halfedge_of(face_of(h0)) = h0;
+    halfedge_of(face_of(h1)) = h1;
+
+    // fix half-edges
+    to_vertex_of(h0) = to_vertex_of(h0_next);
+    face_of(h0_next) = face_of(h1);
+
+    // move to next
+    connect_prev_next(h1_prev, h0_next);
+    connect_prev_next(h0_next, h1);
+    connect_prev_next(h0, h0_next_next);
+
+    // fix boundary state
+    fix_boundary_state_of(face_of(h0));
+    fix_boundary_state_of(face_of(h1));
+}
+
+inline void low_level_api_mutable::halfedge_rotate_prev(halfedge_index h) const
+{
+    assert(m.handle_of(h).prev().prev().prev() != h && "does not work for triangles");
+    assert(!m.handle_of(h).edge().is_boundary() && "does not work on boundaries");
+    assert(m.handle_of(h).vertex_from().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
+
+    auto h0 = h;
+    auto h1 = opposite(h);
+
+    auto h0_prev = prev_halfedge_of(h0);
+    auto h1_next = next_halfedge_of(h1);
+    auto h0_prev_prev = prev_halfedge_of(h0_prev);
+
+    // fix vertex
+    auto &v_out = outgoing_halfedge_of(to_vertex_of(h1));
+    if (v_out == h0)
+        v_out = h1_next;
+
+    // fix faces
+    halfedge_of(face_of(h0)) = h0;
+    halfedge_of(face_of(h1)) = h1;
+
+    // fix half-edge
+    to_vertex_of(h1) = to_vertex_of(h0_prev_prev);
+    face_of(h0_prev) = face_of(h1);
+
+    // move to next
+    connect_prev_next(h0_prev, h1_next);
+    connect_prev_next(h1, h0_prev);
+    connect_prev_next(h0_prev_prev, h0);
+
+    // fix boundary state
+    fix_boundary_state_of(face_of(h0));
+    fix_boundary_state_of(face_of(h1));
+}
+}
diff --git a/src/polymesh/impl/impl_mesh.hh b/src/polymesh/impl/impl_mesh.hh
index 58ec20f..bc83b52 100644
--- a/src/polymesh/impl/impl_mesh.hh
+++ b/src/polymesh/impl/impl_mesh.hh
@@ -6,7 +6,19 @@
 
 namespace polymesh
 {
-inline vertex_index Mesh::add_vertex() { return alloc_vertex(); }
+inline face_index &Mesh::face_of(halfedge_index idx) { return mHalfedgeToFace[(int)idx]; }
+inline vertex_index &Mesh::to_vertex_of(halfedge_index idx) { return mHalfedgeToVertex[(int)idx]; }
+inline halfedge_index &Mesh::next_halfedge_of(halfedge_index idx) { return mHalfedgeToNextHalfedge[(int)idx]; }
+inline halfedge_index &Mesh::prev_halfedge_of(halfedge_index idx) { return mHalfedgeToPrevHalfedge[(int)idx]; }
+inline halfedge_index &Mesh::halfedge_of(face_index idx) { return mFaceToHalfedge[(int)idx]; }
+inline halfedge_index &Mesh::outgoing_halfedge_of(vertex_index idx) { return mVertexToOutgoingHalfedge[(int)idx]; }
+
+inline face_index Mesh::face_of(halfedge_index idx) const { return mHalfedgeToFace[(int)idx]; }
+inline vertex_index Mesh::to_vertex_of(halfedge_index idx) const { return mHalfedgeToVertex[(int)idx]; }
+inline halfedge_index Mesh::next_halfedge_of(halfedge_index idx) const { return mHalfedgeToNextHalfedge[(int)idx]; }
+inline halfedge_index Mesh::prev_halfedge_of(halfedge_index idx) const { return mHalfedgeToPrevHalfedge[(int)idx]; }
+inline halfedge_index Mesh::halfedge_of(face_index idx) const { return mFaceToHalfedge[(int)idx]; }
+inline halfedge_index Mesh::outgoing_halfedge_of(vertex_index idx) const { return mVertexToOutgoingHalfedge[(int)idx]; }
 
 inline vertex_index Mesh::alloc_vertex()
 {
@@ -84,968 +96,6 @@ inline void Mesh::alloc_primitives(int vertices, int faces, int halfedges)
         p->resize(hCnt, false);
 }
 
-inline face_index Mesh::add_face(const vertex_handle *v_handles, int vcnt)
-{
-    mFaceInsertCache.resize(vcnt);
-    for (auto i = 0; i < vcnt; ++i)
-        mFaceInsertCache[i] = add_or_get_halfedge(v_handles[i].idx, v_handles[(i + 1) % vcnt].idx);
-    return add_face(mFaceInsertCache.data(), vcnt);
-}
-
-inline face_index Mesh::add_face(const vertex_index *v_indices, int vcnt)
-{
-    mFaceInsertCache.resize(vcnt);
-    for (auto i = 0; i < vcnt; ++i)
-        mFaceInsertCache[i] = add_or_get_halfedge(v_indices[i], v_indices[(i + 1) % vcnt]);
-    return add_face(mFaceInsertCache.data(), vcnt);
-}
-
-inline face_index Mesh::add_face(const halfedge_handle *half_loop, int vcnt)
-{
-    mFaceInsertCache.resize(vcnt);
-    for (auto i = 0; i < vcnt; ++i)
-        mFaceInsertCache[i] = half_loop[i].idx;
-    return add_face(mFaceInsertCache.data(), vcnt);
-}
-
-inline face_index Mesh::add_face(const halfedge_index *half_loop, int vcnt)
-{
-    assert(vcnt >= 3 && "no support for less-than-triangular faces");
-
-    auto fidx = alloc_face();
-
-    // ensure that half-edges are adjacent at each vertex
-    for (auto i = 0; i < vcnt; ++i)
-    {
-        auto h0 = half_loop[i];
-        auto h1 = half_loop[(i + 1) % vcnt];
-
-        // half-edge must form a chain
-        assert(to_vertex_of(h0) == from_vertex_of(h1) && "half-edges do not form a chain");
-        // half-edge must be free, i.e. allow a new polygon
-        assert(is_free(h0) && "half-edge already contains a face");
-
-        // make them adjacent
-        make_adjacent(h0, h1);
-
-        // link face
-        face_of(h0) = fidx;
-    }
-
-    // fix boundary states
-    for (auto i = 0; i < vcnt; ++i)
-    {
-        auto h = half_loop[i];
-        auto v = to_vertex_of(h);
-        auto f = opposite_face_of(h);
-
-        // fix vertex
-        fix_boundary_state_of(v);
-
-        // fix face
-        if (f.is_valid())
-            fix_boundary_state_of(f);
-    }
-
-    // set up face data
-    halfedge_of(fidx) = half_loop[0];
-
-    // fix new face
-    fix_boundary_state_of(fidx);
-
-    return fidx;
-}
-
-inline edge_index Mesh::add_or_get_edge(vertex_index v_from, vertex_index v_to)
-{
-    assert(v_from != v_to);
-
-    // already exists?
-    auto he = find_halfedge(v_from, v_to);
-    if (he.is_valid())
-        return edge_of(he);
-
-    // allocate new
-    auto e = alloc_edge();
-    auto h_from_to = halfedge_of(e, 0);
-    auto h_to_from = halfedge_of(e, 1);
-
-    // setup data (self-connected edge)
-    to_vertex_of(h_from_to) = v_to;
-    to_vertex_of(h_to_from) = v_from;
-    next_halfedge_of(h_from_to) = h_to_from;
-    next_halfedge_of(h_to_from) = h_from_to;
-
-    // link from vertex
-    if (is_isolated(v_from))
-        outgoing_halfedge_of(v_from) = h_from_to;
-    else
-    {
-        auto from_in = find_free_incident(v_from);
-        assert(from_in.is_valid() && "vertex is already fully connected");
-
-        auto from_out = next_halfedge_of(from_in);
-
-        connect_prev_next(from_in, h_from_to);
-        connect_prev_next(h_to_from, from_out);
-    }
-
-    // link to vertex
-    if (is_isolated(v_to))
-        outgoing_halfedge_of(v_to) = h_to_from;
-    else
-    {
-        auto to_in = find_free_incident(v_to);
-        assert(to_in.is_valid() && "vertex is already fully connected");
-
-        auto to_out = next_halfedge_of(to_in);
-
-        connect_prev_next(to_in, h_to_from);
-        connect_prev_next(h_from_to, to_out);
-    }
-
-    return e;
-}
-
-inline halfedge_index Mesh::add_or_get_halfedge(vertex_index v_from, vertex_index v_to)
-{
-    auto e = add_or_get_edge(v_from, v_to);
-    auto h0 = halfedge_of(e, 0);
-    auto h1 = halfedge_of(e, 1);
-    return to_vertex_of(h0) == v_to ? h0 : h1;
-}
-
-inline edge_index Mesh::add_or_get_edge(halfedge_index h_from, halfedge_index h_to)
-{
-    assert(h_from != h_to);
-
-    auto v_from = to_vertex_of(h_from);
-    auto v_to = to_vertex_of(h_to);
-
-    auto ex_he = find_halfedge(v_from, v_to);
-    if (ex_he.is_valid())
-    {
-        assert(prev_halfedge_of(ex_he) == h_from && prev_halfedge_of(opposite(ex_he)) == h_to);
-
-        // TODO: Maybe try rewriting an existing halfedge that does NOT yet have the right connection.
-        return edge_of(ex_he);
-    }
-
-    assert(is_free(h_from) && is_free(h_to) && "Cannot insert into a face");
-
-    // allocate new
-    auto e = alloc_edge();
-    auto h_from_to = halfedge_of(e, 0);
-    auto h_to_from = halfedge_of(e, 1);
-
-    // setup data (self-connected edge)
-    to_vertex_of(h_from_to) = v_to;
-    to_vertex_of(h_to_from) = v_from;
-
-    // Link from side
-    auto h_from_next = next_halfedge_of(h_from);
-
-    connect_prev_next(h_from, h_from_to);
-    connect_prev_next(h_from_next, h_to_from);
-
-    // Link to side
-    auto h_to_next = next_halfedge_of(h_to);
-
-    connect_prev_next(h_to, h_to_from);
-    connect_prev_next(h_to_next, h_from_to);
-
-    return e;
-}
-
-inline halfedge_index Mesh::add_or_get_halfedge(halfedge_index h_from, halfedge_index h_to)
-{
-    auto e = add_or_get_edge(h_from, h_to);
-    auto h0 = halfedge_of(e, 0);
-    auto h1 = halfedge_of(e, 1);
-    return next_halfedge_of(h_from) == h0 ? h0 : h1;
-}
-
-inline void Mesh::make_adjacent(halfedge_index he_in, halfedge_index he_out)
-{
-    // see http://kaba.hilvi.org/homepage/blog/halfedge/halfedge.htm ::makeAdjacent
-
-    auto he_b = next_halfedge_of(he_in);
-    auto he_d = prev_halfedge_of(he_out);
-
-    // already correct
-    if (he_b == he_out)
-        return;
-
-    // find free half-edge after `out` but before `in`
-    auto he_g = find_free_incident(opposite(he_out), he_in);
-    assert(he_g.is_valid()); // unable to make adjacent
-
-    auto he_h = next_halfedge_of(he_g);
-
-    // properly rewire
-    connect_prev_next(he_in, he_out);
-    connect_prev_next(he_g, he_b);
-    connect_prev_next(he_d, he_h);
-}
-
-inline void Mesh::remove_face(face_index f_idx)
-{
-    assert(!is_removed(f_idx));
-
-    auto he_begin = halfedge_of(f_idx);
-    auto he = he_begin;
-    do
-    {
-        assert(face_of(he) == f_idx);
-
-        // set half-edge face to invalid
-        face_of(he) = face_index::invalid();
-
-        // fix outgoing vertex half-edge
-        // (vertex correctly reports is_boundary)
-        outgoing_halfedge_of(from_vertex_of(he)) = he;
-
-        // fix opposite face half-edge
-        auto ohe = opposite(he);
-        auto of = face_of(ohe);
-        if (of.is_valid())
-            halfedge_of(of) = ohe;
-
-        // advance
-        he = next_halfedge_of(he);
-    } while (he != he_begin);
-
-    // mark removed
-    // (at the end!)
-    set_removed(f_idx);
-
-    // bookkeeping
-    mRemovedFaces++;
-    mCompact = false;
-}
-
-inline void Mesh::remove_edge(edge_index e_idx)
-{
-    auto h_in = halfedge_of(e_idx, 0);
-    auto h_out = halfedge_of(e_idx, 1);
-
-    assert(!is_removed(h_in));
-    assert(!is_removed(h_out));
-
-    auto f0 = face_of(h_in);
-    auto f1 = face_of(h_out);
-
-    // remove adjacent faces
-    if (f0.is_valid())
-        remove_face(f0);
-    if (f1.is_valid())
-        remove_face(f1);
-
-    // rewire vertices
-    auto v_in_to = to_vertex_of(h_in);
-    auto v_out_to = to_vertex_of(h_out);
-
-    auto hi_out_prev = prev_halfedge_of(h_out);
-    auto hi_out_next = next_halfedge_of(h_out);
-
-    auto hi_in_prev = prev_halfedge_of(h_in);
-    auto hi_in_next = next_halfedge_of(h_in);
-
-    // modify vertex if outgoing half-edge is going to be removed
-    auto &v_in_to_out = outgoing_halfedge_of(v_in_to);
-    if (v_in_to_out == h_out)
-    {
-        if (hi_in_next == h_out) // v_in_to becomes isolated
-            v_in_to_out = halfedge_index::invalid();
-        else
-            v_in_to_out = hi_in_next;
-    }
-
-    auto &v_out_to_out = outgoing_halfedge_of(v_out_to);
-    if (v_out_to_out == h_in)
-    {
-        if (hi_out_next == h_in) // v_out_to becomes isolated
-            v_out_to_out = halfedge_index::invalid();
-        else
-            v_out_to_out = hi_out_next;
-    }
-
-    // reqire half-edges
-    connect_prev_next(hi_out_prev, hi_in_next);
-    connect_prev_next(hi_in_prev, hi_out_next);
-
-    // remove half-edges
-    set_removed(e_idx);
-
-    // bookkeeping
-    mRemovedHalfedges++;
-    mRemovedHalfedges++;
-    mCompact = false;
-}
-
-inline void Mesh::remove_vertex(vertex_index v_idx)
-{
-    assert(!is_removed(v_idx));
-
-    // remove all outgoing edges
-    while (!is_isolated(v_idx))
-        remove_edge(edge_of(outgoing_halfedge_of(v_idx)));
-
-    // mark removed
-    set_removed(v_idx);
-
-    // bookkeeping
-    mRemovedVertices++;
-    mCompact = false;
-}
-
-inline void Mesh::fix_boundary_state_of(vertex_index v_idx)
-{
-    assert(!is_isolated(v_idx));
-
-    auto he_begin = outgoing_halfedge_of(v_idx);
-    auto he = he_begin;
-    do
-    {
-        // if half-edge is boundary, set it
-        if (is_free(he))
-        {
-            outgoing_halfedge_of(v_idx) = he;
-            return;
-        }
-
-        // advance
-        he = next_halfedge_of(opposite(he));
-    } while (he != he_begin);
-}
-
-inline void Mesh::fix_boundary_state_of(face_index f_idx)
-{
-    auto he_begin = halfedge_of(f_idx);
-    auto he = he_begin;
-    do
-    {
-        // if half-edge is boundary, set it
-        if (is_free(opposite(he)))
-        {
-            halfedge_of(f_idx) = he;
-            return;
-        }
-
-        // advance
-        he = next_halfedge_of(he);
-    } while (he != he_begin);
-}
-
-inline void Mesh::fix_boundary_state_of_vertices(face_index f_idx)
-{
-    auto he_begin = halfedge_of(f_idx);
-    auto he = he_begin;
-    do
-    {
-        // fix vertex
-        fix_boundary_state_of(to_vertex_of(he));
-
-        // advance
-        he = next_halfedge_of(he);
-    } while (he != he_begin);
-}
-
-inline halfedge_index Mesh::find_free_incident(halfedge_index in_begin, halfedge_index in_end) const
-{
-    assert(to_vertex_of(in_begin) == to_vertex_of(in_end));
-
-    auto he = in_begin;
-    do
-    {
-        assert(to_vertex_of(he) == to_vertex_of(in_end));
-
-        // free? found one!
-        if (is_free(he))
-            return he;
-
-        // next half-edge of vertex
-        he = opposite(next_halfedge_of(he));
-    } while (he != in_end);
-
-    return halfedge_index::invalid();
-}
-
-inline halfedge_index Mesh::find_free_incident(vertex_index v) const
-{
-    auto in_begin = opposite(outgoing_halfedge_of(v));
-    return find_free_incident(in_begin, in_begin);
-}
-
-inline halfedge_index Mesh::find_halfedge(vertex_index from, vertex_index to) const
-{
-    auto he_begin = outgoing_halfedge_of(from);
-    if (!he_begin.is_valid())
-        return halfedge_index::invalid(); // isolated vertex
-
-    auto he = he_begin;
-    do
-    {
-        // found?
-        if (to_vertex_of(he) == to)
-            return he;
-
-        // advance
-        he = next_halfedge_of(opposite(he));
-
-    } while (he != he_begin);
-
-    return halfedge_index::invalid(); // not found
-}
-
-inline bool Mesh::is_boundary(vertex_index idx) const
-{
-    auto oh = outgoing_halfedge_of(idx);
-    return !oh.is_valid() || is_boundary(oh);
-}
-
-inline bool Mesh::is_free(halfedge_index idx) const { return face_of(idx).is_invalid(); }
-inline bool Mesh::is_boundary(halfedge_index idx) const { return is_free(idx); }
-inline bool Mesh::is_boundary(face_index idx) const { return is_free(opposite(halfedge_of(idx))); }
-inline bool Mesh::is_boundary(edge_index idx) const { return is_free(halfedge_of(idx, 0)) || is_free(halfedge_of(idx, 1)); }
-
-inline bool Mesh::is_isolated(vertex_index idx) const { return outgoing_halfedge_of(idx).is_invalid(); }
-inline bool Mesh::is_isolated(edge_index idx) const { return is_free(halfedge_of(idx, 0)) && is_free(halfedge_of(idx, 1)); }
-
-inline bool Mesh::is_removed(vertex_index idx) const { return outgoing_halfedge_of(idx).value == -2; }
-inline bool Mesh::is_removed(face_index idx) const { return halfedge_of(idx).is_invalid(); }
-inline bool Mesh::is_removed(edge_index idx) const { return to_vertex_of(halfedge_of(idx, 0)).is_invalid(); }
-inline bool Mesh::is_removed(halfedge_index idx) const { return to_vertex_of(idx).is_invalid(); }
-
-inline void Mesh::set_removed(vertex_index idx) { outgoing_halfedge_of(idx) = halfedge_index::invalid(); }
-inline void Mesh::set_removed(face_index idx) { halfedge_of(idx) = halfedge_index::invalid(); }
-inline void Mesh::set_removed(edge_index idx)
-{
-    to_vertex_of(halfedge_of(idx, 0)) = vertex_index::invalid();
-    to_vertex_of(halfedge_of(idx, 1)) = vertex_index::invalid();
-}
-
-inline face_index &Mesh::face_of(halfedge_index idx) { return mHalfedgeToFace[(int)idx]; }
-inline vertex_index &Mesh::to_vertex_of(halfedge_index idx) { return mHalfedgeToVertex[(int)idx]; }
-inline halfedge_index &Mesh::next_halfedge_of(halfedge_index idx) { return mHalfedgeToNextHalfedge[(int)idx]; }
-inline halfedge_index &Mesh::prev_halfedge_of(halfedge_index idx) { return mHalfedgeToPrevHalfedge[(int)idx]; }
-inline halfedge_index &Mesh::halfedge_of(face_index idx) { return mFaceToHalfedge[(int)idx]; }
-inline halfedge_index &Mesh::outgoing_halfedge_of(vertex_index idx) { return mVertexToOutgoingHalfedge[(int)idx]; }
-
-inline face_index Mesh::face_of(halfedge_index idx) const { return mHalfedgeToFace[(int)idx]; }
-inline vertex_index Mesh::to_vertex_of(halfedge_index idx) const { return mHalfedgeToVertex[(int)idx]; }
-inline halfedge_index Mesh::next_halfedge_of(halfedge_index idx) const { return mHalfedgeToNextHalfedge[(int)idx]; }
-inline halfedge_index Mesh::prev_halfedge_of(halfedge_index idx) const { return mHalfedgeToPrevHalfedge[(int)idx]; }
-inline halfedge_index Mesh::halfedge_of(face_index idx) const { return mFaceToHalfedge[(int)idx]; }
-inline halfedge_index Mesh::outgoing_halfedge_of(vertex_index idx) const { return mVertexToOutgoingHalfedge[(int)idx]; }
-
-inline halfedge_index Mesh::opposite(halfedge_index he) const { return halfedge_index(he.value ^ 1); }
-inline face_index Mesh::opposite_face_of(halfedge_index he) const { return face_of(opposite(he)); }
-
-inline vertex_index &Mesh::from_vertex_of(halfedge_index idx) { return to_vertex_of(opposite(idx)); }
-inline vertex_index Mesh::from_vertex_of(halfedge_index idx) const { return to_vertex_of(opposite(idx)); }
-
-inline vertex_index Mesh::next_valid_idx_from(vertex_index idx) const
-{
-    auto s = size_all_vertices();
-    auto i = idx;
-    while (i.value < s && is_removed(i))
-        i.value++;
-    return i;
-}
-
-inline vertex_index Mesh::prev_valid_idx_from(vertex_index idx) const
-{
-    auto i = idx;
-    while (i.value >= 0 && is_removed(i))
-        i.value--;
-    return i;
-}
-
-inline edge_index Mesh::next_valid_idx_from(edge_index idx) const
-{
-    auto s = size_all_edges();
-    auto i = idx;
-    while (i.value < s && is_removed(i))
-        i.value++;
-    return i;
-}
-
-inline edge_index Mesh::prev_valid_idx_from(edge_index idx) const
-{
-    auto i = idx;
-    while (i.value >= 0 && is_removed(i))
-        i.value--;
-    return i;
-}
-
-inline face_index Mesh::next_valid_idx_from(face_index idx) const
-{
-    auto s = size_all_faces();
-    auto i = idx;
-    while (i.value < s && is_removed(i))
-        i.value++;
-    return i;
-}
-
-inline face_index Mesh::prev_valid_idx_from(face_index idx) const
-{
-    auto i = idx;
-    while (i.value >= 0 && is_removed(i))
-        i.value--;
-    return i;
-}
-
-inline halfedge_index Mesh::next_valid_idx_from(halfedge_index idx) const
-{
-    auto s = size_all_halfedges();
-    auto i = idx;
-    while (i.value < s && is_removed(i))
-        i.value++;
-    return i;
-}
-
-inline halfedge_index Mesh::prev_valid_idx_from(halfedge_index idx) const
-{
-    auto i = idx;
-    while (i.value >= 0 && is_removed(i))
-        i.value--;
-    return i;
-}
-
-inline void Mesh::connect_prev_next(halfedge_index prev, halfedge_index next)
-{
-    next_halfedge_of(prev) = next;
-    prev_halfedge_of(next) = prev;
-}
-
-inline vertex_index Mesh::face_split(face_index f)
-{
-    // TODO: can be made more performant
-
-    auto h_begin = halfedge_of(f);
-
-    // remove face
-    remove_face(f);
-
-    // add vertex
-    vertex_index vs[3];
-    vs[0] = add_vertex();
-
-    // add triangles
-    auto h = h_begin;
-    do
-    {
-        vs[1] = from_vertex_of(h);
-        vs[2] = to_vertex_of(h);
-
-        add_face(vs, 3);
-
-        // NOTE: add_face inserted a new halfedge
-        h = next_halfedge_of(opposite(next_halfedge_of(h)));
-    } while (h != h_begin);
-
-    return vs[0];
-}
-
-inline vertex_index Mesh::edge_split(edge_index e)
-{
-    auto h0 = halfedge_of(e, 0);
-    auto h1 = halfedge_of(e, 1);
-
-    auto v0 = to_vertex_of(h0);
-    auto v1 = to_vertex_of(h1);
-    auto f0 = face_of(h0);
-    auto f1 = face_of(h1);
-
-    // add vertex
-    auto v = add_vertex();
-
-    // add two new edges
-    auto e1 = alloc_edge();
-    auto e2 = alloc_edge();
-    auto e1h0 = halfedge_of(e1, 0);
-    auto e1h1 = halfedge_of(e1, 1);
-    auto e2h0 = halfedge_of(e2, 0);
-    auto e2h1 = halfedge_of(e2, 1);
-
-    // rewire edges
-    auto h0_prev = prev_halfedge_of(h0);
-    auto h0_next = next_halfedge_of(h0);
-    auto h1_prev = prev_halfedge_of(h1);
-    auto h1_next = next_halfedge_of(h1);
-
-    face_of(e1h0) = f0;
-    face_of(e2h0) = f0;
-    face_of(e1h1) = f1;
-    face_of(e2h1) = f1;
-
-    to_vertex_of(e1h0) = v0;
-    to_vertex_of(e2h0) = v;
-    to_vertex_of(e1h1) = v;
-    to_vertex_of(e2h1) = v1;
-
-    connect_prev_next(h0_prev, e2h0);
-    connect_prev_next(e2h0, e1h0);
-    connect_prev_next(e1h0, h0_next);
-
-    connect_prev_next(h1_prev, e1h1);
-    connect_prev_next(e1h1, e2h1);
-    connect_prev_next(e2h1, h1_next);
-
-    // rewire vertices
-    auto &v0_out = outgoing_halfedge_of(v0);
-    auto &v1_out = outgoing_halfedge_of(v1);
-    if (v0_out == h1)
-        v0_out = e1h1;
-    if (v1_out == h0)
-        v1_out = e2h0;
-
-    // rewire faces
-    if (f0.is_valid())
-    {
-        auto &f0_h = halfedge_of(f0);
-        if (f0_h == h0)
-            f0_h = e1h0;
-    }
-    if (f1.is_valid())
-    {
-        auto &f1_h = halfedge_of(f1);
-        if (f1_h == h1)
-            f1_h = e2h1;
-    }
-
-    // remove edge
-    set_removed(e);
-    mRemovedHalfedges += 2;
-    mCompact = false;
-
-    return v;
-}
-
-inline vertex_index Mesh::halfedge_split(halfedge_index h)
-{
-    // add vertex and edge
-    auto v = add_vertex();
-    auto e = alloc_edge();
-
-    auto h0 = h;
-    auto h1 = opposite(h);
-    auto h2 = halfedge_of(e, 0);
-    auto h3 = halfedge_of(e, 1);
-
-    auto v0 = to_vertex_of(h0);
-    auto v1 = to_vertex_of(h1);
-
-    // rewire edges
-    auto h0_next = next_halfedge_of(h0);
-    auto h1_prev = prev_halfedge_of(h1);
-
-    auto f0 = face_of(h0);
-    auto f1 = face_of(h1);
-
-    face_of(h2) = f0;
-    face_of(h3) = f1;
-
-    to_vertex_of(h0) = v;
-    to_vertex_of(h1) = v1; //< already there
-    to_vertex_of(h2) = v0;
-    to_vertex_of(h3) = v;
-
-    connect_prev_next(h0, h2);
-    connect_prev_next(h2, h0_next);
-
-    connect_prev_next(h1_prev, h3);
-    connect_prev_next(h3, h1);
-
-    // rewire vertices
-    auto &v0_out = outgoing_halfedge_of(v0);
-    if (v0_out == h1)
-        v0_out = h3;
-
-    outgoing_halfedge_of(v) = is_free(h1) ? h1 : h2; // boundary
-
-    // rewire faces
-    // -> already ok
-
-    return v;
-}
-
-inline face_index Mesh::face_fill(halfedge_index h)
-{
-    assert(is_boundary(h));
-
-    auto f = alloc_face();
-
-    halfedge_of(f) = h;
-
-    auto h_begin = h;
-    do
-    {
-        // set face
-        face_of(h) = f;
-
-        // fix face boundary
-        if (is_boundary(opposite(h)))
-            halfedge_of(f) = h;
-
-        // fix adj face boundary
-        auto adj_face = opposite_face_of(h);
-        if (adj_face.is_valid())
-            fix_boundary_state_of(adj_face);
-
-        // advance
-        h = next_halfedge_of(h);
-    } while (h != h_begin);
-
-    // fix vertex boundaries
-    fix_boundary_state_of_vertices(f);
-
-    return f;
-}
-
-inline void Mesh::halfedge_attach(halfedge_index h, vertex_index v)
-{
-    assert(is_isolated(v));
-
-    auto h_next = next_halfedge_of(h);
-    auto v_to = to_vertex_of(h);
-
-    auto f = face_of(h);
-
-    auto e = alloc_edge();
-    auto h0 = halfedge_of(e, 0);
-    auto h1 = halfedge_of(e, 1);
-
-    face_of(h0) = f;
-    to_vertex_of(h0) = v;
-
-    face_of(h1) = f;
-    to_vertex_of(h1) = v_to;
-
-    outgoing_halfedge_of(v) = h1;
-
-    connect_prev_next(h, h0);
-    connect_prev_next(h0, h1);
-    connect_prev_next(h1, h_next);
-}
-
-inline void Mesh::vertex_collapse(vertex_index v)
-{
-    // isolated vertices are just removed
-    if (is_isolated(v))
-    {
-        remove_vertex(v);
-    }
-    // boundary vertices are special
-    else if (is_boundary(v))
-    {
-        assert(0 && "not implemented");
-    }
-    else // interior vertex
-    {
-        auto h_begin = next_halfedge_of(outgoing_halfedge_of(v));
-
-        remove_vertex(v);
-
-        assert(is_boundary(h_begin));
-
-        // TODO: optimize
-        std::vector<halfedge_index> hs;
-        auto h = h_begin;
-        do
-        {
-            // add half-edge ring
-            hs.push_back(h);
-
-            // advance
-            h = next_halfedge_of(h);
-        } while (h != h_begin);
-
-        // add face
-        add_face(hs.data(), (int)hs.size());
-    }
-}
-
-inline void Mesh::halfedge_collapse(halfedge_index h)
-{
-    // TODO: collapse half-edge
-    // preserve adjacent non-triangles
-
-    assert(0 && "not implemented");
-}
-
-inline void Mesh::edge_rotate_next(edge_index e)
-{
-    assert(!handle_of(e).is_boundary() && "does not work on boundaries");
-    assert(handle_of(e).vertexA().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
-    assert(handle_of(e).vertexB().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
-
-    auto h0 = halfedge_of(e, 0);
-    auto h1 = halfedge_of(e, 1);
-
-    auto h0_next = next_halfedge_of(h0);
-    auto h0_prev = prev_halfedge_of(h0);
-    auto h1_next = next_halfedge_of(h1);
-    auto h1_prev = prev_halfedge_of(h1);
-
-    auto h0_next_next = next_halfedge_of(h0_next);
-    auto h1_next_next = next_halfedge_of(h1_next);
-
-    // fix vertices
-    auto &v0_out = outgoing_halfedge_of(to_vertex_of(h0));
-    if (v0_out == h1)
-        v0_out = h0_next;
-    auto &v1_out = outgoing_halfedge_of(to_vertex_of(h1));
-    if (v1_out == h0)
-        v1_out = h1_next;
-
-    // fix faces
-    halfedge_of(face_of(h0)) = h0;
-    halfedge_of(face_of(h1)) = h1;
-
-    // fix half-edges
-    to_vertex_of(h0) = to_vertex_of(h0_next);
-    to_vertex_of(h1) = to_vertex_of(h1_next);
-    face_of(h0_next) = face_of(h1);
-    face_of(h1_next) = face_of(h0);
-
-    // move to next
-    connect_prev_next(h1_prev, h0_next);
-    connect_prev_next(h0_prev, h1_next);
-
-    connect_prev_next(h0_next, h1);
-    connect_prev_next(h1_next, h0);
-
-    connect_prev_next(h0, h0_next_next);
-    connect_prev_next(h1, h1_next_next);
-
-    // fix boundary state
-    fix_boundary_state_of(face_of(h0));
-    fix_boundary_state_of(face_of(h1));
-}
-
-inline void Mesh::edge_rotate_prev(edge_index e)
-{
-    assert(!handle_of(e).is_boundary() && "does not work on boundaries");
-    assert(handle_of(e).vertexA().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
-    assert(handle_of(e).vertexB().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
-
-    auto h0 = halfedge_of(e, 0);
-    auto h1 = halfedge_of(e, 1);
-
-    auto h0_next = next_halfedge_of(h0);
-    auto h0_prev = prev_halfedge_of(h0);
-    auto h1_next = next_halfedge_of(h1);
-    auto h1_prev = prev_halfedge_of(h1);
-
-    auto h0_prev_prev = prev_halfedge_of(h0_prev);
-    auto h1_prev_prev = prev_halfedge_of(h1_prev);
-
-    // fix vertex
-    auto &v0_out = outgoing_halfedge_of(to_vertex_of(h0));
-    if (v0_out == h1)
-        v0_out = h0_next;
-    auto &v1_out = outgoing_halfedge_of(to_vertex_of(h1));
-    if (v1_out == h0)
-        v1_out = h1_next;
-
-    // fix faces
-    halfedge_of(face_of(h0)) = h0;
-    halfedge_of(face_of(h1)) = h1;
-
-    // fix half-edge
-    to_vertex_of(h1) = to_vertex_of(h0_prev_prev);
-    to_vertex_of(h0) = to_vertex_of(h1_prev_prev);
-    face_of(h0_prev) = face_of(h1);
-    face_of(h1_prev) = face_of(h0);
-
-    // move to next
-    connect_prev_next(h0_prev, h1_next);
-    connect_prev_next(h1_prev, h0_next);
-
-    connect_prev_next(h1, h0_prev);
-    connect_prev_next(h0, h1_prev);
-
-    connect_prev_next(h0_prev_prev, h0);
-    connect_prev_next(h1_prev_prev, h1);
-
-    // fix boundary state
-    fix_boundary_state_of(face_of(h0));
-    fix_boundary_state_of(face_of(h1));
-}
-
-inline void Mesh::halfedge_rotate_next(halfedge_index h)
-{
-    assert(handle_of(h).next().next().next() != h && "does not work for triangles");
-    assert(!handle_of(h).edge().is_boundary() && "does not work on boundaries");
-    assert(handle_of(h).vertex_to().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
-
-    auto h0 = h;
-    auto h1 = opposite(h);
-
-    auto h0_next = next_halfedge_of(h0);
-    auto h1_prev = prev_halfedge_of(h1);
-    auto h0_next_next = next_halfedge_of(h0_next);
-
-    // fix vertex
-    auto &v_out = outgoing_halfedge_of(to_vertex_of(h0));
-    if (v_out == h1)
-        v_out = h0_next;
-
-    // fix faces
-    halfedge_of(face_of(h0)) = h0;
-    halfedge_of(face_of(h1)) = h1;
-
-    // fix half-edges
-    to_vertex_of(h0) = to_vertex_of(h0_next);
-    face_of(h0_next) = face_of(h1);
-
-    // move to next
-    connect_prev_next(h1_prev, h0_next);
-    connect_prev_next(h0_next, h1);
-    connect_prev_next(h0, h0_next_next);
-
-    // fix boundary state
-    fix_boundary_state_of(face_of(h0));
-    fix_boundary_state_of(face_of(h1));
-}
-
-inline void Mesh::halfedge_rotate_prev(halfedge_index h)
-{
-    assert(handle_of(h).prev().prev().prev() != h && "does not work for triangles");
-    assert(!handle_of(h).edge().is_boundary() && "does not work on boundaries");
-    assert(handle_of(h).vertex_from().adjacent_vertices().size() > 2 && "does not work on valence <= 2 vertices");
-
-    auto h0 = h;
-    auto h1 = opposite(h);
-
-    auto h0_prev = prev_halfedge_of(h0);
-    auto h1_next = next_halfedge_of(h1);
-    auto h0_prev_prev = prev_halfedge_of(h0_prev);
-
-    // fix vertex
-    auto &v_out = outgoing_halfedge_of(to_vertex_of(h1));
-    if (v_out == h0)
-        v_out = h1_next;
-
-    // fix faces
-    halfedge_of(face_of(h0)) = h0;
-    halfedge_of(face_of(h1)) = h1;
-
-    // fix half-edge
-    to_vertex_of(h1) = to_vertex_of(h0_prev_prev);
-    face_of(h0_prev) = face_of(h1);
-
-    // move to next
-    connect_prev_next(h0_prev, h1_next);
-    connect_prev_next(h1, h0_prev);
-    connect_prev_next(h0_prev_prev, h0);
-
-    // fix boundary state
-    fix_boundary_state_of(face_of(h0));
-    fix_boundary_state_of(face_of(h1));
-}
 
 inline void Mesh::permute_vertices(std::vector<int> const &p)
 {
@@ -1149,6 +199,8 @@ inline void Mesh::compactify()
     if (is_compact())
         return;
 
+    auto ll = low_level_api(this);
+
     // calculate remappings
     int v_cnt = size_all_vertices();
     int f_cnt = size_all_faces();
@@ -1167,25 +219,25 @@ inline void Mesh::compactify()
     std::vector<int> f_old_to_new(f_cnt, -1);
 
     for (auto i = 0; i < v_cnt; ++i)
-        if (!is_removed(vertex_index(i)))
+        if (!ll.is_removed(vertex_index(i)))
         {
             v_old_to_new[i] = (int)v_new_to_old.size();
             v_new_to_old.push_back(i);
         }
 
     for (auto i = 0; i < f_cnt; ++i)
-        if (!is_removed(face_index(i)))
+        if (!ll.is_removed(face_index(i)))
         {
             f_old_to_new[i] = (int)f_new_to_old.size();
             f_new_to_old.push_back(i);
         }
 
     for (auto i = 0; i < e_cnt; ++i)
-        if (!is_removed(edge_index(i)))
+        if (!ll.is_removed(edge_index(i)))
             e_new_to_old.push_back(i);
 
     for (auto i = 0; i < h_cnt; ++i)
-        if (!is_removed(halfedge_index(i)))
+        if (!ll.is_removed(halfedge_index(i)))
         {
             h_old_to_new[i] = (int)h_new_to_old.size();
             h_new_to_old.push_back(i);
diff --git a/src/polymesh/impl/impl_primitive.hh b/src/polymesh/impl/impl_primitive.hh
index f60b142..3e5281d 100644
--- a/src/polymesh/impl/impl_primitive.hh
+++ b/src/polymesh/impl/impl_primitive.hh
@@ -5,22 +5,22 @@
 namespace polymesh
 {
 // primitive::all_size
-inline int primitive<vertex_tag>::all_size(Mesh const &m) { return m.size_all_vertices(); }
-inline int primitive<face_tag>::all_size(Mesh const &m) { return m.size_all_faces(); }
-inline int primitive<edge_tag>::all_size(Mesh const &m) { return m.size_all_edges(); }
-inline int primitive<halfedge_tag>::all_size(Mesh const &m) { return m.size_all_halfedges(); }
+inline int primitive<vertex_tag>::all_size(Mesh const &m) { return low_level_api(m).size_all_vertices(); }
+inline int primitive<face_tag>::all_size(Mesh const &m) { return low_level_api(m).size_all_faces(); }
+inline int primitive<edge_tag>::all_size(Mesh const &m) { return low_level_api(m).size_all_edges(); }
+inline int primitive<halfedge_tag>::all_size(Mesh const &m) { return low_level_api(m).size_all_halfedges(); }
 
 // primitive::valid_size
-inline int primitive<vertex_tag>::valid_size(Mesh const &m) { return m.size_valid_vertices(); }
-inline int primitive<face_tag>::valid_size(Mesh const &m) { return m.size_valid_faces(); }
-inline int primitive<edge_tag>::valid_size(Mesh const &m) { return m.size_valid_edges(); }
-inline int primitive<halfedge_tag>::valid_size(Mesh const &m) { return m.size_valid_halfedges(); }
+inline int primitive<vertex_tag>::valid_size(Mesh const &m) { return low_level_api(m).size_valid_vertices(); }
+inline int primitive<face_tag>::valid_size(Mesh const &m) { return low_level_api(m).size_valid_faces(); }
+inline int primitive<edge_tag>::valid_size(Mesh const &m) { return low_level_api(m).size_valid_edges(); }
+inline int primitive<halfedge_tag>::valid_size(Mesh const &m) { return low_level_api(m).size_valid_halfedges(); }
 
 // primitive::reserve
-inline void primitive<vertex_tag>::reserve(Mesh &m, int capacity) { m.reserve_vertices(capacity); }
-inline void primitive<face_tag>::reserve(Mesh &m, int capacity) { m.reserve_faces(capacity); }
-inline void primitive<edge_tag>::reserve(Mesh &m, int capacity) { m.reserve_edges(capacity); }
-inline void primitive<halfedge_tag>::reserve(Mesh &m, int capacity) { m.reserve_halfedges(capacity); }
+inline void primitive<vertex_tag>::reserve(Mesh &m, int capacity) { low_level_api(m).reserve_vertices(capacity); }
+inline void primitive<face_tag>::reserve(Mesh &m, int capacity) { low_level_api(m).reserve_faces(capacity); }
+inline void primitive<edge_tag>::reserve(Mesh &m, int capacity) { low_level_api(m).reserve_edges(capacity); }
+inline void primitive<halfedge_tag>::reserve(Mesh &m, int capacity) { low_level_api(m).reserve_halfedges(capacity); }
 
 // primitive::all_collection_of
 inline all_vertex_collection primitive<vertex_tag>::all_collection_of(Mesh &m) { return m.all_vertices(); }
diff --git a/src/polymesh/impl/impl_ranges.hh b/src/polymesh/impl/impl_ranges.hh
index f7b145a..da6786f 100644
--- a/src/polymesh/impl/impl_ranges.hh
+++ b/src/polymesh/impl/impl_ranges.hh
@@ -468,43 +468,43 @@ typename primitive<tag>::handle smart_collection<mesh_ptr, tag, iterator>::rando
 template <class iterator>
 void vertex_collection<iterator>::remove(vertex_handle v) const
 {
-    this->mesh->remove_vertex(v.idx);
+    low_level_api(this->mesh).remove_vertex(v.idx);
 }
 
 template <class iterator>
 void face_collection<iterator>::remove(face_handle f) const
 {
-    this->mesh->remove_face(f.idx);
+    low_level_api(this->mesh).remove_face(f.idx);
 }
 
 template <class iterator>
 void edge_collection<iterator>::remove(edge_handle e) const
 {
-    this->mesh->remove_edge(e.idx);
+    low_level_api(this->mesh).remove_edge(e.idx);
 }
 
 template <class iterator>
 void halfedge_collection<iterator>::remove_edge(halfedge_handle h) const
 {
-    this->mesh->remove_edge(this->mesh->edge_of(h.idx));
+    low_level_api(this->mesh).remove_edge(low_level_api(this->mesh).edge_of(h.idx));
 }
 
 template <class iterator>
 void face_collection<iterator>::permute(std::vector<int> const &p) const
 {
-    this->mesh->permute_faces(p);
+    low_level_api(this->mesh).permute_faces(p);
 }
 
 template <class iterator>
 void edge_collection<iterator>::permute(std::vector<int> const &p) const
 {
-    this->mesh->permute_edges(p);
+    low_level_api(this->mesh).permute_edges(p);
 }
 
 template <class iterator>
 void vertex_collection<iterator>::permute(std::vector<int> const &p) const
 {
-    this->mesh->permute_vertices(p);
+    low_level_api(this->mesh).permute_vertices(p);
 }
 
 
@@ -608,19 +608,19 @@ inline all_halfedge_const_collection Mesh::all_halfedges() const
 template <class iterator>
 vertex_handle vertex_collection<iterator>::add() const
 {
-    return this->mesh->handle_of(this->mesh->add_vertex());
+    return this->mesh->handle_of(low_level_api(this->mesh).add_vertex());
 }
 
 template <class iterator>
 face_handle face_collection<iterator>::add(const vertex_handle *v_handles, int vcnt) const
 {
-    return this->mesh->handle_of(this->mesh->add_face(v_handles, vcnt));
+    return this->mesh->handle_of(low_level_api(this->mesh).add_face(v_handles, vcnt));
 }
 
 template <class iterator>
 face_handle face_collection<iterator>::add(const halfedge_handle *half_loop, int vcnt) const
 {
-    return this->mesh->handle_of(this->mesh->add_face(half_loop, vcnt));
+    return this->mesh->handle_of(low_level_api(this->mesh).add_face(half_loop, vcnt));
 }
 
 template <class iterator>
@@ -639,37 +639,37 @@ template <class iterator>
 face_handle face_collection<iterator>::add(vertex_handle v0, vertex_handle v1, vertex_handle v2) const
 {
     halfedge_index hs[3] = {
-        this->mesh->add_or_get_halfedge(v0.idx, v1.idx), //
-        this->mesh->add_or_get_halfedge(v1.idx, v2.idx), //
-        this->mesh->add_or_get_halfedge(v2.idx, v0.idx), //
+        low_level_api(this->mesh).add_or_get_halfedge(v0.idx, v1.idx), //
+        low_level_api(this->mesh).add_or_get_halfedge(v1.idx, v2.idx), //
+        low_level_api(this->mesh).add_or_get_halfedge(v2.idx, v0.idx), //
     };
-    return this->mesh->handle_of(this->mesh->add_face(hs, 3));
+    return this->mesh->handle_of(low_level_api(this->mesh).add_face(hs, 3));
 }
 
 template <class iterator>
 face_handle face_collection<iterator>::add(vertex_handle v0, vertex_handle v1, vertex_handle v2, vertex_handle v3) const
 {
     halfedge_index hs[4] = {
-        this->mesh->add_or_get_halfedge(v0.idx, v1.idx), //
-        this->mesh->add_or_get_halfedge(v1.idx, v2.idx), //
-        this->mesh->add_or_get_halfedge(v2.idx, v3.idx), //
-        this->mesh->add_or_get_halfedge(v3.idx, v0.idx), //
+        low_level_api(this->mesh).add_or_get_halfedge(v0.idx, v1.idx), //
+        low_level_api(this->mesh).add_or_get_halfedge(v1.idx, v2.idx), //
+        low_level_api(this->mesh).add_or_get_halfedge(v2.idx, v3.idx), //
+        low_level_api(this->mesh).add_or_get_halfedge(v3.idx, v0.idx), //
     };
-    return this->mesh->handle_of(this->mesh->add_face(hs, 4));
+    return this->mesh->handle_of(low_level_api(this->mesh).add_face(hs, 4));
 }
 
 template <class iterator>
 face_handle face_collection<iterator>::add(halfedge_handle h0, halfedge_handle h1, halfedge_handle h2) const
 {
     halfedge_index hs[3] = {h0.idx, h1.idx, h2.idx};
-    return this->mesh->handle_of(this->mesh->add_face(hs, 3));
+    return this->mesh->handle_of(low_level_api(this->mesh).add_face(hs, 3));
 }
 
 template <class iterator>
 face_handle face_collection<iterator>::add(halfedge_handle h0, halfedge_handle h1, halfedge_handle h2, halfedge_handle h3) const
 {
     halfedge_index hs[4] = {h0.idx, h1.idx, h2.idx, h3.idx};
-    return this->mesh->handle_of(this->mesh->add_face(hs, 4));
+    return this->mesh->handle_of(low_level_api(this->mesh).add_face(hs, 4));
 }
 
 template <class iterator>
@@ -678,8 +678,8 @@ face_handle face_collection<iterator>::add(const vertex_handle (&v_handles)[N])
 {
     halfedge_index hs[N];
     for (auto i = 0; i < N; ++i)
-        hs[i] = this->mesh->find_halfedge(v_handles[i].idx, v_handles[(i + 1) % N].idx);
-    return this->mesh->handle_of(this->mesh->add_face(hs, N));
+        hs[i] = low_level_api(this->mesh).find_halfedge(v_handles[i].idx, v_handles[(i + 1) % N].idx);
+    return this->mesh->handle_of(low_level_api(this->mesh).add_face(hs, N));
 }
 
 template <class iterator>
@@ -689,108 +689,108 @@ face_handle face_collection<iterator>::add(const halfedge_handle (&half_loop)[N]
     halfedge_index hs[N];
     for (auto i = 0; i < N; ++i)
         hs[i] = half_loop[i].idx;
-    return this->mesh->handle_of(this->mesh->add_face(hs, N));
+    return this->mesh->handle_of(low_level_api(this->mesh).add_face(hs, N));
 }
 
 template <class iterator>
 edge_handle edge_collection<iterator>::add_or_get(vertex_handle v_from, vertex_handle v_to) const
 {
-    return this->mesh->handle_of(this->mesh->add_or_get_edge(v_from.idx, v_to.idx));
+    return this->mesh->handle_of(low_level_api(this->mesh).add_or_get_edge(v_from.idx, v_to.idx));
 }
 
 template <class iterator>
 halfedge_handle halfedge_collection<iterator>::add_or_get(vertex_handle v_from, vertex_handle v_to) const
 {
-    return this->mesh->handle_of(this->mesh->add_or_get_halfedge(v_from.idx, v_to.idx));
+    return this->mesh->handle_of(low_level_api(this->mesh).add_or_get_halfedge(v_from.idx, v_to.idx));
 }
 
 template <class iterator>
 edge_handle edge_collection<iterator>::add_or_get(halfedge_handle h_from, halfedge_handle h_to) const
 {
-    return this->mesh->handle_of(this->mesh->add_or_get_edge(h_from.idx, h_to.idx));
+    return this->mesh->handle_of(low_level_api(this->mesh).add_or_get_edge(h_from.idx, h_to.idx));
 }
 
 template <class iterator>
 halfedge_handle halfedge_collection<iterator>::add_or_get(halfedge_handle h_from, halfedge_handle h_to) const
 {
-    return this->mesh->handle_of(this->mesh->add_or_get_halfedge(h_from.idx, h_to.idx));
+    return this->mesh->handle_of(low_level_api(this->mesh).add_or_get_halfedge(h_from.idx, h_to.idx));
 }
 
 template <class iterator>
 edge_handle edge_collection<iterator>::find(vertex_handle v_from, vertex_handle v_to) const
 {
-    return this->mesh->handle_of(this->mesh->edge_of(this->mesh->find_halfedge(v_from.idx, v_to.idx)));
+    return this->mesh->handle_of(low_level_api(this->mesh).edge_of(low_level_api(this->mesh).find_halfedge(v_from.idx, v_to.idx)));
 }
 
 template <class iterator>
 halfedge_handle halfedge_collection<iterator>::find(vertex_handle v_from, vertex_handle v_to) const
 {
-    return this->mesh->handle_of(this->mesh->find_halfedge(v_from.idx, v_to.idx));
+    return this->mesh->handle_of(low_level_api(this->mesh).find_halfedge(v_from.idx, v_to.idx));
 }
 
 template <class iterator>
 vertex_handle face_collection<iterator>::split(face_handle f) const
 {
-    return this->mesh->handle_of(this->mesh->face_split(f.idx));
+    return this->mesh->handle_of(low_level_api(this->mesh).face_split(f.idx));
 }
 
 template <class iterator>
 face_handle face_collection<iterator>::fill(halfedge_handle h) const
 {
-    return this->mesh->handle_of(this->mesh->face_fill(h.idx));
+    return this->mesh->handle_of(low_level_api(this->mesh).face_fill(h.idx));
 }
 
 template <class iterator>
 vertex_handle edge_collection<iterator>::split(edge_handle e) const
 {
-    return this->mesh->handle_of(this->mesh->edge_split(e.idx));
+    return this->mesh->handle_of(low_level_api(this->mesh).edge_split(e.idx));
 }
 
 template <class iterator>
 void vertex_collection<iterator>::collapse(vertex_handle v) const
 {
-    this->mesh->vertex_collapse(v.idx);
+    low_level_api(this->mesh).vertex_collapse(v.idx);
 }
 
 template <class iterator>
 void halfedge_collection<iterator>::collapse(halfedge_handle h) const
 {
-    this->mesh->halfedge_collapse(h.idx);
+    low_level_api(this->mesh).halfedge_collapse(h.idx);
 }
 
 template <class iterator>
 void edge_collection<iterator>::rotate_next(edge_handle e) const
 {
-    this->mesh->edge_rotate_next(e.idx);
+    low_level_api(this->mesh).edge_rotate_next(e.idx);
 }
 
 template <class iterator>
 void edge_collection<iterator>::rotate_prev(edge_handle e) const
 {
-    this->mesh->edge_rotate_prev(e.idx);
+    low_level_api(this->mesh).edge_rotate_prev(e.idx);
 }
 
 template <class iterator>
 vertex_handle halfedge_collection<iterator>::split(halfedge_handle h) const
 {
-    return this->mesh->handle_of(this->mesh->halfedge_split(h.idx));
+    return this->mesh->handle_of(low_level_api(this->mesh).halfedge_split(h.idx));
 }
 
 template <class iterator>
 void halfedge_collection<iterator>::attach(halfedge_handle h, vertex_handle v) const
 {
-    this->mesh->halfedge_attach(h.idx, v.idx);
+    low_level_api(this->mesh).halfedge_attach(h.idx, v.idx);
 }
 
 template <class iterator>
 void halfedge_collection<iterator>::rotate_next(halfedge_handle h) const
 {
-    this->mesh->halfedge_rotate_next(h.idx);
+    low_level_api(this->mesh).halfedge_rotate_next(h.idx);
 }
 
 template <class iterator>
 void halfedge_collection<iterator>::rotate_prev(halfedge_handle h) const
 {
-    this->mesh->halfedge_rotate_prev(h.idx);
+    low_level_api(this->mesh).halfedge_rotate_prev(h.idx);
 }
 }
diff --git a/src/polymesh/iterators.hh b/src/polymesh/iterators.hh
index dcf97f1..2f18e01 100644
--- a/src/polymesh/iterators.hh
+++ b/src/polymesh/iterators.hh
@@ -20,13 +20,13 @@ struct valid_primitive_iterator
     static const bool is_valid_iterator = true;
 
     valid_primitive_iterator() = default;
-    valid_primitive_iterator(handle_t handle) : handle(handle) { this->handle.idx = handle.mesh->next_valid_idx_from(handle.idx); }
+    valid_primitive_iterator(handle_t handle) : handle(handle) { this->handle.idx = low_level_api(handle.mesh).next_valid_idx_from(handle.idx); }
 
     handle_t operator*() const { return handle; }
     valid_primitive_iterator& operator++()
     {
         ++handle.idx.value;
-        handle.idx = handle.mesh->next_valid_idx_from(handle.idx);
+        handle.idx = low_level_api(handle.mesh).next_valid_idx_from(handle.idx);
         return *this;
     }
     valid_primitive_iterator operator++(int) const
diff --git a/src/polymesh/low_level_api.hh b/src/polymesh/low_level_api.hh
index 4166e6e..c300c5e 100644
--- a/src/polymesh/low_level_api.hh
+++ b/src/polymesh/low_level_api.hh
@@ -1,14 +1,281 @@
 #pragma once
 
+#include <type_traits>
+#include <vector>
+#include "cursors.hh"
+#include "tmp.hh"
+
 namespace polymesh
 {
 class Mesh;
 
-struct low_level_api
+struct low_level_api_const;
+struct low_level_api_mutable;
+
+inline low_level_api_const low_level_api(Mesh const& m);
+inline low_level_api_const low_level_api(Mesh const* m);
+inline low_level_api_mutable low_level_api(Mesh& m);
+inline low_level_api_mutable low_level_api(Mesh* m);
+
+/**
+ * The Low-Level API allows access to the internal datastructure of a Mesh
+ *
+ * CAUTION: should only be used if you know what you do!
+ *
+ * Usage:
+ *      auto ll = low_level_api(myMesh);
+ *      ll.<something>();
+ */
+template <class MeshT>
+struct low_level_api_base
 {
-    Mesh const& m;
+    // primitive access
+public:
+    tmp::ref_if_mut<vertex_index, MeshT> to_vertex_of(halfedge_index idx) const;
+    tmp::ref_if_mut<face_index, MeshT> face_of(halfedge_index idx) const;
+    tmp::ref_if_mut<halfedge_index, MeshT> next_halfedge_of(halfedge_index idx) const;
+    tmp::ref_if_mut<halfedge_index, MeshT> prev_halfedge_of(halfedge_index idx) const;
+    tmp::ref_if_mut<halfedge_index, MeshT> halfedge_of(face_index idx) const;
+    tmp::ref_if_mut<halfedge_index, MeshT> outgoing_halfedge_of(vertex_index idx) const;
+
+    // number of primitives
+public:
+    int size_all_faces() const;
+    int size_all_vertices() const;
+    int size_all_edges() const;
+    int size_all_halfedges() const;
+
+    int size_valid_faces() const;
+    int size_valid_vertices() const;
+    int size_valid_edges() const;
+    int size_valid_halfedges() const;
+
+    // traversal helper
+public:
+    // returns the next valid idx (returns the given one if valid)
+    // NOTE: the result can be invalid if no valid one was found
+    vertex_index next_valid_idx_from(vertex_index idx) const;
+    edge_index next_valid_idx_from(edge_index idx) const;
+    face_index next_valid_idx_from(face_index idx) const;
+    halfedge_index next_valid_idx_from(halfedge_index idx) const;
+    // returns the next valid idx (returns the given one if valid) counting DOWNWARDS
+    vertex_index prev_valid_idx_from(vertex_index idx) const;
+    edge_index prev_valid_idx_from(edge_index idx) const;
+    face_index prev_valid_idx_from(face_index idx) const;
+    halfedge_index prev_valid_idx_from(halfedge_index idx) const;
+
+    // topology helper
+public:
+    /// Returns the opposite of a given valid half-edge
+    halfedge_index opposite(halfedge_index he) const;
+    /// Returns the opposite face of a given valid half-edge
+    face_index opposite_face_of(halfedge_index he) const;
+    /// Returns the from-vertex of a given valid half-edge
+    vertex_index from_vertex_of(halfedge_index idx) const;
+
+    /// finds the next free incoming half-edge around a certain vertex
+    /// starting from in_begin, EXCLUDING in_end (if in_begin == in_end, the whole vertex is searched)
+    /// returns invalid index if no edge is found
+    halfedge_index find_free_incident(halfedge_index in_begin, halfedge_index in_end) const;
+    /// finds a free incident incoming half-edge around a given vertex
+    halfedge_index find_free_incident(vertex_index v) const;
+
+    /// returns half-edge going from `from`, point to `to`
+    /// returns invalid index if not exists
+    halfedge_index find_halfedge(vertex_index from, vertex_index to) const;
+
+    /// returns edge index belonging to a half-edge
+    edge_index edge_of(halfedge_index idx) const { return edge_index(idx.value >> 1); }
+    /// returns a half-edge belonging to an edge
+    halfedge_index halfedge_of(edge_index idx, int i) const { return halfedge_index((idx.value << 1) + i); }
+    /// returns the to-vertex belonging to the halfedge of an edge
+    vertex_index to_vertex_of(edge_index idx, int i) const { return to_vertex_of(halfedge_of(idx, i)); }
+    /// returns the face belonging to the halfedge of an edge
+    face_index face_of(edge_index idx, int i) const { return face_of(halfedge_of(idx, i)); }
+
+    // properties
+public:
+    /// a half-edge is free if it has no adjacent face
+    bool is_free(halfedge_index idx) const;
 
-    low_level_api(Mesh const& m) : m(m) { }
-    low_level_api(Mesh const* m) : m(*m) { }
+    /// a vertex is boundary if any outgoing half-edges is a boundary (isolated vertices are boundary)
+    bool is_boundary(vertex_index idx) const;
+    /// a half-edge is boundary if it has no face
+    bool is_boundary(halfedge_index idx) const;
+    /// an edge is boundary if at least one half-edge is boundary
+    bool is_boundary(edge_index idx) const;
+    /// a face is boundary if at least one half-edge is boundary
+    bool is_boundary(face_index idx) const;
+
+    /// true iff vertex is marked for removal
+    bool is_removed(vertex_index idx) const;
+    /// true iff face is marked for removal
+    bool is_removed(face_index idx) const;
+    /// true iff edge is marked for removal
+    bool is_removed(edge_index idx) const;
+    /// true iff half-edge is marked for removal
+    bool is_removed(halfedge_index idx) const;
+
+    /// vertices without outgoing half-edges are isolated
+    bool is_isolated(vertex_index idx) const;
+    /// edges without faces are isolated
+    bool is_isolated(edge_index idx) const;
+
+protected:
+    MeshT& m;
+
+    low_level_api_base(MeshT& m) : m(m) {}
 };
+
+struct low_level_api_mutable : low_level_api_base<Mesh>
+{
+    // allocation
+public:
+    // reserves a certain number of primitives
+    void reserve_faces(int capacity) const;
+    void reserve_vertices(int capacity) const;
+    void reserve_edges(int capacity) const;
+    void reserve_halfedges(int capacity) const;
+
+    /// Allocates a new vertex
+    vertex_index alloc_vertex() const;
+    /// Allocates a new face
+    face_index alloc_face() const;
+    /// Allocates a new edge
+    edge_index alloc_edge() const;
+    /// Allocates a given amount of vertices, faces, and halfedges
+    /// NOTE: leaves ALL of them in an unspecified state
+    void alloc_primitives(int vertices, int faces, int halfedges) const;
+
+    // adding primitives
+public:
+    /// Adds a single non-connected vertex
+    /// Does NOT invalidate iterators!
+    vertex_index add_vertex() const;
+
+    /// Adds a face consisting of N vertices
+    /// The vertices must already be sorted in CCW order
+    /// (note: trying to add already existing halfedges triggers assertions)
+    face_index add_face(vertex_handle const* v_handles, int vcnt) const;
+    face_index add_face(vertex_index const* v_indices, int vcnt) const;
+    face_index add_face(halfedge_handle const* half_loop, int vcnt) const;
+    face_index add_face(halfedge_index const* half_loop, int vcnt) const;
+
+    /// Adds an edge between two existing, distinct vertices
+    /// if edge already exists, returns it
+    edge_index add_or_get_edge(vertex_index v_from, vertex_index v_to) const;
+    /// Adds an edge between to existing, distinct vertices that are pointed to by two given halfedges
+    /// if edge already exists, returns it
+    edge_index add_or_get_edge(halfedge_index h_from, halfedge_index h_to) const;
+
+    /// same as add_or_get_edge but returns the apattrriate half-edge
+    /// Assures:
+    ///     return_value.from_vertex == v_from
+    ///     return_value.to_vertex == v_to
+    halfedge_index add_or_get_halfedge(vertex_index v_from, vertex_index v_to) const;
+    /// same as add_or_get_edge but returns the apattrriate half-edge
+    /// Assures:
+    ///     return_value.from_vertex == h_from.to_vertex
+    ///     return_value.to_vertex == h_to.to_vertex
+    halfedge_index add_or_get_halfedge(halfedge_index h_from, halfedge_index h_to) const;
+
+    // removing primitives
+public:
+    /// Marks this vertex as removed
+    void set_removed(vertex_index idx) const;
+    /// Marks this face as removed
+    void set_removed(face_index idx) const;
+    /// Marks this edge as removed
+    void set_removed(edge_index idx) const;
+
+    /// removes a face (actually sets the removed status)
+    /// modifies all adjacent vertices so that they correctly report is_boundary true
+    void remove_face(face_index f_idx) const;
+    /// removes both adjacent faces, then removes both half edges
+    void remove_edge(edge_index e_idx) const;
+    /// removes all adjacent edges, then the vertex
+    void remove_vertex(vertex_index v_idx) const;
+
+    // reordering
+public:
+    /// applies an index remapping to all face indices (p[curr_idx] = new_idx)
+    void permute_faces(std::vector<int> const& p) const;
+    /// applies an index remapping to all edge (and half-edge) indices (p[curr_idx] = new_idx)
+    void permute_edges(std::vector<int> const& p) const;
+    /// applies an index remapping to all vertices indices (p[curr_idx] = new_idx)
+    void permute_vertices(std::vector<int> const& p) const;
+
+    // topology modification
+public:
+    /// Makes two half-edges adjacent
+    /// Ensures:
+    ///     * he_in->next == he_out
+    ///     * he_out->prev == he_in
+    /// Requires:
+    ///     * he_in->is_free()
+    ///     * he_out->is_free()
+    /// Only works if a free incident half-edge is available
+    void make_adjacent(halfedge_index he_in, halfedge_index he_out) const;
+
+    /// connect two half-edges in a prev-next relation
+    void connect_prev_next(halfedge_index prev, halfedge_index next) const;
+
+    /// splits a face
+    vertex_index face_split(face_index f) const;
+    /// splits an edge
+    vertex_index edge_split(edge_index e) const;
+    /// splits a half-edge
+    vertex_index halfedge_split(halfedge_index h) const;
+
+    /// fills a face
+    face_index face_fill(halfedge_index h) const;
+
+    /// attaches a given vertex to the to-vertex of a given half-edge
+    void halfedge_attach(halfedge_index h, vertex_index v) const;
+
+    /// collapse a vertex
+    void vertex_collapse(vertex_index v) const;
+    /// collapse a half-edge
+    void halfedge_collapse(halfedge_index h) const;
+
+    /// rotates an edge to next
+    void edge_rotate_next(edge_index e) const;
+    /// rotates an edge to prev
+    void edge_rotate_prev(edge_index e) const;
+    /// rotates a half-edge to next
+    void halfedge_rotate_next(halfedge_index h) const;
+    /// rotates a half-edge to prev
+    void halfedge_rotate_prev(halfedge_index h) const;
+
+    // boundary states
+public:
+    /// choses a new outgoing half-edge for a given vertex, prefers boundary ones
+    /// assumes non-isolated vertex
+    void fix_boundary_state_of(vertex_index v_idx) const;
+    /// choses a new half-edge for a given face, prefers boundary ones
+    void fix_boundary_state_of(face_index f_idx) const;
+    /// choses a new half-edge for all vertices of a face, prefers boundary ones
+    void fix_boundary_state_of_vertices(face_index f_idx) const;
+
+public:
+    using low_level_api_base<Mesh>::low_level_api_base;
+
+    friend low_level_api_mutable low_level_api(Mesh& m);
+    friend low_level_api_mutable low_level_api(Mesh* m);
+};
+
+struct low_level_api_const : low_level_api_base<Mesh const>
+{
+public:
+    using low_level_api_base<Mesh const>::low_level_api_base;
+
+    friend low_level_api_const low_level_api(Mesh const& m);
+    friend low_level_api_const low_level_api(Mesh const* m);
+};
+
+// Convenience functions
+inline low_level_api_const low_level_api(Mesh const& m) { return {m}; }
+inline low_level_api_const low_level_api(Mesh const* m) { return {*m}; }
+inline low_level_api_mutable low_level_api(Mesh& m) { return {m}; }
+inline low_level_api_mutable low_level_api(Mesh* m) { return {*m}; }
 }
diff --git a/src/polymesh/tmp.hh b/src/polymesh/tmp.hh
index 7cd237e..d536a47 100644
--- a/src/polymesh/tmp.hh
+++ b/src/polymesh/tmp.hh
@@ -24,5 +24,30 @@ struct result_of
 
 template <class FuncT, class ArgT>
 using result_type_of = typename result_of<FuncT, ArgT>::type;
+
+template <class T>
+using enable_if_const_t = typename std::enable_if<std::is_const<T>::value, std::nullptr_t>::type;
+template <class T>
+using enable_if_mutable_t = typename std::enable_if<!std::is_const<T>::value, std::nullptr_t>::type;
+
+template <bool Condition, class TrueType, class FalseType>
+struct if_then_else;
+
+template <class TrueType, class FalseType>
+struct if_then_else<true, TrueType, FalseType>
+{
+    using result = TrueType;
+};
+template <class TrueType, class FalseType>
+struct if_then_else<false, TrueType, FalseType>
+{
+    using result = FalseType;
+};
+
+template <class TargetT, class TestT>
+using ref_if_mut = typename if_then_else<std::is_const<TestT>::value, TargetT, typename std::add_lvalue_reference<TargetT>::type>::result;
+
+// std::add_lvalue_reference
+// template <class T>
 }
 }
-- 
GitLab