diff --git a/src/compas/datastructures/volmesh/volmesh.py b/src/compas/datastructures/volmesh/volmesh.py
index 6c23c547d1a..afa03aa8c77 100644
--- a/src/compas/datastructures/volmesh/volmesh.py
+++ b/src/compas/datastructures/volmesh/volmesh.py
@@ -383,7 +383,7 @@ def from_obj(cls, filepath, precision=None):
 
     @classmethod
     def from_vertices_and_cells(cls, vertices, cells):
-        # type: (list[list[float]], list[list[list[int]]]) -> VolMesh
+        # type: (list[list[float]] | dict[int, list[float]], list[list[list[int]]]) -> VolMesh
         """Construct a volmesh object from vertices and cells.
 
         Parameters
@@ -407,11 +407,11 @@ def from_vertices_and_cells(cls, vertices, cells):
         volmesh = cls()
 
         if isinstance(vertices, Mapping):
-            for key, xyz in vertices.items():
+            for key, xyz in vertices.items():  # type: ignore
                 volmesh.add_vertex(key=key, attr_dict=dict(zip(("x", "y", "z"), xyz)))
         else:
-            for x, y, z in iter(vertices):
-                volmesh.add_vertex(x=x, y=y, z=z)
+            for x, y, z in iter(vertices):  # type: ignore
+                volmesh.add_vertex(x=x, y=y, z=z)  # type: ignore
 
         for cell in cells:
             volmesh.add_cell(cell)
@@ -531,9 +531,9 @@ def to_obj(self, filepath, precision=None, **kwargs):
         the faces to the file.
 
         """
-        meshes = [self.cell_to_mesh(cell) for cell in self.cells()]
+        meshes = [self.cell_to_mesh(cell) for cell in self.cells()]  # type: ignore
         obj = OBJ(filepath, precision=precision)
-        obj.write(meshes, **kwargs)
+        obj.write(meshes, **kwargs)  # type: ignore
 
     def to_vertices_and_cells(self):
         # type: () -> tuple[list[list[float]], list[list[list[int]]]]
@@ -659,7 +659,7 @@ def vertex_sample(self, size=1):
         :meth:`edge_sample`, :meth:`face_sample`, :meth:`cell_sample`
 
         """
-        return sample(list(self.vertices()), size)
+        return sample(list(self.vertices()), size)  # type: ignore
 
     def edge_sample(self, size=1):
         # type: (int) -> list[tuple[int, int]]
@@ -680,7 +680,7 @@ def edge_sample(self, size=1):
         :meth:`vertex_sample`, :meth:`face_sample`, :meth:`cell_sample`
 
         """
-        return sample(list(self.edges()), size)
+        return sample(list(self.edges()), size)  # type: ignore
 
     def face_sample(self, size=1):
         # type: (int) -> list[int]
@@ -701,7 +701,7 @@ def face_sample(self, size=1):
         :meth:`vertex_sample`, :meth:`edge_sample`, :meth:`cell_sample`
 
         """
-        return sample(list(self.faces()), size)
+        return sample(list(self.faces()), size)  # type: ignore
 
     def cell_sample(self, size=1):
         # type: (int) -> list[int]
@@ -722,7 +722,7 @@ def cell_sample(self, size=1):
         :meth:`vertex_sample`, :meth:`edge_sample`, :meth:`face_sample`
 
         """
-        return sample(list(self.cells()), size)
+        return sample(list(self.cells()), size)  # type: ignore
 
     def vertex_index(self):
         # type: () -> dict[int, int]
@@ -739,7 +739,7 @@ def vertex_index(self):
         :meth:`index_vertex`
 
         """
-        return {key: index for index, key in enumerate(self.vertices())}
+        return {key: index for index, key in enumerate(self.vertices())}  # type: ignore
 
     def index_vertex(self):
         # type: () -> dict[int, int]
@@ -756,7 +756,7 @@ def index_vertex(self):
         :meth:`vertex_index`
 
         """
-        return dict(enumerate(self.vertices()))
+        return dict(enumerate(self.vertices()))  # type: ignore
 
     def vertex_gkey(self, precision=None):
         # type: (int | None) -> dict[int, str]
@@ -781,7 +781,7 @@ def vertex_gkey(self, precision=None):
         """
         gkey = TOL.geometric_key
         xyz = self.vertex_coordinates
-        return {vertex: gkey(xyz(vertex), precision) for vertex in self.vertices()}
+        return {vertex: gkey(xyz(vertex), precision) for vertex in self.vertices()}  # type: ignore
 
     def gkey_vertex(self, precision=None):
         # type: (int | None) -> dict[str, int]
@@ -806,7 +806,7 @@ def gkey_vertex(self, precision=None):
         """
         gkey = TOL.geometric_key
         xyz = self.vertex_coordinates
-        return {gkey(xyz(vertex), precision): vertex for vertex in self.vertices()}
+        return {gkey(xyz(vertex), precision): vertex for vertex in self.vertices()}  # type: ignore
 
     # --------------------------------------------------------------------------
     # Builders & Modifiers
@@ -1760,8 +1760,10 @@ def vertex_halffaces(self, vertex):
         u = vertex
         faces = []
         for v in self._plane[u]:
-            for face in self._plane[u][v]:
-                if face is not None:
+            for w in self._plane[u][v]:
+                cell = self._plane[u][v][w]
+                if cell is not None:
+                    face = self.cell_halfedge_face(cell, (u, v))
                     faces.append(face)
         return faces
 
@@ -1789,7 +1791,8 @@ def vertex_cells(self, vertex):
             for w in self._plane[u][v]:
                 cell = self._plane[u][v][w]
                 if cell is not None:
-                    cells.append(cell)
+                    if cell not in cells:
+                        cells.append(cell)
         return cells
 
     def is_vertex_on_boundary(self, vertex):
@@ -3108,36 +3111,36 @@ def halfface_opposite_halfface(self, halfface):
         nbr = self._plane[w][v][u]
         return None if nbr is None else self._cell[nbr][w][v]
 
-    def halfface_adjacent_halfface(self, halfface, halfedge):
-        """Return the halfface adjacent to the halfface across the halfedge.
-
-        Parameters
-        ----------
-        halfface : int
-            The identifier of the halfface.
-        halfedge : tuple[int, int]
-            The identifier of the halfedge.
-
-        Returns
-        -------
-        int | None
-            The identifier of the adjacent half-face, or None if `halfedge` is on the boundary.
-
-        See Also
-        --------
-        :meth:`halfface_opposite_halfface`
-
-        Notes
-        -----
-        The adjacent face belongs a to one of the cell neighbors over faces of the initial cell.
-        A face and its adjacent face share two common vertices.
-
-        """
-        u, v = halfedge
-        cell = self.halfface_cell(halfface)
-        nbr_halfface = self._cell[cell][v][u]
-        nbr_cell = self._plane[u][v][nbr_halfface]
-        return None if nbr_cell is None else self._cell[nbr_cell][v][u]
+    # def halfface_adjacent_halfface(self, halfface, halfedge):
+    #     """Return the halfface adjacent to the halfface across the halfedge.
+
+    #     Parameters
+    #     ----------
+    #     halfface : int
+    #         The identifier of the halfface.
+    #     halfedge : tuple[int, int]
+    #         The identifier of the halfedge.
+
+    #     Returns
+    #     -------
+    #     int | None
+    #         The identifier of the adjacent half-face, or None if `halfedge` is on the boundary.
+
+    #     See Also
+    #     --------
+    #     :meth:`halfface_opposite_halfface`
+
+    #     Notes
+    #     -----
+    #     The adjacent face belongs to one of the cell neighbors over faces of the initial cell.
+    #     A face and its adjacent face share two common vertices.
+
+    #     """
+    #     u, v = halfedge
+    #     cell = self.halfface_cell(halfface)
+    #     nbr_halfface = self._cell[cell][v][u]
+    #     nbr_cell = self._plane[u][v][nbr_halfface]
+    #     return None if nbr_cell is None else self._cell[nbr_cell][v][u]
 
     def halfface_vertex_ancestor(self, halfface, vertex):
         """Return the vertex before the specified vertex in a specific face.
@@ -3282,8 +3285,8 @@ def is_halfface_on_boundary(self, halfface):
         :meth:`is_vertex_on_boundary`, :meth:`is_edge_on_boundary`, :meth:`is_cell_on_boundary`
 
         """
-        u, v = self._halfface[halfface][:2]
-        return self._plane[v][u][halfface] is None
+        u, v, w = self._halfface[halfface][:3]
+        return self._plane[w][v][u] is None
 
     # --------------------------------------------------------------------------
     # Face Geometry
diff --git a/tests/compas/datastructures/test_volmesh.py b/tests/compas/datastructures/test_volmesh.py
index 60795d929f2..c7f3e92d2f5 100644
--- a/tests/compas/datastructures/test_volmesh.py
+++ b/tests/compas/datastructures/test_volmesh.py
@@ -72,6 +72,111 @@ def test_volmesh_data():
 # Samples
 # ==============================================================================
 
+# ==============================================================================
+# Topology
+# ==============================================================================
+
+
+@pytest.mark.parametrize(
+    "nx,ny,nz",
+    [
+        [1, 1, 1],
+        [2, 2, 2],
+        [3, 3, 3],
+    ],
+)
+def test_vertex_neighbours(nx, ny, nz):
+    volmesh = VolMesh.from_meshgrid(1, 1, 1, nx, ny, nz)
+
+    for vertex in volmesh.vertices():
+        count = len(volmesh.vertex_neighbors(vertex))
+
+        if volmesh.is_vertex_on_boundary(vertex):
+            assert 2 < count < 6
+        else:
+            assert count == 6
+
+
+@pytest.mark.parametrize(
+    "nx,ny,nz",
+    [
+        [1, 1, 1],
+        [2, 2, 2],
+        [3, 3, 3],
+    ],
+)
+def test_vertex_cells(nx, ny, nz):
+    volmesh = VolMesh.from_meshgrid(1, 1, 1, nx, ny, nz)
+
+    for vertex in volmesh.vertices():
+        nbrs = len(volmesh.vertex_neighbors(vertex))
+        cells = len(volmesh.vertex_cells(vertex))
+
+        if nbrs == 6:
+            assert cells == 8
+        elif nbrs == 5:
+            assert cells == 4
+        elif nbrs == 4:
+            assert cells == 2
+        elif nbrs == 3:
+            assert cells == 1
+
+
+@pytest.mark.parametrize(
+    "nx,ny,nz",
+    [
+        [1, 1, 1],
+        [2, 2, 2],
+        [3, 3, 3],
+    ],
+)
+def test_edge_cells(nx, ny, nz):
+    volmesh = VolMesh.from_meshgrid(1, 1, 1, nx, ny, nz)
+
+    for edge in volmesh.edges():
+        cells = len(volmesh.edge_cells(edge))
+
+        if volmesh.is_edge_on_boundary(edge):
+            assert 0 < cells < 3
+        else:
+            assert cells == 4
+
+
+@pytest.mark.parametrize(
+    "nx,ny,nz",
+    [
+        [1, 1, 1],
+        [2, 2, 2],
+        [3, 3, 3],
+    ],
+)
+def test_edge_halffaces(nx, ny, nz):
+    volmesh = VolMesh.from_meshgrid(1, 1, 1, nx, ny, nz)
+
+    for edge in volmesh.edges():
+        faces = len(volmesh.edge_halffaces(edge))
+
+        if volmesh.is_edge_on_boundary(edge):
+            assert 0 < faces < 3
+        else:
+            assert faces == 4
+
+
+@pytest.mark.parametrize(
+    "nx,ny,nz",
+    [
+        [1, 1, 1],
+        [2, 2, 2],
+        [3, 3, 3],
+    ],
+)
+def test_halffaces_on_boundary(nx, ny, nz):
+    volmesh = VolMesh.from_meshgrid(1, 1, 1, nx, ny, nz)
+
+    count = sum(volmesh.is_halfface_on_boundary(face) for face in volmesh.halffaces())
+    assert count == 2 * nx * ny + 2 * ny * nz + 2 * nx * nz
+
+
 # ==============================================================================
 # Vertex Attributes
 # ==============================================================================