Merge pull request #1602 from paulromano/composite-surfaces

Implement several composite surfaces

docs/source/pythonapi/model.rst

@@ -17,6 +17,20 @@ Convenience Functions
    openmc.model.subdivide
    openmc.model.pin
 
+Composite Surfaces
+------------------
+
+.. autosummary::
+   :toctree: generated
+   :nosignatures:
+   :template: myclass.rst
+
+   openmc.model.RectangularParallelepiped
+   openmc.model.RightCircularCylinder
+   openmc.model.XConeOneSided
+   openmc.model.YConeOneSided
+   openmc.model.ZConeOneSided
+
 TRISO Fuel Modeling
 -------------------
 

openmc/model/__init__.py

@@ -1,3 +1,4 @@
 from .triso import *
 from .model import *
 from .funcs import *
+from .surface_composite import *

openmc/model/surface_composite.py

@@ -0,0 +1,328 @@
+from abc import ABC, abstractmethod
+from copy import copy
+
+import openmc
+from openmc.checkvalue import check_greater_than, check_value
+
+
+class CompositeSurface(ABC):
+    """Multiple primitive surfaces combined into a composite surface"""
+
+    def translate(self, vector, inplace=False):
+        surf = self if inplace else copy(self)
+        for name in self._surface_names:
+            s = getattr(surf, name)
+            setattr(surf, name, s.translate(vector, inplace))
+        return surf
+
+    def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False):
+        surf = copy(self)
+        for name in self._surface_names:
+            s = getattr(surf, name)
+            setattr(surf, name, s.rotate(rotation, pivot, order, inplace))
+        return surf
+
+    @property
+    def boundary_type(self):
+        return getattr(self, self._surface_names[0]).boundary_type
+
+    @boundary_type.setter
+    def boundary_type(self, boundary_type):
+        # Set boundary type on underlying surfaces, but not for ambiguity plane
+        # on one-sided cones
+        for name in self._surface_names:
+            if name != 'plane':
+                getattr(self, name).boundary_type = boundary_type
+
+    def __repr__(self):
+        return "<{} at 0x{:x}>".format(type(self).__name__, id(self))
+
+    @property
+    @abstractmethod
+    def _surface_names(self):
+        """Iterable of attribute names corresponding to underlying surfaces."""
+
+    @abstractmethod
+    def __pos__(self):
+        """Return the positive half-space of the composite surface."""
+
+    @abstractmethod
+    def __neg__(self):
+        """Return the negative half-space of the composite surface."""
+
+
+class RightCircularCylinder(CompositeSurface):
+    """Right circular cylinder composite surface
+
+    A right circular cylinder is composed of a cylinder and two planar surface
+    perpendicular to the axis of the cylinder. 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
+    right circular cylinder.
+
+    .. versionadded:: 0.12
+
+    Parameters
+    ----------
+    center_base : iterable of float
+        Cartesian coordinate of the center of the base of the cylinder
+    height : float
+        Height of the cylinder
+    radius : float
+        Radius of the cylinder
+    axis : {'x', 'y', 'z'}
+        Axis of the cylinder
+    **kwargs
+        Keyword arguments passed to underlying cylinder and plane classes
+
+    Attributes
+    ----------
+    cyl : openmc.Cylinder
+        Underlying cylinder surface
+    bottom : openmc.Plane
+        Bottom planar surface of the cylinder
+    top : openmc.Plane
+        Top planar surface of the cylinder
+
+    """
+    _surface_names = ('cyl', 'bottom', 'top')
+
+    def __init__(self, center_base, height, radius, axis='z', **kwargs):
+        cx, cy, cz = center_base
+        check_greater_than('cylinder height', height, 0.0)
+        check_greater_than('cylinder radius', radius, 0.0)
+        check_value('cylinder axis', axis, ('x', 'y', 'z'))
+        if axis == 'x':
+            self.cyl = openmc.XCylinder(y0=cy, z0=cz, r=radius, **kwargs)
+            self.bottom = openmc.XPlane(x0=cx, **kwargs)
+            self.top = openmc.XPlane(x0=cx + height, **kwargs)
+        elif axis == 'y':
+            self.cyl = openmc.YCylinder(x0=cx, z0=cz, r=radius, **kwargs)
+            self.bottom = openmc.YPlane(y0=cy, **kwargs)
+            self.top = openmc.YPlane(y0=cy + height, **kwargs)
+        elif axis == 'z':
+            self.cyl = openmc.ZCylinder(x0=cx, y0=cy, r=radius, **kwargs)
+            self.bottom = openmc.ZPlane(z0=cz, **kwargs)
+            self.top = openmc.ZPlane(z0=cz + height, **kwargs)
+
+    def __neg__(self):
+        return -self.cyl & +self.bottom & -self.top
+
+    def __pos__(self):
+        return +self.cyl | -self.bottom | +self.top
+
+
+class RectangularParallelepiped(CompositeSurface):
+    """Rectangular parallelpiped composite surface
+
+    A rectangular parallelpiped is composed of six planar surfaces. 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 rectangular parallelpiped.
+
+    .. versionadded:: 0.12
+
+    Parameters
+    ----------
+    xmin, xmax : float
+        Minimum and maximum x coordinates of the parallelepiped
+    ymin, ymax : float
+        Minimum and maximum y coordinates of the parallelepiped
+    zmin, zmax : float
+        Minimum and maximum z coordinates of the parallelepiped
+    **kwargs
+        Keyword arguments passed to underlying plane classes
+
+    Attributes
+    ----------
+    xmin, xmax : openmc.XPlane
+        Sides of the parallelepiped
+    ymin, ymax : openmc.YPlane
+        Sides of the parallelepiped
+    zmin, zmax : openmc.ZPlane
+        Sides of the parallelepiped
+
+    """
+    _surface_names = ('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax')
+
+    def __init__(self, xmin, xmax, ymin, ymax, zmin, zmax, **kwargs):
+        if xmin >= xmax:
+            raise ValueError('xmin must be less than xmax')
+        if ymin >= ymax:
+            raise ValueError('ymin must be less than ymax')
+        if zmin >= zmax:
+            raise ValueError('zmin must be less than zmax')
+        self.xmin = openmc.XPlane(x0=xmin, **kwargs)
+        self.xmax = openmc.XPlane(x0=xmax, **kwargs)
+        self.ymin = openmc.YPlane(y0=ymin, **kwargs)
+        self.ymax = openmc.YPlane(y0=ymax, **kwargs)
+        self.zmin = openmc.ZPlane(z0=zmin, **kwargs)
+        self.zmax = openmc.ZPlane(z0=zmax, **kwargs)
+
+    def __neg__(self):
+        return +self.xmin & -self.xmax & +self.ymin & -self.ymax & +self.zmin & -self.zmax
+
+    def __pos__(self):
+        return -self.xmin | +self.ymax | -self.ymin | +self.ymax | -self.zmin | +self.zmax
+
+
+class XConeOneSided(CompositeSurface):
+    """One-sided cone parallel the x-axis
+
+    A one-sided cone is composed of a normal cone surface and an "ambiguity"
+    surface that eliminates the ambiguity as to which region of space is
+    included. 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 cone.
+
+    .. versionadded:: 0.12
+
+    Parameters
+    ----------
+    x0 : float, optional
+        x-coordinate of the apex. Defaults to 0.
+    y0 : float, optional
+        y-coordinate of the apex. Defaults to 0.
+    z0 : float, optional
+        z-coordinate of the apex. Defaults to 0.
+    r2 : float, optional
+        Parameter related to the aperature. Defaults to 1.
+    up : bool
+        Whether to select the side of the cone that extends to infinity in the
+        positive direction of the coordinate axis (the positive half-space of
+        the ambiguity plane)
+    **kwargs
+        Keyword arguments passed to underlying plane classes
+
+    Attributes
+    ----------
+    cone : openmc.XCone
+        Regular two-sided cone
+    plane : openmc.XPlane
+        Ambiguity surface
+    up : bool
+        Whether to select the side of the cone that extends to infinity in the
+        positive direction of the coordinate axis (the positive half-space of
+        the ambiguity plane)
+
+    """
+    _surface_names = ('cone', 'plane')
+
+    def __init__(self, x0=0., y0=0., z0=0., r2=1., up=True, **kwargs):
+        check_greater_than('cone R^2', r2, 0.0)
+        self.cone = openmc.XCone(x0, y0, z0, r2, **kwargs)
+        self.plane = openmc.XPlane(x0)
+        self.up = up
+
+    def __neg__(self):
+        return -self.cone & (+self.plane if self.up else -self.plane)
+
+    def __pos__(self):
+        if self.up:
+            return (+self.cone & +self.plane) | -self.plane
+        else:
+            return (+self.cone & -self.plane) | +self.plane
+
+
+class YConeOneSided(CompositeSurface):
+    """One-sided cone parallel the y-axis
+
+    A one-sided cone is composed of a normal cone surface and an "ambiguity"
+    surface that eliminates the ambiguity as to which region of space is
+    included. 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 cone.
+
+    .. versionadded:: 0.12
+
+    Parameters
+    ----------
+    x0 : float, optional
+        x-coordinate of the apex. Defaults to 0.
+    y0 : float, optional
+        y-coordinate of the apex. Defaults to 0.
+    z0 : float, optional
+        z-coordinate of the apex. Defaults to 0.
+    r2 : float, optional
+        Parameter related to the aperature. Defaults to 1.
+    up : bool
+        Whether to select the side of the cone that extends to infinity in the
+        positive direction of the coordinate axis (the positive half-space of
+        the ambiguity plane)
+    **kwargs
+        Keyword arguments passed to underlying plane classes
+
+    Attributes
+    ----------
+    cone : openmc.YCone
+        Regular two-sided cone
+    plane : openmc.YPlane
+        Ambiguity surface
+    up : bool
+        Whether to select the side of the cone that extends to infinity in the
+        positive direction of the coordinate axis (the positive half-space of
+        the ambiguity plane)
+
+    """
+    _surface_names = ('cone', 'plane')
+
+    def __init__(self, x0=0., y0=0., z0=0., r2=1., up=True, **kwargs):
+        check_greater_than('cone R^2', r2, 0.0)
+        self.cone = openmc.YCone(x0, y0, z0, r2, **kwargs)
+        self.plane = openmc.YPlane(y0)
+        self.up = up
+
+    __neg__ = XConeOneSided.__neg__
+    __pos__ = XConeOneSided.__pos__
+
+
+class ZConeOneSided(CompositeSurface):
+    """One-sided cone parallel the z-axis
+
+    A one-sided cone is composed of a normal cone surface and an "ambiguity"
+    surface that eliminates the ambiguity as to which region of space is
+    included. 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 cone.
+
+    .. versionadded:: 0.12
+
+    Parameters
+    ----------
+    x0 : float, optional
+        x-coordinate of the apex. Defaults to 0.
+    y0 : float, optional
+        y-coordinate of the apex. Defaults to 0.
+    z0 : float, optional
+        z-coordinate of the apex. Defaults to 0.
+    r2 : float, optional
+        Parameter related to the aperature. Defaults to 1.
+    up : bool
+        Whether to select the side of the cone that extends to infinity in the
+        positive direction of the coordinate axis (the positive half-space of
+        the ambiguity plane)
+    **kwargs
+        Keyword arguments passed to underlying plane classes
+
+    Attributes
+    ----------
+    cone : openmc.ZCone
+        Regular two-sided cone
+    plane : openmc.ZPlane
+        Ambiguity surface
+    up : bool
+        Whether to select the side of the cone that extends to infinity in the
+        positive direction of the coordinate axis (the positive half-space of
+        the ambiguity plane)
+
+    """
+    _surface_names = ('cone', 'plane')
+
+    def __init__(self, x0=0., y0=0., z0=0., r2=1., up=True, **kwargs):
+        check_greater_than('cone R^2', r2, 0.0)
+        self.cone = openmc.ZCone(x0, y0, z0, r2, **kwargs)
+        self.plane = openmc.ZPlane(z0)
+        self.up = up
+
+    __neg__ = XConeOneSided.__neg__
+    __pos__ = XConeOneSided.__pos__

tests/unit_tests/test_surface_composite.py

@@ -0,0 +1,152 @@
+from random import uniform
+
+import numpy as np
+import openmc
+import pytest
+
+
+def test_rectangular_parallelepiped():
+    xmin = uniform(-5., 5.)
+    xmax = xmin + uniform(0., 5.)
+    ymin = uniform(-5., 5.)
+    ymax = ymin + uniform(0., 5.)
+    zmin = uniform(-5., 5.)
+    zmax = zmin + uniform(0., 5.)
+    s = openmc.model.RectangularParallelepiped(xmin, xmax, ymin, ymax, zmin, zmax)
+    assert isinstance(s.xmin, openmc.XPlane)
+    assert isinstance(s.xmax, openmc.XPlane)
+    assert isinstance(s.ymin, openmc.YPlane)
+    assert isinstance(s.ymax, openmc.YPlane)
+    assert isinstance(s.zmin, openmc.ZPlane)
+    assert isinstance(s.zmax, openmc.ZPlane)
+
+    # Make sure boundary condition propagates
+    s.boundary_type = 'reflective'
+    assert s.boundary_type == 'reflective'
+    for axis in 'xyz':
+        assert getattr(s, '{}min'.format(axis)).boundary_type == 'reflective'
+        assert getattr(s, '{}max'.format(axis)).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 == pytest.approx((xmax, ymax, zmax))
+    assert ll == pytest.approx((xmin, ymin, zmin))
+
+    # __contains__ on associated half-spaces
+    assert (xmin - 0.1, 0., 0.) in +s
+    assert (xmin - 0.1, 0., 0.) not in -s
+    dx, dy, dz = xmax - xmin, ymax - ymin, zmax - zmin
+    assert (xmin + dx/2, ymin + dy/2, zmin + dz/2) in -s
+    assert (xmin + dx/2, ymin + dy/2, zmin + dz/2) 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)
+
+
+@pytest.mark.parametrize(
+    "axis, indices", [
+        ("X", [0, 1, 2]),
+        ("Y", [1, 2, 0]),
+        ("Z", [2, 0, 1]),
+    ]
+)
+def test_right_circular_cylinder(axis, indices):
+    x, y, z = 1.0, -2.5, 3.0
+    h, r = 5.0, 3.0
+    s = openmc.model.RightCircularCylinder((x, y, z), h, r, axis=axis.lower())
+    assert isinstance(s.cyl, getattr(openmc, axis + "Cylinder"))
+    assert isinstance(s.top, getattr(openmc, axis + "Plane"))
+    assert isinstance(s.bottom, getattr(openmc, axis + "Plane"))
+
+    # Make sure boundary condition propagates
+    s.boundary_type = 'reflective'
+    assert s.boundary_type == 'reflective'
+    assert s.cyl.boundary_type == 'reflective'
+    assert s.bottom.boundary_type == 'reflective'
+    assert s.top.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 == pytest.approx((x, y, z) + np.roll([0, -r, -r], indices[0]))
+    assert ur == pytest.approx((x, y, z) + np.roll([h, r, r], indices[0]))
+
+    # __contains__ on associated half-spaces
+    point_pos = (x, y, z) + np.roll([h/2, r+1, r+1], indices[0])
+    assert point_pos in +s
+    assert point_pos not in -s
+    point_neg = (x, y, z) + np.roll([h/2, 0, 0], indices[0])
+    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)
+
+
+@pytest.mark.parametrize(
+    "axis, point_pos, point_neg, ll_true", [
+        ("X", (8., 0., 0.), (12., 0., 0.), (10., -np.inf, -np.inf)),
+        ("Y", (10., -2., 0.), (10., 2., 0.), (-np.inf, 0., -np.inf)),
+        ("Z", (10., 0., -3.), (10., 0., 3.), (-np.inf, -np.inf, 0.))
+    ]
+)
+def test_cone_one_sided(axis, point_pos, point_neg, ll_true):
+    cone_oneside = getattr(openmc.model, axis + "ConeOneSided")
+    cone_twoside = getattr(openmc, axis + "Cone")
+    plane = getattr(openmc, axis + "Plane")
+
+    x, y, z = 10., 0., 0.
+    r2 = 4.
+    s = cone_oneside(x, y, z, r2, True)
+    assert isinstance(s.cone, cone_twoside)
+    assert isinstance(s.plane, plane)
+    assert s.up
+
+    # Make sure boundary condition propagates
+    s.boundary_type = 'reflective'
+    assert s.boundary_type == 'reflective'
+    assert s.cone.boundary_type == 'reflective'
+    assert s.plane.boundary_type == 'transmission'
+
+    # 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 np.all(np.isinf(ur))
+    assert ll == pytest.approx(ll_true)
+
+    # __contains__ on associated half-spaces
+    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)