diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 3a17b6169..5962897de 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -135,6 +135,11 @@ class CylinderSector(CompositeSurface): 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 @@ class CylinderSector(CompositeSurface): **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 @@ class CylinderSector(CompositeSurface): 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 diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py index da7ffbcc4..a04fc73c5 100644 --- a/tests/unit_tests/test_surface_composite.py +++ b/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))