From 37ab9c75605746221cebd76b82e4e6a8a174fc71 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 28 Feb 2020 11:38:20 -0600 Subject: [PATCH 1/8] Add composite surfaces --- openmc/model/__init__.py | 1 + openmc/model/surface_composite.py | 164 ++++++++++++++++++++++++++++++ 2 files changed, 165 insertions(+) create mode 100644 openmc/model/surface_composite.py 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..66afd4b1ce --- /dev/null +++ b/openmc/model/surface_composite.py @@ -0,0 +1,164 @@ +from copy import copy + +import openmc + + +class CompositeMixin: + def evaluate(self, point): + raise NotImplementedError('Composite surfaces do not have a surface equation.') + + def _get_base_coeffs(self): + raise NotImplementedError('Composite surfaces do not have base coefficients.') + + 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): + surf = copy(self) + for name in self._surface_names: + s = getattr(surf, name) + setattr(surf, name, s.rotate(rotation)) + return surf + + def __repr__(self): + return "<{} at 0x{:x}>".format(type(self).__name__, id(self)) + + +class RightCircularCylinder(CompositeMixin, openmc.Surface): + _surface_names = ('cyl', 'bottom', 'top') + + def __init__(self, center_base, height, radius, axis='z', + boundary_type='transmission'): + kwargs = {'boundary_type': boundary_type} + cx, cy, cz = center_base + 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(CompositeMixin, openmc.Surface): + _surface_names = ('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax') + + def __init__(self, xmin, xmax, ymin, ymax, zmin, zmax, boundary_type='transmission'): + kwargs = {'boundary_type': boundary_type} + 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 Box(CompositeMixin, openmc.Surface): + _surface_names = ('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax') + + def __init__(self, v, a1, a2, a3, boundary_type='transmission'): + kwargs = {'boundary_type': boundary_type} + vx, vy, vz = v + a1x, a1y, a1z = a1 + a2x, a2y, a2z = a2 + a3x, a3y, a3z = a3 + + # Only support boxes with axis-aligned vectors + if any(x != 0.0 for x in (a1y, a1z, a2x, a2z, a3x, a3y)): + raise NotImplementedError('Box macrobody with non-axis-aligned ' + 'vector not supported.') + + # Determine each side of the box + if a1x > 0: + xmin, xmax = vx, vx + a1x + else: + xmin, xmax = vx + a1x, vx + if a2y > 0: + ymin, ymax = vy, vy + a2y + else: + ymin, ymax = vy + a2y, vy + if a3z > 0: + zmin, zmax = vz, vz + a3z + else: + zmin, zmax = vz + a3z, vz + + # Create surfaces + 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.xmax | + -self.ymin | +self.ymax | + -self.zmin | +self.zmax) + + +class XConeOneSided(CompositeMixin, openmc.Surface): + _surface_names = ('cone', 'plane') + + def __init__(self, x0=0., y0=0., z0=0., r2=1., up=True, **kwargs): + 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(CompositeMixin, openmc.Surface): + _surface_names = ('cone', 'plane') + + def __init__(self, x0=0., y0=0., z0=0., r2=1., up=True, **kwargs): + 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(CompositeMixin, openmc.Surface): + _surface_names = ('cone', 'plane') + + def __init__(self, x0=0., y0=0., z0=0., r2=1., up=True, **kwargs): + self.cone = openmc.ZCone(x0, y0, z0, r2, **kwargs) + self.plane = openmc.ZPlane(z0) + self.up = up + + __neg__ = XConeOneSided.__neg__ + __pos__ = XConeOneSided.__pos__ From df7b52530621d65015e947d1ce032db12b1476ec Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 2 Jul 2020 07:34:36 -0500 Subject: [PATCH 2/8] Add docstrings to composite surface classes --- docs/source/pythonapi/model.rst | 15 ++ openmc/model/surface_composite.py | 224 ++++++++++++++++++++++++++++-- 2 files changed, 229 insertions(+), 10 deletions(-) diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 1091d7caef..97704e76da 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -17,6 +17,21 @@ Convenience Functions openmc.model.subdivide openmc.model.pin +Composite Surfaces +------------------ + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + openmc.model.Box + openmc.model.RectangularParallelepiped + openmc.model.RightCircularCylinder + openmc.model.XConeOneSided + openmc.model.YConeOneSided + openmc.model.ZConeOneSided + TRISO Fuel Modeling ------------------- diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 66afd4b1ce..e730390e9b 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -4,6 +4,8 @@ import openmc class CompositeMixin: + """Mixin class for composite surfaces""" + def evaluate(self, point): raise NotImplementedError('Composite surfaces do not have a surface equation.') @@ -17,7 +19,7 @@ class CompositeMixin: setattr(surf, name, s.translate(vector, inplace)) return surf - def rotate(self, rotation): + 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) @@ -29,11 +31,42 @@ class CompositeMixin: class RightCircularCylinder(CompositeMixin, openmc.Surface): - _surface_names = ('cyl', 'bottom', 'top') + """Right circular cylinder composite surface - def __init__(self, center_base, height, radius, axis='z', - boundary_type='transmission'): - kwargs = {'boundary_type': boundary_type} + 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 if axis == 'x': self.cyl = openmc.XCylinder(y0=cy, z0=cz, r=radius, **kwargs) @@ -56,10 +89,39 @@ class RightCircularCylinder(CompositeMixin, openmc.Surface): class RectangularParallelepiped(CompositeMixin, openmc.Surface): - _surface_names = ('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax') + """Rectangular parallelpiped composite surface - def __init__(self, xmin, xmax, ymin, ymax, zmin, zmax, boundary_type='transmission'): - kwargs = {'boundary_type': boundary_type} + 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): self.xmin = openmc.XPlane(x0=xmin, **kwargs) self.xmax = openmc.XPlane(x0=xmax, **kwargs) self.ymin = openmc.YPlane(y0=ymin, **kwargs) @@ -75,10 +137,41 @@ class RectangularParallelepiped(CompositeMixin, openmc.Surface): class Box(CompositeMixin, openmc.Surface): + """Arbitrarily oriented orthogonal box + + An arbitrarily oriented orthogonal box 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 box. + + .. versionadded:: 0.12 + + Parameters + ---------- + v : iterable of float + Cartesian coordinates of a corner of the box + a1 : iterable of float + Vector of first side from the specified corner coordinates + a2 : iterable of float + Vector of second side from the specified corner coordinates + a3 : iterable of float + Vector of third side from the specified corner coordinates + **kwargs + Keyword arguments passed to underlying plane classes + + Attributes + ---------- + xmin, xmax : openmc.XPlane + Sides of the box + ymin, ymax : openmc.YPlane + Sides of the box + zmin, zmax : openmc.ZPlane + Sides of the box + + """ _surface_names = ('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax') - def __init__(self, v, a1, a2, a3, boundary_type='transmission'): - kwargs = {'boundary_type': boundary_type} + def __init__(self, v, a1, a2, a3, **kwargs): vx, vy, vz = v a1x, a1y, a1z = a1 a2x, a2y, a2z = a2 @@ -123,6 +216,43 @@ class Box(CompositeMixin, openmc.Surface): class XConeOneSided(CompositeMixin, openmc.Surface): + """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 + **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 + + """ _surface_names = ('cone', 'plane') def __init__(self, x0=0., y0=0., z0=0., r2=1., up=True, **kwargs): @@ -141,6 +271,43 @@ class XConeOneSided(CompositeMixin, openmc.Surface): class YConeOneSided(CompositeMixin, openmc.Surface): + """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 + **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 + + """ _surface_names = ('cone', 'plane') def __init__(self, x0=0., y0=0., z0=0., r2=1., up=True, **kwargs): @@ -153,6 +320,43 @@ class YConeOneSided(CompositeMixin, openmc.Surface): class ZConeOneSided(CompositeMixin, openmc.Surface): + """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 + **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 + + """ _surface_names = ('cone', 'plane') def __init__(self, x0=0., y0=0., z0=0., r2=1., up=True, **kwargs): From f24a0586249fb32ddd3c2c1097a783b27614c331 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 2 Jul 2020 16:13:04 -0500 Subject: [PATCH 3/8] Fix rotate and add boundary_type property on composite surfaces --- openmc/model/surface_composite.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index e730390e9b..b3883d0c63 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -23,9 +23,21 @@ class CompositeMixin: surf = copy(self) for name in self._surface_names: s = getattr(surf, name) - setattr(surf, name, s.rotate(rotation)) + 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)) @@ -64,7 +76,7 @@ class RightCircularCylinder(CompositeMixin, openmc.Surface): Top planar surface of the cylinder """ - surface_names = ('cyl', 'bottom', 'top') + _surface_names = ('cyl', 'bottom', 'top') def __init__(self, center_base, height, radius, axis='z', **kwargs): cx, cy, cz = center_base @@ -119,7 +131,7 @@ class RectangularParallelepiped(CompositeMixin, openmc.Surface): Sides of the parallelepiped """ - surface_names = ('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax') + _surface_names = ('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax') def __init__(self, xmin, xmax, ymin, ymax, zmin, zmax, **kwargs): self.xmin = openmc.XPlane(x0=xmin, **kwargs) From 675b4f539c24e88a421a6b7e6e68d3fa23b74eac Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 2 Jul 2020 16:14:19 -0500 Subject: [PATCH 4/8] Add tests for composite surfaces --- tests/unit_tests/test_surface_composite.py | 221 +++++++++++++++++++++ 1 file changed, 221 insertions(+) create mode 100644 tests/unit_tests/test_surface_composite.py diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py new file mode 100644 index 0000000000..f7eec7a6d2 --- /dev/null +++ b/tests/unit_tests/test_surface_composite.py @@ -0,0 +1,221 @@ +from random import uniform + +import numpy as np +import openmc +import pytest + + +def test_box(): + x, y, z = 1.0, -2.5, 3.0 + dx, dy, dz = 4.0, 3.0, 6.0 + s = openmc.model.Box((x, y, z), (dx, 0, 0), (0, dy, 0), (0, 0, dz)) + 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((x + dx, y + dy, z + dz)) + assert ll == pytest.approx((x, y, z)) + + # __contains__ on associated half-spaces + assert (x - 1, 0, 0) in +s + assert (x - 1, 0, 0) not in -s + assert (x + dx/2, y + dy/2, z + dz/2) in -s + assert (x + dx/2, y + dy/2, z + dz/2) not in +s + + # evaluate method + with pytest.raises(NotImplementedError): + s.evaluate((0., 0., 0.)) + + # 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) + + # rotate method + s_rot = s.rotate((0., 0., 90.)) + ll, ur = (-s_rot).bounding_box + assert ll == pytest.approx((-(y + dy), x, z)) + assert ur == pytest.approx((-y, x + dx, z + dz)) + + # Make sure repr works + repr(s) + + # non-axis-aligned vectors not implemented + with pytest.raises(NotImplementedError): + openmc.model.Box((0, 0, 0), (1, 1, 1), (1, 1, 1), (1, 1, 1)) + + +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 + + # evaluate method + with pytest.raises(NotImplementedError): + s.evaluate((0., 0., 0.)) + + # 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 + + # evaluate method + with pytest.raises(NotImplementedError): + s.evaluate((0., 0., 0.)) + + # 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 + + # evaluate method + with pytest.raises(NotImplementedError): + s.evaluate((0., 0., 0.)) + + # 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) From 9bbd63185b686ab26847c873f28be60b8349887a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 6 Jul 2020 21:46:43 -0500 Subject: [PATCH 5/8] Remove Surface as parent class of composite surfaces --- openmc/model/surface_composite.py | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index b3883d0c63..914ffc3389 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -1,10 +1,11 @@ +from abc import ABC, abstractmethod from copy import copy import openmc -class CompositeMixin: - """Mixin class for composite surfaces""" +class CompositeSurface(ABC): + """Multiple primitive surfaces combined into a composite surface""" def evaluate(self, point): raise NotImplementedError('Composite surfaces do not have a surface equation.') @@ -41,8 +42,21 @@ class CompositeMixin: 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.""" -class RightCircularCylinder(CompositeMixin, openmc.Surface): + @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 @@ -100,7 +114,7 @@ class RightCircularCylinder(CompositeMixin, openmc.Surface): return +self.cyl | -self.bottom | +self.top -class RectangularParallelepiped(CompositeMixin, openmc.Surface): +class RectangularParallelepiped(CompositeSurface): """Rectangular parallelpiped composite surface A rectangular parallelpiped is composed of six planar surfaces. This class @@ -148,7 +162,7 @@ class RectangularParallelepiped(CompositeMixin, openmc.Surface): return -self.xmin | +self.ymax | -self.ymin | +self.ymax | -self.zmin | +self.zmax -class Box(CompositeMixin, openmc.Surface): +class Box(CompositeSurface): """Arbitrarily oriented orthogonal box An arbitrarily oriented orthogonal box is composed of six planar surfaces. @@ -227,7 +241,7 @@ class Box(CompositeMixin, openmc.Surface): -self.zmin | +self.zmax) -class XConeOneSided(CompositeMixin, openmc.Surface): +class XConeOneSided(CompositeSurface): """One-sided cone parallel the x-axis A one-sided cone is composed of a normal cone surface and an "ambiguity" @@ -282,7 +296,7 @@ class XConeOneSided(CompositeMixin, openmc.Surface): return (+self.cone & -self.plane) | +self.plane -class YConeOneSided(CompositeMixin, openmc.Surface): +class YConeOneSided(CompositeSurface): """One-sided cone parallel the y-axis A one-sided cone is composed of a normal cone surface and an "ambiguity" @@ -331,7 +345,7 @@ class YConeOneSided(CompositeMixin, openmc.Surface): __pos__ = XConeOneSided.__pos__ -class ZConeOneSided(CompositeMixin, openmc.Surface): +class ZConeOneSided(CompositeSurface): """One-sided cone parallel the z-axis A one-sided cone is composed of a normal cone surface and an "ambiguity" From f3a88bc3f95a25aee22e127c18b7b02101aac69c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Jul 2020 09:28:21 -0500 Subject: [PATCH 6/8] Address @eepeterson comments on #1602 --- openmc/model/surface_composite.py | 24 +++++++++++----------- tests/unit_tests/test_surface_composite.py | 12 ----------- 2 files changed, 12 insertions(+), 24 deletions(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 914ffc3389..8c4b99ae40 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -7,12 +7,6 @@ import openmc class CompositeSurface(ABC): """Multiple primitive surfaces combined into a composite surface""" - def evaluate(self, point): - raise NotImplementedError('Composite surfaces do not have a surface equation.') - - def _get_base_coeffs(self): - raise NotImplementedError('Composite surfaces do not have base coefficients.') - def translate(self, vector, inplace=False): surf = self if inplace else copy(self) for name in self._surface_names: @@ -264,7 +258,8 @@ class XConeOneSided(CompositeSurface): 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 + positive direction of the coordinate axis (the positive half-space of + the ambiguity plane) **kwargs Keyword arguments passed to underlying plane classes @@ -276,7 +271,8 @@ class XConeOneSided(CompositeSurface): Ambiguity surface up : bool Whether to select the side of the cone that extends to infinity in the - positive direction of the coordinate axis + positive direction of the coordinate axis (the positive half-space of + the ambiguity plane) """ _surface_names = ('cone', 'plane') @@ -319,7 +315,8 @@ class YConeOneSided(CompositeSurface): 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 + positive direction of the coordinate axis (the positive half-space of + the ambiguity plane) **kwargs Keyword arguments passed to underlying plane classes @@ -331,7 +328,8 @@ class YConeOneSided(CompositeSurface): Ambiguity surface up : bool Whether to select the side of the cone that extends to infinity in the - positive direction of the coordinate axis + positive direction of the coordinate axis (the positive half-space of + the ambiguity plane) """ _surface_names = ('cone', 'plane') @@ -368,7 +366,8 @@ class ZConeOneSided(CompositeSurface): 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 + positive direction of the coordinate axis (the positive half-space of + the ambiguity plane) **kwargs Keyword arguments passed to underlying plane classes @@ -380,7 +379,8 @@ class ZConeOneSided(CompositeSurface): Ambiguity surface up : bool Whether to select the side of the cone that extends to infinity in the - positive direction of the coordinate axis + positive direction of the coordinate axis (the positive half-space of + the ambiguity plane) """ _surface_names = ('cone', 'plane') diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py index f7eec7a6d2..28ded03cba 100644 --- a/tests/unit_tests/test_surface_composite.py +++ b/tests/unit_tests/test_surface_composite.py @@ -99,10 +99,6 @@ def test_rectangular_parallelepiped(): 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 - # evaluate method - with pytest.raises(NotImplementedError): - s.evaluate((0., 0., 0.)) - # translate method t = uniform(-5.0, 5.0) s_t = s.translate((t, t, t)) @@ -152,10 +148,6 @@ def test_right_circular_cylinder(axis, indices): assert point_neg in -s assert point_neg not in +s - # evaluate method - with pytest.raises(NotImplementedError): - s.evaluate((0., 0., 0.)) - # translate method t = uniform(-5.0, 5.0) s_t = s.translate((t, t, t)) @@ -206,10 +198,6 @@ def test_cone_one_sided(axis, point_pos, point_neg, ll_true): assert point_neg in -s assert point_neg not in +s - # evaluate method - with pytest.raises(NotImplementedError): - s.evaluate((0., 0., 0.)) - # translate method t = uniform(-5.0, 5.0) s_t = s.translate((t, t, t)) From 0ce9490187f9c81fe194a2d2afdda7b4d377b145 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Jul 2020 21:21:51 -0500 Subject: [PATCH 7/8] Remove (incomplete) Box composite surface --- docs/source/pythonapi/model.rst | 1 - openmc/model/surface_composite.py | 79 ---------------------- tests/unit_tests/test_surface_composite.py | 57 ---------------- 3 files changed, 137 deletions(-) diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 97704e76da..a6c89be7cd 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -25,7 +25,6 @@ Composite Surfaces :nosignatures: :template: myclass.rst - openmc.model.Box openmc.model.RectangularParallelepiped openmc.model.RightCircularCylinder openmc.model.XConeOneSided diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 8c4b99ae40..7f444d4b10 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -156,85 +156,6 @@ class RectangularParallelepiped(CompositeSurface): return -self.xmin | +self.ymax | -self.ymin | +self.ymax | -self.zmin | +self.zmax -class Box(CompositeSurface): - """Arbitrarily oriented orthogonal box - - An arbitrarily oriented orthogonal box 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 box. - - .. versionadded:: 0.12 - - Parameters - ---------- - v : iterable of float - Cartesian coordinates of a corner of the box - a1 : iterable of float - Vector of first side from the specified corner coordinates - a2 : iterable of float - Vector of second side from the specified corner coordinates - a3 : iterable of float - Vector of third side from the specified corner coordinates - **kwargs - Keyword arguments passed to underlying plane classes - - Attributes - ---------- - xmin, xmax : openmc.XPlane - Sides of the box - ymin, ymax : openmc.YPlane - Sides of the box - zmin, zmax : openmc.ZPlane - Sides of the box - - """ - _surface_names = ('xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax') - - def __init__(self, v, a1, a2, a3, **kwargs): - vx, vy, vz = v - a1x, a1y, a1z = a1 - a2x, a2y, a2z = a2 - a3x, a3y, a3z = a3 - - # Only support boxes with axis-aligned vectors - if any(x != 0.0 for x in (a1y, a1z, a2x, a2z, a3x, a3y)): - raise NotImplementedError('Box macrobody with non-axis-aligned ' - 'vector not supported.') - - # Determine each side of the box - if a1x > 0: - xmin, xmax = vx, vx + a1x - else: - xmin, xmax = vx + a1x, vx - if a2y > 0: - ymin, ymax = vy, vy + a2y - else: - ymin, ymax = vy + a2y, vy - if a3z > 0: - zmin, zmax = vz, vz + a3z - else: - zmin, zmax = vz + a3z, vz - - # Create surfaces - 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.xmax | - -self.ymin | +self.ymax | - -self.zmin | +self.zmax) - - class XConeOneSided(CompositeSurface): """One-sided cone parallel the x-axis diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py index 28ded03cba..a1a218311d 100644 --- a/tests/unit_tests/test_surface_composite.py +++ b/tests/unit_tests/test_surface_composite.py @@ -5,63 +5,6 @@ import openmc import pytest -def test_box(): - x, y, z = 1.0, -2.5, 3.0 - dx, dy, dz = 4.0, 3.0, 6.0 - s = openmc.model.Box((x, y, z), (dx, 0, 0), (0, dy, 0), (0, 0, dz)) - 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((x + dx, y + dy, z + dz)) - assert ll == pytest.approx((x, y, z)) - - # __contains__ on associated half-spaces - assert (x - 1, 0, 0) in +s - assert (x - 1, 0, 0) not in -s - assert (x + dx/2, y + dy/2, z + dz/2) in -s - assert (x + dx/2, y + dy/2, z + dz/2) not in +s - - # evaluate method - with pytest.raises(NotImplementedError): - s.evaluate((0., 0., 0.)) - - # 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) - - # rotate method - s_rot = s.rotate((0., 0., 90.)) - ll, ur = (-s_rot).bounding_box - assert ll == pytest.approx((-(y + dy), x, z)) - assert ur == pytest.approx((-y, x + dx, z + dz)) - - # Make sure repr works - repr(s) - - # non-axis-aligned vectors not implemented - with pytest.raises(NotImplementedError): - openmc.model.Box((0, 0, 0), (1, 1, 1), (1, 1, 1), (1, 1, 1)) - - def test_rectangular_parallelepiped(): xmin = uniform(-5., 5.) xmax = xmin + uniform(0., 5.) From f3dd04ebcc5eae7efd7994eb95f2af4773ed0979 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 9 Jul 2020 16:07:56 -0500 Subject: [PATCH 8/8] Add some value checks on parameters for composite surfaces --- openmc/model/surface_composite.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 7f444d4b10..9f9cd39c60 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -2,6 +2,7 @@ from abc import ABC, abstractmethod from copy import copy import openmc +from openmc.checkvalue import check_greater_than, check_value class CompositeSurface(ABC): @@ -88,6 +89,9 @@ class RightCircularCylinder(CompositeSurface): 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) @@ -142,6 +146,12 @@ class RectangularParallelepiped(CompositeSurface): _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) @@ -199,6 +209,7 @@ class XConeOneSided(CompositeSurface): _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 @@ -256,6 +267,7 @@ class YConeOneSided(CompositeSurface): _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 @@ -307,6 +319,7 @@ class ZConeOneSided(CompositeSurface): _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