diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml
index 8aacfd5..28f3ce9 100644
--- a/.github/workflows/run_tests.yml
+++ b/.github/workflows/run_tests.yml
@@ -13,7 +13,7 @@ jobs:
       - uses: actions/checkout@v3
 
       - name: Install dependencies
-        run: python3 -m pip install --break-system-packages pyvista
+        run: python3 -m pip install --break-system-packages pyvista pytest
 
       - name: Run unit tests
         run: python3 -m pytest test/
diff --git a/febug/__init__.py b/febug/__init__.py
index e433eca..6cb70e6 100644
--- a/febug/__init__.py
+++ b/febug/__init__.py
@@ -9,6 +9,7 @@
     plot_streamlines_from_source,
     plot_entity_indices, plot_entity_indices_global,
     plot_entity_indices_original,
+    plot_geometry_indices,
     plot_mesh_quality, plot_function_dofs,
     plot_point_cloud)
 
diff --git a/febug/plot.py b/febug/plot.py
index e24a826..01d358c 100644
--- a/febug/plot.py
+++ b/febug/plot.py
@@ -473,6 +473,38 @@ def get_original_indices(
         local_to_label_mapping=get_original_indices)
 
 
+def plot_geometry_indices(mesh: dolfinx.mesh.Mesh,
+                          plotter: pyvista.Plotter=None):
+    if plotter is None:
+        plotter = pyvista.Plotter()
+
+    plot_mesh(mesh, tdim=0, plotter=plotter, show_owners=False)
+
+    size_local = mesh.geometry.index_map().size_local
+    num_ghosts = mesh.geometry.index_map().num_ghosts
+    entities = np.arange(size_local, dtype=np.int32)
+    ghosts = np.arange(size_local, size_local + num_ghosts, dtype=np.int32)
+
+    if size_local > 0:
+        x = mesh.geometry.x[entities]
+        x_polydata = pyvista.PolyData(x)
+        labels = np.arange(size_local)
+        x_polydata["labels"] = [f"{i}" for i in labels]
+        plotter.add_point_labels(x_polydata, "labels", **entity_label_args,
+                                 point_color="grey")
+
+    if num_ghosts > 0:
+        x_ghost = mesh.geometry.x[ghosts]
+        x_ghost_polydata = pyvista.PolyData(x_ghost)
+        ghost_labels = np.arange(size_local, size_local + num_ghosts)
+        x_ghost_polydata["labels"] = [f"{i}" for i in ghost_labels]
+        plotter.add_point_labels(
+            x_ghost_polydata, "labels", **entity_label_args,
+            point_color="pink")
+
+    return plotter
+
+
 def plot_point_cloud(xp: np.ndarray[float],
                      plotter: pyvista.Plotter=None):
     if plotter is None:
diff --git a/setup.py b/setup.py
index ca555ab..a744909 100644
--- a/setup.py
+++ b/setup.py
@@ -13,9 +13,9 @@
     print("Python 3.10 or higher required, please upgrade.")
     sys.exit(1)
 
-VERSION = "0.8.0"
+VERSION = "0.9.0"
 
-REQUIREMENTS = ["pyvista", "fenics-dolfinx>0.8.0"]
+REQUIREMENTS = ["pyvista", "fenics-dolfinx>0.9.0"]
 
 
 class CMakeExtension(Extension):
diff --git a/test/test_plot.py b/test/test_plot.py
index 009ee92..572dbf5 100644
--- a/test/test_plot.py
+++ b/test/test_plot.py
@@ -165,6 +165,13 @@ def test_plot_entity_indices_global(mesh):
         febug.plot_entity_indices_global(mesh, d)
 
 
+@pytest.mark.parametrize("mesh", all_meshes)
+def test_plot_geometry_indices(mesh):
+    for d in range(mesh.topology.dim+1):
+        mesh.topology.create_connectivity(d, mesh.topology.dim)
+        febug.plot_geometry_indices(mesh)
+
+
 @pytest.mark.parametrize(
     "mesh,tdim", [pytest.param(msh, tdim, marks=pytest.mark.xfail(
         tdim not in (0, msh.topology.dim),