diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 1091d7caef..a6c89be7cd 100644 --- a/docs/source/pythonapi/model.rst +++ b/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 ------------------- diff --git a/openmc/model/__init__.py b/openmc/model/__init__.py index 9fa999dd4e..bd43b14ce8 100644 --- a/openmc/model/__init__.py +++ b/openmc/model/__init__.py @@ -1,3 +1,4 @@ from .triso import * from .model import * from .funcs import * +from .surface_composite import * diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py new file mode 100644 index 0000000000..9f9cd39c60 --- /dev/null +++ b/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__ diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py new file mode 100644 index 0000000000..a1a218311d --- /dev/null +++ b/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)