diff --git a/src/pyransame/util.py b/src/pyransame/util.py
index 13d074f..3529f54 100644
--- a/src/pyransame/util.py
+++ b/src/pyransame/util.py
@@ -3,6 +3,8 @@
 from typing import Tuple
 
 import numpy as np
+import pyvista as pv
+from vtk import mutable
 
 import pyransame
 
@@ -351,3 +353,41 @@ def _generate_points_in_hexagonal_prism(points: np.ndarray, n: int = 1) -> np.nd
             ] = _generate_points_in_hexahedron(points[hexahedron1, :], n=count)
 
     return out
+
+
+def _generate_points_in_polyhedron(cell: pv.Cell, n: int = 1) -> np.ndarray:
+    faces = cell.faces
+
+    # version >0.39 pyvista could be used in the future
+    para_center = [0.0, 0.0, 0.0]
+    sub_id = cell.GetParametricCenter(para_center)
+    # EvaluateLocation requires mutable sub_id
+    sub_id = mutable(sub_id)
+    # center and weights are returned from EvaluateLocation
+    cell_center = [0.0, 0.0, 0.0]
+    weights = [0.0] * cell.n_points
+    cell.EvaluateLocation(sub_id, para_center, cell_center, weights)
+
+    tetras = []
+
+    for face in faces:
+        face_points = face.points
+        ntri = face_points.shape[0] - 2
+        for i in range(ntri):
+            tetra = face_points[[0, i + 1, i + 2], :]
+            tetra = np.append(tetra, [cell_center], axis=0)
+            tetras.append(tetra)
+
+    areas = np.array([_area_tetra(tetra) for tetra in tetras])
+
+    p = areas / areas.sum()
+    out = np.empty((n, 3))
+
+    chosen_cells, unique_counts, point_indices = _random_cells(len(tetras), n, p)
+
+    for i, (chosen_cell, count) in enumerate(zip(chosen_cells, unique_counts)):
+        out[point_indices[i] : point_indices[i + 1], :] = _generate_points_in_tetra(
+            tetras[chosen_cell], n=count
+        )
+
+    return out
diff --git a/src/pyransame/volume.py b/src/pyransame/volume.py
index 27684ae..07b5da8 100644
--- a/src/pyransame/volume.py
+++ b/src/pyransame/volume.py
@@ -120,6 +120,10 @@ def random_volume_points(
             points[
                 point_indices[i] : point_indices[i + 1], :
             ] = util._generate_points_in_hexagonal_prism(c.points, count)
+        elif c.type == CellType.POLYHEDRON:
+            points[
+                point_indices[i] : point_indices[i + 1], :
+            ] = util._generate_points_in_polyhedron(c, count)
         else:
             raise NotImplementedError(
                 f"Random generation for {c.type.name} not yet supported"
diff --git a/tests/test_random_volume_points.py b/tests/test_random_volume_points.py
index c89c889..fe40597 100644
--- a/tests/test_random_volume_points.py
+++ b/tests/test_random_volume_points.py
@@ -68,6 +68,27 @@ def make_hexagonal_prism():
     return pv.UnstructuredGrid(cells, celltypes, p)
 
 
+def make_polyhedron():
+    dodecahedron_poly = pv.Dodecahedron()
+
+    faces = []
+    for cell in dodecahedron_poly.cell:
+        faces.append(cell.point_ids)
+
+    faces = []
+    faces.append(dodecahedron_poly.n_cells)
+    for cell in dodecahedron_poly.cell:
+        point_ids = cell.point_ids
+        faces.append(len(point_ids))
+        [faces.append(id) for id in point_ids]
+
+    faces.insert(0, len(faces))
+
+    return pv.UnstructuredGrid(
+        faces, [pv.CellType.POLYHEDRON], dodecahedron_poly.points
+    )
+
+
 def test_cell_types():
     mesh = pv.UniformGrid(dimensions=(4, 4, 4))
     assert mesh.get_cell(0).type == pv.CellType.VOXEL
@@ -97,6 +118,10 @@ def test_cell_types():
     assert mesh.get_cell(0).type == pv.CellType.HEXAGONAL_PRISM
     pyransame.random_volume_points(mesh, 20)
 
+    mesh = make_polyhedron()
+    assert mesh.get_cell(0).type == pv.CellType.POLYHEDRON
+    pyransame.random_volume_points(mesh, 20)
+
 
 def test_mixed_types():
     # as long as there are 3D cells, we should be able to sample even if there are other cell types
diff --git a/tests/test_util.py b/tests/test_util.py
index 0128797..fd02697 100644
--- a/tests/test_util.py
+++ b/tests/test_util.py
@@ -7,6 +7,7 @@
 from hypothesis import assume, given, settings
 from hypothesis import strategies as st
 from hypothesis.extra.numpy import arrays
+from vtk import VTK_POLYHEDRON, vtkIdList, vtkPoints, vtkUnstructuredGrid
 
 from pyransame.util import (
     _generate_points_in_hexagonal_prism,
@@ -14,6 +15,7 @@
     _generate_points_in_pentagonal_prism,
     _generate_points_in_pixel,
     _generate_points_in_polygon,
+    _generate_points_in_polyhedron,
     _generate_points_in_pyramid,
     _generate_points_in_quad,
     _generate_points_in_tetra,
@@ -301,3 +303,28 @@ def test_uniformity_hexagonal_prism():
     center = np.array([0.0, 0.0, 0.5])
     points = _generate_points_in_hexagonal_prism(p, 2000000)
     assert np.allclose(points.mean(axis=0), center, rtol=1e-3, atol=1e-3)
+
+
+def test_uniformity_polyhedron():
+    dodecahedron_poly = pv.Dodecahedron()
+
+    faces = []
+    for cell in dodecahedron_poly.cell:
+        faces.append(cell.point_ids)
+
+    faces = []
+    faces.append(dodecahedron_poly.n_cells)
+    for cell in dodecahedron_poly.cell:
+        point_ids = cell.point_ids
+        faces.append(len(point_ids))
+        [faces.append(id) for id in point_ids]
+
+    faces.insert(0, len(faces))
+
+    mesh = pv.UnstructuredGrid(
+        faces, [pv.CellType.POLYHEDRON], dodecahedron_poly.points
+    )
+
+    points = _generate_points_in_polyhedron(mesh.get_cell(0), 2000000)
+
+    assert np.allclose(points.mean(axis=0), np.array(mesh.center), rtol=1e-3, atol=1e-3)