Add OrthogonalBox composite surface (openmc-dev#3118)

docs/source/pythonapi/model.rst

@@ -26,6 +26,7 @@ Composite Surfaces
    openmc.model.CylinderSector
    openmc.model.HexagonalPrism
    openmc.model.IsogonalOctagon
+   openmc.model.OrthogonalBox
    openmc.model.Polygon
    openmc.model.RectangularParallelepiped
    openmc.model.RectangularPrism

openmc/model/surface_composite.py

@@ -645,6 +645,86 @@ def __pos__(self):
         return +self.xmax | -self.xmin | +self.ymax | -self.ymin | +self.zmax | -self.zmin
 
 
+class OrthogonalBox(CompositeSurface):
+    """Arbitrarily oriented orthogonal box
+
+    This composite surface is composed of four or six planar surfaces that form
+    an arbitrarily oriented orthogonal box when combined.
+
+    Parameters
+    ----------
+    v : iterable of float
+        (x,y,z) coordinates of a corner of the box
+    a1 : iterable of float
+        Vector of first side starting from ``v``
+    a2 : iterable of float
+        Vector of second side starting from ``v``
+    a3 : iterable of float, optional
+        Vector of third side starting from ``v``. When not specified, it is
+        assumed that the box will be infinite along the vector normal to the
+        plane specified by ``a1`` and ``a2``.
+    **kwargs
+        Keyword arguments passed to underlying plane classes
+
+    Attributes
+    ----------
+    ax1_min, ax1_max : openmc.Plane
+        Planes representing minimum and maximum along first axis
+    ax2_min, ax2_max : openmc.Plane
+        Planes representing minimum and maximum along second axis
+    ax3_min, ax3_max : openmc.Plane
+        Planes representing minimum and maximum along third axis
+
+    """
+    _surface_names = ('ax1_min', 'ax1_max', 'ax2_min', 'ax2_max', 'ax3_min', 'ax3_max')
+
+    def __init__(self, v, a1, a2, a3=None, **kwargs):
+        v = np.array(v)
+        a1 = np.array(a1)
+        a2 = np.array(a2)
+        if has_a3 := a3 is not None:
+            a3 = np.array(a3)
+        else:
+            a3 = np.cross(a1, a2)  # normal to plane specified by a1 and a2
+
+        # Generate corners of box
+        p1 = v
+        p2 = v + a1
+        p3 = v + a2
+        p4 = v + a3
+        p5 = v + a1 + a2
+        p6 = v + a2 + a3
+        p7 = v + a1 + a3
+
+        # Generate 6 planes of box
+        self.ax1_min = openmc.Plane.from_points(p1, p3, p4, **kwargs)
+        self.ax1_max = openmc.Plane.from_points(p2, p5, p7, **kwargs)
+        self.ax2_min = openmc.Plane.from_points(p1, p4, p2, **kwargs)
+        self.ax2_max = openmc.Plane.from_points(p3, p6, p5, **kwargs)
+        if has_a3:
+            self.ax3_min = openmc.Plane.from_points(p1, p2, p3, **kwargs)
+            self.ax3_max = openmc.Plane.from_points(p4, p7, p6, **kwargs)
+
+        # Make sure a point inside the box produces the correct senses. If not,
+        # flip the plane coefficients so it does.
+        mid_point = v + (a1 + a2 + a3)/2
+        nums = (1, 2, 3) if has_a3 else (1, 2)
+        for num in nums:
+            min_surf = getattr(self, f'ax{num}_min')
+            max_surf = getattr(self, f'ax{num}_max')
+            if mid_point in -min_surf:
+                min_surf.flip_normal()
+            if mid_point in +max_surf:
+                max_surf.flip_normal()
+
+    def __neg__(self):
+        region = (+self.ax1_min & -self.ax1_max &
+                  +self.ax2_min & -self.ax2_max)
+        if hasattr(self, 'ax3_min'):
+            region &= (+self.ax3_min & -self.ax3_max)
+        return region
+
+
 class XConeOneSided(CompositeSurface):
     """One-sided cone parallel the x-axis
 

openmc/surface.py

@@ -802,6 +802,13 @@ def from_points(cls, p1, p2, p3, **kwargs):
         d = np.dot(n, p1)
         return cls(a=a, b=b, c=c, d=d, **kwargs)
 
+    def flip_normal(self):
+        """Modify plane coefficients to reverse the normal vector."""
+        self.a = -self.a
+        self.b = -self.b
+        self.c = -self.c
+        self.d = -self.d
+
 
 class XPlane(PlaneMixin, Surface):
     """A plane perpendicular to the x axis of the form :math:`x - x_0 = 0`
@@ -1762,7 +1769,7 @@ class Cone(QuadricMixin, Surface):
         z-coordinate of the apex in [cm]. Defaults to 0.
     r2 : float, optional
         Parameter related to the aperture [:math:`\\rm cm^2`].
-        It can be interpreted as the increase in the radius squared per cm along 
+        It can be interpreted as the increase in the radius squared per cm along
         the cone's axis of revolution.
     dx : float, optional
         x-component of the vector representing the axis of the cone.
@@ -1920,7 +1927,7 @@ class XCone(QuadricMixin, Surface):
         z-coordinate of the apex in [cm]. Defaults to 0.
     r2 : float, optional
         Parameter related to the aperture [:math:`\\rm cm^2`].
-        It can be interpreted as the increase in the radius squared per cm along 
+        It can be interpreted as the increase in the radius squared per cm along
         the cone's axis of revolution.
     boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
@@ -2021,7 +2028,7 @@ class YCone(QuadricMixin, Surface):
         z-coordinate of the apex in [cm]. Defaults to 0.
     r2 : float, optional
         Parameter related to the aperture [:math:`\\rm cm^2`].
-        It can be interpreted as the increase in the radius squared per cm along 
+        It can be interpreted as the increase in the radius squared per cm along
         the cone's axis of revolution.
     boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the

tests/unit_tests/test_surface_composite.py

@@ -498,3 +498,57 @@ def test_cruciform_prism(axis):
         openmc.model.CruciformPrism([1.0, 0.5, 2.0, 3.0])
     with pytest.raises(ValueError):
         openmc.model.CruciformPrism([3.0, 2.0, 0.5, 1.0])
+
+
+def test_box():
+    v = (-1.0, -1.0, -2.5)
+    a1 = (2.0, -1.0, 0.0)
+    a2 = (1.0, 2.0, 0.0)
+    a3 = (0.0, 0.0, 5.0)
+    s = openmc.model.OrthogonalBox(v, a1, a2, a3)
+    for num in (1, 2, 3):
+        assert isinstance(getattr(s, f'ax{num}_min'), openmc.Plane)
+        assert isinstance(getattr(s, f'ax{num}_max'), openmc.Plane)
+
+    # Make sure boundary condition propagates
+    s.boundary_type = 'reflective'
+    assert s.boundary_type == 'reflective'
+    for num in (1, 2, 3):
+        assert getattr(s, f'ax{num}_min').boundary_type == 'reflective'
+        assert getattr(s, f'ax{num}_max').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 ll[2] == pytest.approx(-2.5)
+    assert ur[2] == pytest.approx(2.5)
+
+    # __contains__ on associated half-spaces
+    assert (0., 0., 0.) in -s
+    assert (-2., 0., 0.) not in -s
+    assert (0., 0.9, 0.) in -s
+    assert (0., 0., -3.) in +s
+    assert (0., 0., 3.) in +s
+
+    # translate method
+    s_t = s.translate((1., 1., 0.))
+    assert (-0.01, 0., 0.) in +s_t
+    assert (0.01, 0., 0.) in -s_t
+
+    # Make sure repr works
+    repr(s)
+
+    # Version with infinite 3rd dimension
+    s = openmc.model.OrthogonalBox(v, a1, a2)
+    assert not hasattr(s, 'ax3_min')
+    assert not hasattr(s, 'ax3_max')
+    ll, ur = (-s).bounding_box
+    assert np.all(np.isinf(ll))
+    assert np.all(np.isinf(ur))
+    assert (0., 0., 0.) in -s
+    assert (-2., 0., 0.) not in -s
+    assert (0., 0.9, 0.) in -s
+    assert (0., 0., -3.) not in +s
+    assert (0., 0., 3.) not in +s