From 463299d04aa15c87baedaf33bf84c07b77090eed Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 8 Apr 2024 09:45:01 -0400 Subject: [PATCH] Polygon fix to better handle colinear points (#2935) Co-authored-by: Paul Romano --- openmc/model/surface_composite.py | 26 +++++++++++++--------- tests/unit_tests/test_surface_composite.py | 26 ++++++++++++++++++++++ 2 files changed, 41 insertions(+), 11 deletions(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 8f48b2f7117..5543a88e827 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -635,7 +635,7 @@ class XConeOneSided(CompositeSurface): z-coordinate of the apex. Defaults to 0. r2 : float, optional Parameter related to the aperture [:math:`\\rm cm^2`]. - It can be interpreted as the increase in the radius squared per cm along + It can be interpreted as the increase in the radius squared per cm along the cone's axis of revolution. up : bool Whether to select the side of the cone that extends to infinity in the @@ -695,7 +695,7 @@ class YConeOneSided(CompositeSurface): z-coordinate of the apex. Defaults to 0. r2 : float, optional Parameter related to the aperture [:math:`\\rm cm^2`]. - It can be interpreted as the increase in the radius squared per cm along + It can be interpreted as the increase in the radius squared per cm along the cone's axis of revolution. up : bool Whether to select the side of the cone that extends to infinity in the @@ -749,7 +749,7 @@ class ZConeOneSided(CompositeSurface): z-coordinate of the apex. Defaults to 0. r2 : float, optional Parameter related to the aperture [:math:`\\rm cm^2`]. - It can be interpreted as the increase in the radius squared per cm along + It can be interpreted as the increase in the radius squared per cm along the cone's axis of revolution. up : bool Whether to select the side of the cone that extends to infinity in the @@ -1029,9 +1029,9 @@ def _constrain_triangulation(self, points, depth=0): ------- None """ - # Only attempt the triangulation up to 3 times. - if depth > 2: - raise RuntimeError('Could not create a valid triangulation after 3' + # Only attempt the triangulation up to 5 times. + if depth > 4: + raise RuntimeError('Could not create a valid triangulation after 5' ' attempts') tri = Delaunay(points, qhull_options='QJ') @@ -1039,7 +1039,7 @@ def _constrain_triangulation(self, points, depth=0): # included in the triangulation, break it into two line segments. n = len(points) new_pts = [] - for i, j in zip(range(n), range(1, n +1)): + for i, j in zip(range(n), range(1, n + 1)): # If both vertices of any edge are not found in any simplex, insert # a new point between them. if not any([i in s and j % n in s for s in tri.simplices]): @@ -1077,7 +1077,8 @@ def _group_simplices(self, neighbor_map, group=None): return group # If group is empty, grab the next simplex in the dictionary and recurse if group is None: - sidx = next(iter(neighbor_map)) + # Start with smallest neighbor lists + sidx = sorted(neighbor_map.items(), key=lambda item: len(item[1]))[0][0] return self._group_simplices(neighbor_map, group=[sidx]) # Otherwise use the last simplex in the group else: @@ -1087,13 +1088,16 @@ def _group_simplices(self, neighbor_map, group=None): # For each neighbor check if it is part of the same convex # hull as the rest of the group. If yes, recurse. If no, continue on. for n in neighbors: - if n in group or neighbor_map.get(n, None) is None: + if n in group or neighbor_map.get(n) is None: continue test_group = group + [n] test_point_idx = np.unique(self._tri.simplices[test_group, :]) test_points = self.points[test_point_idx] - # If test_points are convex keep adding to this group - if len(test_points) == len(ConvexHull(test_points).vertices): + test_hull = ConvexHull(test_points, qhull_options='Qc') + pts_on_hull = len(test_hull.vertices) + len(test_hull.coplanar) + # If test_points are convex (including coplanar) keep adding to + # this group + if len(test_points) == pts_on_hull: group = self._group_simplices(neighbor_map, group=test_group) return group diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py index 73519c38368..19221212e06 100644 --- a/tests/unit_tests/test_surface_composite.py +++ b/tests/unit_tests/test_surface_composite.py @@ -399,6 +399,32 @@ def test_polygon(): with pytest.raises(ValueError): openmc.model.Polygon(rz_points) + # Test "M" shaped polygon + points = np.array([[8.5151581, -17.988337], + [10.381711000000001, -17.988337], + [12.744357, -24.288728000000003], + [15.119406000000001, -17.988337], + [16.985959, -17.988337], + [16.985959, -27.246687], + [15.764328, -27.246687], + [15.764328, -19.116951], + [13.376877, -25.466951], + [12.118039, -25.466951], + [9.7305877, -19.116951], + [9.7305877, -27.246687], + [8.5151581, -27.246687]]) + + # Test points inside and outside by using offset method + m_polygon = openmc.model.Polygon(points, basis='xz') + inner_pts = m_polygon.offset(-0.1).points + assert all([(pt[0], 0, pt[1]) in -m_polygon for pt in inner_pts]) + outer_pts = m_polygon.offset(0.1).points + assert all([(pt[0], 0, pt[1]) in +m_polygon for pt in outer_pts]) + + # Offset of -0.2 will cause self-intersection + with pytest.raises(ValueError): + m_polygon.offset(-0.2) + @pytest.mark.parametrize("axis", ["x", "y", "z"]) def test_cruciform_prism(axis):