diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 002fad9d9e..354a628832 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -659,21 +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) - self._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 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 @@ -723,7 +712,7 @@ class Polygon(CompositeSurface): @property def points(self): - return self._points + return self._tri.points @property def basis(self): @@ -738,7 +727,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 @@ -761,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 ---------- @@ -773,15 +763,145 @@ 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) + + if len(points) != len(np.unique(points, axis=0)): + raise ValueError('Duplicate points were detected in the Polygon input') + + # 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 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: + points = points[::-1, :] - # If signed area is positive the curve is oriented counter-clockwise - if np.sum(ypts*(np.roll(xpts, 1) - np.roll(xpts, -1))) > 0: - return points + # 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) - return points[::-1, :] + # 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, 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 + # 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 == 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 + 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 + + return points + + def _constrain_triangulation(self, points, depth=0): + """Generate a constrained triangulation by ensuring all edges of the + Polygon are contained within the simplices. + + Parameters + ---------- + points : np.ndarray (Nx2) + An Nx2 array of coordinate pairs describing the vertices. These + points represent a planar straight line graph. + + Returns + ------- + 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 + # 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, depth=depth + 1) def _group_simplices(self, neighbor_map, group=None): """Generate a convex grouping of simplices. @@ -819,7 +939,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 +1022,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 +1051,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 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)