diff --git a/applications/incremental_delaunay/internals/mod.rs b/applications/incremental_delaunay/internals/mod.rs
index 9bc29550a..9cd7c1a81 100644
--- a/applications/incremental_delaunay/internals/mod.rs
+++ b/applications/incremental_delaunay/internals/mod.rs
@@ -135,6 +135,7 @@ pub fn delaunay_box_3d<T: CoordsFloat>(
                 T::from(ly).unwrap(),
                 T::from(lz).unwrap(),
             ])
+            // .enable_vertex_id_cache(true)
             .split_cells(true)
             .build()
             .unwrap()
diff --git a/honeycomb-core/src/cmap/builder/io.rs b/honeycomb-core/src/cmap/builder/io.rs
index d2e132bd5..18188dade 100644
--- a/honeycomb-core/src/cmap/builder/io.rs
+++ b/honeycomb-core/src/cmap/builder/io.rs
@@ -206,6 +206,7 @@ pub fn build_2d_from_cmap_file<T: CoordsFloat>(
 #[allow(clippy::too_many_lines)]
 pub fn build_3d_from_cmap_file<T: CoordsFloat>(
     f: CMapFile,
+    enable_vid_cache: bool,
     manager: AttrStorageManager, // FIXME: find a cleaner solution to populate the manager
 ) -> Result<CMap3<T>, BuilderError> {
     if f.meta.1 != 3 {
@@ -214,7 +215,7 @@ pub fn build_3d_from_cmap_file<T: CoordsFloat>(
             "mismatch between requested dimension and header",
         ));
     }
-    let map = CMap3::new_with_undefined_attributes(f.meta.2, manager);
+    let map = CMap3::new_with_undefined_attributes(f.meta.2, enable_vid_cache, manager);
 
     // putting it in a scope to drop the data
     let betas = f.betas.lines().collect::<Vec<_>>();
diff --git a/honeycomb-core/src/cmap/builder/structure.rs b/honeycomb-core/src/cmap/builder/structure.rs
index 0a64e04ed..46dc583d8 100644
--- a/honeycomb-core/src/cmap/builder/structure.rs
+++ b/honeycomb-core/src/cmap/builder/structure.rs
@@ -67,6 +67,7 @@ pub enum BuilderError {
 pub struct CMapBuilder<const D: usize> {
     builder_kind: BuilderType,
     attributes: AttrStorageManager,
+    enable_vid_cache: bool,
 }
 
 enum BuilderType {
@@ -101,9 +102,12 @@ impl<T: CoordsFloat> Builder<T> for CMapBuilder<3> {
 
     fn build(self) -> Result<Self::MapType, BuilderError> {
         match self.builder_kind {
-            BuilderType::CMap(cfile) => super::io::build_3d_from_cmap_file(cfile, self.attributes),
+            BuilderType::CMap(cfile) => {
+                super::io::build_3d_from_cmap_file(cfile, self.enable_vid_cache, self.attributes)
+            }
             BuilderType::FreeDarts(n_darts) => Ok(CMap3::new_with_undefined_attributes(
                 n_darts,
+                self.enable_vid_cache,
                 self.attributes,
             )),
             BuilderType::Vtk(_vfile) => unimplemented!(),
@@ -119,6 +123,7 @@ impl<const D: usize> CMapBuilder<D> {
         Self {
             builder_kind: BuilderType::FreeDarts(n_darts),
             attributes: other.attributes,
+            enable_vid_cache: false,
         }
     }
 
@@ -128,6 +133,7 @@ impl<const D: usize> CMapBuilder<D> {
         Self {
             builder_kind: BuilderType::FreeDarts(n_darts),
             attributes: AttrStorageManager::default(),
+            enable_vid_cache: false,
         }
     }
 
@@ -147,6 +153,7 @@ impl<const D: usize> CMapBuilder<D> {
         Self {
             builder_kind: BuilderType::CMap(cmap_file),
             attributes: AttrStorageManager::default(),
+            enable_vid_cache: false,
         }
     }
 
@@ -163,6 +170,7 @@ impl<const D: usize> CMapBuilder<D> {
         Self {
             builder_kind: BuilderType::Vtk(vtk_file),
             attributes: AttrStorageManager::default(),
+            enable_vid_cache: false,
         }
     }
 
@@ -184,6 +192,18 @@ impl<const D: usize> CMapBuilder<D> {
         self
     }
 
+    /// Enable usage of an internal vertex ID cache.
+    ///
+    /// By default, vertex IDs are recomputed at each call of the `vertex_id(_tx)` methods. By
+    /// enabling this cache, vertex IDs associated to darts are instead stored in a dedicated
+    /// collection, and updated on (un)sews. This can be useful for algorithm which frequently
+    /// use this data.
+    #[must_use = "unused builder object"]
+    pub fn enable_vertex_id_cache(mut self, enable: bool) -> Self {
+        self.enable_vid_cache = enable;
+        self
+    }
+
     #[allow(clippy::missing_errors_doc)]
     /// Consumes the builder and produce a combinatorial map object.
     ///
diff --git a/honeycomb-core/src/cmap/dim3/basic_ops.rs b/honeycomb-core/src/cmap/dim3/basic_ops.rs
index 23da81e54..1c1941ae3 100644
--- a/honeycomb-core/src/cmap/dim3/basic_ops.rs
+++ b/honeycomb-core/src/cmap/dim3/basic_ops.rs
@@ -181,35 +181,39 @@ impl<T: CoordsFloat> CMap3<T> {
         t: &mut Transaction,
         dart_id: DartIdType,
     ) -> Result<VertexIdType, StmError> {
-        AUXILIARIES.with(|cell| {
-            let (pending, marked) = &mut *cell.borrow_mut();
-            // clear from previous computations
-            pending.clear();
-            marked.clear();
-            // initialize
-            pending.push_front(dart_id);
-            marked.insert(NULL_DART_ID); // we don't want to include the null dart in the orbit
-
-            let mut min = dart_id;
-
-            while let Some(d) = pending.pop_front() {
-                if marked.insert(d) {
-                    min = min.min(d);
-                    let (b0, b2, b3) = (
-                        self.beta_tx::<0>(t, d)?, // ?
-                        self.beta_tx::<2>(t, d)?,
-                        self.beta_tx::<3>(t, d)?,
-                    );
-                    pending.push_back(self.beta_tx::<1>(t, b3)?);
-                    pending.push_back(self.beta_tx::<3>(t, b2)?);
-                    pending.push_back(self.beta_tx::<1>(t, b2)?);
-                    pending.push_back(self.beta_tx::<3>(t, b0)?); // ?
-                    pending.push_back(self.beta_tx::<2>(t, b0)?); // ?
+        if let Some(ref vids) = self.vid_cache {
+            vids[dart_id as usize].read(t)
+        } else {
+            AUXILIARIES.with(|cell| {
+                let (pending, marked) = &mut *cell.borrow_mut();
+                // clear from previous computations
+                pending.clear();
+                marked.clear();
+                // initialize
+                pending.push_front(dart_id);
+                marked.insert(NULL_DART_ID); // we don't want to include the null dart in the orbit
+
+                let mut min = dart_id;
+
+                while let Some(d) = pending.pop_front() {
+                    if marked.insert(d) {
+                        min = min.min(d);
+                        let (b0, b2, b3) = (
+                            self.beta_tx::<0>(t, d)?,
+                            self.beta_tx::<2>(t, d)?,
+                            self.beta_tx::<3>(t, d)?,
+                        );
+                        pending.push_back(self.beta_tx::<1>(t, b3)?);
+                        pending.push_back(self.beta_tx::<3>(t, b2)?);
+                        pending.push_back(self.beta_tx::<1>(t, b2)?);
+                        pending.push_back(self.beta_tx::<3>(t, b0)?);
+                        pending.push_back(self.beta_tx::<2>(t, b0)?);
+                    }
                 }
-            }
 
-            Ok(min)
-        })
+                Ok(min)
+            })
+        }
     }
 
     /// Compute the ID of the edge a given dart is part of.
diff --git a/honeycomb-core/src/cmap/dim3/links/one.rs b/honeycomb-core/src/cmap/dim3/links/one.rs
index 547c4cfc2..3dcacb621 100644
--- a/honeycomb-core/src/cmap/dim3/links/one.rs
+++ b/honeycomb-core/src/cmap/dim3/links/one.rs
@@ -1,6 +1,6 @@
 //! 1D link implementations
 
-use crate::cmap::{CMap3, DartIdType, LinkError, NULL_DART_ID};
+use crate::cmap::{CMap3, DartIdType, LinkError, NULL_DART_ID, OrbitPolicy};
 use crate::geometry::CoordsFloat;
 use crate::stm::{Transaction, TransactionClosureResult, abort, atomically_with_err};
 
@@ -18,6 +18,24 @@ impl<T: CoordsFloat> CMap3<T> {
         if b3_ld != NULL_DART_ID && b3_rd != NULL_DART_ID {
             self.betas.one_link_core(t, b3_rd, b3_ld)?;
         }
+        if let Some(ref vids) = self.vid_cache {
+            let b2_ld = self.beta_tx::<2>(t, ld)?;
+            let ll = b2_ld.max(b3_ld);
+            if ll != NULL_DART_ID {
+                let lvid = vids[ll as usize].read(t)?;
+                let rvid = vids[rd as usize].read(t)?;
+                if lvid != rvid {
+                    let new_vid = lvid.min(rvid);
+                    let mut darts = Vec::with_capacity(16);
+                    for d in self.orbit_tx(t, OrbitPolicy::Vertex, ll) {
+                        darts.push(d?);
+                    }
+                    for d in darts {
+                        vids[d as usize].write(t, new_vid)?;
+                    }
+                }
+            }
+        }
         Ok(())
     }
 
@@ -48,6 +66,38 @@ impl<T: CoordsFloat> CMap3<T> {
             }
             self.betas.one_unlink_core(t, b3_rd)?;
         }
+        if let Some(ref vids) = self.vid_cache {
+            let b2_ld = self.beta_tx::<2>(t, ld)?;
+            let ll = b2_ld.max(b3_ld);
+            if ll != NULL_DART_ID {
+                let mut l_orbit = Vec::with_capacity(16);
+                let mut lvid = ll;
+                for d in self.orbit_tx(t, OrbitPolicy::Vertex, ll) {
+                    let d = d?;
+                    l_orbit.push(d);
+                    if d < lvid {
+                        lvid = d;
+                    }
+                }
+                let mut r_orbit = Vec::with_capacity(16);
+                let mut rvid = rd;
+                for d in self.orbit_tx(t, OrbitPolicy::Vertex, rd) {
+                    let d = d?;
+                    r_orbit.push(d);
+                    if d < rvid {
+                        rvid = d;
+                    }
+                }
+                if lvid != rvid {
+                    for d in l_orbit {
+                        vids[d as usize].write(t, lvid)?;
+                    }
+                    for d in r_orbit {
+                        vids[d as usize].write(t, rvid)?;
+                    }
+                }
+            }
+        }
         Ok(())
     }
 
diff --git a/honeycomb-core/src/cmap/dim3/links/three.rs b/honeycomb-core/src/cmap/dim3/links/three.rs
index ca259d7a5..b6f67d17f 100644
--- a/honeycomb-core/src/cmap/dim3/links/three.rs
+++ b/honeycomb-core/src/cmap/dim3/links/three.rs
@@ -1,6 +1,6 @@
 //! 3D link implementations
 
-use crate::cmap::{CMap3, DartIdType, LinkError, NULL_DART_ID};
+use crate::cmap::{CMap3, DartIdType, LinkError, NULL_DART_ID, OrbitPolicy};
 use crate::geometry::CoordsFloat;
 use crate::stm::{Transaction, TransactionClosureResult, abort, atomically_with_err};
 
@@ -14,15 +14,19 @@ impl<T: CoordsFloat> CMap3<T> {
         rd: DartIdType,
     ) -> TransactionClosureResult<(), LinkError> {
         self.betas.three_link_core(t, ld, rd)?;
+        let mut pairs = Vec::with_capacity(16);
         let (mut lside, mut rside) = (self.beta_tx::<1>(t, ld)?, self.beta_tx::<0>(t, rd)?);
+        pairs.push((lside.max(self.beta_tx::<2>(t, ld)?), rd));
         // while we haven't completed the loop, or reached an end
         while lside != ld && lside != NULL_DART_ID {
+            let (b1l, b2l) = (self.beta_tx::<1>(t, lside)?, self.beta_tx::<2>(t, lside)?);
+            pairs.push((b1l.max(b2l), rside));
             if rside == NULL_DART_ID {
                 // (*)
                 abort(LinkError::AsymmetricalFaces(ld, rd))?;
             }
             self.betas.three_link_core(t, lside, rside)?;
-            (lside, rside) = (self.beta_tx::<1>(t, lside)?, self.beta_tx::<0>(t, rside)?);
+            (lside, rside) = (b1l, self.beta_tx::<0>(t, rside)?);
         }
         // the face was open, so we need to cover the other direction
         // for meshes, we should be working on complete faces at all times,
@@ -34,12 +38,21 @@ impl<T: CoordsFloat> CMap3<T> {
             }
             (lside, rside) = (self.beta_tx::<0>(t, ld)?, self.beta_tx::<1>(t, rd)?);
             while lside != NULL_DART_ID {
+                let (b1l, b2l) = (self.beta_tx::<1>(t, lside)?, self.beta_tx::<2>(t, lside)?);
+                pairs.push((b1l.max(b2l), rside));
                 if rside == NULL_DART_ID {
                     // (*)
                     abort(LinkError::AsymmetricalFaces(ld, rd))?;
                 }
                 self.betas.three_link_core(t, lside, rside)?;
-                (lside, rside) = (self.beta_tx::<0>(t, lside)?, self.beta_tx::<1>(t, rside)?);
+                let b1r = self.beta_tx::<1>(t, rside)?;
+                if b1r == NULL_DART_ID {
+                    let b2r = self.beta_tx::<2>(t, rside)?;
+                    if b2r != NULL_DART_ID {
+                        pairs.push((lside, b2r));
+                    }
+                }
+                (lside, rside) = (self.beta_tx::<0>(t, lside)?, b1r);
             }
         }
         // (*): if we land on NULL on one side, the other side should be NULL as well
@@ -47,6 +60,27 @@ impl<T: CoordsFloat> CMap3<T> {
         //      - we're trying to sew open faces with a different number of darts
         //      - we're trying to sew open faces that are offset by one (or more) dart(s)
         //      in both case, this is way too clunky to be considered valid
+
+        if let Some(ref vids) = self.vid_cache {
+            for (ll, rr) in pairs {
+                if ll == NULL_DART_ID || rr == NULL_DART_ID {
+                    continue;
+                }
+                let lvid = vids[ll as usize].read(t)?;
+                let rvid = vids[rr as usize].read(t)?;
+                if lvid != rvid {
+                    let new_vid = lvid.min(rvid);
+                    let mut darts = Vec::with_capacity(16);
+                    for d in self.orbit_tx(t, OrbitPolicy::Vertex, ll) {
+                        darts.push(d?);
+                    }
+                    for d in darts {
+                        vids[d as usize].write(t, new_vid)?;
+                    }
+                }
+            }
+        }
+
         Ok(())
     }
 
@@ -67,15 +101,19 @@ impl<T: CoordsFloat> CMap3<T> {
         let rd = self.beta_tx::<3>(t, ld)?;
 
         self.betas.three_unlink_core(t, ld)?;
+        let mut pairs = Vec::with_capacity(16);
         let (mut lside, mut rside) = (self.beta_tx::<1>(t, ld)?, self.beta_tx::<0>(t, rd)?);
+        pairs.push((lside.max(self.beta_tx::<2>(t, ld)?), rd));
         // while we haven't completed the loop, or reached an end
         while lside != ld && lside != NULL_DART_ID {
+            let (b1l, b2l) = (self.beta_tx::<1>(t, lside)?, self.beta_tx::<2>(t, lside)?);
+            pairs.push((b1l.max(b2l), rside));
             if lside != self.beta_tx::<3>(t, rside)? {
                 // (*); FIXME: add dedicated err ~LinkError::DivergentStructures ?
                 abort(LinkError::AsymmetricalFaces(ld, rd))?;
             }
             self.betas.three_unlink_core(t, lside)?;
-            (lside, rside) = (self.beta_tx::<1>(t, lside)?, self.beta_tx::<0>(t, rside)?);
+            (lside, rside) = (b1l, self.beta_tx::<0>(t, rside)?);
         }
         // the face was open, so we need to cover the other direction
         // for meshes, we should be working on complete faces at all times,
@@ -87,18 +125,62 @@ impl<T: CoordsFloat> CMap3<T> {
             }
             (lside, rside) = (self.beta_tx::<0>(t, ld)?, self.beta_tx::<1>(t, rd)?);
             while lside != NULL_DART_ID {
+                let (b1l, b2l) = (self.beta_tx::<1>(t, lside)?, self.beta_tx::<2>(t, lside)?);
+                pairs.push((b1l.max(b2l), rside));
                 if lside != self.beta_tx::<3>(t, rside)? {
                     // (*); FIXME: add dedicated err ~LinkError::DivergentStructures ?
                     abort(LinkError::AsymmetricalFaces(ld, rd))?;
                 }
                 assert_eq!(lside, self.beta_tx::<3>(t, rside)?); // (*)
                 self.betas.three_unlink_core(t, lside)?;
+                let b1r = self.beta_tx::<1>(t, rside)?;
+                if b1r == NULL_DART_ID {
+                    let b2r = self.beta_tx::<2>(t, rside)?;
+                    if b2r != NULL_DART_ID {
+                        pairs.push((lside, b2r));
+                    }
+                }
                 (lside, rside) = (self.beta_tx::<0>(t, lside)?, self.beta_tx::<1>(t, rside)?);
             }
         }
         // (*) : this can be changed, but the idea here is to ensure we're unlinking the expected
         //       construct
         // (**): if we land on NULL on one side, the other side should be NULL as well
+
+        if let Some(ref vids) = self.vid_cache {
+            for (ll, rr) in pairs {
+                if ll == NULL_DART_ID || rr == NULL_DART_ID {
+                    continue;
+                }
+                let mut l_orbit = Vec::with_capacity(16);
+                let mut lvid = ll;
+                for d in self.orbit_tx(t, OrbitPolicy::Vertex, ll) {
+                    let d = d?;
+                    l_orbit.push(d);
+                    if d < lvid {
+                        lvid = d;
+                    }
+                }
+                let mut r_orbit = Vec::with_capacity(16);
+                let mut rvid = rr;
+                for d in self.orbit_tx(t, OrbitPolicy::Vertex, rr) {
+                    let d = d?;
+                    r_orbit.push(d);
+                    if d < rvid {
+                        rvid = d;
+                    }
+                }
+                if lvid != rvid {
+                    for d in l_orbit {
+                        vids[d as usize].write(t, lvid)?;
+                    }
+                    for d in r_orbit {
+                        vids[d as usize].write(t, rvid)?;
+                    }
+                }
+            }
+        }
+
         Ok(())
     }
 
diff --git a/honeycomb-core/src/cmap/dim3/links/two.rs b/honeycomb-core/src/cmap/dim3/links/two.rs
index 83921a40b..e020774d8 100644
--- a/honeycomb-core/src/cmap/dim3/links/two.rs
+++ b/honeycomb-core/src/cmap/dim3/links/two.rs
@@ -1,6 +1,6 @@
 //! 2D link implementations
 
-use crate::cmap::{CMap3, DartIdType, LinkError};
+use crate::cmap::{CMap3, DartIdType, LinkError, NULL_DART_ID, OrbitPolicy};
 use crate::geometry::CoordsFloat;
 use crate::stm::{Transaction, TransactionClosureResult, atomically_with_err};
 
@@ -10,10 +10,56 @@ impl<T: CoordsFloat> CMap3<T> {
     pub(crate) fn two_link(
         &self,
         t: &mut Transaction,
-        lhs_dart_id: DartIdType,
-        rhs_dart_id: DartIdType,
+        ld: DartIdType,
+        rd: DartIdType,
     ) -> TransactionClosureResult<(), LinkError> {
-        self.betas.two_link_core(t, lhs_dart_id, rhs_dart_id)
+        self.betas.two_link_core(t, ld, rd)?;
+        if let Some(ref vids) = self.vid_cache {
+            // first vertex
+            {
+                let b1l = self.beta_tx::<1>(t, ld)?;
+                let b3l = self.beta_tx::<3>(t, ld)?;
+                let d = b1l.max(b3l);
+                // one or both is not zero -> two orbits were merged
+                if d != NULL_DART_ID {
+                    let ll = rd;
+                    let rr = d;
+                    let lvid = vids[ll as usize].read(t)?;
+                    let rvid = vids[rr as usize].read(t)?;
+                    if lvid != rvid {
+                        let new_vid = lvid.min(rvid);
+                        let mut darts = Vec::with_capacity(16);
+                        for d in self.orbit_tx(t, OrbitPolicy::Vertex, ll) {
+                            darts.push(d?);
+                        }
+                        for d in darts {
+                            vids[d as usize].write(t, new_vid)?;
+                        }
+                    }
+                }
+            }
+            // second vertex
+            let b1r = self.beta_tx::<1>(t, rd)?;
+            let b3r = self.beta_tx::<3>(t, rd)?;
+            let d = b1r.max(b3r);
+            if d != NULL_DART_ID {
+                let ll = ld;
+                let rr = d;
+                let lvid = vids[ll as usize].read(t)?;
+                let rvid = vids[rr as usize].read(t)?;
+                if lvid != rvid {
+                    let new_vid = lvid.min(rvid);
+                    let mut darts = Vec::with_capacity(16);
+                    for d in self.orbit_tx(t, OrbitPolicy::Vertex, ll) {
+                        darts.push(d?);
+                    }
+                    for d in darts {
+                        vids[d as usize].write(t, new_vid)?;
+                    }
+                }
+            }
+        }
+        Ok(())
     }
 
     /// 2-link operation.
@@ -22,7 +68,7 @@ impl<T: CoordsFloat> CMap3<T> {
         lhs_dart_id: DartIdType,
         rhs_dart_id: DartIdType,
     ) -> Result<(), LinkError> {
-        atomically_with_err(|t| self.betas.two_link_core(t, lhs_dart_id, rhs_dart_id))
+        atomically_with_err(|t| self.two_link(t, lhs_dart_id, rhs_dart_id))
     }
 }
 
@@ -32,13 +78,87 @@ impl<T: CoordsFloat> CMap3<T> {
     pub(crate) fn two_unlink(
         &self,
         t: &mut Transaction,
-        lhs_dart_id: DartIdType,
+        ld: DartIdType,
     ) -> TransactionClosureResult<(), LinkError> {
-        self.betas.two_unlink_core(t, lhs_dart_id)
+        let rd = self.beta_tx::<2>(t, ld)?;
+        self.betas.two_unlink_core(t, ld)?;
+        if let Some(ref vids) = self.vid_cache {
+            let mut l_orbit = Vec::with_capacity(16);
+            let mut r_orbit = Vec::with_capacity(16);
+            // first vertex
+            {
+                let ll = ld;
+                let b1r = self.beta_tx::<1>(t, rd)?;
+                let b3r = self.beta_tx::<3>(t, rd)?;
+                let rr = b1r.max(b3r);
+                if rr != NULL_DART_ID {
+                    let mut lvid = ll;
+                    for d in self.orbit_tx(t, OrbitPolicy::Vertex, ll) {
+                        let d = d?;
+                        l_orbit.push(d);
+                        if d < lvid {
+                            lvid = d;
+                        }
+                    }
+                    let mut rvid = rr;
+                    for d in self.orbit_tx(t, OrbitPolicy::Vertex, rr) {
+                        let d = d?;
+                        r_orbit.push(d);
+                        if d < rvid {
+                            rvid = d;
+                        }
+                    }
+                    if lvid != rvid {
+                        for &d in &l_orbit {
+                            vids[d as usize].write(t, lvid)?;
+                        }
+                        for &d in &r_orbit {
+                            vids[d as usize].write(t, rvid)?;
+                        }
+                    }
+                }
+            }
+            l_orbit.clear();
+            r_orbit.clear();
+            // second vertex
+            {
+                let b1l = self.beta_tx::<1>(t, ld)?;
+                let b3l = self.beta_tx::<3>(t, ld)?;
+                let ll = b1l.max(b3l);
+                let rr = rd;
+                if ll != NULL_DART_ID {
+                    let mut lvid = ll;
+                    for d in self.orbit_tx(t, OrbitPolicy::Vertex, ll) {
+                        let d = d?;
+                        l_orbit.push(d);
+                        if d < lvid {
+                            lvid = d;
+                        }
+                    }
+                    let mut rvid = rr;
+                    for d in self.orbit_tx(t, OrbitPolicy::Vertex, rr) {
+                        let d = d?;
+                        r_orbit.push(d);
+                        if d < rvid {
+                            rvid = d;
+                        }
+                    }
+                    if lvid != rvid {
+                        for d in l_orbit {
+                            vids[d as usize].write(t, lvid)?;
+                        }
+                        for d in r_orbit {
+                            vids[d as usize].write(t, rvid)?;
+                        }
+                    }
+                }
+            }
+        }
+        Ok(())
     }
 
     /// 2-unlink operation.
     pub(crate) fn force_two_unlink(&self, lhs_dart_id: DartIdType) -> Result<(), LinkError> {
-        atomically_with_err(|t| self.betas.two_unlink_core(t, lhs_dart_id))
+        atomically_with_err(|t| self.two_unlink(t, lhs_dart_id))
     }
 }
diff --git a/honeycomb-core/src/cmap/dim3/sews/three.rs b/honeycomb-core/src/cmap/dim3/sews/three.rs
index e923b86af..95696610a 100644
--- a/honeycomb-core/src/cmap/dim3/sews/three.rs
+++ b/honeycomb-core/src/cmap/dim3/sews/three.rs
@@ -19,22 +19,33 @@ impl<T: CoordsFloat> CMap3<T> {
     ) -> TransactionClosureResult<(), SewError> {
         // using these custom orbits, I can get both dart of all sides, directly ordered
         // for the merges
-        let l_face = self
-            .orbit(OrbitPolicy::Custom(&[1, 0]), ld)
-            .min()
-            .expect("E: unreachable");
-        let r_face = self
-            .orbit(OrbitPolicy::Custom(&[0, 1]), rd)
-            .min()
-            .expect("E: unreachable");
+        let mut l_face = ld;
+        for d in self.orbit_tx(t, OrbitPolicy::Custom(&[1, 0]), ld) {
+            let d = d?;
+            if d < l_face {
+                l_face = d;
+            }
+        }
+        let mut r_face = rd;
+        for d in self.orbit_tx(t, OrbitPolicy::Custom(&[1, 0]), rd) {
+            let d = d?;
+            if d < r_face {
+                r_face = d;
+            }
+        }
         let mut edges: Vec<(EdgeIdType, EdgeIdType)> = Vec::with_capacity(10);
         let mut vertices: Vec<(VertexIdType, VertexIdType)> = Vec::with_capacity(10);
 
         // read edge + vertex on the b1ld side. if b0ld == NULL, we need to read the left vertex
-        for (l, r) in self
-            .orbit(OrbitPolicy::Custom(&[1, 0]), ld)
-            .zip(self.orbit(OrbitPolicy::Custom(&[0, 1]), rd))
-        {
+        let mut lside = Vec::with_capacity(16);
+        let mut rside = Vec::with_capacity(16);
+        for l in self.orbit_tx(t, OrbitPolicy::Custom(&[1, 0]), ld) {
+            lside.push(l?);
+        }
+        for r in self.orbit_tx(t, OrbitPolicy::Custom(&[0, 1]), rd) {
+            rside.push(r?);
+        }
+        for (l, r) in lside.into_iter().zip(rside.into_iter()) {
             edges.push((self.edge_id_tx(t, l)?, self.edge_id_tx(t, r)?));
             let b1l = self.beta_tx::<1>(t, l)?;
             let b2l = self.beta_tx::<2>(t, l)?;
@@ -161,14 +172,20 @@ impl<T: CoordsFloat> CMap3<T> {
         try_or_coerce!(self.unlink::<3>(t, ld), SewError);
 
         // faces
-        let l_face = self
-            .orbit(OrbitPolicy::Custom(&[1, 0]), ld)
-            .min()
-            .expect("E: unreachable");
-        let r_face = self
-            .orbit(OrbitPolicy::Custom(&[0, 1]), rd)
-            .min()
-            .expect("E: unreachable");
+        let mut l_face = ld;
+        for d in self.orbit_tx(t, OrbitPolicy::Custom(&[1, 0]), ld) {
+            let d = d?;
+            if d < l_face {
+                l_face = d;
+            }
+        }
+        let mut r_face = rd;
+        for d in self.orbit_tx(t, OrbitPolicy::Custom(&[1, 0]), rd) {
+            let d = d?;
+            if d < r_face {
+                r_face = d;
+            }
+        }
         try_or_coerce!(
             self.attributes.split_attributes(
                 t,
@@ -180,10 +197,15 @@ impl<T: CoordsFloat> CMap3<T> {
             SewError
         );
 
-        for (l, r) in self
-            .orbit(OrbitPolicy::Custom(&[1, 0]), ld)
-            .zip(self.orbit(OrbitPolicy::Custom(&[0, 1]), rd))
-        {
+        let mut lside = Vec::with_capacity(16);
+        let mut rside = Vec::with_capacity(16);
+        for l in self.orbit_tx(t, OrbitPolicy::Custom(&[1, 0]), ld) {
+            lside.push(l?);
+        }
+        for r in self.orbit_tx(t, OrbitPolicy::Custom(&[0, 1]), rd) {
+            rside.push(r?);
+        }
+        for (l, r) in lside.into_iter().zip(rside.into_iter()) {
             // edge
             let (eid_l, eid_r) = (self.edge_id_tx(t, l)?, self.edge_id_tx(t, r)?);
             try_or_coerce!(
diff --git a/honeycomb-core/src/cmap/dim3/sews/two.rs b/honeycomb-core/src/cmap/dim3/sews/two.rs
index c5d100bb0..43aa70c14 100644
--- a/honeycomb-core/src/cmap/dim3/sews/two.rs
+++ b/honeycomb-core/src/cmap/dim3/sews/two.rs
@@ -17,8 +17,8 @@ impl<T: CoordsFloat> CMap3<T> {
         ld: DartIdType,
         rd: DartIdType,
     ) -> TransactionClosureResult<(), SewError> {
-        let b1ld = self.betas[(1, ld)].read(t)?;
-        let b1rd = self.betas[(1, rd)].read(t)?;
+        let b1ld = self.betas[(1, ld)].read(t)?.max(self.beta_tx::<3>(t, ld)?);
+        let b1rd = self.betas[(1, rd)].read(t)?.max(self.beta_tx::<3>(t, rd)?);
         // match (is lhs 1-free, is rhs 1-free)
         match (b1ld == NULL_DART_ID, b1rd == NULL_DART_ID) {
             // trivial case, no update needed
@@ -199,14 +199,14 @@ impl<T: CoordsFloat> CMap3<T> {
         ld: DartIdType,
     ) -> TransactionClosureResult<(), SewError> {
         let rd = self.betas[(2, ld)].read(t)?;
-        let b1ld = self.betas[(1, ld)].read(t)?;
-        let b1rd = self.betas[(1, rd)].read(t)?;
+        let b1ld = self.betas[(1, ld)].read(t)?.max(self.beta_tx::<3>(t, ld)?);
+        let b1rd = self.betas[(1, rd)].read(t)?.max(self.beta_tx::<3>(t, rd)?);
         // match (is lhs 1-free, is rhs 1-free)
         match (b1ld == NULL_DART_ID, b1rd == NULL_DART_ID) {
             (true, true) => {
                 // fetch IDs before topology update
                 // update the topology
-                try_or_coerce!(self.betas.two_unlink_core(t, ld), SewError);
+                try_or_coerce!(self.two_unlink(t, ld), SewError);
                 // split attributes from the old ID to the new ones
                 let (eid_newl, eid_newr) = (self.edge_id_tx(t, ld)?, self.edge_id_tx(t, rd)?);
                 try_or_coerce!(
@@ -222,7 +222,7 @@ impl<T: CoordsFloat> CMap3<T> {
             }
             (true, false) => {
                 // update the topology
-                try_or_coerce!(self.betas.two_unlink_core(t, ld), SewError);
+                try_or_coerce!(self.two_unlink(t, ld), SewError);
                 // split vertex & attributes from the old ID to the new ones
                 let (eid_newl, eid_newr) = (self.edge_id_tx(t, ld)?, self.edge_id_tx(t, rd)?);
                 try_or_coerce!(
@@ -256,7 +256,7 @@ impl<T: CoordsFloat> CMap3<T> {
             }
             (false, true) => {
                 // update the topology
-                try_or_coerce!(self.betas.two_unlink_core(t, ld), SewError);
+                try_or_coerce!(self.two_unlink(t, ld), SewError);
                 // split vertex & attributes from the old ID to the new ones
                 let (eid_newl, eid_newr) = (self.edge_id_tx(t, ld)?, self.edge_id_tx(t, rd)?);
                 try_or_coerce!(
@@ -290,7 +290,7 @@ impl<T: CoordsFloat> CMap3<T> {
             }
             (false, false) => {
                 // update the topology
-                try_or_coerce!(self.betas.two_unlink_core(t, ld), SewError);
+                try_or_coerce!(self.two_unlink(t, ld), SewError);
                 // split vertices & attributes from the old ID to the new ones
                 let (eid_newl, eid_newr) = (self.edge_id_tx(t, ld)?, self.edge_id_tx(t, rd)?);
                 try_or_coerce!(
diff --git a/honeycomb-core/src/cmap/dim3/structure.rs b/honeycomb-core/src/cmap/dim3/structure.rs
index d8909a191..2082eeaf3 100644
--- a/honeycomb-core/src/cmap/dim3/structure.rs
+++ b/honeycomb-core/src/cmap/dim3/structure.rs
@@ -3,13 +3,14 @@
 //! This module contains the main structure definition ([`CMap3`]) as well as its constructor
 //! implementation.
 
+use fast_stm::TVar;
 #[cfg(feature = "par-internals")]
 use rayon::prelude::*;
 
 use crate::{
     attributes::{AttrSparseVec, AttrStorageManager, UnknownAttributeStorage},
     cmap::{
-        DartIdType, DartReleaseError, DartReservationError,
+        DartIdType, DartReleaseError, DartReservationError, VertexIdType,
         components::{betas::BetaFunctions, unused::UnusedDarts},
     },
     geometry::{CoordsFloat, Vertex3},
@@ -29,6 +30,7 @@ pub struct CMap3<T: CoordsFloat> {
     pub(super) unused_darts: UnusedDarts,
     /// Array representation of the beta functions
     pub(super) betas: BetaFunctions<CMAP3_BETA>,
+    pub(super) vid_cache: Option<Vec<TVar<VertexIdType>>>,
 }
 
 unsafe impl<T: CoordsFloat> Send for CMap3<T> {}
@@ -79,12 +81,21 @@ impl<T: CoordsFloat> CMap3<T> {
     /// Creates a new 3D combinatorial map.
     #[allow(unused)]
     #[must_use = "unused return value"]
-    pub(crate) fn new(n_darts: usize) -> Self {
+    pub(crate) fn new(n_darts: usize, enable_vid_cache: bool) -> Self {
         Self {
             attributes: AttrStorageManager::default(),
             vertices: AttrSparseVec::new(n_darts + 1),
             unused_darts: UnusedDarts::new(n_darts + 1),
             betas: BetaFunctions::new(n_darts + 1),
+            vid_cache: if enable_vid_cache {
+                Some(
+                    (0..n_darts as VertexIdType + 1)
+                        .map(|v| TVar::new(v))
+                        .collect(),
+                )
+            } else {
+                None
+            },
         }
     }
 
@@ -95,6 +106,7 @@ impl<T: CoordsFloat> CMap3<T> {
     #[must_use = "unused return value"]
     pub(crate) fn new_with_undefined_attributes(
         n_darts: usize,
+        enable_vid_cache: bool,
         mut attr_storage_manager: AttrStorageManager,
     ) -> Self {
         // extend all storages to the expected length: n_darts + 1 (for the null dart)
@@ -104,6 +116,15 @@ impl<T: CoordsFloat> CMap3<T> {
             vertices: AttrSparseVec::new(n_darts + 1),
             unused_darts: UnusedDarts::new(n_darts + 1),
             betas: BetaFunctions::new(n_darts + 1),
+            vid_cache: if enable_vid_cache {
+                Some(
+                    (0..n_darts as VertexIdType + 1)
+                        .map(|v| TVar::new(v))
+                        .collect(),
+                )
+            } else {
+                None
+            },
         }
     }
 }
@@ -156,6 +177,9 @@ impl<T: CoordsFloat> CMap3<T> {
         self.unused_darts.extend_with(n_darts, unused);
         self.vertices.extend(n_darts);
         self.attributes.extend_storages(n_darts);
+        if let Some(ref mut vids) = self.vid_cache {
+            vids.extend((new_id..new_id + n_darts as VertexIdType).map(|v| TVar::new(v)));
+        }
         new_id
     }
 
@@ -305,6 +329,9 @@ impl<T: CoordsFloat> CMap3<T> {
         }
         self.attributes.clear_attribute_values(t, dart_id)?;
         self.vertices.clear_slot(t, dart_id)?;
+        if let Some(ref vids) = self.vid_cache {
+            vids[dart_id as usize].write(t, dart_id)?;
+        }
         Ok(self.unused_darts[dart_id].replace(t, true)?) // Ok(_?) necessary for err type coercion
     }
 }
diff --git a/honeycomb-core/src/cmap/dim3/tests.rs b/honeycomb-core/src/cmap/dim3/tests.rs
index 9b9994a5d..4e6dd8eda 100644
--- a/honeycomb-core/src/cmap/dim3/tests.rs
+++ b/honeycomb-core/src/cmap/dim3/tests.rs
@@ -231,7 +231,7 @@ fn atomically_rebuild_edge(map: &CMap3<f64>, dart: DartIdType) {
 #[test]
 fn example_test_txtional() {
     // Build a tetrahedron (A)
-    let mut map: CMap3<f64> = CMap3::new(12); // 3*4 darts
+    let mut map: CMap3<f64> = CMap3::new(12, true); // 3*4 darts
 
     // face z- (base)
     let res = atomically_with_err(|t| {
@@ -473,7 +473,7 @@ fn reserve_darts() {
 
 #[test]
 fn remove_vertex_twice() {
-    let map: CMap3<f64> = CMap3::new(4);
+    let map: CMap3<f64> = CMap3::new(4, false);
     assert!(map.force_write_vertex(1, (1.0, 1.0, 1.0)).is_none());
     assert_eq!(map.force_remove_vertex(1), Some(Vertex3(1.0, 1.0, 1.0)));
     assert!(map.force_remove_vertex(1).is_none());
@@ -484,7 +484,7 @@ fn remove_dart_twice() {
     // in its default state, all darts are:
     // - used
     // - free
-    let map: CMap3<f64> = CMap3::new(4);
+    let map: CMap3<f64> = CMap3::new(4, false);
     // set dart 1 as unused
     assert!(!map.release_dart(1).unwrap());
     // set dart 1 as unused, again
@@ -495,7 +495,7 @@ fn remove_dart_twice() {
 
 #[test]
 fn one_sew() {
-    let map: CMap3<f64> = CMap3::new(8);
+    let map: CMap3<f64> = CMap3::new(8, false);
     // map.force_link::<1>(1, 2);
     map.force_link::<1>(2, 3).unwrap();
     map.force_link::<1>(3, 4).unwrap();
@@ -530,7 +530,7 @@ fn one_sew() {
 
 #[test]
 fn three_sew() {
-    let map: CMap3<f64> = CMap3::new(8);
+    let map: CMap3<f64> = CMap3::new(8, true);
     map.force_link::<1>(1, 2).unwrap();
     map.force_link::<1>(2, 3).unwrap();
     map.force_link::<1>(3, 4).unwrap();
@@ -558,7 +558,7 @@ fn three_sew() {
 
 #[test]
 fn two_sew_bad_orientation() {
-    let map: CMap3<f64> = CMap3::new(8);
+    let map: CMap3<f64> = CMap3::new(8, false);
     map.force_link::<1>(1, 2).unwrap();
     map.force_link::<1>(2, 3).unwrap();
     map.force_link::<1>(3, 4).unwrap();
@@ -583,7 +583,7 @@ fn two_sew_bad_orientation() {
 
 #[test]
 fn three_sew_bad_orientation() {
-    let map: CMap3<f64> = CMap3::new(8);
+    let map: CMap3<f64> = CMap3::new(8, true);
     map.force_link::<1>(1, 2).unwrap();
     map.force_link::<1>(2, 3).unwrap();
     map.force_link::<1>(3, 4).unwrap();
@@ -612,7 +612,7 @@ fn three_sew_bad_orientation() {
 fn sew_ordering() {
     loom::model(|| {
         // setup the map
-        let map: CMap3<f64> = CMap3::new(5);
+        let map: CMap3<f64> = CMap3::new(5, false);
         map.force_link::<2>(1, 2).unwrap();
         map.force_link::<1>(4, 5).unwrap();
         map.force_write_vertex(2, Vertex3(1.0, 1.0, 1.0));
@@ -652,7 +652,7 @@ fn sew_ordering() {
 fn sew_ordering_with_txtions() {
     loom::model(|| {
         // setup the map
-        let map: CMap3<f64> = CMap3::new(5);
+        let map: CMap3<f64> = CMap3::new(5, false);
         map.force_link::<2>(1, 2).unwrap();
         map.force_link::<2>(3, 4).unwrap();
         // only one vertex is defined
@@ -762,7 +762,7 @@ impl AttributeBind for Weight {
 fn unsew_ordering() {
     loom::model(|| {
         // setup the map FIXME: use the builder
-        let mut map: CMap3<f64> = CMap3::new(5);
+        let mut map: CMap3<f64> = CMap3::new(5, false);
         map.attributes.add_storage::<Weight>(6);
 
         map.force_link::<2>(1, 2).unwrap();
@@ -805,7 +805,7 @@ fn unsew_ordering() {
 fn unsew_ordering_with_txtions() {
     loom::model(|| {
         // setup the map FIXME: use the builder
-        let mut map: CMap3<f64> = CMap3::new(5);
+        let mut map: CMap3<f64> = CMap3::new(5, false);
         map.attributes.add_storage::<Weight>(6);
 
         let res = atomically_with_err(|t| {
diff --git a/honeycomb-core/src/cmap/dim3/utils.rs b/honeycomb-core/src/cmap/dim3/utils.rs
index e27439dd4..e01005275 100644
--- a/honeycomb-core/src/cmap/dim3/utils.rs
+++ b/honeycomb-core/src/cmap/dim3/utils.rs
@@ -2,7 +2,11 @@
 //!
 //! This module contains utility code for the [`CMap3`] structure.
 
-use crate::cmap::{CMap3, DartIdType};
+use std::collections::HashSet;
+
+use rayon::iter::{IntoParallelIterator, ParallelIterator};
+
+use crate::cmap::{CMap3, DartIdType, OrbitPolicy};
 use crate::geometry::CoordsFloat;
 use crate::stm::atomically;
 
@@ -39,4 +43,18 @@ impl<T: CoordsFloat> CMap3<T> {
             Ok(())
         });
     }
+
+    pub fn update_vertex_id_cache(&self) {
+        if let Some(ref vids) = self.vid_cache {
+            (1..self.n_darts() as DartIdType)
+                .into_par_iter()
+                .for_each(|d| {
+                    let min = self
+                        .orbit(OrbitPolicy::Vertex, d)
+                        .min()
+                        .expect("E: unreachable");
+                    atomically(|t| vids[d as usize].write(t, min));
+                });
+        }
+    }
 }
diff --git a/honeycomb-kernels/src/grid_generation/internals.rs b/honeycomb-kernels/src/grid_generation/internals.rs
index 335309034..234d6ae0f 100644
--- a/honeycomb-kernels/src/grid_generation/internals.rs
+++ b/honeycomb-kernels/src/grid_generation/internals.rs
@@ -239,11 +239,13 @@ pub(crate) fn build_3d_grid<T: CoordsFloat>(
     origin: Vertex3<T>,
     n_cells_per_axis: [usize; 3],
     lengths: [T; 3],
+    enable_vid_cache: bool,
 ) -> CMap3<T> {
     let [n_square_x, n_square_y, n_square_z] = n_cells_per_axis;
     let n_darts = 24 * n_square_x * n_square_y * n_square_z;
 
     let map: CMap3<T> = CMapBuilder::<3>::from_n_darts_and_attributes(n_darts, builder)
+        .enable_vertex_id_cache(enable_vid_cache)
         .build()
         .unwrap();
 
@@ -255,6 +257,8 @@ pub(crate) fn build_3d_grid<T: CoordsFloat>(
             map.set_betas(dart, images);
         });
 
+    map.update_vertex_id_cache();
+
     // place vertices
     (1..=n_darts as DartIdType)
         .into_par_iter()
@@ -430,11 +434,13 @@ pub(crate) fn build_3d_tetgrid<T: CoordsFloat>(
     origin: Vertex3<T>,
     n_cells_per_axis: [usize; 3],
     lengths: [T; 3],
+    enable_vid_cache: bool,
 ) -> CMap3<T> {
     let [n_square_x, n_square_y, n_square_z] = n_cells_per_axis;
     let n_darts = 60 * n_square_x * n_square_y * n_square_z;
 
     let map: CMap3<T> = CMapBuilder::<3>::from_n_darts_and_attributes(n_darts, builder)
+        .enable_vertex_id_cache(enable_vid_cache)
         .build()
         .unwrap();
 
@@ -446,6 +452,8 @@ pub(crate) fn build_3d_tetgrid<T: CoordsFloat>(
             map.set_betas(dart, images);
         });
 
+    map.update_vertex_id_cache();
+
     // place vertices
     (1..=n_darts as DartIdType)
         .into_par_iter()
diff --git a/honeycomb-kernels/src/grid_generation/mod.rs b/honeycomb-kernels/src/grid_generation/mod.rs
index aaf87568f..bf6004a28 100644
--- a/honeycomb-kernels/src/grid_generation/mod.rs
+++ b/honeycomb-kernels/src/grid_generation/mod.rs
@@ -56,6 +56,7 @@ pub struct GridBuilder<const D: usize, T: CoordsFloat> {
     pub(crate) len_per_cell: Option<[T; D]>,
     pub(crate) lens: Option<[T; D]>,
     pub(crate) split_cells: bool,
+    pub(crate) enable_vid_cache: bool,
 }
 
 impl<const D: usize, T: CoordsFloat> Default for GridBuilder<D, T> {
@@ -67,6 +68,7 @@ impl<const D: usize, T: CoordsFloat> Default for GridBuilder<D, T> {
             len_per_cell: None,
             lens: None,
             split_cells: false,
+            enable_vid_cache: false,
         }
     }
 }
@@ -111,6 +113,15 @@ impl<const D: usize, T: CoordsFloat> GridBuilder<D, T> {
         self
     }
 
+    /// Enable vertex ID cache
+    ///
+    /// See [`CMapBuilder`] for more information.
+    #[must_use = "unused builder object"]
+    pub fn enable_vertex_id_cache(mut self, enable: bool) -> Self {
+        self.enable_vid_cache = enable;
+        self
+    }
+
     /// Add the attribute `A` to the attributes the created map will contain.
     ///
     /// # Usage
@@ -391,12 +402,14 @@ impl<T: CoordsFloat> GridBuilder<3, T> {
     /// This method may panic if type casting goes wrong during parameters parsing.
     pub fn build(self) -> Result<CMap3<T>, GridBuilderError> {
         let split = self.split_cells;
-        self.parse_3d().map(|(builder, origin, ns, lens)| {
+        let enable = self.enable_vid_cache;
+        let map = self.parse_3d().map(|(builder, origin, ns, lens)| {
             if split {
-                internals::build_3d_tetgrid(builder, origin, ns, lens)
+                internals::build_3d_tetgrid(builder, origin, ns, lens, enable)
             } else {
-                internals::build_3d_grid(builder, origin, ns, lens)
+                internals::build_3d_grid(builder, origin, ns, lens, enable)
             }
-        })
+        });
+        map
     }
 }