Merge pull request #2362 from eepeterson/polygon_fix

try to constrain triangulation

Author: paulromano

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 @@ def boundary_type(self, boundary_type):
 
     @property
     def points(self):
-        return self._points
+        return self._tri.points
 
     @property
     def basis(self):
@@ -738,7 +727,7 @@ def _normals(self):
         # 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 @@ def regions(self):
     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 @@ def _make_ccw(self, points):
         -------
         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, :]
+
+        # 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, 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
 
-        # 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
+                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[::-1, :]
+        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 @@ def _group_simplices(self, neighbor_map, group=None):
                     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 @@ def _decompose_polygon_into_convex_sets(self):
 
         # 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 @@ def _decompose_polygon_into_convex_sets(self):
             # 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

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)