From afbbc3d8457854114fb75d73e884251c8522e410 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Tue, 24 Jan 2023 23:00:27 -0500 Subject: [PATCH 1/8] try to constrain triangulation --- openmc/model/surface_composite.py | 54 +++++++++++++++++++++++++------ 1 file changed, 45 insertions(+), 9 deletions(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 002fad9d9e..e44b81ac44 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -670,10 +670,11 @@ class Polygon(CompositeSurface): check_length('points', points, 3) # Order the points counter-clockwise (necessary for offset method) - self._points = self._make_ccw(points) + points = self._make_ccw(points) - # Create a triangulation of the points. - self._tri = Delaunay(self._points, qhull_options='QJ') + # Create a constrained triangulation of the points. + # The constrained triangulation is set to _tri attribute + self._constrain_triangulation(points) # Decompose the polygon into groups of simplices forming convex subsets # and get the sets of (surface, operator) pairs defining the polygon @@ -723,7 +724,7 @@ class Polygon(CompositeSurface): @property def points(self): - return self._points + return self._tri.points @property def basis(self): @@ -738,7 +739,7 @@ class Polygon(CompositeSurface): # Get the unit vectors that point from one point in the polygon to the # next given that they are ordered counterclockwise and that the final # point is connected to the first point - tangents = np.diff(self._points, axis=0, append=[self._points[0, :]]) + tangents = np.diff(self.points, axis=0, append=[self.points[0, :]]) tangents /= np.linalg.norm(tangents, axis=-1, keepdims=True) # Rotate the tangent vectors clockwise by 90 degrees, which for a # counter-clockwise ordered polygon will produce the outward normal @@ -783,6 +784,41 @@ class Polygon(CompositeSurface): return points[::-1, :] + def _constrain_triangulation(self, points): + """Generate a constrained triangulation. + + Parameters + ---------- + points : np.ndarray (Nx2) + An Nx2 array of coordinate pairs describing the vertices. These + points represent a planar straight line graph. + + Returns + ------- + None + """ + + tri = Delaunay(points, qhull_options='QJ') + # Loop through the boundary edges of the polygon. If an edge is not + # 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)): + # 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]): + newpt = (points[i, :] + points[j % n, :]) / 2 + new_pts.append((j, newpt)) + + # If all the edges are included in the triangulation set it, otherwise + # try again with additional points inserted on offending edges + if not new_pts: + self._tri = tri + else: + for i, pt in new_pts[::-1]: + points = np.insert(points, i, pt, axis=0) + self._constrain_triangulation(points) + def _group_simplices(self, neighbor_map, group=None): """Generate a convex grouping of simplices. @@ -819,7 +855,7 @@ class Polygon(CompositeSurface): continue test_group = group + [n] test_point_idx = np.unique(self._tri.simplices[test_group, :]) - test_points = self._tri.points[test_point_idx] + 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): group = self._group_simplices(neighbor_map, group=test_group) @@ -902,8 +938,8 @@ class Polygon(CompositeSurface): # Get centroids of all the simplices and determine if they are inside # the polygon defined by input vertices or not. - centroids = np.mean(self._points[self._tri.simplices], axis=1) - in_polygon = Path(self._points).contains_points(centroids) + centroids = np.mean(self.points[self._tri.simplices], axis=1) + in_polygon = Path(self.points).contains_points(centroids) self._in_polygon = in_polygon # Build a map with keys of simplex indices inside the polygon whose @@ -931,7 +967,7 @@ class Polygon(CompositeSurface): # generate the convex hull and find the resulting surfaces and # unary operators that represent this convex subset of the polygon. idx = np.unique(self._tri.simplices[group, :]) - qhull = ConvexHull(self._tri.points[idx, :]) + qhull = ConvexHull(self.points[idx, :]) surf_ops = self._get_convex_hull_surfs(qhull) surfsets.append(surf_ops) return surfsets From 4fade673a74681e51e0ed198ac0f4f7bcd24fc26 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Wed, 25 Jan 2023 16:27:12 -0500 Subject: [PATCH 2/8] refactor to single validate_points method --- openmc/model/surface_composite.py | 63 +++++++++++++++++++------------ 1 file changed, 39 insertions(+), 24 deletions(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index e44b81ac44..04e0c4d80d 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -659,22 +659,10 @@ class Polygon(CompositeSurface): def __init__(self, points, basis='rz'): check_value('basis', basis, ('xy', 'yz', 'xz', 'rz')) self._basis = basis - points = np.asarray(points, dtype=float) - check_iterable_type('points', points, float, min_depth=2, max_depth=2) - check_length('points', points[0, :], 2, 2) - # If the last point is the same as the first, remove it and make sure - # there are still at least 3 points for a valid polygon. - if np.allclose(points[0, :], points[-1, :]): - points = points[:-1, :] - check_length('points', points, 3) - - # Order the points counter-clockwise (necessary for offset method) - points = self._make_ccw(points) - - # Create a constrained triangulation of the points. - # The constrained triangulation is set to _tri attribute - self._constrain_triangulation(points) + # Create a constrained triangulation of the validated points. + # The constrained triangulation is set to the _tri attribute + self._constrain_triangulation(self._validate_points(points)) # Decompose the polygon into groups of simplices forming convex subsets # and get the sets of (surface, operator) pairs defining the polygon @@ -762,8 +750,9 @@ class Polygon(CompositeSurface): def region(self): return self._region - def _make_ccw(self, points): - """Order a set of points counter-clockwise. + def _validate_points(self, points): + """Ensure the closed path defined by points does not intersect and is + oriented counter-clockwise. Parameters ---------- @@ -774,18 +763,40 @@ class Polygon(CompositeSurface): ------- ordered_points : the input points ordered counter-clockwise """ + points = np.asarray(points, dtype=float) + check_iterable_type('points', points, float, min_depth=2, max_depth=2) + check_length('points', points[0, :], 2, 2) + + # If the last point is the same as the first, remove it and make sure + # there are still at least 3 points for a valid polygon. + if np.allclose(points[0, :], points[-1, :]): + points = points[:-1, :] + check_length('points', points, 3) + + # Check if polygon is self-intersecting by comparing edges pairwise + n = len(points) + is_self_intersecting = False + for i in range(n): + p1x, p1y = points[i, :] + p2x, p2y = points[(i + 1) % n, :] + for j in range(i, n): + p3x, p3y = points[j, :] + p4x, p4y = points[(j + 1) % n, :] + + + # Order the points counter-clockwise (necessary for offset method) # Calculates twice the signed area of the polygon using the "Shoelace # Formula" https://en.wikipedia.org/wiki/Shoelace_formula - xpts, ypts = points.T - # If signed area is positive the curve is oriented counter-clockwise + xpts, ypts = points.T if np.sum(ypts*(np.roll(xpts, 1) - np.roll(xpts, -1))) > 0: return points return points[::-1, :] - def _constrain_triangulation(self, points): - """Generate a constrained triangulation. + def _constrain_triangulation(self, points, depth=0): + """Generate a constrained triangulation by ensuring all edges of the + Polygon are contained within the simplices. Parameters ---------- @@ -797,6 +808,10 @@ class Polygon(CompositeSurface): ------- None """ + # Only attempt the triangulation up to 3 times. + if depth > 2: + raise RuntimeError('Could not create a valid triangulation after 3' + ' attempts') tri = Delaunay(points, qhull_options='QJ') # Loop through the boundary edges of the polygon. If an edge is not @@ -805,19 +820,19 @@ class Polygon(CompositeSurface): new_pts = [] 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 + # a new point between them. if not any([i in s and j % n in s for s in tri.simplices]): newpt = (points[i, :] + points[j % n, :]) / 2 new_pts.append((j, newpt)) # If all the edges are included in the triangulation set it, otherwise - # try again with additional points inserted on offending edges + # try again with additional points inserted on offending edges. if not new_pts: self._tri = tri else: for i, pt in new_pts[::-1]: points = np.insert(points, i, pt, axis=0) - self._constrain_triangulation(points) + self._constrain_triangulation(points, depth=depth + 1) def _group_simplices(self, neighbor_map, group=None): """Generate a convex grouping of simplices. From e20e648690b33923e280476f97efd041f63e1435 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 30 Jan 2023 19:40:00 -0500 Subject: [PATCH 3/8] check for duplicate points --- openmc/model/surface_composite.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 04e0c4d80d..a570ec15d1 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -773,13 +773,16 @@ class Polygon(CompositeSurface): points = points[:-1, :] check_length('points', points, 3) + if len(points) != len(np.unique(points, axis=0)): + raise ValueError('Duplicate points were detected in the Polygon input') + # Check if polygon is self-intersecting by comparing edges pairwise n = len(points) is_self_intersecting = False for i in range(n): p1x, p1y = points[i, :] p2x, p2y = points[(i + 1) % n, :] - for j in range(i, n): + for j in range(i + 1, n): p3x, p3y = points[j, :] p4x, p4y = points[(j + 1) % n, :] From 78b0fbaa27b98410ecd0a5103a6aa27fdf1705e4 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Sun, 5 Feb 2023 22:02:14 -0500 Subject: [PATCH 4/8] finished polygon self-intersection check --- openmc/model/surface_composite.py | 91 ++++++++++++++++++++++++++----- 1 file changed, 76 insertions(+), 15 deletions(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index a570ec15d1..047cfe7fcd 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -776,26 +776,87 @@ class Polygon(CompositeSurface): if len(points) != len(np.unique(points, axis=0)): raise ValueError('Duplicate points were detected in the Polygon input') - # Check if polygon is self-intersecting by comparing edges pairwise - n = len(points) - is_self_intersecting = False - for i in range(n): - p1x, p1y = points[i, :] - p2x, p2y = points[(i + 1) % n, :] - for j in range(i + 1, n): - p3x, p3y = points[j, :] - p4x, p4y = points[(j + 1) % n, :] - - # Order the points counter-clockwise (necessary for offset method) # Calculates twice the signed area of the polygon using the "Shoelace # Formula" https://en.wikipedia.org/wiki/Shoelace_formula - # If signed area is positive the curve is oriented counter-clockwise + # If signed area is positive the curve is oriented counter-clockwise. + # If the signed area is negative the curve is oriented clockwise. xpts, ypts = points.T - if np.sum(ypts*(np.roll(xpts, 1) - np.roll(xpts, -1))) > 0: - return points + if np.sum(ypts*(np.roll(xpts, 1) - np.roll(xpts, -1))) < 0: + points = points[::-1, :] - return points[::-1, :] + + # Check if polygon is self-intersecting by comparing edges pairwise + n = len(points) + for i in range(n): + p0 = points[i, :] + p1 = points[(i + 1) % n, :] + for j in range(i + 1, n): + p2 = points[j, :] + p3 = points[(j + 1) % n, :] + # Compute orientation of p0 wrt p2->p3 line segment + cp0 = np.cross(p3-p0, p2-p0) + # Compute orientation of p1 wrt p2->p3 line segment + cp1 = np.cross(p3-p1, p2-p1) + # Compute orientation of p2 wrt p0->p1 line segment + cp2 = np.cross(p1-p2, p0-p2) + # Compute orientation of p3 wrt p0->p1 line segment + cp3 = np.cross(p1-p3, p0-p3) + + # Group cross products in an array and find out how many are 0 + cross_products = np.array([[cp0, cp1], [cp2, cp3]]) + cps_near_zero = np.isclose(cross_products, 0).astype(int) + num_zeros = np.sum(cps_near_zero) + + # Topologies of 2 finite line segments categorized by the number + # of zero-valued cross products: + # + # 0: No 3 points lie on the same line + # 1: 1 point lies on the same line defined by the other line + # segment, but is not coincident with either of the points + # 2: 2 points are coincident, but the line segments are not + # collinear which guarantees no intersection + # 3: not possible + # 4: Both line segments are collinear, simply need to check if + # they overlap or not + + if num_zeros == 0: + # If the orientations of p0 and p1 have opposite signs + # and the orientations of p2 and p3 have opposite signs + # then there is an intersection. + if all(np.prod(cross_products, axis=-1) < 0): + raise ValueError('Polygon cannot be self-intersecting') + continue + + elif num_zeros == 1: + # determine which line segment has 2 out of the 3 collinear + # points + idx = np.argwhere(np.sum(cps_near_zero, axis=-1) == 0) + if np.prod(cross_products[idx, :]) < 0: + raise ValueError('Polygon cannot be self-intersecting') + continue + + elif num_zeros == 2: + continue + + elif num_zeros == 4: + # Determine number of unique points, x span and y span for + # both line segments + #unique_pts = np.unique(np.vstack((p0, p1, p2, p3)), axis=0) + xmin1, xmax1 = min(p0[0], p1[0]), max(p0[0], p1[0]) + ymin1, ymax1 = min(p0[1], p1[1]), max(p0[1], p1[1]) + xmin2, xmax2 = min(p2[0], p3[0]), max(p2[0], p3[0]) + ymin2, ymax2 = min(p2[1], p3[1]), max(p2[1], p3[1]) + xlap = xmin1 < xmax2 and xmin2 < xmax1 + ylap = ymin1 < ymax2 and ymin2 < ymax1 + if xlap or ylap: + raise ValueError('Polygon cannot be self-intersecting') + continue + + else: + warnings.warn('Unclear if Polygon is self-intersecting') + + return points def _constrain_triangulation(self, points, depth=0): """Generate a constrained triangulation by ensuring all edges of the From c07386fc62bcfea422a03889821b9b3f053d546f Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Sun, 5 Feb 2023 22:15:10 -0500 Subject: [PATCH 5/8] added more tests for Polygon class --- tests/unit_tests/test_surface_composite.py | 54 ++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py index dbfaa62b17..3f94893edb 100644 --- a/tests/unit_tests/test_surface_composite.py +++ b/tests/unit_tests/test_surface_composite.py @@ -339,3 +339,57 @@ def test_polygon(): offset_star = star_poly.offset(.6) assert (0, 0, 0) in -offset_star assert any([(0, 0, 0) in reg for reg in offset_star.regions]) + + # check invalid Polygon input points + # duplicate points not just at start and end + rz_points = np.array([[6.88, 3.02], + [6.88, 2.72], + [6.88, 3.02], + [7.63, 0.0], + [5.75, 0.0], + [5.75, 1.22], + [7.63, 0.0], + [6.30, 1.22], + [6.30, 3.02], + [6.88, 3.02]]) + with pytest.raises(ValueError): + openmc.model.Polygon(rz_points) + + # segment traces back on previous segment + rz_points = np.array([[6.88, 3.02], + [6.88, 2.72], + [6.88, 2.32], + [6.88, 2.52], + [7.63, 0.0], + [5.75, 0.0], + [6.75, 0.0], + [5.75, 1.22], + [6.30, 1.22], + [6.30, 3.02], + [6.88, 3.02]]) + with pytest.raises(ValueError): + openmc.model.Polygon(rz_points) + + # segments intersect (line-line) + rz_points = np.array([[6.88, 3.02], + [5.88, 2.32], + [7.63, 0.0], + [5.75, 0.0], + [5.75, 1.22], + [6.30, 1.22], + [6.30, 3.02], + [6.88, 3.02]]) + with pytest.raises(ValueError): + openmc.model.Polygon(rz_points) + + # segments intersect (line-point) + rz_points = np.array([[6.88, 3.02], + [6.3, 2.32], + [7.63, 0.0], + [5.75, 0.0], + [5.75, 1.22], + [6.30, 1.22], + [6.30, 3.02], + [6.88, 3.02]]) + with pytest.raises(ValueError): + openmc.model.Polygon(rz_points) From 637080a0c2b33f93f26c6647977aef6d99daf849 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Mon, 6 Feb 2023 15:52:46 -0500 Subject: [PATCH 6/8] address comments from @hassec review --- openmc/model/surface_composite.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 047cfe7fcd..044ecf61e6 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -816,9 +816,14 @@ class Polygon(CompositeSurface): # segment, but is not coincident with either of the points # 2: 2 points are coincident, but the line segments are not # collinear which guarantees no intersection - # 3: not possible + # 3: not possible, except maybe floating point issues? # 4: Both line segments are collinear, simply need to check if # they overlap or not + # adapted from algorithm linked below and modified to only + # consider intersections on the interior of line segments as + # proper intersections: i.e. segments sharing end points do not + # count as intersections. + # https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/ if num_zeros == 0: # If the orientations of p0 and p1 have opposite signs @@ -839,7 +844,12 @@ class Polygon(CompositeSurface): elif num_zeros == 2: continue - elif num_zeros == 4: + elif num_zeros == 3: + warnings.warn('Unclear if Polygon is self-intersecting') + continue + + else: + # All 4 cross products are zero # Determine number of unique points, x span and y span for # both line segments #unique_pts = np.unique(np.vstack((p0, p1, p2, p3)), axis=0) @@ -853,9 +863,6 @@ class Polygon(CompositeSurface): raise ValueError('Polygon cannot be self-intersecting') continue - else: - warnings.warn('Unclear if Polygon is self-intersecting') - return points def _constrain_triangulation(self, points, depth=0): From 1a2d1237da0b46f8cc0dddeaccee3c7f1036381d Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Thu, 16 Feb 2023 12:27:15 -0500 Subject: [PATCH 7/8] Update openmc/model/surface_composite.py Co-authored-by: Jonathan Shimwell --- openmc/model/surface_composite.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 044ecf61e6..09c2f4fee9 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -785,7 +785,6 @@ class Polygon(CompositeSurface): if np.sum(ypts*(np.roll(xpts, 1) - np.roll(xpts, -1))) < 0: points = points[::-1, :] - # Check if polygon is self-intersecting by comparing edges pairwise n = len(points) for i in range(n): From 2f5428f739808abc67674fc9fc126e197f3201a9 Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Fri, 24 Mar 2023 12:13:00 -0400 Subject: [PATCH 8/8] Update openmc/model/surface_composite.py Co-authored-by: Patrick Shriwise --- openmc/model/surface_composite.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 09c2f4fee9..354a628832 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -851,7 +851,6 @@ class Polygon(CompositeSurface): # All 4 cross products are zero # Determine number of unique points, x span and y span for # both line segments - #unique_pts = np.unique(np.vstack((p0, p1, p2, p3)), axis=0) xmin1, xmax1 = min(p0[0], p1[0]), max(p0[0], p1[0]) ymin1, ymax1 = min(p0[1], p1[1]), max(p0[1], p1[1]) xmin2, xmax2 = min(p2[0], p3[0]), max(p2[0], p3[0])