Implement CruciformPrism composite surface

docs/source/pythonapi/model.rst

@@ -24,14 +24,15 @@ Composite Surfaces
    :nosignatures:
    :template: myclass.rst
 
+   openmc.model.CruciformPrism
    openmc.model.CylinderSector
    openmc.model.IsogonalOctagon
+   openmc.model.Polygon
    openmc.model.RectangularParallelepiped
    openmc.model.RightCircularCylinder
    openmc.model.XConeOneSided
    openmc.model.YConeOneSided
    openmc.model.ZConeOneSided
-   openmc.model.Polygon
 
 TRISO Fuel Modeling
 -------------------

openmc/model/surface_composite.py

@@ -630,7 +630,7 @@ class ZConeOneSided(CompositeSurface):
 
 
 class Polygon(CompositeSurface):
-    """Create a polygon composite surface from a path of closed points.
+    """Polygon formed from a path of closed points.
 
     .. versionadded:: 0.13.3
 
@@ -1079,3 +1079,75 @@ class Polygon(CompositeSurface):
         disp_vec = distance / costheta * unit_nvec
 
         return type(self)(self.points + disp_vec, basis=self.basis)
+
+
+class CruciformPrism(CompositeSurface):
+    """Generalized cruciform prism
+
+    This surface represents a prism parallel to an axis formed by planes at
+    multiple distances from the center. Equivalent to the 'gcross' derived
+    surface in Serpent.
+
+    .. versionadded:: 0.13.4
+
+    Parameters
+    ----------
+    distances : iterable of float
+        A monotonically increasing (or decreasing) iterable of distances in [cm]
+        that form the planes of the generalized cruciform.
+    center : iterable of float
+        The center of the prism in the two non-parallel axes (e.g., (x, y) when
+        axis is 'z') in [cm]
+    axis : {'x', 'y', 'z'}
+        Axis to which the prism is parallel
+    **kwargs
+        Keyword arguments passed to underlying plane classes
+
+    """
+
+    def __init__(self, distances, center=(0., 0.), axis='z', **kwargs):
+        x0, y0 = center
+        self.distances = distances
+
+        if axis == 'x':
+            cls_horizontal = openmc.YPlane
+            cls_vertical = openmc.ZPlane
+        elif axis == 'y':
+            cls_horizontal = openmc.XPlane
+            cls_vertical = openmc.ZPlane
+        elif axis == 'z':
+            cls_horizontal = openmc.XPlane
+            cls_vertical = openmc.YPlane
+        else:
+            raise ValueError("axis must be 'x', 'y', or 'z'")
+
+        # Create each planar surface
+        surfnames = []
+        for i, d in enumerate(distances):
+            setattr(self, f'hmin{i}', cls_horizontal(x0 - d, **kwargs))
+            setattr(self, f'hmax{i}', cls_horizontal(x0 + d, **kwargs))
+            setattr(self, f'vmin{i}', cls_vertical(y0 - d, **kwargs))
+            setattr(self, f'vmax{i}', cls_vertical(y0 + d, **kwargs))
+            surfnames.extend([f'hmin{i}', f'hmax{i}', f'vmin{i}', f'vmax{i}'])
+
+        # Set _surfnames to satisfy CompositeSurface protocol
+        self._surfnames = tuple(surfnames)
+
+    @property
+    def _surface_names(self):
+        return self._surfnames
+
+    def __neg__(self):
+        n = len(self.distances)
+        regions = []
+        for i in range(n):
+            regions.append(
+                +getattr(self, f'hmin{i}') &
+                -getattr(self, f'hmax{i}') &
+                +getattr(self, f'vmin{n-1-i}') &
+                -getattr(self, f'vmax{n-1-i}')
+            )
+        return openmc.Union(regions)
+
+    def __pos__(self):
+        return ~(-self)

tests/unit_tests/test_surface_composite.py

@@ -315,6 +315,7 @@ def test_isogonal_octagon(axis, plane_tb, plane_lr, axis_idx):
     # Make sure repr works
     repr(s)
 
+
 def test_polygon():
     # define a 5 pointed star centered on 1, 1
     star = np.array([[1.        , 2.        ],
@@ -393,3 +394,64 @@ def test_polygon():
                           [6.88, 3.02]])
     with pytest.raises(ValueError):
         openmc.model.Polygon(rz_points)
+
+
+@pytest.mark.parametrize("axis", ["x", "y", "z"])
+def test_cruciform_prism(axis):
+    center = x0, y0 = (3., 4.)
+    distances = [2., 3., 5.]
+    s = openmc.model.CruciformPrism(distances, center, axis=axis)
+
+    if axis == 'x':
+        i1, i2 = 1, 2
+    elif axis == 'y':
+        i1, i2 = 0, 2
+    elif axis == 'z':
+        i1, i2 = 0, 1
+    plane_cls = (openmc.XPlane, openmc.YPlane, openmc.ZPlane)
+
+    # Check type of surfaces
+    for i in range(3):
+        assert isinstance(getattr(s, f'hmin{i}'), plane_cls[i1])
+        assert isinstance(getattr(s, f'hmax{i}'), plane_cls[i1])
+        assert isinstance(getattr(s, f'vmin{i}'), plane_cls[i2])
+        assert isinstance(getattr(s, f'vmax{i}'), plane_cls[i2])
+
+    # Make sure boundary condition propagates
+    s.boundary_type = 'reflective'
+    for i in range(3):
+        assert getattr(s, f'hmin{i}').boundary_type == 'reflective'
+        assert getattr(s, f'hmax{i}').boundary_type == 'reflective'
+        assert getattr(s, f'vmin{i}').boundary_type == 'reflective'
+        assert getattr(s, f'vmax{i}').boundary_type == 'reflective'
+
+    # Check bounding box
+    ll, ur = (+s).bounding_box
+    assert np.all(np.isinf(ll))
+    assert np.all(np.isinf(ur))
+    ll, ur = (-s).bounding_box
+    assert ur[i1] == pytest.approx(x0 + distances[-1])
+    assert ur[i2] == pytest.approx(y0 + distances[-1])
+    assert ll[i1] == pytest.approx(x0 - distances[-1])
+    assert ll[i2] == pytest.approx(y0 - distances[-1])
+
+    # __contains__ on associated half-spaces
+    point_pos, point_neg = np.zeros(3), np.zeros(3)
+    point_pos[i1] = x0 + 3.1
+    point_pos[i2] = y0 + 2.05
+    point_neg[i1] = x0 + 3.5
+    point_neg[i2] = y0 + 1.99
+    assert point_pos in +s
+    assert point_pos not in -s
+    assert point_neg in -s
+    assert point_neg not in +s
+
+    # translate method
+    t = uniform(-5.0, 5.0)
+    s_t = s.translate((t, t, t))
+    ll_t, ur_t = (-s_t).bounding_box
+    assert ur_t == pytest.approx(ur + t)
+    assert ll_t == pytest.approx(ll + t)
+
+    # Make sure repr works
+    repr(s)