Handle reflex angles in CylinderSector (#3303)

Author: paulromano

openmc/model/surface_composite.py

@@ -135,6 +135,11 @@ def __init__(self,
         if theta2 <= theta1:
             raise ValueError('theta2 must be greater than theta1.')
 
+        # Determine whether the angle between theta1 and theta2 is a reflex
+        # angle, in which case we need to use a union between the planar
+        # half-spaces
+        self._reflex = (theta2 - theta1 > 180.0)
+
         phi1 = pi / 180 * theta1
         phi2 = pi / 180 * theta2
 
@@ -169,6 +174,9 @@ def __init__(self,
                                                **kwargs)
         self.plane2 = openmc.Plane.from_points(p1, p2_plane2, p3_plane2,
                                                **kwargs)
+        if axis == 'y':
+            self.plane1.flip_normal()
+            self.plane2.flip_normal()
 
     @classmethod
     def from_theta_alpha(cls,
@@ -223,17 +231,11 @@ def from_theta_alpha(cls,
         return cls(r1, r2, theta1, theta2, center=center, axis=axis, **kwargs)
 
     def __neg__(self):
-        if isinstance(self.inner_cyl, openmc.YCylinder):
-            return -self.outer_cyl & +self.inner_cyl & +self.plane1 & -self.plane2
+        if self._reflex:
+            return -self.outer_cyl & +self.inner_cyl & (-self.plane1 | +self.plane2)
         else:
             return -self.outer_cyl & +self.inner_cyl & -self.plane1 & +self.plane2
 
-    def __pos__(self):
-        if isinstance(self.inner_cyl, openmc.YCylinder):
-            return +self.outer_cyl | -self.inner_cyl | -self.plane1 | +self.plane2
-        else:
-            return +self.outer_cyl | -self.inner_cyl | +self.plane1 | -self.plane2
-
 
 class IsogonalOctagon(CompositeSurface):
     r"""Infinite isogonal octagon composite surface

tests/unit_tests/test_surface_composite.py

@@ -207,6 +207,19 @@ def test_cylinder_sector(axis, indices, center):
     assert point_neg[indices] in -s
     assert point_neg[indices] not in +s
 
+    # Check __contains__ for sector with reflex angle
+    s_reflex = openmc.model.CylinderSector(
+        r1, r2, 0., 270., center=center, axis=axis.lower())
+    points = [
+        np.array([c1 + r1 + d, c2 + 0.01, 0.]),
+        np.array([c1, c2 + r1 + d, 0.]),
+        np.array([c1 - r1 - d, c2, 0.]),
+        np.array([c1 - 0.01, c2 - r1 - d, 0.])
+    ]
+    for point_neg in points:
+        assert point_neg[indices] in -s_reflex
+        assert point_neg[indices] not in +s_reflex
+
     # translate method
     t = uniform(-5.0, 5.0)
     s_t = s.translate((t, t, t))