From f6152d82505492c2a2fb14b4b1b7f56c7df29912 Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Thu, 13 Oct 2022 16:28:20 -0400
Subject: [PATCH 01/15] working on other basis options

---
 openmc/model/surface_composite.py | 286 +++++++++++++++++++++++++++++-
 1 file changed, 285 insertions(+), 1 deletion(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 4c76c7786..b8f952f19 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -1,7 +1,9 @@
 from abc import ABC, abstractmethod
 from copy import copy
-from math import sqrt, pi, sin, cos
+from math import sqrt, pi, sin, cos, isclose
 import numpy as np
+from scipy.spatial import ConvexHull, Delaunay
+from matplotlib.path import Path
 
 import openmc
 from openmc.checkvalue import check_greater_than, check_value
@@ -621,3 +623,285 @@ class ZConeOneSided(CompositeSurface):
 
     __neg__ = XConeOneSided.__neg__
     __pos__ = XConeOneSided.__pos__
+
+
+class Polygon(CompositeSurface):
+    """Create a polygon composite surface from connected points.
+
+    Parameters
+    ----------
+    points : np.ndarray (Nx2)
+        Points defining the vertices of the polygon.
+    basis : str, {'rz', 'xy', 'yz', 'xz'}, optional
+        2D basis set for the polygon.
+
+    Attributes
+    ----------
+    """
+    _basis_surface_map = {}
+    _basis_surface_map['xy'] = 2*(openmc.XPlane, openmc.YPlane)
+    _basis_surface_map['yz'] = 2*(openmc.YPlane, openmc.ZPlane)
+    _basis_surface_map['xz'] = 2*(openmc.XPlane, openmc.ZPlane)
+
+    def __init__(self, points, basis='rz'):
+        # If the last point is the same as the first remove it and order the
+        # vertices in a counter-clockwise sense.
+        if np.allclose(points[0, :], points[-1, :]):
+            points = points[:-1, :]
+        self._points = self.make_ccw(np.asarray(points, dtype=float))
+        self._basis = basis
+
+        # Create a triangulation and convex hull of the points. The
+        # Polygon region will be primarily defined by the intersection
+        # of the convex hull with the intersection of the complements of the
+        # simplices outside the polygon, but inside the convex hull.
+        self._tri = Delaunay(self._points, qhull_options='QJ')
+        self._convex_hull = ConvexHull(self._points)
+        self._convex_hull_surfs = self.get_convex_hull_surfs()
+
+        # Get centroids of all the simplices and determine if they are inside
+        # the polygon defined by input vertices or not. If they are not, they
+        # are added to the convex_subsets
+        centroids = np.mean(self._points[self._tri.simplices], axis=1)
+        path = Path(self._points)
+        in_polygon = np.array([path.contains_point(c) for c in centroids])
+        self._in_polygon = in_polygon
+        # Loop through simplices that aren't in the polygon and add them to the
+        # surfaces that need to be removed
+        self._surfs_to_remove = []
+        for simplex in self._tri.simplices[~in_polygon, :]:
+            qhull = ConvexHull(self._points[simplex, :])
+            idx = self.get_ordered_simplex_indices(qhull)
+            eqns = qhull.equations[idx, :]
+            self._surfs_to_remove.append(self.get_convex_hull_surfs(eqns))
+
+        # Set surface names as required by CompositeSurface protocol
+        surfnames = []
+        for i, (surf, _) in enumerate(self._convex_hull_surfs):
+            setattr(self, f'hull_surf_{i}', surf)
+            surfnames.append(f'hull_surf_{i}')
+
+        i = 0
+        for surfs_ops in self._surfs_to_remove:
+            for surf, _ in surfs_ops:
+                setattr(self, f'aux_surf_{i}', surf)
+                surfnames.append(f'aux_surf_{i}')
+                i += 1
+
+        self._surfnames = tuple(surfnames)
+
+    def __neg__(self):
+        # inside convex surface and outside all convex surfaces formed from
+        # concave points
+        return self.region
+
+    def __pos__(self):
+        # outside convex hull or inside one of the convex shapes formed from
+        # concave points
+        return ~self.region
+
+    @property
+    def _surface_names(self):
+        return self._surfnames
+
+    @property
+    def points(self):
+        return self._points
+
+    @property
+    def basis(self):
+        return self._basis
+
+    @property
+    def tangents(self):
+        return np.diff(self._points, axis=0, append=[self._points[0, :]])
+
+    @property
+    def normals(self):
+        rotation = np.array([[0, 1], [-1, 0]])
+        tangents = self.tangents
+        tangents /= np.linalg.norm(tangents, axis=-1, keepdims=True)
+        return rotation.dot(tangents.T).T
+
+    @property
+    def hull_points(self):
+        return self._points[self._convex_hull.vertices]
+
+    @property
+    def hull_tangents(self):
+        pts = self._points[self._convex_hull.vertices]
+        return np.diff(pts, axis=0, append=[pts[0, :]])
+
+    @property
+    def hull_equations(self):
+        idx = self.get_ordered_simplex_indices()
+        return self._convex_hull.equations[idx, :]
+
+    @property
+    def convex_hull_surfs(self):
+        return self._convex_hull_surfs
+
+    @property
+    def hull_region(self):
+        surfs_ops = self.convex_hull_surfs
+        regions = [getattr(surf, op)() for surf, op in surfs_ops]
+        return openmc.Intersection(regions)
+
+    @property
+    def region(self):
+        hull_reg = self.hull_region
+        surfs_ops_sets = self._surfs_to_remove
+        complements = []
+        for surfs_ops in surfs_ops_sets:
+            regions = [getattr(surf, op)() for surf, op in surfs_ops]
+            complements.append(~openmc.Intersection(regions))
+        return hull_reg & openmc.Intersection(complements)
+
+    def offset(self, distance):
+        """Offset this polygon by a set distance
+
+        Parameters
+        ----------
+        distance : float
+            The distance to offset the polygon by. Positive is outward
+            (expanding) and negative is inward (shrinking).
+
+
+        Returns
+        -------
+        offset_polygon : openmc.model.Polygon
+        """
+        points = self.points
+        normals = self.normals
+        normals = np.insert(normals, 0, normals[-1, :], axis=0)
+        ndotv1 = np.sum(normals[:-1, :]*(points + distance*normals[:-1, :]),
+                        axis=-1, keepdims=True)
+        ndotv2 = np.sum(normals[1:, :]*(points + distance*normals[1:, :]),
+                        axis=-1, keepdims=True)
+
+        new_points = np.empty_like(points)
+        denom = normals[:-1, 0]*normals[1:, 1] - normals[:-1, 1]*normals[1:, 0]
+        new_points[:, 0] = normals[1:, 1]*ndotv1.T - normals[:-1, 1]*ndotv2.T
+        new_points[:, 1] = -normals[1:, 0]*ndotv1.T + normals[:-1, 0]*ndotv2.T
+        new_points /= denom[:, None]
+
+        return type(self)(new_points, basis=self.basis)
+
+    def get_ordered_simplex_indices(self, qhull=None):
+        """Return simplex indices in same order as ConvexHull.vertices
+
+        Parameters
+        ----------
+        qhull : scipy.spatial.ConvexHull, optional
+            A ConvexHull object.
+
+        Returns
+        -------
+        idxlist : np.array of ordered simplex indices
+        """
+        qhull = self._convex_hull if qhull is None else qhull
+        idxlist = []
+        verts = qhull.vertices
+        nverts = len(verts)
+        simplices = qhull.simplices
+        for i in range(nverts):
+            next_i = (i + 1) % nverts
+            tmp_verts = [verts[i], verts[next_i]]
+            for j, simplex in enumerate(simplices):
+                if all(idx in simplex for idx in tmp_verts):
+                    idxlist.append(j)
+        return np.array(idxlist)
+
+    def get_convex_hull_surfs(self, hull_equations=None):
+        """Generate a list of surfaces given by a set of linear equations
+
+        Parameters
+        ----------
+        hull_equations : np.ndarray, optional
+            An Nx3 array where N is the number of facets (or sides) that represent
+            the equations and each row is given by (nx, ny, c) where (nx, ny) is
+            the unit vector normal to the facet and c is a constant such that the
+            surface described by the equation is n dot x + c = 0.
+
+        Returns
+        -------
+        surfs_ops : list of (surface, operator) tuples
+
+        """
+        hull_equations = self.hull_equations if hull_equations is None else \
+            hull_equations
+        # Collect surface/operator pairs such that the intersection of the
+        # regions defined by these pairs is the inside of the polygon.
+        surfs_ops = []
+        # hull facet equation: dx*x + dy*y + c = 0
+        for dx, dy, c in hull_equations:
+            # default to negative halfspace operator for inside the polygon
+            op = '__neg__'
+            # Check if the facet is horizontal
+            if isclose(dx, 0):
+                if self.basis in ('xz', 'yz', 'rz'):
+                    surf = openmc.ZPlane(z0=-c/dy)
+                else:
+                    surf = openmc.YPlane(y0=-c/dy)
+                # if (0, 1).(dx, dy) < 0 we want positive halfspace instead
+                if dy < 0:
+                    op = '__pos__'
+            # Check if the facet is vertical
+            elif isclose(dy, 0):
+                if self.basis in ('xy', 'xz'):
+                    surf = openmc.XPlane(x0=-c/dx)
+                elif self.basis == 'yz':
+                    surf = openmc.YPlane(y0=-c/dx)
+                else:
+                    surf = openmc.ZCylinder(r=-c/dx)
+                # if (1, 0).(dx, dy) < 0 we want positive halfspace instead
+                if dx < 0:
+                    op = '__pos__'
+            # Otherwise the facet is at an angle
+            else:
+                y0 = -c/dy
+                r2 = dy**2 / dx**2
+                # Check if the *slope* of the facet is positive
+                if dy / dx < 0:
+                    if self.basis == 'rz':
+                        surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=True)
+                    else:
+                        raise NotImplementedError
+                    # if (1, -1).(dx, dy) < 0 we want positive halfspace instead
+                    if dx - dy < 0:
+                        op = '__pos__'
+                else:
+                    if self.basis == 'rz':
+                        surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=False)
+                    else:
+                        raise NotImplementedError
+                    # if (1, 1).(dx, dy) < 0 we want positive halfspace instead
+                    if dx + dy < 0:
+                        op = '__pos__'
+
+            surfs_ops.append((surf, op))
+
+        return surfs_ops
+
+    def make_ccw(self, verts):
+        """Determine whether the vertices are ordered counter-clockwise
+
+        Parameters
+        ----------
+        verts : np.ndarray (Nx2)
+            An Nx2 array of coordinate pairs describing the vertices.
+
+
+        Returns
+        -------
+        bool : True if vertices are in counter-clockwise order. False otherwise.
+
+        """
+        vector = np.empty(verts.shape[0])
+        vector[:-1] = (verts[1:, 0] - verts[:-1, 0])*(verts[1:, 1] + verts[:-1, 1])
+        vector[-1] = (verts[0, 0] - verts[-1, 0])*(verts[0, 1] + verts[-1, 1])
+
+        if np.sum(vector) < 0:
+            return verts
+
+        return verts[::-1, :]

From aa7aab1768b28a8d87ef4a7b762a04598f145155 Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Fri, 14 Oct 2022 15:38:27 -0400
Subject: [PATCH 02/15] method to break polygon into convex sets complete

---
 openmc/model/surface_composite.py | 311 +++++++++++++++++-------------
 1 file changed, 181 insertions(+), 130 deletions(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index b8f952f19..71f95c01c 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -638,17 +638,13 @@ class Polygon(CompositeSurface):
     Attributes
     ----------
     """
-    _basis_surface_map = {}
-    _basis_surface_map['xy'] = 2*(openmc.XPlane, openmc.YPlane)
-    _basis_surface_map['yz'] = 2*(openmc.YPlane, openmc.ZPlane)
-    _basis_surface_map['xz'] = 2*(openmc.XPlane, openmc.ZPlane)
 
     def __init__(self, points, basis='rz'):
         # If the last point is the same as the first remove it and order the
         # vertices in a counter-clockwise sense.
         if np.allclose(points[0, :], points[-1, :]):
             points = points[:-1, :]
-        self._points = self.make_ccw(np.asarray(points, dtype=float))
+        self._points = make_ccw(np.asarray(points, dtype=float))
         self._basis = basis
 
         # Create a triangulation and convex hull of the points. The
@@ -656,8 +652,6 @@ class Polygon(CompositeSurface):
         # of the convex hull with the intersection of the complements of the
         # simplices outside the polygon, but inside the convex hull.
         self._tri = Delaunay(self._points, qhull_options='QJ')
-        self._convex_hull = ConvexHull(self._points)
-        self._convex_hull_surfs = self.get_convex_hull_surfs()
 
         # Get centroids of all the simplices and determine if they are inside
         # the polygon defined by input vertices or not. If they are not, they
@@ -666,26 +660,35 @@ class Polygon(CompositeSurface):
         path = Path(self._points)
         in_polygon = np.array([path.contains_point(c) for c in centroids])
         self._in_polygon = in_polygon
-        # Loop through simplices that aren't in the polygon and add them to the
-        # surfaces that need to be removed
-        self._surfs_to_remove = []
-        for simplex in self._tri.simplices[~in_polygon, :]:
-            qhull = ConvexHull(self._points[simplex, :])
-            idx = self.get_ordered_simplex_indices(qhull)
-            eqns = qhull.equations[idx, :]
-            self._surfs_to_remove.append(self.get_convex_hull_surfs(eqns))
+
+        # ndict maps simplex indices to a list of their neighbors inside the
+        # polygon
+        ndict = {}
+        for i, nlist in enumerate(self._tri.neighbors):
+            if not in_polygon[i]:
+                continue
+            ndict[i] = [n for n in nlist if in_polygon[n] and n >=0]
+        #for key, value in ndict.items():
+        #    print(key, ' : ', value)
+
+        groups = group_wrapper(self._tri, ndict)
+        self._groups = groups
+        for g in groups:
+            print(g)
+
+        self._surfsets = []
+        for pts in get_ordered_points(self._tri, groups):
+            qhull = ConvexHull(pts)
+            surf_ops = get_convex_hull_surfs(qhull)
+            self._surfsets.append(surf_ops)
 
         # Set surface names as required by CompositeSurface protocol
         surfnames = []
-        for i, (surf, _) in enumerate(self._convex_hull_surfs):
-            setattr(self, f'hull_surf_{i}', surf)
-            surfnames.append(f'hull_surf_{i}')
-
         i = 0
-        for surfs_ops in self._surfs_to_remove:
+        for surfs_ops in self._surfsets:
             for surf, _ in surfs_ops:
-                setattr(self, f'aux_surf_{i}', surf)
-                surfnames.append(f'aux_surf_{i}')
+                setattr(self, f'surface_{i}', surf)
+                surfnames.append(f'surface_{i}')
                 i += 1
 
         self._surfnames = tuple(surfnames)
@@ -747,15 +750,18 @@ class Polygon(CompositeSurface):
         regions = [getattr(surf, op)() for surf, op in surfs_ops]
         return openmc.Intersection(regions)
 
+    @property
+    def regions(self):
+        regions = []
+        for surfs_ops in self._surfsets:
+            reg = openmc.Intersection([getattr(surf, op)() for surf, op in
+                                       surfs_ops])
+            regions.append(reg)
+        return regions
+
     @property
     def region(self):
-        hull_reg = self.hull_region
-        surfs_ops_sets = self._surfs_to_remove
-        complements = []
-        for surfs_ops in surfs_ops_sets:
-            regions = [getattr(surf, op)() for surf, op in surfs_ops]
-            complements.append(~openmc.Intersection(regions))
-        return hull_reg & openmc.Intersection(complements)
+        return openmc.Union(self.regions)
 
     def offset(self, distance):
         """Offset this polygon by a set distance
@@ -787,121 +793,166 @@ class Polygon(CompositeSurface):
 
         return type(self)(new_points, basis=self.basis)
 
-    def get_ordered_simplex_indices(self, qhull=None):
-        """Return simplex indices in same order as ConvexHull.vertices
+def get_ordered_simplex_indices(qhull):
+    """Return simplex indices in same order as ConvexHull.vertices
 
-        Parameters
-        ----------
-        qhull : scipy.spatial.ConvexHull, optional
-            A ConvexHull object.
+    Parameters
+    ----------
+    qhull : scipy.spatial.ConvexHull, optional
+        A ConvexHull object.
 
-        Returns
-        -------
-        idxlist : np.array of ordered simplex indices
-        """
-        qhull = self._convex_hull if qhull is None else qhull
-        idxlist = []
-        verts = qhull.vertices
-        nverts = len(verts)
-        simplices = qhull.simplices
-        for i in range(nverts):
-            next_i = (i + 1) % nverts
-            tmp_verts = [verts[i], verts[next_i]]
-            for j, simplex in enumerate(simplices):
-                if all(idx in simplex for idx in tmp_verts):
-                    idxlist.append(j)
-        return np.array(idxlist)
+    Returns
+    -------
+    idxlist : np.array of ordered simplex indices
+    """
+    idxlist = []
+    verts = qhull.vertices
+    nverts = len(verts)
+    simplices = qhull.simplices
+    for i in range(nverts):
+        next_i = (i + 1) % nverts
+        tmp_verts = [verts[i], verts[next_i]]
+        for j, simplex in enumerate(simplices):
+            if all(idx in simplex for idx in tmp_verts):
+                idxlist.append(j)
+    return np.array(idxlist)
 
-    def get_convex_hull_surfs(self, hull_equations=None):
-        """Generate a list of surfaces given by a set of linear equations
+def get_convex_hull_surfs(qhull, basis='rz'):
+    """Generate a list of surfaces given by a set of linear equations
 
-        Parameters
-        ----------
-        hull_equations : np.ndarray, optional
-            An Nx3 array where N is the number of facets (or sides) that represent
-            the equations and each row is given by (nx, ny, c) where (nx, ny) is
-            the unit vector normal to the facet and c is a constant such that the
-            surface described by the equation is n dot x + c = 0.
+    Parameters
+    ----------
+    hull_equations : np.ndarray, optional
+        An Nx3 array where N is the number of facets (or sides) that represent
+        the equations and each row is given by (nx, ny, c) where (nx, ny) is
+        the unit vector normal to the facet and c is a constant such that the
+        surface described by the equation is n dot x + c = 0.
 
-        Returns
-        -------
-        surfs_ops : list of (surface, operator) tuples
+    Returns
+    -------
+    surfs_ops : list of (surface, operator) tuples
 
-        """
-        hull_equations = self.hull_equations if hull_equations is None else \
-            hull_equations
-        # Collect surface/operator pairs such that the intersection of the
-        # regions defined by these pairs is the inside of the polygon.
-        surfs_ops = []
-        # hull facet equation: dx*x + dy*y + c = 0
-        for dx, dy, c in hull_equations:
-            # default to negative halfspace operator for inside the polygon
-            op = '__neg__'
-            # Check if the facet is horizontal
-            if isclose(dx, 0):
-                if self.basis in ('xz', 'yz', 'rz'):
-                    surf = openmc.ZPlane(z0=-c/dy)
-                else:
-                    surf = openmc.YPlane(y0=-c/dy)
-                # if (0, 1).(dx, dy) < 0 we want positive halfspace instead
-                if dy < 0:
-                    op = '__pos__'
-            # Check if the facet is vertical
-            elif isclose(dy, 0):
-                if self.basis in ('xy', 'xz'):
-                    surf = openmc.XPlane(x0=-c/dx)
-                elif self.basis == 'yz':
-                    surf = openmc.YPlane(y0=-c/dx)
-                else:
-                    surf = openmc.ZCylinder(r=-c/dx)
-                # if (1, 0).(dx, dy) < 0 we want positive halfspace instead
-                if dx < 0:
-                    op = '__pos__'
-            # Otherwise the facet is at an angle
+    """
+    idx = get_ordered_simplex_indices(qhull)
+    hull_equations = qhull.equations[idx, :]
+    # Collect surface/operator pairs such that the intersection of the
+    # regions defined by these pairs is the inside of the polygon.
+    surfs_ops = []
+    # hull facet equation: dx*x + dy*y + c = 0
+    for dx, dy, c in hull_equations:
+        # default to negative halfspace operator for inside the polygon
+        op = '__neg__'
+        # Check if the facet is horizontal
+        if isclose(dx, 0):
+            if basis in ('xz', 'yz', 'rz'):
+                surf = openmc.ZPlane(z0=-c/dy)
             else:
-                y0 = -c/dy
-                r2 = dy**2 / dx**2
-                # Check if the *slope* of the facet is positive
-                if dy / dx < 0:
-                    if self.basis == 'rz':
-                        surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=True)
-                    else:
-                        raise NotImplementedError
-                    # if (1, -1).(dx, dy) < 0 we want positive halfspace instead
-                    if dx - dy < 0:
-                        op = '__pos__'
+                surf = openmc.YPlane(y0=-c/dy)
+            # if (0, 1).(dx, dy) < 0 we want positive halfspace instead
+            if dy < 0:
+                op = '__pos__'
+        # Check if the facet is vertical
+        elif isclose(dy, 0):
+            if basis in ('xy', 'xz'):
+                surf = openmc.XPlane(x0=-c/dx)
+            elif basis == 'yz':
+                surf = openmc.YPlane(y0=-c/dx)
+            else:
+                surf = openmc.ZCylinder(r=-c/dx)
+            # if (1, 0).(dx, dy) < 0 we want positive halfspace instead
+            if dx < 0:
+                op = '__pos__'
+        # Otherwise the facet is at an angle
+        else:
+            y0 = -c/dy
+            r2 = dy**2 / dx**2
+            # Check if the *slope* of the facet is positive
+            if dy / dx < 0:
+                if basis == 'rz':
+                    surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=True)
                 else:
-                    if self.basis == 'rz':
-                        surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=False)
-                    else:
-                        raise NotImplementedError
-                    # if (1, 1).(dx, dy) < 0 we want positive halfspace instead
-                    if dx + dy < 0:
-                        op = '__pos__'
+                    raise NotImplementedError
+                # if (1, -1).(dx, dy) < 0 we want positive halfspace instead
+                if dx - dy < 0:
+                    op = '__pos__'
+            else:
+                if basis == 'rz':
+                    surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=False)
+                else:
+                    raise NotImplementedError
+                # if (1, 1).(dx, dy) < 0 we want positive halfspace instead
+                if dx + dy < 0:
+                    op = '__pos__'
 
-            surfs_ops.append((surf, op))
+        surfs_ops.append((surf, op))
 
-        return surfs_ops
+    return surfs_ops
 
-    def make_ccw(self, verts):
-        """Determine whether the vertices are ordered counter-clockwise
+def make_ccw(verts):
+    """Determine whether the vertices are ordered counter-clockwise
 
-        Parameters
-        ----------
-        verts : np.ndarray (Nx2)
-            An Nx2 array of coordinate pairs describing the vertices.
+    Parameters
+    ----------
+    verts : np.ndarray (Nx2)
+        An Nx2 array of coordinate pairs describing the vertices.
 
 
-        Returns
-        -------
-        bool : True if vertices are in counter-clockwise order. False otherwise.
+    Returns
+    -------
+    bool : True if vertices are in counter-clockwise order. False otherwise.
 
-        """
-        vector = np.empty(verts.shape[0])
-        vector[:-1] = (verts[1:, 0] - verts[:-1, 0])*(verts[1:, 1] + verts[:-1, 1])
-        vector[-1] = (verts[0, 0] - verts[-1, 0])*(verts[0, 1] + verts[-1, 1])
+    """
+    vector = np.empty(verts.shape[0])
+    vector[:-1] = (verts[1:, 0] - verts[:-1, 0])*(verts[1:, 1] + verts[:-1, 1])
+    vector[-1] = (verts[0, 0] - verts[-1, 0])*(verts[0, 1] + verts[-1, 1])
 
-        if np.sum(vector) < 0:
-            return verts
+    if np.sum(vector) < 0:
+        return verts
 
-        return verts[::-1, :]
+    return verts[::-1, :]
+
+
+def get_ordered_points(tri, groups):
+    points = []
+    for g in groups:
+        idx = np.unique(tri.simplices[g, :])
+        qhull = ConvexHull(tri.points[idx, :])
+        points.append(qhull.points[qhull.vertices, :])
+    return points
+
+def group_wrapper(tri, simp_dict):
+    groups = []
+    while simp_dict:
+        groups.append(group_simplices(tri, simp_dict))
+    return groups
+
+def group_simplices(tri, simp_dict, group=None):
+    """Generate a list of convex subsets"""
+    # If dictionary is empty there's nothing left to do
+    if not simp_dict:
+        return group
+    # If group is empty grab the next simplex in the dictionary and recurse
+    if group is None:
+        sidx = next(iter(simp_dict))
+        return group_simplices(tri, simp_dict, group=[sidx])
+    # Otherwise use the last simplex in the group 
+    else:
+        # Remove current simplex from dictionary since it is in a group
+        sidx = group[-1]
+        neighbors = simp_dict.pop(sidx, [])
+        # 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 simp_dict.get(n, None) is None:
+                continue
+            test_group = group + [n]
+            #print('group :', group)
+            #print('test_group :', test_group)
+            test_point_idx = np.unique(tri.simplices[test_group, :])
+            test_points = tri.points[test_point_idx]
+            if is_convex(test_points):
+                group = group_simplices(tri, simp_dict, group=test_group)
+        return group
+
+def is_convex(points):
+    return len(points) == len(ConvexHull(points).vertices)

From f6f2ed7b171c763d089550fec7ad7c80bef7c876 Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Sat, 15 Oct 2022 17:56:54 -0400
Subject: [PATCH 03/15] got other basis sets working

---
 openmc/model/surface_composite.py | 75 ++++++++++++++++---------------
 1 file changed, 38 insertions(+), 37 deletions(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 71f95c01c..b0c83b205 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -626,17 +626,21 @@ class ZConeOneSided(CompositeSurface):
 
 
 class Polygon(CompositeSurface):
-    """Create a polygon composite surface from connected points.
+    """Create a polygon composite surface from a path of closed points.
 
     Parameters
     ----------
-    points : np.ndarray (Nx2)
-        Points defining the vertices of the polygon.
+    points : np.ndarray
+        An Nx2 array of points defining the vertices of the polygon.
     basis : str, {'rz', 'xy', 'yz', 'xz'}, optional
         2D basis set for the polygon.
 
     Attributes
     ----------
+    points : np.ndarray
+        An Nx2 array of points defining the vertices of the polygon.
+    basis : str, {'rz', 'xy', 'yz', 'xz'}
+        2D basis set for the polygon.
     """
 
     def __init__(self, points, basis='rz'):
@@ -644,24 +648,24 @@ class Polygon(CompositeSurface):
         # vertices in a counter-clockwise sense.
         if np.allclose(points[0, :], points[-1, :]):
             points = points[:-1, :]
+
+        # TODO check if path is self intersecting, throw error if yes
         self._points = make_ccw(np.asarray(points, dtype=float))
+        check_value('basis', basis, ('xy', 'yz', 'xz', 'rz'))
         self._basis = basis
 
-        # Create a triangulation and convex hull of the points. The
-        # Polygon region will be primarily defined by the intersection
-        # of the convex hull with the intersection of the complements of the
-        # simplices outside the polygon, but inside the convex hull.
+        # Create a triangulation of the points.
         self._tri = Delaunay(self._points, qhull_options='QJ')
 
         # Get centroids of all the simplices and determine if they are inside
-        # the polygon defined by input vertices or not. If they are not, they
-        # are added to the convex_subsets
+        # the polygon defined by input vertices or not.
         centroids = np.mean(self._points[self._tri.simplices], axis=1)
         path = Path(self._points)
         in_polygon = np.array([path.contains_point(c) for c in centroids])
         self._in_polygon = in_polygon
 
-        # ndict maps simplex indices to a list of their neighbors inside the
+        # Build a map with keys of simplex indices inside the polygon whose
+        # values are lists of that simplex's neighbors inside the
         # polygon
         ndict = {}
         for i, nlist in enumerate(self._tri.neighbors):
@@ -679,7 +683,7 @@ class Polygon(CompositeSurface):
         self._surfsets = []
         for pts in get_ordered_points(self._tri, groups):
             qhull = ConvexHull(pts)
-            surf_ops = get_convex_hull_surfs(qhull)
+            surf_ops = get_convex_hull_surfs(qhull, basis=self.basis)
             self._surfsets.append(surf_ops)
 
         # Set surface names as required by CompositeSurface protocol
@@ -809,9 +813,8 @@ def get_ordered_simplex_indices(qhull):
     verts = qhull.vertices
     nverts = len(verts)
     simplices = qhull.simplices
-    for i in range(nverts):
-        next_i = (i + 1) % nverts
-        tmp_verts = [verts[i], verts[next_i]]
+    for i, vert in enumerate(verts):
+        tmp_verts = [vert, verts[(i + 1) % nverts]]
         for j, simplex in enumerate(simplices):
             if all(idx in simplex for idx in tmp_verts):
                 idxlist.append(j)
@@ -835,13 +838,12 @@ def get_convex_hull_surfs(qhull, basis='rz'):
     """
     idx = get_ordered_simplex_indices(qhull)
     hull_equations = qhull.equations[idx, :]
+    check_value('basis', basis, ('xy', 'yz', 'xz', 'rz'))
     # Collect surface/operator pairs such that the intersection of the
     # regions defined by these pairs is the inside of the polygon.
     surfs_ops = []
     # hull facet equation: dx*x + dy*y + c = 0
     for dx, dy, c in hull_equations:
-        # default to negative halfspace operator for inside the polygon
-        op = '__neg__'
         # Check if the facet is horizontal
         if isclose(dx, 0):
             if basis in ('xz', 'yz', 'rz'):
@@ -849,8 +851,7 @@ def get_convex_hull_surfs(qhull, basis='rz'):
             else:
                 surf = openmc.YPlane(y0=-c/dy)
             # if (0, 1).(dx, dy) < 0 we want positive halfspace instead
-            if dy < 0:
-                op = '__pos__'
+            op = '__pos__' if dy < 0 else '__neg__'
         # Check if the facet is vertical
         elif isclose(dy, 0):
             if basis in ('xy', 'xz'):
@@ -860,29 +861,29 @@ def get_convex_hull_surfs(qhull, basis='rz'):
             else:
                 surf = openmc.ZCylinder(r=-c/dx)
             # if (1, 0).(dx, dy) < 0 we want positive halfspace instead
-            if dx < 0:
-                op = '__pos__'
+            op = '__pos__' if dx < 0 else '__neg__'
         # Otherwise the facet is at an angle
         else:
-            y0 = -c/dy
-            r2 = dy**2 / dx**2
-            # Check if the *slope* of the facet is positive
-            if dy / dx < 0:
-                if basis == 'rz':
-                    surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=True)
-                else:
-                    raise NotImplementedError
-                # if (1, -1).(dx, dy) < 0 we want positive halfspace instead
-                if dx - dy < 0:
-                    op = '__pos__'
+            op = '__neg__'
+            if basis == 'xy':
+                surf = openmc.Plane(a=dx, b=dy, d=-c)
+            elif basis == 'yz':
+                surf = openmc.Plane(b=dx, c=dy, d=-c)
+            elif basis == 'xz':
+                surf = openmc.Plane(a=dx, c=dy, d=-c)
             else:
-                if basis == 'rz':
-                    surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=False)
-                else:
-                    raise NotImplementedError
-                # if (1, 1).(dx, dy) < 0 we want positive halfspace instead
-                if dx + dy < 0:
+                y0 = -c/dy
+                r2 = dy**2 / dx**2
+                # Check if the *slope* of the facet is positive. If dy/dx < 0
+                # then we want up to be True for the one-sided cones.
+                up = dy / dx < 0
+                surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=up)
+                # if (1, -1).(dx, dy) < 0 for up cones we want positive halfspace
+                # if (1, 1).(dx, dy) < 0 for down cones we want positive halfspace
+                if (up and dx - dy < 0) or (not up and dx + dy < 0):
                     op = '__pos__'
+                else:
+                    op = '__neg__'
 
         surfs_ops.append((surf, op))
 

From 3fe2bf4cb650dab3c4ab76da258d587923520c41 Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Sat, 15 Oct 2022 18:22:47 -0400
Subject: [PATCH 04/15] cleaning up unnecessary methods

---
 openmc/model/surface_composite.py | 66 +++++++------------------------
 1 file changed, 15 insertions(+), 51 deletions(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index b0c83b205..72c625ba9 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -660,26 +660,24 @@ 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)
-        path = Path(self._points)
-        in_polygon = np.array([path.contains_point(c) for c in centroids])
+        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
-        # values are lists of that simplex's neighbors inside the
+        # values are lists of that simplex's neighbors also inside the
         # polygon
         ndict = {}
         for i, nlist in enumerate(self._tri.neighbors):
             if not in_polygon[i]:
                 continue
             ndict[i] = [n for n in nlist if in_polygon[n] and n >=0]
-        #for key, value in ndict.items():
-        #    print(key, ' : ', value)
 
+        # Get the groups of simplices forming convex polygons whose union
+        # comprises the full input polygon.
         groups = group_wrapper(self._tri, ndict)
         self._groups = groups
-        for g in groups:
-            print(g)
 
+        # Get the sets of surface, operator pairs defining the polygon
         self._surfsets = []
         for pts in get_ordered_points(self._tri, groups):
             qhull = ConvexHull(pts)
@@ -698,14 +696,10 @@ class Polygon(CompositeSurface):
         self._surfnames = tuple(surfnames)
 
     def __neg__(self):
-        # inside convex surface and outside all convex surfaces formed from
-        # concave points
-        return self.region
+        return self._region
 
     def __pos__(self):
-        # outside convex hull or inside one of the convex shapes formed from
-        # concave points
-        return ~self.region
+        return ~self._region
 
     @property
     def _surface_names(self):
@@ -720,52 +714,22 @@ class Polygon(CompositeSurface):
         return self._basis
 
     @property
-    def tangents(self):
-        return np.diff(self._points, axis=0, append=[self._points[0, :]])
-
-    @property
-    def normals(self):
+    def _normals(self):
         rotation = np.array([[0, 1], [-1, 0]])
-        tangents = self.tangents
+        tangents = np.diff(self._points, axis=0, append=[self._points[0, :]])
         tangents /= np.linalg.norm(tangents, axis=-1, keepdims=True)
         return rotation.dot(tangents.T).T
 
     @property
-    def hull_points(self):
-        return self._points[self._convex_hull.vertices]
-
-    @property
-    def hull_tangents(self):
-        pts = self._points[self._convex_hull.vertices]
-        return np.diff(pts, axis=0, append=[pts[0, :]])
-
-    @property
-    def hull_equations(self):
-        idx = self.get_ordered_simplex_indices()
-        return self._convex_hull.equations[idx, :]
-
-    @property
-    def convex_hull_surfs(self):
-        return self._convex_hull_surfs
-
-    @property
-    def hull_region(self):
-        surfs_ops = self.convex_hull_surfs
-        regions = [getattr(surf, op)() for surf, op in surfs_ops]
-        return openmc.Intersection(regions)
-
-    @property
-    def regions(self):
+    def _regions(self):
         regions = []
         for surfs_ops in self._surfsets:
-            reg = openmc.Intersection([getattr(surf, op)() for surf, op in
-                                       surfs_ops])
-            regions.append(reg)
-        return regions
+            regions.append([getattr(surf, op)() for surf, op in surfs_ops])
+        return [openmc.Intersection(regs) for regs in regions]
 
     @property
-    def region(self):
-        return openmc.Union(self.regions)
+    def _region(self):
+        return openmc.Union(self._regions)
 
     def offset(self, distance):
         """Offset this polygon by a set distance
@@ -782,7 +746,7 @@ class Polygon(CompositeSurface):
         offset_polygon : openmc.model.Polygon
         """
         points = self.points
-        normals = self.normals
+        normals = self._normals
         normals = np.insert(normals, 0, normals[-1, :], axis=0)
         ndotv1 = np.sum(normals[:-1, :]*(points + distance*normals[:-1, :]),
                         axis=-1, keepdims=True)

From a26dcc65b728b297d2e1cf09b1069784ad78902f Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Sun, 16 Oct 2022 21:38:29 -0400
Subject: [PATCH 05/15] cleaning up/consolidating functions

---
 openmc/model/surface_composite.py | 181 ++++++++++++++++--------------
 1 file changed, 95 insertions(+), 86 deletions(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 72c625ba9..2036950f5 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -657,38 +657,17 @@ class Polygon(CompositeSurface):
         # Create a triangulation of the points.
         self._tri = Delaunay(self._points, qhull_options='QJ')
 
-        # 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)
-        self._in_polygon = in_polygon
+        # Decompose the polygon into groups of simplices forming convex subsets
+        self._groups = self._decompose_polygon_into_convex_sets()
 
-        # Build a map with keys of simplex indices inside the polygon whose
-        # values are lists of that simplex's neighbors also inside the
-        # polygon
-        ndict = {}
-        for i, nlist in enumerate(self._tri.neighbors):
-            if not in_polygon[i]:
-                continue
-            ndict[i] = [n for n in nlist if in_polygon[n] and n >=0]
-
-        # Get the groups of simplices forming convex polygons whose union
-        # comprises the full input polygon.
-        groups = group_wrapper(self._tri, ndict)
-        self._groups = groups
-
-        # Get the sets of surface, operator pairs defining the polygon
-        self._surfsets = []
-        for pts in get_ordered_points(self._tri, groups):
-            qhull = ConvexHull(pts)
-            surf_ops = get_convex_hull_surfs(qhull, basis=self.basis)
-            self._surfsets.append(surf_ops)
+        # Get the sets of (surface, operator) pairs defining the polygon
+        self._surfsets = self._get_surfsets()
 
         # Set surface names as required by CompositeSurface protocol
         surfnames = []
         i = 0
-        for surfs_ops in self._surfsets:
-            for surf, _ in surfs_ops:
+        for surfset in self._surfsets:
+            for surf, op in surfset:
                 setattr(self, f'surface_{i}', surf)
                 surfnames.append(f'surface_{i}')
                 i += 1
@@ -715,6 +694,7 @@ class Polygon(CompositeSurface):
 
     @property
     def _normals(self):
+        """Generate the outward normal unit vectors for the polygon."""
         rotation = np.array([[0, 1], [-1, 0]])
         tangents = np.diff(self._points, axis=0, append=[self._points[0, :]])
         tangents /= np.linalg.norm(tangents, axis=-1, keepdims=True)
@@ -722,6 +702,7 @@ class Polygon(CompositeSurface):
 
     @property
     def _regions(self):
+        """Generate a list of regions whose union represents the polygon."""
         regions = []
         for surfs_ops in self._surfsets:
             regions.append([getattr(surf, op)() for surf, op in surfs_ops])
@@ -731,6 +712,57 @@ class Polygon(CompositeSurface):
     def _region(self):
         return openmc.Union(self._regions)
 
+    def _decompose_polygon_into_convex_sets(self):
+        """Decompose the Polygon into a set of convex polygons.
+
+        Returns
+        -------
+        list of sets of simplices
+        """
+
+        # 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)
+        self._in_polygon = in_polygon
+
+        # Build a map with keys of simplex indices inside the polygon whose
+        # values are lists of that simplex's neighbors also inside the
+        # polygon
+        neighbor_map = {}
+        for i, nlist in enumerate(self._tri.neighbors):
+            if not in_polygon[i]:
+                continue
+            neighbor_map[i] = [n for n in nlist if in_polygon[n] and n >=0]
+
+        # Get the groups of simplices forming convex polygons whose union
+        # comprises the full input polygon. While there are still simplices
+        # left in the neighbor map, group them together into convex sets.
+        groups = []
+        while neighbor_map:
+            groups.append(group_simplices(self._tri, neighbor_map))
+
+        return groups
+
+    def _get_surfsets(self):
+        """Generate lists of surface, operator pairs for the sub-polygons.
+
+        Returns
+        -------
+        surfsets : a list of lists of surface, operator pairs
+        """
+
+        surfsets = []
+        for g in self._groups:
+            # Find all the unique points in the convex group of simplices,
+            # 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[g, :])
+            qhull = ConvexHull(self._tri.points[idx, :])
+            surf_ops = get_convex_hull_surfs(qhull, basis=self.basis)
+            surfsets.append(surf_ops)
+        return surfsets
+
     def offset(self, distance):
         """Offset this polygon by a set distance
 
@@ -745,6 +777,8 @@ class Polygon(CompositeSurface):
         -------
         offset_polygon : openmc.model.Polygon
         """
+        # Get the points of the polygon and outward normals such that
+        # normals[i] corresponds to the edge between points[i-1] and points[i]
         points = self.points
         normals = self._normals
         normals = np.insert(normals, 0, normals[-1, :], axis=0)
@@ -761,53 +795,28 @@ class Polygon(CompositeSurface):
 
         return type(self)(new_points, basis=self.basis)
 
-def get_ordered_simplex_indices(qhull):
-    """Return simplex indices in same order as ConvexHull.vertices
-
-    Parameters
-    ----------
-    qhull : scipy.spatial.ConvexHull, optional
-        A ConvexHull object.
-
-    Returns
-    -------
-    idxlist : np.array of ordered simplex indices
-    """
-    idxlist = []
-    verts = qhull.vertices
-    nverts = len(verts)
-    simplices = qhull.simplices
-    for i, vert in enumerate(verts):
-        tmp_verts = [vert, verts[(i + 1) % nverts]]
-        for j, simplex in enumerate(simplices):
-            if all(idx in simplex for idx in tmp_verts):
-                idxlist.append(j)
-    return np.array(idxlist)
 
 def get_convex_hull_surfs(qhull, basis='rz'):
     """Generate a list of surfaces given by a set of linear equations
 
     Parameters
     ----------
-    hull_equations : np.ndarray, optional
-        An Nx3 array where N is the number of facets (or sides) that represent
-        the equations and each row is given by (nx, ny, c) where (nx, ny) is
-        the unit vector normal to the facet and c is a constant such that the
-        surface described by the equation is n dot x + c = 0.
+    qhull : scipy.spatial.ConvexHull
+        A ConvexHull object representing the sub-region of the polygon.
+    basis : str, {'rz', 'xy', 'yz', 'xz'}, optional
+        2D basis set for the polygon.
 
     Returns
     -------
     surfs_ops : list of (surface, operator) tuples
 
     """
-    idx = get_ordered_simplex_indices(qhull)
-    hull_equations = qhull.equations[idx, :]
     check_value('basis', basis, ('xy', 'yz', 'xz', 'rz'))
     # Collect surface/operator pairs such that the intersection of the
     # regions defined by these pairs is the inside of the polygon.
     surfs_ops = []
     # hull facet equation: dx*x + dy*y + c = 0
-    for dx, dy, c in hull_equations:
+    for dx, dy, c in qhull.equations:
         # Check if the facet is horizontal
         if isclose(dx, 0):
             if basis in ('xz', 'yz', 'rz'):
@@ -853,6 +862,7 @@ def get_convex_hull_surfs(qhull, basis='rz'):
 
     return surfs_ops
 
+
 def make_ccw(verts):
     """Determine whether the vertices are ordered counter-clockwise
 
@@ -877,47 +887,46 @@ def make_ccw(verts):
     return verts[::-1, :]
 
 
-def get_ordered_points(tri, groups):
-    points = []
-    for g in groups:
-        idx = np.unique(tri.simplices[g, :])
-        qhull = ConvexHull(tri.points[idx, :])
-        points.append(qhull.points[qhull.vertices, :])
-    return points
+def group_simplices(tri, neighbor_map, group=None):
+    """Generate a list of convex groups of simplices.
 
-def group_wrapper(tri, simp_dict):
-    groups = []
-    while simp_dict:
-        groups.append(group_simplices(tri, simp_dict))
-    return groups
+    Parameters
+    ----------
+    tri : scipy.spatial.Delaunay
+        A Delaunay triangulation of points generated from the polygon.
+    neighbor_map : dict
+        A map whose keys are simplex indices for simplices inside the polygon
+        and whose values are a list of simplex indices that neighbor this
+        simplex and are also inside the polygon.
+    group : list
+        A list of simplex indices that comprise the current convex group.
 
-def group_simplices(tri, simp_dict, group=None):
-    """Generate a list of convex subsets"""
-    # If dictionary is empty there's nothing left to do
-    if not simp_dict:
+    Returns
+    -------
+    group : list
+        The list of simplex indices that comprise the complete convex group.
+    """
+    # If neighbor_map is empty there's nothing left to do
+    if not neighbor_map:
         return group
-    # If group is empty grab the next simplex in the dictionary and recurse
+    # If group is empty, grab the next simplex in the dictionary and recurse
     if group is None:
-        sidx = next(iter(simp_dict))
-        return group_simplices(tri, simp_dict, group=[sidx])
+        sidx = next(iter(neighbor_map))
+        return group_simplices(tri, neighbor_map, group=[sidx])
     # Otherwise use the last simplex in the group 
     else:
-        # Remove current simplex from dictionary since it is in a group
         sidx = group[-1]
-        neighbors = simp_dict.pop(sidx, [])
+        # Remove current simplex from dictionary since it is in a group
+        neighbors = neighbor_map.pop(sidx, [])
         # 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 simp_dict.get(n, None) is None:
+            if n in group or neighbor_map.get(n, None) is None:
                 continue
             test_group = group + [n]
-            #print('group :', group)
-            #print('test_group :', test_group)
             test_point_idx = np.unique(tri.simplices[test_group, :])
             test_points = tri.points[test_point_idx]
-            if is_convex(test_points):
-                group = group_simplices(tri, simp_dict, group=test_group)
+            # If test_points are convex keep adding to this group
+            if len(test_points) == len(ConvexHull(test_points).vertices):
+                group = group_simplices(tri, neighbor_map, group=test_group)
         return group
-
-def is_convex(points):
-    return len(points) == len(ConvexHull(points).vertices)

From bf58a46a9cc79e6d52cb84f3fbc9ba6ddc93077c Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Tue, 18 Oct 2022 14:56:26 -0400
Subject: [PATCH 06/15] cleaning up and moving all methods to instance methods
 rather than module

---
 openmc/model/surface_composite.py          | 319 ++++++++++-----------
 tests/unit_tests/test_surface_composite.py |   4 +
 2 files changed, 151 insertions(+), 172 deletions(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 2036950f5..258a7bae7 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -6,7 +6,8 @@ from scipy.spatial import ConvexHull, Delaunay
 from matplotlib.path import Path
 
 import openmc
-from openmc.checkvalue import check_greater_than, check_value
+from openmc.checkvalue import (check_greater_than, check_value,
+                               check_iterable_type, check_length)
 
 
 class CompositeSurface(ABC):
@@ -641,27 +642,34 @@ class Polygon(CompositeSurface):
         An Nx2 array of points defining the vertices of the polygon.
     basis : str, {'rz', 'xy', 'yz', 'xz'}
         2D basis set for the polygon.
+    regions : list of openmc.Region
+        A list of openmc.Region objects, one for each of the convex polygons
+        formed during the decomposition of the input polygon.
+    region : openmc.Union
+        The union of all the regions comprising the polygon.
     """
 
     def __init__(self, points, basis='rz'):
-        # If the last point is the same as the first remove it and order the
-        # vertices in a counter-clockwise sense.
-        if np.allclose(points[0, :], points[-1, :]):
-            points = points[:-1, :]
-
-        # TODO check if path is self intersecting, throw error if yes
-        self._points = make_ccw(np.asarray(points, dtype=float))
         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)
+
+        self._points = points
 
         # Create a triangulation of the points.
         self._tri = Delaunay(self._points, qhull_options='QJ')
 
         # Decompose the polygon into groups of simplices forming convex subsets
-        self._groups = self._decompose_polygon_into_convex_sets()
-
-        # Get the sets of (surface, operator) pairs defining the polygon
-        self._surfsets = self._get_surfsets()
+        # and get the sets of (surface, operator) pairs defining the polygon
+        self._surfsets = self._decompose_polygon_into_convex_sets()
 
         # Set surface names as required by CompositeSurface protocol
         surfnames = []
@@ -671,9 +679,17 @@ class Polygon(CompositeSurface):
                 setattr(self, f'surface_{i}', surf)
                 surfnames.append(f'surface_{i}')
                 i += 1
-
         self._surfnames = tuple(surfnames)
 
+        # Generate a list of regions whose union represents the polygon.
+        regions = []
+        for surfs_ops in self._surfsets:
+            regions.append([getattr(surf, op)() for surf, op in surfs_ops])
+        self._regions = [openmc.Intersection(regs) for regs in regions]
+
+        # Create the union of all the convex subsets
+        self._region = openmc.Union(self._regions)
+
     def __neg__(self):
         return self._region
 
@@ -701,23 +717,125 @@ class Polygon(CompositeSurface):
         return rotation.dot(tangents.T).T
 
     @property
-    def _regions(self):
-        """Generate a list of regions whose union represents the polygon."""
-        regions = []
-        for surfs_ops in self._surfsets:
-            regions.append([getattr(surf, op)() for surf, op in surfs_ops])
-        return [openmc.Intersection(regs) for regs in regions]
+    def regions(self):
+        return self._regions
 
     @property
-    def _region(self):
-        return openmc.Union(self._regions)
+    def region(self):
+        return self._region
+
+    def _group_simplices(self, neighbor_map, group=None):
+        """Generate a convex grouping of simplices.
+
+        Parameters
+        ----------
+        neighbor_map : dict
+            A map whose keys are simplex indices for simplices inside the polygon
+            and whose values are a list of simplex indices that neighbor this
+            simplex and are also inside the polygon.
+        group : list
+            A list of simplex indices that comprise the current convex group.
+
+        Returns
+        -------
+        group : list
+            The list of simplex indices that comprise the complete convex group.
+        """
+        # If neighbor_map is empty there's nothing left to do
+        if not neighbor_map:
+            return group
+        # If group is empty, grab the next simplex in the dictionary and recurse
+        if group is None:
+            sidx = next(iter(neighbor_map))
+            return self._group_simplices(neighbor_map, group=[sidx])
+        # Otherwise use the last simplex in the group 
+        else:
+            sidx = group[-1]
+            # Remove current simplex from dictionary since it is in a group
+            neighbors = neighbor_map.pop(sidx, [])
+            # 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:
+                    continue
+                test_group = group + [n]
+                test_point_idx = np.unique(self._tri.simplices[test_group, :])
+                test_points = self._tri.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)
+            return group
+
+    def _get_convex_hull_surfs(self, qhull):
+        """Generate a list of surfaces given by a set of linear equations
+
+        Parameters
+        ----------
+        qhull : scipy.spatial.ConvexHull
+            A ConvexHull object representing the sub-region of the polygon.
+
+        Returns
+        -------
+        surfs_ops : list of (surface, operator) tuples
+
+        """
+        basis = self.basis
+        # Collect surface/operator pairs such that the intersection of the
+        # regions defined by these pairs is the inside of the polygon.
+        surfs_ops = []
+        # hull facet equation: dx*x + dy*y + c = 0
+        for dx, dy, c in qhull.equations:
+            # Check if the facet is horizontal
+            if isclose(dx, 0):
+                if basis in ('xz', 'yz', 'rz'):
+                    surf = openmc.ZPlane(z0=-c/dy)
+                else:
+                    surf = openmc.YPlane(y0=-c/dy)
+                # if (0, 1).(dx, dy) < 0 we want positive halfspace instead
+                op = '__pos__' if dy < 0 else '__neg__'
+            # Check if the facet is vertical
+            elif isclose(dy, 0):
+                if basis in ('xy', 'xz'):
+                    surf = openmc.XPlane(x0=-c/dx)
+                elif basis == 'yz':
+                    surf = openmc.YPlane(y0=-c/dx)
+                else:
+                    surf = openmc.ZCylinder(r=-c/dx)
+                # if (1, 0).(dx, dy) < 0 we want positive halfspace instead
+                op = '__pos__' if dx < 0 else '__neg__'
+            # Otherwise the facet is at an angle
+            else:
+                op = '__neg__'
+                if basis == 'xy':
+                    surf = openmc.Plane(a=dx, b=dy, d=-c)
+                elif basis == 'yz':
+                    surf = openmc.Plane(b=dx, c=dy, d=-c)
+                elif basis == 'xz':
+                    surf = openmc.Plane(a=dx, c=dy, d=-c)
+                else:
+                    y0 = -c/dy
+                    r2 = dy**2 / dx**2
+                    # Check if the *slope* of the facet is positive. If dy/dx < 0
+                    # then we want up to be True for the one-sided cones.
+                    up = dy / dx < 0
+                    surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=up)
+                    # if (1, -1).(dx, dy) < 0 for up cones we want positive halfspace
+                    # if (1, 1).(dx, dy) < 0 for down cones we want positive halfspace
+                    if (up and dx - dy < 0) or (not up and dx + dy < 0):
+                        op = '__pos__'
+                    else:
+                        op = '__neg__'
+
+            surfs_ops.append((surf, op))
+
+        return surfs_ops
 
     def _decompose_polygon_into_convex_sets(self):
         """Decompose the Polygon into a set of convex polygons.
 
         Returns
         -------
-        list of sets of simplices
+        surfsets : a list of lists of surface, operator pairs
         """
 
         # Get centroids of all the simplices and determine if they are inside
@@ -740,26 +858,19 @@ class Polygon(CompositeSurface):
         # left in the neighbor map, group them together into convex sets.
         groups = []
         while neighbor_map:
-            groups.append(group_simplices(self._tri, neighbor_map))
-
-        return groups
-
-    def _get_surfsets(self):
-        """Generate lists of surface, operator pairs for the sub-polygons.
-
-        Returns
-        -------
-        surfsets : a list of lists of surface, operator pairs
-        """
+            groups.append(self._group_simplices(neighbor_map))
+        self._groups = groups
 
+        # Generate lists of (surface, operator) pairs for each convex
+        # sub-region.
         surfsets = []
-        for g in self._groups:
+        for group in groups:
             # Find all the unique points in the convex group of simplices,
             # 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[g, :])
+            idx = np.unique(self._tri.simplices[group, :])
             qhull = ConvexHull(self._tri.points[idx, :])
-            surf_ops = get_convex_hull_surfs(qhull, basis=self.basis)
+            surf_ops = self._get_convex_hull_surfs(qhull)
             surfsets.append(surf_ops)
         return surfsets
 
@@ -794,139 +905,3 @@ class Polygon(CompositeSurface):
         new_points /= denom[:, None]
 
         return type(self)(new_points, basis=self.basis)
-
-
-def get_convex_hull_surfs(qhull, basis='rz'):
-    """Generate a list of surfaces given by a set of linear equations
-
-    Parameters
-    ----------
-    qhull : scipy.spatial.ConvexHull
-        A ConvexHull object representing the sub-region of the polygon.
-    basis : str, {'rz', 'xy', 'yz', 'xz'}, optional
-        2D basis set for the polygon.
-
-    Returns
-    -------
-    surfs_ops : list of (surface, operator) tuples
-
-    """
-    check_value('basis', basis, ('xy', 'yz', 'xz', 'rz'))
-    # Collect surface/operator pairs such that the intersection of the
-    # regions defined by these pairs is the inside of the polygon.
-    surfs_ops = []
-    # hull facet equation: dx*x + dy*y + c = 0
-    for dx, dy, c in qhull.equations:
-        # Check if the facet is horizontal
-        if isclose(dx, 0):
-            if basis in ('xz', 'yz', 'rz'):
-                surf = openmc.ZPlane(z0=-c/dy)
-            else:
-                surf = openmc.YPlane(y0=-c/dy)
-            # if (0, 1).(dx, dy) < 0 we want positive halfspace instead
-            op = '__pos__' if dy < 0 else '__neg__'
-        # Check if the facet is vertical
-        elif isclose(dy, 0):
-            if basis in ('xy', 'xz'):
-                surf = openmc.XPlane(x0=-c/dx)
-            elif basis == 'yz':
-                surf = openmc.YPlane(y0=-c/dx)
-            else:
-                surf = openmc.ZCylinder(r=-c/dx)
-            # if (1, 0).(dx, dy) < 0 we want positive halfspace instead
-            op = '__pos__' if dx < 0 else '__neg__'
-        # Otherwise the facet is at an angle
-        else:
-            op = '__neg__'
-            if basis == 'xy':
-                surf = openmc.Plane(a=dx, b=dy, d=-c)
-            elif basis == 'yz':
-                surf = openmc.Plane(b=dx, c=dy, d=-c)
-            elif basis == 'xz':
-                surf = openmc.Plane(a=dx, c=dy, d=-c)
-            else:
-                y0 = -c/dy
-                r2 = dy**2 / dx**2
-                # Check if the *slope* of the facet is positive. If dy/dx < 0
-                # then we want up to be True for the one-sided cones.
-                up = dy / dx < 0
-                surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=up)
-                # if (1, -1).(dx, dy) < 0 for up cones we want positive halfspace
-                # if (1, 1).(dx, dy) < 0 for down cones we want positive halfspace
-                if (up and dx - dy < 0) or (not up and dx + dy < 0):
-                    op = '__pos__'
-                else:
-                    op = '__neg__'
-
-        surfs_ops.append((surf, op))
-
-    return surfs_ops
-
-
-def make_ccw(verts):
-    """Determine whether the vertices are ordered counter-clockwise
-
-    Parameters
-    ----------
-    verts : np.ndarray (Nx2)
-        An Nx2 array of coordinate pairs describing the vertices.
-
-
-    Returns
-    -------
-    bool : True if vertices are in counter-clockwise order. False otherwise.
-
-    """
-    vector = np.empty(verts.shape[0])
-    vector[:-1] = (verts[1:, 0] - verts[:-1, 0])*(verts[1:, 1] + verts[:-1, 1])
-    vector[-1] = (verts[0, 0] - verts[-1, 0])*(verts[0, 1] + verts[-1, 1])
-
-    if np.sum(vector) < 0:
-        return verts
-
-    return verts[::-1, :]
-
-
-def group_simplices(tri, neighbor_map, group=None):
-    """Generate a list of convex groups of simplices.
-
-    Parameters
-    ----------
-    tri : scipy.spatial.Delaunay
-        A Delaunay triangulation of points generated from the polygon.
-    neighbor_map : dict
-        A map whose keys are simplex indices for simplices inside the polygon
-        and whose values are a list of simplex indices that neighbor this
-        simplex and are also inside the polygon.
-    group : list
-        A list of simplex indices that comprise the current convex group.
-
-    Returns
-    -------
-    group : list
-        The list of simplex indices that comprise the complete convex group.
-    """
-    # If neighbor_map is empty there's nothing left to do
-    if not neighbor_map:
-        return group
-    # If group is empty, grab the next simplex in the dictionary and recurse
-    if group is None:
-        sidx = next(iter(neighbor_map))
-        return group_simplices(tri, neighbor_map, group=[sidx])
-    # Otherwise use the last simplex in the group 
-    else:
-        sidx = group[-1]
-        # Remove current simplex from dictionary since it is in a group
-        neighbors = neighbor_map.pop(sidx, [])
-        # 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:
-                continue
-            test_group = group + [n]
-            test_point_idx = np.unique(tri.simplices[test_group, :])
-            test_points = tri.points[test_point_idx]
-            # If test_points are convex keep adding to this group
-            if len(test_points) == len(ConvexHull(test_points).vertices):
-                group = group_simplices(tri, 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 4d680ae7c..0c618250e 100644
--- a/tests/unit_tests/test_surface_composite.py
+++ b/tests/unit_tests/test_surface_composite.py
@@ -314,3 +314,7 @@ def test_isogonal_octagon(axis, plane_tb, plane_lr, axis_idx):
 
     # Make sure repr works
     repr(s)
+
+@pytest.mark.parametrize("basis", [("xz",), ("yz",), ("xy",), ("rz",)])
+def test_polygon(basis=basis):
+

From 13b1bb88d030d2f79ea3e6ebc7ed18b812c332b4 Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Tue, 18 Oct 2022 17:01:22 -0400
Subject: [PATCH 07/15] added unit test

---
 tests/unit_tests/test_surface_composite.py | 24 ++++++++++++++++++++--
 1 file changed, 22 insertions(+), 2 deletions(-)

diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py
index 0c618250e..2a77ac5f9 100644
--- a/tests/unit_tests/test_surface_composite.py
+++ b/tests/unit_tests/test_surface_composite.py
@@ -315,6 +315,26 @@ def test_isogonal_octagon(axis, plane_tb, plane_lr, axis_idx):
     # Make sure repr works
     repr(s)
 
-@pytest.mark.parametrize("basis", [("xz",), ("yz",), ("xy",), ("rz",)])
-def test_polygon(basis=basis):
+def test_polygon():
+    # define a 5 pointed star centered on 1, 1
+    star = np.array([[1.        , 2.        ],
+                     [0.70610737, 1.4045085 ],
+                     [0.04894348, 1.30901699],
+                     [0.52447174, 0.8454915 ],
+                     [0.41221475, 0.19098301],
+                     [1.        , 0.5       ],
+                     [1.58778525, 0.19098301],
+                     [1.47552826, 0.8454915 ],
+                     [1.95105652, 1.30901699],
+                     [1.29389263, 1.4045085 ],
+                     [1.        , 2.        ]])
+    points_in = [(1, 1, 0), (0, 1, 1), (1, 0, 1), (.707, .707, 1)]
+    for i, basis in enumerate(('xy', 'yz', 'xz', 'rz')):
+        star_poly = openmc.model.Polygon(star, basis=basis)
+        assert points_in[i] in -star_poly
+        assert points_in[i] not in +star_poly
+        assert (0, 0, 0) not in -star_poly
+        if basis != 'rz':
+            assert(0, 0, 0) in -star_poly.offset(.6)
+
 

From 299525984fc04e0479ffa6ca1d945ca5337f5852 Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Tue, 18 Oct 2022 20:56:58 -0400
Subject: [PATCH 08/15] simplified offset method

---
 openmc/model/surface_composite.py | 50 ++++++++++++++++++++-----------
 1 file changed, 32 insertions(+), 18 deletions(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 258a7bae7..82ebbfa03 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -662,7 +662,8 @@ class Polygon(CompositeSurface):
             points = points[:-1, :]
         check_length('points', points, 3)
 
-        self._points = points
+        # 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')
@@ -724,6 +725,29 @@ class Polygon(CompositeSurface):
     def region(self):
         return self._region
 
+    def _make_ccw(self, points):
+        """Order a set of points counter-clockwise.
+
+        Parameters
+        ----------
+        points : np.ndarray (Nx2)
+            An Nx2 array of coordinate pairs describing the vertices.
+
+        Returns
+        -------
+        ordered_points : the input points ordered counter-clockwise
+        """
+        vector = np.empty(points.shape[0])
+        this_x, this_y = points[:-1, 0], points[:-1, 1]
+        next_x, next_y = points[1:, 0], points[1:, 1]
+        vector[:-1] = (next_x - this_x)*(next_y + this_y)
+        vector[-1] = (this_x[0] - next_x[-1])*(this_y[0] + next_y[-1])
+
+        if np.sum(vector) < 0:
+            return points
+
+        return points[::-1, :]
+
     def _group_simplices(self, neighbor_map, group=None):
         """Generate a convex grouping of simplices.
 
@@ -883,25 +907,15 @@ class Polygon(CompositeSurface):
             The distance to offset the polygon by. Positive is outward
             (expanding) and negative is inward (shrinking).
 
-
         Returns
         -------
         offset_polygon : openmc.model.Polygon
         """
-        # Get the points of the polygon and outward normals such that
-        # normals[i] corresponds to the edge between points[i-1] and points[i]
-        points = self.points
-        normals = self._normals
-        normals = np.insert(normals, 0, normals[-1, :], axis=0)
-        ndotv1 = np.sum(normals[:-1, :]*(points + distance*normals[:-1, :]),
-                        axis=-1, keepdims=True)
-        ndotv2 = np.sum(normals[1:, :]*(points + distance*normals[1:, :]),
-                        axis=-1, keepdims=True)
+        normals = np.insert(self._normals, 0, self._normals[-1, :], axis=0)
+        cos2theta = np.sum(normals[1:, :]*normals[:-1, :], axis=-1, keepdims=True)
+        costheta = np.cos(np.arccos(cos2theta) / 2)
+        nvec = (normals[1:, :] + normals[:-1, :])
+        unit_nvec = nvec / np.linalg.norm(nvec, axis=-1, keepdims=True)
+        disp_vec = distance / costheta * unit_nvec
 
-        new_points = np.empty_like(points)
-        denom = normals[:-1, 0]*normals[1:, 1] - normals[:-1, 1]*normals[1:, 0]
-        new_points[:, 0] = normals[1:, 1]*ndotv1.T - normals[:-1, 1]*ndotv2.T
-        new_points[:, 1] = -normals[1:, 0]*ndotv1.T + normals[:-1, 0]*ndotv2.T
-        new_points /= denom[:, None]
-
-        return type(self)(new_points, basis=self.basis)
+        return type(self)(self.points + disp_vec, basis=self.basis)

From fdc3cfcfbba9fe4f48221850757b2c86f415b0d3 Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Thu, 20 Oct 2022 12:46:29 -0400
Subject: [PATCH 09/15] fixed proper boundary_type setting?

---
 openmc/model/surface_composite.py | 26 ++++++++++++++++++++------
 1 file changed, 20 insertions(+), 6 deletions(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 82ebbfa03..aeb77c868 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -676,16 +676,18 @@ class Polygon(CompositeSurface):
         surfnames = []
         i = 0
         for surfset in self._surfsets:
-            for surf, op in surfset:
-                setattr(self, f'surface_{i}', surf)
-                surfnames.append(f'surface_{i}')
-                i += 1
+            print(surfset)
+            for surf, op, on_boundary in surfset:
+                if on_boundary:
+                    setattr(self, f'surface_{i}', surf)
+                    surfnames.append(f'surface_{i}')
+                    i += 1
         self._surfnames = tuple(surfnames)
 
         # Generate a list of regions whose union represents the polygon.
         regions = []
         for surfs_ops in self._surfsets:
-            regions.append([getattr(surf, op)() for surf, op in surfs_ops])
+            regions.append([getattr(surf, op)() for surf, op, _ in surfs_ops])
         self._regions = [openmc.Intersection(regs) for regs in regions]
 
         # Create the union of all the convex subsets
@@ -717,6 +719,14 @@ class Polygon(CompositeSurface):
         tangents /= np.linalg.norm(tangents, axis=-1, keepdims=True)
         return rotation.dot(tangents.T).T
 
+    @property
+    def _equations(self):
+        normals = self._normals
+        equations = np.empty((normals.shape[0], 3))
+        equations[:, 0:2] = normals
+        equations[:, 2] = -np.sum(normals*self.points, axis=-1)
+        return equations
+
     @property
     def regions(self):
         return self._regions
@@ -804,11 +814,15 @@ class Polygon(CompositeSurface):
 
         """
         basis = self.basis
+        boundary_eqns = self._equations
         # Collect surface/operator pairs such that the intersection of the
         # regions defined by these pairs is the inside of the polygon.
         surfs_ops = []
         # hull facet equation: dx*x + dy*y + c = 0
         for dx, dy, c in qhull.equations:
+            # check if this facet is on the boundary of the polygon
+            facet_eq = np.array([dx, dy, c])
+            on_boundary = any([np.allclose(facet_eq, eq) for eq in boundary_eqns])
             # Check if the facet is horizontal
             if isclose(dx, 0):
                 if basis in ('xz', 'yz', 'rz'):
@@ -850,7 +864,7 @@ class Polygon(CompositeSurface):
                     else:
                         op = '__neg__'
 
-            surfs_ops.append((surf, op))
+            surfs_ops.append((surf, op, on_boundary))
 
         return surfs_ops
 

From a587f5e1c9966190d1e0f984e41d28145ae3ec2b Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Thu, 20 Oct 2022 12:50:33 -0400
Subject: [PATCH 10/15] forgot to remove print

---
 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 aeb77c868..8c47ac130 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -676,7 +676,6 @@ class Polygon(CompositeSurface):
         surfnames = []
         i = 0
         for surfset in self._surfsets:
-            print(surfset)
             for surf, op, on_boundary in surfset:
                 if on_boundary:
                     setattr(self, f'surface_{i}', surf)

From 8ddb045c7d6501db3a85acadd1d9c46f972d6aa2 Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Thu, 20 Oct 2022 15:28:26 -0400
Subject: [PATCH 11/15] added warning about non-transmissive boundary_type

---
 openmc/model/surface_composite.py | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 8c47ac130..922e60f0a 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -4,6 +4,7 @@ from math import sqrt, pi, sin, cos, isclose
 import numpy as np
 from scipy.spatial import ConvexHull, Delaunay
 from matplotlib.path import Path
+import warnings
 
 import openmc
 from openmc.checkvalue import (check_greater_than, check_value,
@@ -702,6 +703,18 @@ class Polygon(CompositeSurface):
     def _surface_names(self):
         return self._surfnames
 
+    @CompositeSurface.boundary_type.setter
+    def boundary_type(self, boundary_type):
+        if boundary_type != 'transmission':
+            warnings.warn("Setting boundary_type to a value other than "
+                          "'transmission' on Polygon composite surfaces can "
+                          "result in unintended behavior. Please use the "
+                          "regions property of the Polygon to generate "
+                          "individual openmc.Cell objects to avoid unwanted "
+                          "behavior.")
+        for name in self._surface_names:
+            getattr(self, name).boundary_type = boundary_type
+
     @property
     def points(self):
         return self._points

From b0d71b31466b4c1d835e28d11b48b28d0d730bc0 Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Tue, 8 Nov 2022 21:21:23 -0500
Subject: [PATCH 12/15] Apply suggestions from @paulromano code review

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
---
 openmc/model/surface_composite.py          | 10 +++++-----
 tests/unit_tests/test_surface_composite.py |  2 +-
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 922e60f0a..3ddaa6653 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -634,17 +634,17 @@ class Polygon(CompositeSurface):
     ----------
     points : np.ndarray
         An Nx2 array of points defining the vertices of the polygon.
-    basis : str, {'rz', 'xy', 'yz', 'xz'}, optional
+    basis : {'rz', 'xy', 'yz', 'xz'}, optional
         2D basis set for the polygon.
 
     Attributes
     ----------
     points : np.ndarray
         An Nx2 array of points defining the vertices of the polygon.
-    basis : str, {'rz', 'xy', 'yz', 'xz'}
+    basis : {'rz', 'xy', 'yz', 'xz'}
         2D basis set for the polygon.
     regions : list of openmc.Region
-        A list of openmc.Region objects, one for each of the convex polygons
+        A list of :class:`openmc.Region` objects, one for each of the convex polygons
         formed during the decomposition of the input polygon.
     region : openmc.Union
         The union of all the regions comprising the polygon.
@@ -726,7 +726,7 @@ class Polygon(CompositeSurface):
     @property
     def _normals(self):
         """Generate the outward normal unit vectors for the polygon."""
-        rotation = np.array([[0, 1], [-1, 0]])
+        rotation = np.array([[0., 1.], [-1., 0.]])
         tangents = np.diff(self._points, axis=0, append=[self._points[0, :]])
         tangents /= np.linalg.norm(tangents, axis=-1, keepdims=True)
         return rotation.dot(tangents.T).T
@@ -735,7 +735,7 @@ class Polygon(CompositeSurface):
     def _equations(self):
         normals = self._normals
         equations = np.empty((normals.shape[0], 3))
-        equations[:, 0:2] = normals
+        equations[:, :2] = normals
         equations[:, 2] = -np.sum(normals*self.points, axis=-1)
         return equations
 
diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py
index 2a77ac5f9..14579df22 100644
--- a/tests/unit_tests/test_surface_composite.py
+++ b/tests/unit_tests/test_surface_composite.py
@@ -335,6 +335,6 @@ def test_polygon():
         assert points_in[i] not in +star_poly
         assert (0, 0, 0) not in -star_poly
         if basis != 'rz':
-            assert(0, 0, 0) in -star_poly.offset(.6)
+            assert (0, 0, 0) in -star_poly.offset(.6)
 
 

From d314bedbd7951be8e8844e439c0349de9b621040 Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Tue, 8 Nov 2022 21:40:16 -0500
Subject: [PATCH 13/15] improving documentation and style

---
 docs/source/pythonapi/model.rst   |  1 +
 openmc/model/surface_composite.py | 23 ++++++++++++++---------
 2 files changed, 15 insertions(+), 9 deletions(-)

diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst
index 636935c6b..cdabee396 100644
--- a/docs/source/pythonapi/model.rst
+++ b/docs/source/pythonapi/model.rst
@@ -31,6 +31,7 @@ Composite Surfaces
    openmc.model.XConeOneSided
    openmc.model.YConeOneSided
    openmc.model.ZConeOneSided
+   openmc.model.Polygon
 
 TRISO Fuel Modeling
 -------------------
diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 3ddaa6653..5229e6f3d 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -1,10 +1,12 @@
 from abc import ABC, abstractmethod
 from copy import copy
 from math import sqrt, pi, sin, cos, isclose
+import warnings
+import operator
+
 import numpy as np
 from scipy.spatial import ConvexHull, Delaunay
 from matplotlib.path import Path
-import warnings
 
 import openmc
 from openmc.checkvalue import (check_greater_than, check_value,
@@ -635,7 +637,11 @@ class Polygon(CompositeSurface):
     points : np.ndarray
         An Nx2 array of points defining the vertices of the polygon.
     basis : {'rz', 'xy', 'yz', 'xz'}, optional
-        2D basis set for the polygon.
+        2D basis set for the polygon. The polygon is two dimensional and has
+        infinite extent in the third (unspecified) dimension. For example, the
+        'xy' basis produces a polygon with infinite extent in the +/- z
+        direction. For the 'rz' basis the phi extent is infinite, thus forming
+        an axisymmetric surface.
 
     Attributes
     ----------
@@ -687,7 +693,7 @@ class Polygon(CompositeSurface):
         # Generate a list of regions whose union represents the polygon.
         regions = []
         for surfs_ops in self._surfsets:
-            regions.append([getattr(surf, op)() for surf, op, _ in surfs_ops])
+            regions.append([op(surf) for surf, op, _ in surfs_ops])
         self._regions = [openmc.Intersection(regs) for regs in regions]
 
         # Create the union of all the convex subsets
@@ -842,7 +848,7 @@ class Polygon(CompositeSurface):
                 else:
                     surf = openmc.YPlane(y0=-c/dy)
                 # if (0, 1).(dx, dy) < 0 we want positive halfspace instead
-                op = '__pos__' if dy < 0 else '__neg__'
+                op = operator.pos if dy < 0 else operator.neg
             # Check if the facet is vertical
             elif isclose(dy, 0):
                 if basis in ('xy', 'xz'):
@@ -852,10 +858,10 @@ class Polygon(CompositeSurface):
                 else:
                     surf = openmc.ZCylinder(r=-c/dx)
                 # if (1, 0).(dx, dy) < 0 we want positive halfspace instead
-                op = '__pos__' if dx < 0 else '__neg__'
+                op = operator.pos if dx < 0 else operator.neg
             # Otherwise the facet is at an angle
             else:
-                op = '__neg__'
+                op = operator.neg
                 if basis == 'xy':
                     surf = openmc.Plane(a=dx, b=dy, d=-c)
                 elif basis == 'yz':
@@ -871,10 +877,9 @@ class Polygon(CompositeSurface):
                     surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=up)
                     # if (1, -1).(dx, dy) < 0 for up cones we want positive halfspace
                     # if (1, 1).(dx, dy) < 0 for down cones we want positive halfspace
+                    # otherwise we keep the negative halfspace operator
                     if (up and dx - dy < 0) or (not up and dx + dy < 0):
-                        op = '__pos__'
-                    else:
-                        op = '__neg__'
+                        op = operator.pos
 
             surfs_ops.append((surf, op, on_boundary))
 

From 2adbe971942f134babcfae92090f8fc9bc374fad Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Thu, 10 Nov 2022 11:14:20 -0500
Subject: [PATCH 14/15] updated algorithms and documentation

---
 openmc/model/surface_composite.py          | 19 +++++++++++++------
 tests/unit_tests/test_surface_composite.py |  7 ++++---
 2 files changed, 17 insertions(+), 9 deletions(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 5229e6f3d..5d1553a38 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -732,9 +732,17 @@ class Polygon(CompositeSurface):
     @property
     def _normals(self):
         """Generate the outward normal unit vectors for the polygon."""
+        # Rotation matrix for 90 degree clockwise rotation (-90 degrees about z
+        # axis for an 'xy' basis).
         rotation = np.array([[0., 1.], [-1., 0.]])
+        # 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.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
+        # vectors.
         return rotation.dot(tangents.T).T
 
     @property
@@ -765,13 +773,12 @@ class Polygon(CompositeSurface):
         -------
         ordered_points : the input points ordered counter-clockwise
         """
-        vector = np.empty(points.shape[0])
-        this_x, this_y = points[:-1, 0], points[:-1, 1]
-        next_x, next_y = points[1:, 0], points[1:, 1]
-        vector[:-1] = (next_x - this_x)*(next_y + this_y)
-        vector[-1] = (this_x[0] - next_x[-1])*(this_y[0] + next_y[-1])
+        # Calculates twice the signed area of the polygon using the "Shoelace
+        # Formula" https://en.wikipedia.org/wiki/Shoelace_formula
+        xpts, ypts = points.T
 
-        if np.sum(vector) < 0:
+        # 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
 
         return points[::-1, :]
diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py
index 14579df22..dbfaa62b1 100644
--- a/tests/unit_tests/test_surface_composite.py
+++ b/tests/unit_tests/test_surface_composite.py
@@ -332,9 +332,10 @@ def test_polygon():
     for i, basis in enumerate(('xy', 'yz', 'xz', 'rz')):
         star_poly = openmc.model.Polygon(star, basis=basis)
         assert points_in[i] in -star_poly
+        assert any([points_in[i] in reg for reg in star_poly.regions])
         assert points_in[i] not in +star_poly
         assert (0, 0, 0) not in -star_poly
         if basis != 'rz':
-            assert (0, 0, 0) in -star_poly.offset(.6)
-
-
+            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])

From 2283e08fdcccbc4b2f49f6fc212f2d424fb419b7 Mon Sep 17 00:00:00 2001
From: Ethan Peterson <ethan.peterson@mit.edu>
Date: Thu, 10 Nov 2022 13:00:11 -0500
Subject: [PATCH 15/15] add abs_tol to prevent cones with infinite or 0 slope

---
 openmc/model/surface_composite.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 5d1553a38..002fad9d9 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -849,7 +849,7 @@ class Polygon(CompositeSurface):
             facet_eq = np.array([dx, dy, c])
             on_boundary = any([np.allclose(facet_eq, eq) for eq in boundary_eqns])
             # Check if the facet is horizontal
-            if isclose(dx, 0):
+            if isclose(dx, 0, abs_tol=1e-8):
                 if basis in ('xz', 'yz', 'rz'):
                     surf = openmc.ZPlane(z0=-c/dy)
                 else:
@@ -857,7 +857,7 @@ class Polygon(CompositeSurface):
                 # if (0, 1).(dx, dy) < 0 we want positive halfspace instead
                 op = operator.pos if dy < 0 else operator.neg
             # Check if the facet is vertical
-            elif isclose(dy, 0):
+            elif isclose(dy, 0, abs_tol=1e-8):
                 if basis in ('xy', 'xz'):
                     surf = openmc.XPlane(x0=-c/dx)
                 elif basis == 'yz':