Merge pull request #2023 from yardasol/octagon

`IsogonalOctagon` composite surface

docs/source/pythonapi/model.rst

@@ -25,6 +25,7 @@ Composite Surfaces
    :nosignatures:
    :template: myclass.rst
 
+   openmc.model.IsogonalOctagon
    openmc.model.RectangularParallelepiped
    openmc.model.RightCircularCylinder
    openmc.model.XConeOneSided

openmc/model/surface_composite.py

@@ -1,10 +1,11 @@
 from abc import ABC, abstractmethod
 from copy import copy
+from math import sqrt
+import numpy as np
 
 import openmc
 from openmc.checkvalue import check_greater_than, check_value
 
-
 class CompositeSurface(ABC):
     """Multiple primitive surfaces combined into a composite surface"""
 
@@ -51,6 +52,137 @@ class CompositeSurface(ABC):
         """Return the negative half-space of the composite surface."""
 
 
+class IsogonalOctagon(CompositeSurface):
+    """Infinite isogonal octagon composite surface
+
+    An isogonal octagon is composed of eight planar surfaces. The prism is
+    parallel to the x, y, or z axis. The remaining two axes (y and z, z and x,
+    or x and y) serve as a basis for constructing the surfaces. Two surfaces
+    are parallel to the first basis axis, two surfaces are parallel
+    to the second basis axis, and the remaining four surfaces intersect both
+    basis axes at 45 degree angles.
+
+    This class acts as a proper surface, meaning that unary `+` and `-`
+    operators applied to it will produce a half-space. The negative side is
+    defined to be the region inside of the octogonal prism.
+
+    Parameters
+    ----------
+    center : iterable of float
+        Coordinate for the central axis of the octagon in the
+        (y, z), (z, x), or (x, y) basis.
+    r1 : float
+        Half-width of octagon across its basis axis-parallel sides in units
+        of cm. Must be less than :math:`r_2\sqrt{2}`.
+    r2 : float
+        Half-width of octagon across its basis axis intersecting sides in
+        units of cm. Must be less than than :math:`r_1\sqrt{2}`.
+    axis : {'x', 'y', 'z'}
+        Central axis of octagon. Defaults to 'z'
+    **kwargs
+        Keyword arguments passed to underlying plane classes
+
+    Attributes
+    ----------
+    top : openmc.ZPlane, openmc.XPlane, or openmc.YPlane
+        Top planar surface of octagon
+    bottom : openmc.ZPlane, openmc.XPlane, or openmc.YPlane
+        Bottom planar surface of octagon
+    right : openmc.YPlane, openmc.ZPlane, or openmc.XPlane
+        Right planar surface of octagon
+    left : openmc.YPlane, openmc.ZPlane, or openmc.XPlane
+        Left planar surface of octagon
+    upper_right : openmc.Plane
+        Upper right planar surface of octagon
+    lower_right : openmc.Plane
+        Lower right planar surface of octagon
+    lower_left : openmc.Plane
+        Lower left planar surface of octagon
+    upper_left : openmc.Plane
+        Upper left planar surface of octagon
+
+    """
+
+    _surface_names = ('top', 'bottom',
+                      'upper_right', 'lower_left',
+                      'right', 'left',
+                      'lower_right', 'upper_left')
+
+    def __init__(self, center, r1, r2, axis='z', **kwargs):
+        c1, c2 = center
+
+        # Coords for axis-perpendicular planes
+        ctop = c1 + r1
+        cbottom = c1 - r1
+
+        cright = c2 + r1
+        cleft = c2 - r1
+
+        # Side lengths
+        if r2 > r1 * sqrt(2):
+            raise ValueError(f'r2 is greater than sqrt(2) * r1. Octagon' + \
+                             ' may be erroneous.')
+        if r1 > r2 * sqrt(2):
+            raise ValueError(f'r1 is greater than sqrt(2) * r2. Octagon' + \
+                             ' may be erroneous.')
+
+        L_basis_ax = (r2 * sqrt(2) - r1)
+
+        # Coords for quadrant planes
+        p1_ur = np.array([L_basis_ax, r1, 0.])
+        p2_ur = np.array([r1, L_basis_ax, 0.])
+        p3_ur = np.array([r1, L_basis_ax, 1.])
+
+        p1_lr = np.array([r1, -L_basis_ax, 0.])
+        p2_lr = np.array([L_basis_ax, -r1, 0.])
+        p3_lr = np.array([L_basis_ax, -r1, 1.])
+
+        points = [p1_ur, p2_ur, p3_ur, p1_lr, p2_lr, p3_lr]
+
+        # Orientation specific variables
+        if axis == 'z':
+            coord_map = [0, 1, 2]
+            self.top = openmc.YPlane(ctop, **kwargs)
+            self.bottom = openmc.YPlane(cbottom, **kwargs)
+            self.right = openmc.XPlane(cright, **kwargs)
+            self.left = openmc.XPlane(cleft, **kwargs)
+        elif axis == 'y':
+            coord_map = [1, 2, 0]
+            self.top = openmc.XPlane(ctop, **kwargs)
+            self.bottom = openmc.XPlane(cbottom, **kwargs)
+            self.right = openmc.ZPlane(cright, **kwargs)
+            self.left = openmc.ZPlane(cleft, **kwargs)
+        elif axis == 'x':
+            coord_map = [2, 0, 1]
+            self.top = openmc.ZPlane(ctop, **kwargs)
+            self.bottom = openmc.ZPlane(cbottom, **kwargs)
+            self.right = openmc.YPlane(cright, **kwargs)
+            self.left = openmc.YPlane(cleft, **kwargs)
+
+        # Put our coordinates in (x,y,z) order
+        for p in points:
+            p[:] = p[coord_map]
+
+        self.upper_right = openmc.Plane.from_points(p1_ur, p2_ur, p3_ur,
+                                                    **kwargs)
+        self.lower_right = openmc.Plane.from_points(p1_lr, p2_lr, p3_lr,
+                                                    **kwargs)
+        self.lower_left = openmc.Plane.from_points(-p1_ur, -p2_ur, -p3_ur,
+                                                   **kwargs)
+        self.upper_left = openmc.Plane.from_points(-p1_lr, -p2_lr, -p3_lr,
+                                                   **kwargs)
+
+    def __neg__(self):
+        return -self.top & +self.bottom & -self.right & +self.left & \
+            +self.upper_right & +self.lower_right & -self.lower_left & \
+            -self.upper_left
+
+    def __pos__(self):
+        return +self.top | -self.bottom | +self.right | -self.left | \
+            -self.upper_right | -self.lower_right | +self.lower_left | \
+            +self.upper_left
+
+
 class RightCircularCylinder(CompositeSurface):
     """Right circular cylinder composite surface
 

tests/unit_tests/test_surface_composite.py

@@ -150,3 +150,80 @@ def test_cone_one_sided(axis, point_pos, point_neg, ll_true):
 
     # Make sure repr works
     repr(s)
+
+
+@pytest.mark.parametrize(
+    "axis, plane_tb, plane_lr, axis_idx", [
+        ("x", "Z", "Y", 0),
+        ("y", "X", "Z", 1),
+        ("z", "Y", "X", 2),
+    ]
+)
+def test_isogonal_octagon(axis, plane_tb, plane_lr, axis_idx):
+    center = np.array([0., 0.])
+    point_pos = np.array([0.8, 0.8])
+    point_neg = np.array([0.7, 0.7])
+    r1 = 1.
+    r2 = 1.
+    plane_top_bottom = getattr(openmc, plane_tb + "Plane")
+    plane_left_right = getattr(openmc, plane_lr + "Plane")
+    s = openmc.model.IsogonalOctagon(center, r1, r2, axis=axis)
+    assert isinstance(s.top, plane_top_bottom)
+    assert isinstance(s.bottom, plane_top_bottom)
+    assert isinstance(s.right, plane_left_right)
+    assert isinstance(s.left, plane_left_right)
+    assert isinstance(s.upper_right, openmc.Plane)
+    assert isinstance(s.lower_right, openmc.Plane)
+    assert isinstance(s.upper_left, openmc.Plane)
+    assert isinstance(s.lower_left, openmc.Plane)
+
+    # Make sure boundary condition propagates
+    s.boundary_type = 'reflective'
+    assert s.boundary_type == 'reflective'
+    assert s.top.boundary_type == 'reflective'
+    assert s.bottom.boundary_type == 'reflective'
+    assert s.right.boundary_type == 'reflective'
+    assert s.left.boundary_type == 'reflective'
+    assert s.upper_right.boundary_type == 'reflective'
+    assert s.lower_right.boundary_type == 'reflective'
+    assert s.lower_left.boundary_type == 'reflective'
+    assert s.upper_left.boundary_type == 'reflective'
+
+    # Check bounding box
+    center = np.insert(center, axis_idx, np.inf)
+    xmax, ymax, zmax = center + r1
+    coord_min = center - r1
+    coord_min[axis_idx] *= -1
+    xmin, ymin, zmin = coord_min
+    ll, ur = (+s).bounding_box
+    assert np.all(np.isinf(ll))
+    assert np.all(np.isinf(ur))
+    ll, ur = (-s).bounding_box
+    assert ur == pytest.approx((xmax, ymax, zmax))
+    assert ll == pytest.approx((xmin, ymin, zmin))
+
+    # __contains__ on associated half-spaces
+    point_pos = np.insert(point_pos, axis_idx, 0)
+    point_neg = np.insert(point_neg, axis_idx, 0)
+    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)
+
+    # Check invalid r1, r2 combinations
+    with pytest.raises(ValueError):
+        openmc.model.IsogonalOctagon(center, r1=1.0, r2=10.)
+    with pytest.raises(ValueError):
+        openmc.model.IsogonalOctagon(center, r1=10., r2=1.)
+
+    # Make sure repr works
+    repr(s)
+
+