Skip to content

Commit

Permalink
add polyhedron
Browse files Browse the repository at this point in the history
  • Loading branch information
MatthewFlamm authored Jun 14, 2023
1 parent fe5c4de commit eb2e278
Show file tree
Hide file tree
Showing 4 changed files with 96 additions and 0 deletions.
40 changes: 40 additions & 0 deletions src/pyransame/util.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,8 @@
from typing import Tuple

import numpy as np
import pyvista as pv
from vtk import mutable

import pyransame

Expand Down Expand Up @@ -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
4 changes: 4 additions & 0 deletions src/pyransame/volume.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
25 changes: 25 additions & 0 deletions tests/test_random_volume_points.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
27 changes: 27 additions & 0 deletions tests/test_util.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,15 @@
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,
_generate_points_in_hexahedron,
_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,
Expand Down Expand Up @@ -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)

0 comments on commit eb2e278

Please sign in to comment.