diff --git a/.gitignore b/.gitignore
index 4f946e3a..9f39ec24 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,6 @@
 *.py[cod]
 
-# C extensions
+# C extensions 
 *.so
 *.lo
 *.la
diff --git a/src/tyssue/behaviors/sheet/basic_events.py b/src/tyssue/behaviors/sheet/basic_events.py
index b47ee9a7..d472dd34 100644
--- a/src/tyssue/behaviors/sheet/basic_events.py
+++ b/src/tyssue/behaviors/sheet/basic_events.py
@@ -7,7 +7,8 @@
 import logging
 
 from ...geometry.sheet_geometry import SheetGeometry
-from ...topology.sheet_topology import cell_division
+from ...topology.sheet_topology import cell_division, type1_transition as T1_transition
+from ...topology.base_topology import drop_face
 from ...utils.decorators import face_lookup
 from .actions import (
     decrease,
@@ -248,3 +249,41 @@ def contraction_line_tension(sheet, manager, **kwargs):
         isotropic=True,
         limit=100,
     )
+
+
+def T1Swap(sheet,manager, geom, t1_threshold, multiplier):
+    """
+    A behaviour function of the T1 transition that performs a T1 swap on an edge that is shorter than T1_threshold.
+    """
+    # First, we need to get the joint index over free and east edges,
+    # that is, the indices of edges that is spanning the entire graph without double edges
+    sheet.get_extra_indices()
+    edge_dataframe = sheet.edge_df.loc[sheet.sgle_edges]
+    short_edges = edge_dataframe.loc[(edge_dataframe['length'] < t1_threshold)]
+    # take advantage of the unique ID to help us tracking the edges,
+    # we need to reindex the df in type 1 transition function, otherwise there will be inconsistency between
+    # dataframes and will cause out of index error.
+    unique_list = short_edges['unique_id'].tolist()
+    for ID in unique_list:
+        idx = sheet.idx_lookup(ID,'edge')
+        print(f'Performed T1 swap on edge {idx}')
+        T1_transition(sheet, idx, do_reindex=True, remove_tri_faces=False, multiplier=multiplier)
+        sheet.reset_index()     # removes disconnected vertices and faces
+        geom.update_all(sheet)
+    manager.append(T1Swap, geom= geom, t1_threshold = t1_threshold, multiplier = multiplier)
+
+
+def T2Swap(sheet, manager, crit_area):
+    """
+    A behaviour function of the T2 transition that should be added to the manager during simulation.
+    It removes the face if it is triangular and its area is smaller than crit_area.
+    """
+    face_dataframe = sheet.face_df
+    Face_list = face_dataframe.loc[(face_dataframe['num_sides'] < 4) & (face_dataframe['area'] < crit_area)].index.tolist()
+    for face in Face_list:
+        drop_face(sheet, face)
+    manager.append(T2Swap, crit_area = crit_area)
+
+
+
+
diff --git a/src/tyssue/behaviors/sheet/bilayer_dummy_set.py b/src/tyssue/behaviors/sheet/bilayer_dummy_set.py
new file mode 100644
index 00000000..6c3bf449
--- /dev/null
+++ b/src/tyssue/behaviors/sheet/bilayer_dummy_set.py
@@ -0,0 +1,33 @@
+""" The function in this script auto-controls the edges are active or not for a bilayer tissue sheet."""
+
+def bilayer_dummy_set(sheet):
+    """
+    Set edges as active or dummy for bilayer tissue sheet.
+
+    Parameters
+    ----------
+    sheet : tyssue.Sheet
+        The tissue sheet instance containing face_df and edge_df DataFrames.
+    """
+
+    # Iterate through each edge in the edge DataFrame
+    for i in sheet.edge_df.index:
+        # Check if the edge has an opposite edge (i.e., it's internal)
+        if sheet.edge_df.loc[i, 'opposite'] != -1:
+            # Get the associated cell (face) for this edge
+            associated_cell = sheet.edge_df.loc[i, 'face']
+            # Get the opposite edge index
+            opposite_edge = sheet.edge_df.loc[i, 'opposite']
+            # Get the opposite cell (face) for the opposite edge
+            opposite_cell = sheet.edge_df.loc[opposite_edge, 'face']
+            # If both associated and opposite cells are of class 'STB', set edge as dummy (inactive)
+            if (sheet.face_df.loc[associated_cell, 'cell_class'] == 'STB' and
+                    sheet.face_df.loc[opposite_cell, 'cell_class'] == 'STB'):
+                sheet.edge_df.loc[i, 'is_active'] = 0
+                sheet.edge_df.loc[opposite_edge, 'is_active'] = 0
+            else:
+                # Otherwise, set edge as active
+                sheet.edge_df.loc[i, 'is_active'] = 1
+        else:
+            # For boundary edges, set as active
+            sheet.edge_df.loc[i, 'is_active'] = 1
diff --git a/src/tyssue/behaviors/sheet/cell_activity_events.py b/src/tyssue/behaviors/sheet/cell_activity_events.py
new file mode 100644
index 00000000..f6be305e
--- /dev/null
+++ b/src/tyssue/behaviors/sheet/cell_activity_events.py
@@ -0,0 +1,125 @@
+"""
+This file contains all behaviour functions that models different cellular activities.
+
+"""
+
+import numpy as np
+import pandas as pd
+from ...geometry.planar_geometry import PlanarGeometry
+from ...topology.sheet_topology import cell_division, type1_transition, remove_face
+
+"""cell proliferation behaviour function
+
+A cell proliferation behaviour function does two thing on the selected cell.
+First, it compares the current cell area with an area threshold to see if the cell large enough for division.
+Secondly: if the cell is large enough, a cell division is performed on the cell; alternatively, if the cell area is 
+smaller than the threshold value, then the target area is added by "growth_rate * dt" to expand the cell.
+"""
+# Note: still needs to improve the function, so it controls cell class change.
+def proliferation(sheet, manager, geom, unique_id, crit_area, growth_rate, dt):
+    idx = sheet.idx_lookup(unique_id,'face') # get the current face index from unique_id.
+    if sheet.face_df.loc[idx, "area"] > crit_area:
+        # restore prefered_area
+        sheet.face_df.loc[idx, "prefered_area"] = 1.0
+        # Do division
+        daughter = cell_division(sheet, idx, geom)
+        # Update the topology
+        sheet.reset_index(order=True)
+        # update geometry
+        geom.update_all(sheet)
+        print(f"cell n°{daughter} is born")
+    else:
+        sheet.face_df.loc[idx, "prefered_area"] *= (1 + dt * growth_rate)
+        manager.append(proliferation, geom = geom, unique_id = unique_id, crit_area = crit_area, growth_rate = growth_rate, dt = dt)
+
+
+"""cell fusion behaviour function
+
+A cell fusion behaviour function is used when a CT is fusing into the STB layer. The cell class of the selection cell 
+should become "STB" at the end of the function.
+First, the edges between the selected cell and its STB neighbours are disabled for edge tension term (coefficient = 0).
+
+Secondly, we need to perform a T1 swap on the edge that connects a boundary vertex and the mutual vertex shared between all 
+STB neighbours and the fusing cell; in this way, the newly fused cell would have an edge that is "open" to "outside".
+
+Thirdly, new dynamic parameters need to be updated to ensure consistent physics rule.
+"""
+def fusion(sheet, manager, geom, unique_id):
+    sheet.get_extra_indices()
+    idx = sheet.idx_lookup('face', unique_id)
+    # store the face index of STB neighbours
+    STB_neighbours = sheet.get_neighbors(idx)
+    STB_neighbours = list(STB_neighbours.intersection(set(sheet.face_df.loc[sheet.face_df['cell_class'] == 'STB'].index)))
+    # Find the edges associated with the fusing face index, filter out the boundary edges.
+    internal_edges = sheet.edge_df[(sheet.edge_df['face'] == idx) & (sheet.edge_df['opposite'] != -1)]
+    # We have to update both the internal arrowed edge and its opposite arrowed edge.
+    for ie in internal_edges.index:
+        sheet.edge_df.loc[ie,'is_active'] = 0
+        sheet.edge_df.loc[sheet.edge_df.loc[ie,'opposite'],'is_active'] = 0
+
+    # Find all the boundary vertices in STB neighbours, based on the opposite == -1 value.
+    STB_boundary_edges = sheet.edge_df[
+        (sheet.edge_df['face'].isin(STB_neighbours)) &
+        (sheet.edge_df['opposite'] == -1)
+        ]
+    STB_boundary_verts = pd.unique(STB_boundary_edges[['srce','trgt']].values.ravel())
+    # Next, extract all the vertices belong to the fusing face, the edge connects a boundary vertex and the fusing face
+    # is the edge that we should perform T1 swap on. We can utilise the variable internal_edges.
+    fusing_face_verts = internal_edges['srce']
+    # We only need to looping over the STB neighbours edges, then find edges connects a boundary vertex to a fusing face vertex.
+    STB_edges = sheet.edge_df[sheet.edge_df['face'].isin(STB_neighbours)]
+    matching_edge = STB_edges[
+        (
+                (STB_edges['srce'].isin(fusing_face_verts)) & (STB_edges['trgt'].isin(STB_boundary_verts))
+        ) |
+        (
+                (STB_edges['trgt'].isin(fusing_face_verts)) & (STB_edges['srce'].isin(STB_boundary_verts))
+        )
+    ]
+    if matching_edge.empty:
+        raise ValueError("No matching edge found between fusing face vertices and STB boundary vertices.")
+    first_edge_index = matching_edge.index[0]
+    New_boundary_edge = type1_transition(sheet,first_edge_index,do_reindex=False, remove_tri_faces=False, multiplier=1.5)
+    # Then make the new boundary edge to be active in tension (was dummy before T1).
+    sheet.edge_df.loc[New_boundary_edge,'is_active'] = 1
+    geom.update_all(sheet)
+
+
+"""cell extrusion behaviour function
+
+A cell extrusion behaviour function has two components.
+
+The first component is to let the STB to detach from the CT layer. This is done by a series of T1 transition on shared
+edges between the selected STB and CTs: Detach STB.
+
+The second component is to let the STB to shed from the layer. The shedding is modelled by cell removal but keep the
+vertices shared between STB units that are still in the system: STB removal.
+"""
+def extrude(sheet, manager, geom, unique_id):
+    idx = sheet.idx_lookup('face', unique_id)
+    remove.face(sheet,idx)
+    geom.update_all(sheet)
+
+def detach(sheet, manager, geom, unique_id):
+    sheet.get_extra_indices()
+    idx = sheet.idx_lookup('face', unique_id)
+    # Identify CT neighbours (non-STB)
+    CT_neighbours = sheet.get_neighbors(idx)
+    CT_neighbours = list(CT_neighbours.intersection(
+        set(sheet.face_df.loc[sheet.face_df['cell_class'] != 'STB'].index)
+    ))
+    # Find internal edges of the detaching face
+    internal_edges = sheet.edge_df[(sheet.edge_df['face'] == idx) & (sheet.edge_df['opposite'] != -1)]
+    # Filter internal edges that are shared with CT neighbours
+    shared_edges = internal_edges[
+        sheet.edge_df.loc[internal_edges['opposite'], 'face'].isin(CT_neighbours).values
+    ]
+    if shared_edges.empty:
+        raise ValueError("No mutual edge found between detaching face and CT neighbours.")
+    # Perform T1 transitions on each shared edge and update geometry
+    for edge_idx in shared_edges.index:
+        new_edge_idx = type1_transition(sheet, edge_idx, do_reindex=False, remove_tri_faces=False, multiplier=1.5)
+        sheet.edge_df.loc[new_edge_idx, 'is_active'] = 1
+        geom.update_all(sheet)
+    # Optionally extrude the detached face
+    manager.append(extrude, geom = geom, unique_id = unique_id)
diff --git a/src/tyssue/behaviors/sheet/cell_class_events.py b/src/tyssue/behaviors/sheet/cell_class_events.py
new file mode 100644
index 00000000..fb2bbe7a
--- /dev/null
+++ b/src/tyssue/behaviors/sheet/cell_class_events.py
@@ -0,0 +1,67 @@
+"""
+Event module for cell class transition rules, modify the details accordingly to your model.
+=======================
+
+"""
+
+import numpy as np
+from ...geometry.planar_geometry import PlanarGeometry
+from ...topology.sheet_topology import cell_division
+
+def cell_cycle_transition(sheet, manager, dt, p_recruit=0.1, G2_duration=0.4, G1_duration=0.11):
+    """
+    Controls cell class state transitions for cell cycle based on timers and probabilities.
+
+    Parameters
+    ----------
+    sheet: tyssue.Sheet
+        The tissue sheet.
+    manager: EventManager
+        The event manager scheduling the behaviour.
+    face_id: Integer
+        ID of the cell being controlled.
+    p_recruit: float
+        Probability for an 'S' cell to be recruited to 'G2'.
+    dt: float
+        Time step increment.
+    G2_duration: float
+        Fixed duration cells stay in G2 phase.
+    G1_duration: float
+        Fixed duration cells stay in G1 phase.
+    """
+    # Generate two df that are cells that need to change its cell class or stay in its current class.
+    cells_stay_in_class = sheet.face_df.loc[sheet.face_df['timer'] > 0].index.tolist()
+    cells_to_change = sheet.face_df.loc[sheet.face_df['timer'] <= 0]
+
+    # For cells that still needs to elapse the timer, just reduce the timer by dt
+    for cell in cells_stay_in_class:
+        sheet.face_df.loc[cell,'timer'] -= dt
+
+    # Then we change the cell type based on the cell cycle diagram.
+    G1_cells = cells_to_change.loc[cells_to_change['cell_class'] == 'G1'].index.tolist()
+    S_cells = cells_to_change.loc[cells_to_change['cell_class'] == 'S'].index.tolist()
+    G2_cells = cells_to_change.loc[cells_to_change['cell_class'] == 'G2'].index.tolist()
+    M_cells = cells_to_change.loc[cells_to_change['cell_class'] == 'M'].index.tolist()
+
+    # For cells in G1_cells indices, we simply change the cell class to S
+    sheet.face_df.loc[G1_cells, 'cell_class'] = 'S'
+
+    # For cells in S_cells , we need to loop and decide based on the random number generated for if it moves into G2.
+    for cell in S_cells:
+        if np.random.rand() < p_recruit:
+            sheet.face_df.loc[cell, 'cell_class'] = 'G2'
+            sheet.face_df.loc[cell, 'timer'] = G2_duration
+
+    # For cells in G2_cells, we move them into M class
+    sheet.face_df.loc[G2_cells, 'cell_class'] = 'M'
+
+    # For cells in M_cells, they need to be undergone cell division
+    for cell in M_cells:
+        daughter = cell_division(sheet, mother=cell, geom = PlanarGeometry )
+        # Set parent and daughter to G1 with G1 timer, note that variable daughter is the index of the new row already.
+        sheet.face_df.loc[cell, 'cell_class'] = 'G1'
+        sheet.face_df.loc[daughter, 'cell_class'] = 'G1'
+        sheet.face_df.loc[cell, 'timer'] = G1_duration
+        sheet.face_df.loc[daughter, 'timer'] = G1_duration
+
+    manager.append(cell_cycle_transition, dt = dt, p_recruit = p_recruit,G2_duration = G2_duration, G1_duration = G1_duration)
diff --git a/src/tyssue/core/objects.py b/src/tyssue/core/objects.py
index ac07db0a..625dcf2d 100644
--- a/src/tyssue/core/objects.py
+++ b/src/tyssue/core/objects.py
@@ -570,7 +570,8 @@ def get_orbits(self, center, periph):
         return orbits
 
     def idx_lookup(self, elem_id, element):
-        """returns the current index of the element with the `"id"` column equal to `elem_id`
+        """returns the current index of the element with the `"id"` column equal to `elem_id`;
+        thse stable IDs are called "ID" in face data frame, but called "unique_id" in vertex and edge data frames
 
         Parameters
         ----------
@@ -579,7 +580,10 @@ def idx_lookup(self, elem_id, element):
         element : {"vert"|"edge"|"face"|"cell"}
           the corresponding dataset.
         """
-        df = self.datasets[element]["id"]
+        if element in ['vert', 'edge','face']:
+            df = self.datasets[element]['unique_id']
+        else:
+            df = self.datasets[element]["id"]
         idx = df[df == elem_id].index
         if len(idx):
             return idx[0]
diff --git a/src/tyssue/draw/bilayer_drawing_tool.py b/src/tyssue/draw/bilayer_drawing_tool.py
new file mode 100644
index 00000000..c65589c3
--- /dev/null
+++ b/src/tyssue/draw/bilayer_drawing_tool.py
@@ -0,0 +1,40 @@
+# This script contains a function that auto updates the drawing specifications for a bilayer tissue sheet.
+
+import numpy as np
+
+def bilayer_draw_spec_update(sheet, specs):
+    """
+    Update drawing specifications for bilayer tissue sheet.
+
+    Parameters
+    ----------
+    sheet : tyssue.Sheet
+        The tissue sheet instance containing face_df and edge_df DataFrames.
+    specs : dict
+        The current drawing specifications dictionary to be updated.
+    """
+
+    # --- FACE COLOR UPDATE ---
+    # Use NumPy vectorization to assign colors:
+    # If 'cell_class' == 'STB', then color = 0.7
+    # Else, then color = 0.1
+    sheet.face_df['color'] = np.where(
+        sheet.face_df['cell_class'] == 'STB', 0.7, 0.1
+    )
+
+    # Update the specs dictionary with the new face colors
+    specs['face']['color'] = sheet.face_df['color'].to_numpy()
+
+    # Set transparency (alpha) for faces
+    specs['face']['alpha'] = 0.2
+
+    # --- EDGE WIDTH UPDATE ---
+    # Use NumPy vectorization to assign edge widths:
+    # If 'is_active' == 0, then width = 2
+    # Else, then width = 0.5
+    sheet.edge_df['width'] = np.where(
+        sheet.edge_df['is_active'] == 0, 2, 0.5
+    )
+
+    # Update the specs dictionary with the new edge widths
+    specs['edge']['width'] = sheet.edge_df['width'].to_numpy()
diff --git a/src/tyssue/draw/plt_draw.py b/src/tyssue/draw/plt_draw.py
index d48c3f0f..01c384ef 100644
--- a/src/tyssue/draw/plt_draw.py
+++ b/src/tyssue/draw/plt_draw.py
@@ -129,7 +129,39 @@ def create_gif(
         plt.close(fig)
 
     try:
-        subprocess.run(["convert", (graph_dir / "movie_*.png").as_posix(), output])
+        # Expand pixels manually (cross-platform safe)
+        pngs = sorted(graph_dir.glob("movie_*.png"))
+        png_paths = [str(p) for p in pngs]
+
+        # Detect the correct ImageMagick executable:
+        # IM6: convert
+        # IM7: magick convert
+        def find_imagemagick():
+            # Try IM7 first
+            try:
+                subprocess.run(["magick", "-version"], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
+                return ["magick", "convert"]
+            except Exception:
+                pass
+
+            # Try IM6
+            try:
+                subprocess.run(["convert", "-version"], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
+                return ["convert"]
+            except Exception:
+                pass
+
+            raise RuntimeError("ImageMagick executable not found on this system.")
+
+        im_cmd = find_imagemagick()
+
+        # Run ImageMagick safely
+        try:
+            subprocess.run(im_cmd + png_paths + [output], check=True)
+        except Exception as e:
+            print("Converting didn't work. Make sure ImageMagick is correctly installed.")
+            raise e
+
     except Exception as e:
         print(
             "Converting didn't work, make sure imagemagick is available on your system"
diff --git a/src/tyssue/generation/hexagonal_grids.py b/src/tyssue/generation/hexagonal_grids.py
index 05e82ea7..d226a3a1 100644
--- a/src/tyssue/generation/hexagonal_grids.py
+++ b/src/tyssue/generation/hexagonal_grids.py
@@ -13,7 +13,7 @@
 
 def hexa_grid2d(nx, ny, distx, disty, noise=None):
     """Creates an hexagonal grid of points"""
-    cy, cx = np.mgrid[0:ny+2, 0:nx+2]
+    cy, cx = np.mgrid[0:ny, 0:nx]
     cx = cx.astype(float)
     cy = cy.astype(float)
     cx[::2, :] += 0.5
diff --git a/src/tyssue/solvers/viscous.py b/src/tyssue/solvers/viscous.py
index fcb4d351..86513001 100644
--- a/src/tyssue/solvers/viscous.py
+++ b/src/tyssue/solvers/viscous.py
@@ -122,6 +122,8 @@ def solve(self, tf, dt, on_topo_change=None, topo_change_args=()):
             self.prev_t = t
             if self.manager is not None:
                 self.manager.execute(self.eptm)
+                self.eptm.reset_index()
+                self.eptm.reset_topo()
                 self.geom.update_all(self.eptm)
                 self.manager.update()
 
@@ -152,6 +154,26 @@ def ode_func(self, t, pos):
             / self.eptm.vert_df.loc[self.eptm.active_verts, "viscosity"].values[:, None]
         ).ravel()
 
+    def single_step_movement(self, dt):
+        """Solves and returns the final movement vector for a single step of the Euler solver defined above.
+
+        Parameters
+        ----------
+        dt : float, time step
+        on_topo_change : function, optional, default None
+             function of `self.eptm`
+        topo_change_args : tuple, arguments passed to `on_topo_change`
+
+        """
+        self.eptm.settings["dt"] = dt
+        for t in np.arange(self.prev_t, dt, dt):
+            pos = self.current_pos
+            dot_r = self.ode_func(t, pos)
+            if self.bounds is not None:
+                dot_r = np.clip(dot_r, *self.bounds)
+        movement = dot_r * dt
+        return movement, dot_r
+
 
 class IVPSolver:
     def __init__(self, *args, **kwargs):
diff --git a/src/tyssue/topology/base_topology.py b/src/tyssue/topology/base_topology.py
index 920f527d..efd023f0 100644
--- a/src/tyssue/topology/base_topology.py
+++ b/src/tyssue/topology/base_topology.py
@@ -14,6 +14,7 @@ def split_vert(sheet, vert, face, to_rewire, epsilon, recenter=False):
     """Creates a new vertex and moves it towards the center of face.
 
     The edges in to_rewire will be connected to the new vertex.
+    This is the low-level function called by topologies when splitting vertices and rewiring edges.
 
     Parameters
     ----------
@@ -48,7 +49,10 @@ def split_vert(sheet, vert, face, to_rewire, epsilon, recenter=False):
     else:
         sheet.vert_df.loc[new_vert, sheet.coords] += shift
 
-    # rewire
+    # rewire the edges.
+    # Updates the rows in the original edge_df at the same indices as to_rewire.
+    # The replacement is done such that any occurrence of 'vert' in the 'srce' column or 'trgt' column of to_rewire
+    # is replaced with new_vert.
     sheet.edge_df.loc[to_rewire.index] = to_rewire.replace(
         {"srce": vert, "trgt": vert}, new_vert
     )
@@ -144,23 +148,26 @@ def close_face(eptm, face):
     faces. Returns the index of the new edge if created, otherwise None
     """
     logger.debug(f"closing face {face}")
+    # Collects all edges of the face.
     face_edges = eptm.edge_df[eptm.edge_df["face"] == face]
     srces = set(face_edges["srce"])
     trgts = set(face_edges["trgt"])
-
+    # If the set of sources equals the set of targets, every vertex has both incoming and outgoing edges, then the loop is complete.
     if srces == trgts:
         logger.debug("Face %d already closed", face)
         return None
     try:
-        (single_srce,) = srces.difference(trgts)
-        (single_trgt,) = trgts.difference(srces)
+        (single_srce,) = srces.difference(trgts)    # vertices that only appear as sources (no incoming edge).
+        (single_trgt,) = trgts.difference(srces)    # vertices that only appear as targets (no outgoing edge).
     except ValueError as err:
         print("Closing only possible with exactly two dangling vertices")
         raise err
-
+    # If there’s exactly one of each, the face is missing a single edge.
+    # Duplicates one existing edge row
     eptm.edge_df = pd.concat([eptm.edge_df, face_edges.iloc[0:1]], ignore_index=True)
     eptm.edge_df.index.name = "edge"
     new_edge = eptm.edge_df.index[-1]
+    # Reassigns its srce and trgt to connect the dangling vertices
     eptm.edge_df.loc[new_edge, ["srce", "trgt"]] = single_trgt, single_srce
     return new_edge
 
@@ -183,6 +190,14 @@ def drop_two_sided_faces(eptm):
     eptm.face_df.drop(two_sided, axis=0, inplace=True)
 
 
+def drop_face(sheet, face, **kwargs):
+    """
+    Removes the face indexed by "face" and all associated edges
+    """
+    edge = sheet.edge_df.loc[(sheet.edge_df['face'] == face)].index
+    print(f"Dropping face '{face}'")
+    sheet.remove(edge, **kwargs)
+
 def remove_face(sheet, face):
     """Removes a face from the mesh.
 
diff --git a/src/tyssue/topology/sheet_topology.py b/src/tyssue/topology/sheet_topology.py
index 372533bc..2845f913 100644
--- a/src/tyssue/topology/sheet_topology.py
+++ b/src/tyssue/topology/sheet_topology.py
@@ -35,14 +35,15 @@ def split_vert(
     if face is None:
         face = np.random.choice(sheet.edge_df[sheet.edge_df["srce"] == vert]["face"])
 
-    face_edges = sheet.edge_df.query(f"face == {face}")
-    (prev_v,) = face_edges[face_edges["trgt"] == vert]["srce"]
-    (next_v,) = face_edges[face_edges["srce"] == vert]["trgt"]
+    face_edges = sheet.edge_df.query(f"face == {face}")     # A filerted view of the edge_df, contains only the edges belonging to the given face.
+    (prev_v,) = face_edges[face_edges["trgt"] == vert]["srce"]  # The vertex connects into vert (the sources of the edges that whose target is vert).
+    (next_v,) = face_edges[face_edges["srce"] == vert]["trgt"]  # The vertex connects out of vert (the targets of the edges whose source is vert).
+    # A filtered view of the edge_df, contains all edges that touch either prev_v or next_v, regardless of face.
     connected = sheet.edge_df[
         sheet.edge_df["trgt"].isin((next_v, prev_v))
         | sheet.edge_df["srce"].isin((next_v, prev_v))
     ]
-
+    # pass the filtered subset of edge_df as connected to rewire.
     base_split_vert(sheet, vert, face, connected, epsilon, recenter)
     new_edges = []
     for face_ in connected["face"]:
@@ -57,7 +58,7 @@ def split_vert(
     return new_edges
 
 
-def type1_transition(sheet, edge01, *, remove_tri_faces=True, multiplier=1.5):
+def type1_transition(sheet, edge01, *, do_reindex =True, remove_tri_faces=True, multiplier=1.5):
     """Performs a type 1 transition around the edge edge01
 
     See ../../doc/illus/t1_transition.png for a sketch of the definition
@@ -70,9 +71,8 @@ def type1_transition(sheet, edge01, *, remove_tri_faces=True, multiplier=1.5):
     sheet : a `Sheet` instance
     edge_01 : int
        index of the edge around which the transition takes place
-    epsilon : float, optional, deprecated
-       default 0.1, the initial length of the new edge, in case "threshold_length"
-       is not in the sheet.settings
+    do_reindex : bool, optional
+        whether or not to reindex  the sheet.
     remove_tri_faces : bool, optional
        if True (the default), will remove triangular cells
        after the T1 transition is performed
@@ -80,13 +80,12 @@ def type1_transition(sheet, edge01, *, remove_tri_faces=True, multiplier=1.5):
        default 1.5, the multiplier to the threshold length, so that the
        length of the new edge is set to multiplier * threshold_length
 
-
     """
 
     srce, trgt, face = sheet.edge_df.loc[edge01, ["srce", "trgt", "face"]].astype(int)
 
     vert = min(srce, trgt)  # find the vertex that won't be reindexed
-    ret_code = collapse_edge(sheet, edge01, reindex=True, allow_two_sided=True)
+    ret_code = collapse_edge(sheet, edge01, reindex=do_reindex, allow_two_sided=True)
     if ret_code < 0:
         warnings.warn(f"Collapse of edge {edge01} failed")
         return ret_code
@@ -96,7 +95,7 @@ def type1_transition(sheet, edge01, *, remove_tri_faces=True, multiplier=1.5):
         vert,
         face,
         multiplier=multiplier,
-        reindex=True,
+        reindex=do_reindex,
         recenter=True,
     )
 
@@ -208,7 +207,7 @@ def face_division(sheet, mother, vert_a, vert_b):
     sheet.face_df = pd.concat([sheet.face_df, face_cols], ignore_index=True)
     sheet.face_df.index.name = "face"
     daughter = int(sheet.face_df.index[-1])
-
+    sheet.face_df.loc[daughter, 'unique_id'] = sheet.face_df['unique_id'].max() + 1
     edge_cols = sheet.edge_df[sheet.edge_df["face"] == mother].iloc[0:1]
 
     sheet.edge_df = pd.concat([sheet.edge_df, edge_cols, edge_cols], ignore_index=True)
diff --git a/tests/topology/test_bulk_topology.py b/tests/topology/test_bulk_topology.py
index a598cfa9..f6b4c06d 100644
--- a/tests/topology/test_bulk_topology.py
+++ b/tests/topology/test_bulk_topology.py
@@ -136,8 +136,8 @@ def test_IH_transition():
     assert eptm.Nv == Nv + 1
 
     invalid = eptm.get_invalid()
-    assert np.alltrue(1 - invalid)
-    assert np.alltrue(eptm.edge_df["sub_vol"] > 0)
+    assert np.all(1 - invalid)
+    assert np.all(eptm.edge_df["sub_vol"] > 0)
     assert (
         eptm.face_df[eptm.face_df.segment == "apical"].shape[0] == eptm.cell_df.shape[0]
     )
@@ -156,7 +156,7 @@ def test_split_vert():
     BulkGeometry.update_all(eptm)
 
     invalid = eptm.get_invalid()
-    assert np.alltrue(1 - invalid)
+    assert np.all(1 - invalid)
 
 
 def test_HI_transition():
@@ -179,8 +179,8 @@ def test_HI_transition():
     assert eptm.Nv == Nv
 
     invalid = eptm.get_invalid()
-    assert np.alltrue(1 - invalid)
-    assert np.alltrue(eptm.edge_df["sub_vol"] > 0)
+    assert np.all(1 - invalid)
+    assert np.all(eptm.edge_df["sub_vol"] > 0)
 
 
 def test_find_transitions():