Polygon fix to better handle colinear points (#2935)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

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,17 +1029,17 @@ 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')
         # 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)):
+        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
 

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):