diff --git a/tests/topology/test_topology.py b/tests/topology/test_topology.py
index bb2f907d..fa1db17e 100644
--- a/tests/topology/test_topology.py
+++ b/tests/topology/test_topology.py
@@ -16,7 +16,7 @@
     condition_4ii,
     remove_face,
 )
-from tyssue.topology.sheet_topology import cell_division, split_vert, type1_transition
+from tyssue.topology.sheet_topology import cell_division, split_vert, type1_transition, drop_face
 
 
 def test_condition4i():
@@ -40,7 +40,6 @@ def test_condition4ii():
 
 
 def test_division():
-
     h5store = os.path.join(stores_dir, "small_hexagonal.hf5")
 
     datasets = load_datasets(h5store, data_names=["face", "vert", "edge"])
@@ -58,7 +57,6 @@ def test_division():
 
 
 def test_t1_transition():
-
     h5store = os.path.join(stores_dir, "small_hexagonal.hf5")
     datasets = load_datasets(h5store, data_names=["face", "vert", "edge"])
     specs = cylindrical_sheet()
@@ -91,7 +89,6 @@ def test_t1_at_border():
 
 
 def test_split_vert():
-
     datasets, specs = three_faces_sheet()
     sheet = Sheet("3cells_2D", datasets, specs)
     geom.update_all(sheet)
@@ -137,10 +134,29 @@ def test_merge_border_edges():
     sheet.get_opposite()
     sheet.sanitize(trim_borders=True)
     assert (
-        sheet.edge_df[sheet.edge_df["opposite"] < 0]
-        .groupby("face")["opposite"]
-        .sum()
-        .min()
-        == -1
+            sheet.edge_df[sheet.edge_df["opposite"] < 0]
+            .groupby("face")["opposite"]
+            .sum()
+            .min()
+            == -1
     )
     assert not set(sheet.vert_df.index).difference(sheet.edge_df.srce)
+
+
+def test_drop_face():
+    init_sheet = Sheet.planar_sheet_2d("planar", 6, 7, 1, 1)
+    init_sheet.sanitize(trim_borders=True, order_edges=True)
+    geom = PlanarGeometry
+    sheet = init_sheet.copy(deep_copy=True)
+
+    drop_face(sheet, 12, geom)
+    drop_face(sheet, 11, geom)
+
+    assert sheet.Nf == init_sheet.Nf - 2
+    assert sheet.Nv == init_sheet.Nv - 2
+    assert sheet.Ne == init_sheet.Ne - 14
+
+    sheet.get_opposite()
+    init_sheet.get_opposite()
+    assert sheet.edge_df[sheet.edge_df['opposite'] == -1].shape[0] == \
+           init_sheet.edge_df[init_sheet.edge_df['opposite'] == -1].shape[0] + 8
diff --git a/tyssue/topology/sheet_topology.py b/tyssue/topology/sheet_topology.py
index a32732dd..372533bc 100644
--- a/tyssue/topology/sheet_topology.py
+++ b/tyssue/topology/sheet_topology.py
@@ -241,6 +241,17 @@ def face_division(sheet, mother, vert_a, vert_b):
     return daughter
 
 
+def drop_face(eptm, face, geom, **kwargs):
+    """
+    Removes the face indexed by "face" and all associated edges to allow holes
+    """
+    edge = eptm.edge_df.loc[(eptm.edge_df['face'] == face)].index
+
+    eptm.remove(edge, **kwargs)
+    eptm.sanitize(trim_borders = True)
+    geom.update_all(eptm)
+
+
 def resolve_t1s(sheet, geom, model, solver, max_iter=60):
 
     l_th = sheet.settings["threshold_length"]