diff --git a/Credits.rst b/Credits.rst index a2376e64..2da7c61f 100644 --- a/Credits.rst +++ b/Credits.rst @@ -112,6 +112,7 @@ Jen Bradley * Added shape families for Archimedean, Catalan, and Johnson solids. * Added shape family for prisms and antiprisms. * Added shape family for equilateral pyramids and dipyramids. +* Added edges, edge_vectors, and num_edges methods. Domagoj Fijan * Rewrote point in polygon check to use NumPy vectorized operations. diff --git a/coxeter/shapes/convex_polyhedron.py b/coxeter/shapes/convex_polyhedron.py index bb3472e0..2b3bcac6 100644 --- a/coxeter/shapes/convex_polyhedron.py +++ b/coxeter/shapes/convex_polyhedron.py @@ -126,6 +126,12 @@ def asphericity(self): """float: Get the asphericity as defined in :cite:`Irrgang2017`.""" return self.mean_curvature * self.surface_area / (3 * self.volume) + @property + def num_edges(self): + """int: Get the number of edges.""" + # Calculate number of edges from Euler Characteristic + return self.num_vertices + self.num_faces - 2 + def is_inside(self, points): """Determine whether points are contained in this polyhedron. diff --git a/coxeter/shapes/polyhedron.py b/coxeter/shapes/polyhedron.py index 22cce8f2..3eafe706 100644 --- a/coxeter/shapes/polyhedron.py +++ b/coxeter/shapes/polyhedron.py @@ -4,6 +4,7 @@ """Defines a polyhedron.""" import warnings +from functools import cached_property import numpy as np import rowan @@ -356,9 +357,45 @@ def vertices(self): @property def faces(self): - """list(:class:`numpy.ndarray`): Get the polyhedron's faces.""" + """list(:class:`numpy.ndarray`): Get the polyhedron's faces. + + Results returned as vertex index lists. + """ return self._faces + @cached_property + def edges(self): + """:class:`numpy.ndarray`: Get the polyhedron's edges. + + Results returned as vertex index pairs, with each edge of the polyhedron + included exactly once. Edge (i,j) pairs are ordered by vertex index with i= 0) + + # Check that all items in the first column are greater than those in the second. + assert np.all(np.diff(poly.edges, axis=1) > 0) + + # Check the second column is in ascending order for each unique item in the first. + # For example, [[0,1],[0,3],[1,2]] is permitted but [[0,1],[0,3],[0,2]] is not. + edges = poly.edges + unique_values = unique_values = np.unique(edges[:, 0]) + assert all( + [ + np.all(np.diff(edges[edges[:, 0] == value, 1]) >= 0) + for value in unique_values + ] + ) + + # Check that there are no duplicate edges. This also double-checks the sorting + assert np.all(np.unique(poly.edges, axis=1) == poly.edges) + + # Check that the edges are immutable + try: + poly.edges[1] = [99, 99] + # If the assignment works, catch that: + assert poly.edges[1] != [99, 99] + except ValueError as ve: + assert "read-only" in str(ve) + + +def test_edge_lengths(): + known_shapes = { + "Tetrahedron": np.sqrt(2) * np.cbrt(3), + "Cube": 1, + "Octahedron": np.power(2, 5 / 6) * np.cbrt(3 / 8), + "Dodecahedron": np.power(2, 2 / 3) * np.cbrt(1 / (15 + np.sqrt(245))), + "Icosahedron": np.cbrt(9 / 5 - 3 / 5 * np.sqrt(5)), + } + for name, edgelength in known_shapes.items(): + poly = PlatonicFamily.get_shape(name) + # Check that edge lengths are correct + veclens = np.linalg.norm( + poly.vertices[poly.edges[:, 1]] - poly.vertices[poly.edges[:, 0]], axis=1 + ) + assert np.allclose(veclens, edgelength) + assert np.allclose(veclens, np.linalg.norm(poly.edge_vectors, axis=1)) + + +def test_num_edges_archimedean(): + known_shapes = { + "Cuboctahedron": 24, + "Icosidodecahedron": 60, + "Truncated Tetrahedron": 18, + "Truncated Octahedron": 36, + "Truncated Cube": 36, + "Truncated Icosahedron": 90, + "Truncated Dodecahedron": 90, + "Rhombicuboctahedron": 48, + "Rhombicosidodecahedron": 120, + "Truncated Cuboctahedron": 72, + "Truncated Icosidodecahedron": 180, + "Snub Cuboctahedron": 60, + "Snub Icosidodecahedron": 150, + } + for name, num_edges in known_shapes.items(): + poly = ArchimedeanFamily.get_shape(name) + assert poly.num_edges == num_edges + + +@given( + EllipsoidSurfaceStrategy, +) +def test_num_edges_polyhedron(points): + hull = ConvexHull(points) + poly = ConvexPolyhedron(points[hull.vertices]) + ppoly = Polyhedron(poly.vertices, poly.faces) + + # Calculate correct number of edges from euler characteristic + euler_characteristic_edge_count = ppoly.num_vertices + ppoly.num_faces - 2 + assert ppoly.num_edges == euler_characteristic_edge_count + + def test_curvature(): """Regression test against values computed with older method.""" # The shapes in the PlatonicFamily are normalized to unit volume.