Add pentagonal prism

src/pyransame/util.py

@@ -160,7 +160,7 @@ def _area_tetra(points: np.ndarray) -> float:
     )
 
 
-def _area_pyramid(points: np.ndarray):
+def _area_pyramid(points: np.ndarray) -> float:
     tetra0 = [0, 1, 2, 4]
     tetra1 = [0, 2, 3, 4]
 
@@ -229,6 +229,13 @@ def _generate_points_in_pixel(
     return a + np.atleast_2d(r[:, 0]).T * v0 + np.atleast_2d(r[:, 1]).T * v1
 
 
+def _area_wedge(points: np.ndarray) -> float:
+    tetra = [0, 2, 1, 4]
+    pyramid = [0, 2, 5, 3, 4]
+
+    return _area_tetra(points[tetra, :]) + _area_pyramid(points[pyramid, :])
+
+
 def _generate_points_in_wedge(points: np.ndarray, n: int = 1) -> np.ndarray:
     tetra = [0, 2, 1, 4]
     pyramid = [0, 2, 5, 3, 4]
@@ -256,6 +263,17 @@ def _generate_points_in_wedge(points: np.ndarray, n: int = 1) -> np.ndarray:
     return out
 
 
+def _area_hexahedron(points: np.ndarray) -> float:
+    tetras = [
+        [0, 1, 4, 3],
+        [3, 7, 6, 4],
+        [1, 5, 4, 6],
+        [2, 3, 6, 1],
+        [3, 1, 6, 4],
+    ]
+    return np.array([_area_tetra(points[tetra, :]) for tetra in tetras]).sum()
+
+
 def _generate_points_in_hexahedron(points: np.ndarray, n: int = 1) -> np.ndarray:
     tetras = [
         [0, 1, 4, 3],
@@ -278,3 +296,29 @@ def _generate_points_in_hexahedron(points: np.ndarray, n: int = 1) -> np.ndarray
         )
 
     return out
+
+
+def _generate_points_in_pentagonal_prism(points: np.ndarray, n: int = 1) -> np.ndarray:
+    wedge = [2, 4, 3, 7, 9, 8]
+    hexahedron = [0, 1, 2, 4, 5, 6, 7, 9]
+
+    areas = np.array(
+        [_area_wedge(points[wedge, :]), _area_hexahedron(points[hexahedron, :])]
+    )
+
+    p = areas / areas.sum()
+    out = np.empty((n, 3))
+
+    chosen_cells, unique_counts, point_indices = _random_cells(2, n, p)
+
+    for i, (chosen_cell, count) in enumerate(zip(chosen_cells, unique_counts)):
+        if chosen_cell == 0:
+            out[point_indices[i] : point_indices[i + 1], :] = _generate_points_in_wedge(
+                points[wedge, :], n=count
+            )
+        else:
+            out[
+                point_indices[i] : point_indices[i + 1], :
+            ] = _generate_points_in_hexahedron(points[hexahedron, :], n=count)
+
+    return out

src/pyransame/volume.py

@@ -22,6 +22,7 @@ def random_volume_points(
     Supported cell types:
 
     - Hexahedron
+    - Pentagonal Prism
     - Pyramid
     - Tetrahedron
     - Voxel
@@ -110,6 +111,10 @@ def random_volume_points(
             points[
                 point_indices[i] : point_indices[i + 1], :
             ] = util._generate_points_in_hexahedron(c.points, count)
+        elif c.type == CellType.PENTAGONAL_PRISM:
+            points[
+                point_indices[i] : point_indices[i + 1], :
+            ] = util._generate_points_in_pentagonal_prism(c.points, count)
         else:
             raise NotImplementedError(
                 f"Random generation for {c.type.name} not yet supported"

tests/test_random_volume_points.py

@@ -40,6 +40,20 @@ def make_hexahedron():
     return pv.UnstructuredGrid(cells, celltypes, points)
 
 
+def make_pentagonal_prism():
+    angles = np.arange(5) * 2 * np.pi / 5  # angles of regular pentagon in radians
+    p = np.zeros(shape=(10, 3))
+    np.sin(angles, out=p[0:5, 0])
+    np.cos(angles, out=p[0:5, 1])
+
+    np.sin(angles, out=p[5:, 0])
+    np.cos(angles, out=p[5:, 1])
+    p[5:, 2] = 1.0
+    celltypes = [pv.CellType.PENTAGONAL_PRISM]
+    cells = [10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
+    return pv.UnstructuredGrid(cells, celltypes, p)
+
+
 def test_cell_types():
     mesh = pv.UniformGrid(dimensions=(4, 4, 4))
     assert mesh.get_cell(0).type == pv.CellType.VOXEL
@@ -61,6 +75,10 @@ def test_cell_types():
     assert mesh.get_cell(0).type == pv.CellType.HEXAHEDRON
     pyransame.random_volume_points(mesh, 20)
 
+    mesh = make_pentagonal_prism()
+    assert mesh.get_cell(0).type == pv.CellType.PENTAGONAL_PRISM
+    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

tests/test_util.py

@@ -10,6 +10,7 @@
 
 from pyransame.util import (
     _generate_points_in_hexahedron,
+    _generate_points_in_pentagonal_prism,
     _generate_points_in_pixel,
     _generate_points_in_polygon,
     _generate_points_in_pyramid,
@@ -269,3 +270,18 @@ def test_uniformity_hexahedra():
 
     points = _generate_points_in_hexahedron(p, 2000000)
     assert np.allclose(points.mean(axis=0), center, rtol=1e-3, atol=1e-3)
+
+
+def test_uniformity_pentagonal_prism():
+    angles = np.arange(5) * 2 * np.pi / 5  # angles of regular pentagon in radians
+    p = np.zeros(shape=(10, 3))
+    np.sin(angles, out=p[0:5, 0])
+    np.cos(angles, out=p[0:5, 1])
+
+    np.sin(angles, out=p[5:, 0])
+    np.cos(angles, out=p[5:, 1])
+    p[5:, 2] = 1.0
+
+    center = np.array([0.0, 0.0, 0.5])
+    points = _generate_points_in_pentagonal_prism(p, 2000000)
+    assert np.allclose(points.mean(axis=0), center, rtol=1e-3, atol=1e-3)