diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 145330127..5dd3d3293 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -25,6 +25,7 @@ Composite Surfaces :nosignatures: :template: myclass.rst + openmc.model.CylinderSector openmc.model.IsogonalOctagon openmc.model.RectangularParallelepiped openmc.model.RightCircularCylinder diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index a34ce98cc..74f64c54a 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -1,11 +1,12 @@ from abc import ABC, abstractmethod from copy import copy -from math import sqrt +from math import sqrt, pi, sin, cos 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""" @@ -52,6 +53,164 @@ class CompositeSurface(ABC): """Return the negative half-space of the composite surface.""" +class CylinderSector(CompositeSurface): + """Infinite cylindrical sector composite surface. + + A cylinder sector is composed of two cylindrical and two planar surfaces. + The cylindrical surfaces are concentric, and the planar surfaces intersect + the central axis of the cylindrical 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 cylinder sector. + + Parameters + ---------- + r1 : float + Inner radius of sector. Must be less than r2. + r2 : float + Outer radius of sector. Must be greater than r1. + theta1 : float + Clockwise-most bound of sector in degrees. Assumed to be in the + counterclockwise direction with respect to the first basis axis + (+y, +z, or +x). Must be less than :attr:`theta2`. + theta2 : float + Counterclockwise-most bound of sector in degrees. Assumed to be in the + counterclockwise direction with respect to the first basis axis + (+y, +z, or +x). Must be greater than :attr:`theta1`. + center : iterable of float + Coordinate for central axes of cylinders in the (y, z), (z, x), or (x, y) + basis. Defaults to (0,0). + axis : {'x', 'y', 'z'} + Central axis of the cylinders defining the inner and outer surfaces of + the sector. Defaults to 'z'. + **kwargs : dict + Keyword arguments passed to the :class:`Cylinder` and + :class:`Plane` constructors. + + Attributes + ---------- + outer_cyl : openmc.ZCylinder, openmc.YCylinder, or openmc.XCylinder + Outer cylinder surface. + inner_cyl : openmc.ZCylinder, openmc.YCylinder, or openmc.XCylinder + Inner cylinder surface. + plane1 : openmc.Plane + Plane at angle :math:`\\theta_1` relative to the first basis axis. + plane2 : openmc.Plane + Plane at angle :math:`\\theta_2` relative to the first basis axis. + + """ + + _surface_names = ('outer_cyl', 'inner_cyl', 'plane1', 'plane2') + + def __init__(self, + r1, + r2, + theta1, + theta2, + center=(0.,0.), + axis='z', + **kwargs): + + if r2 <= r1: + raise ValueError(f'r2 must be greater than r1.') + + if theta2 <= theta1: + raise ValueError(f'theta2 must be greater than theta1.') + + phi1 = pi / 180 * theta1 + phi2 = pi / 180 * theta2 + + # Coords for axis-perpendicular planes + p1 = np.array([0., 0., 1.]) + + p2_plane1 = np.array([r1 * cos(phi1), r1 * sin(phi1), 0.]) + p3_plane1 = np.array([r2 * cos(phi1), r2 * sin(phi1), 0.]) + + p2_plane2 = np.array([r1 * cos(phi2), r1 * sin(phi2), 0.]) + p3_plane2 = np.array([r2 * cos(phi2), r2 * sin(phi2), 0.]) + + points = [p1, p2_plane1, p3_plane1, p2_plane2, p3_plane2] + if axis == 'z': + coord_map = [0, 1, 2] + self.inner_cyl = openmc.ZCylinder(*center, r1, **kwargs) + self.outer_cyl = openmc.ZCylinder(*center, r2, **kwargs) + elif axis == 'y': + coord_map = [1, 2, 0] + self.inner_cyl = openmc.YCylinder(*center, r1, **kwargs) + self.outer_cyl = openmc.YCylinder(*center, r2, **kwargs) + elif axis == 'x': + coord_map = [2, 0, 1] + self.inner_cyl = openmc.XCylinder(*center, r1, **kwargs) + self.outer_cyl = openmc.XCylinder(*center, r2, **kwargs) + + for p in points: + p[:] = p[coord_map] + + self.plane1 = openmc.Plane.from_points(p1, p2_plane1, p3_plane1, + **kwargs) + self.plane2 = openmc.Plane.from_points(p1, p2_plane2, p3_plane2, + **kwargs) + + @classmethod + def from_theta_alpha(cls, + r1, + r2, + theta, + alpha, + center = (0.,0.), + axis='z', + **kwargs): + r"""Alternate constructor for :class:`CylinderSector`. Returns a + :class:`CylinderSector` object based on a central angle :math:`\theta` + and an angular offset :math:`\alpha`. Note that + :math:`\theta_1 = \alpha` and :math:`\theta_2 = \alpha + \theta`. + + Parameters + ---------- + r1 : float + Inner radius of sector. Must be less than r2. + r2 : float + Outer radius of sector. Must be greater than r1. + theta : float + Central angle, :math:`\theta`, of the sector in degrees. Must be + greater that 0 and less than 360. + alpha : float + Angular offset, :math:`\alpha`, of sector in degrees. + The offset is in the counter-clockwise direction + with respect to the first basis axis (+y, +z, or +x). Note that + negative values translate to an offset in the clockwise direction. + center : iterable of float + Coordinate for central axes of cylinders in the (y, z), (z, x), or (x, y) + basis. Defaults to (0,0). + axis : {'x', 'y', 'z'} + Central axis of the cylinders defining the inner and outer surfaces + of the sector. Defaults to 'z'. + **kwargs : dict + Keyword arguments passed to the :class:`Cylinder` and + :class:`Plane` constructors. + + Returns + ------- + CylinderSector + CylinderSector with the given central angle at the given + offset. + """ + if theta >= 360. or theta <= 0: + raise ValueError(f'theta must be less than 360 and greater than 0.') + + theta1 = alpha + theta2 = alpha + theta + + return cls(r1, r2, theta1, theta2, center=center, axis=axis, **kwargs) + + def __neg__(self): + return -self.outer_cyl & +self.inner_cyl & -self.plane1 & +self.plane2 + + def __pos__(self): + return +self.outer_cyl | -self.inner_cyl | +self.plane1 | -self.plane2 + + class IsogonalOctagon(CompositeSurface): """Infinite isogonal octagon composite surface @@ -120,10 +279,10 @@ class IsogonalOctagon(CompositeSurface): # Side lengths if r2 > r1 * sqrt(2): - raise ValueError(f'r2 is greater than sqrt(2) * r1. Octagon' + \ + 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' + \ + raise ValueError(f'r1 is greater than sqrt(2) * r2. Octagon' + ' may be erroneous.') L_basis_ax = (r2 * sqrt(2) - r1) diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py index 14a7c8d3d..4d680ae7c 100644 --- a/tests/unit_tests/test_surface_composite.py +++ b/tests/unit_tests/test_surface_composite.py @@ -12,7 +12,8 @@ def test_rectangular_parallelepiped(): ymax = ymin + uniform(0., 5.) zmin = uniform(-5., 5.) zmax = zmin + uniform(0., 5.) - s = openmc.model.RectangularParallelepiped(xmin, xmax, ymin, ymax, zmin, zmax) + 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) @@ -152,6 +153,94 @@ def test_cone_one_sided(axis, point_pos, point_neg, ll_true): repr(s) +@pytest.mark.parametrize( + "axis, indices", [ + ("X", [0, 1, 2]), + ("Y", [1, 2, 0]), + ("Z", [2, 0, 1]), + ] +) +def test_cylinder_sector(axis, indices): + r1, r2 = 0.5, 1.5 + d = (r2 - r1) / 2 + phi1 = -60. + phi2 = 60 + s = openmc.model.CylinderSector(r1, r2, phi1, phi2, + axis=axis.lower()) + assert isinstance(s.outer_cyl, getattr(openmc, axis + "Cylinder")) + assert isinstance(s.inner_cyl, getattr(openmc, axis + "Cylinder")) + assert isinstance(s.plane1, openmc.Plane) + assert isinstance(s.plane2, openmc.Plane) + + # Make sure boundary condition propagates + s.boundary_type = 'reflective' + assert s.boundary_type == 'reflective' + assert s.outer_cyl.boundary_type == 'reflective' + assert s.inner_cyl.boundary_type == 'reflective' + assert s.plane1.boundary_type == 'reflective' + assert s.plane2.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(np.roll([-np.inf, -r2, -r2], indices[0])) + assert ur == pytest.approx(np.roll([np.inf, r2, r2], indices[0])) + + # __contains__ on associated half-spaces + point_pos = np.roll([0, r2 + 1, 0], indices[0]) + assert point_pos in +s + assert point_pos not in -s + point_neg = np.roll([0, r1 + d, r1 + d], 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) + + # Check invalid r1, r2 combinations + with pytest.raises(ValueError): + openmc.model.CylinderSector(r2, r1, phi1, phi2) + + # Check invalid angles + with pytest.raises(ValueError): + openmc.model.CylinderSector(r1, r2, phi2, phi1) + + # Make sure repr works + repr(s) + + +def test_cylinder_sector_from_theta_alpha(): + r1, r2 = 0.5, 1.5 + d = (r2 - r1) / 2 + theta = 120. + alpha = -60. + theta1 = alpha + theta2 = alpha + theta + s = openmc.model.CylinderSector(r1, r2, theta1, theta2) + s_alt = openmc.model.CylinderSector.from_theta_alpha(r1, + r2, + theta, + alpha) + + # Check that the angles are correct + assert s.plane1.coefficients == s_alt.plane1.coefficients + assert s.plane2.coefficients == s_alt.plane2.coefficients + assert s.inner_cyl.coefficients == s_alt.inner_cyl.coefficients + assert s.outer_cyl.coefficients == s_alt.outer_cyl.coefficients + + # Check invalid sector width + with pytest.raises(ValueError): + openmc.model.CylinderSector.from_theta_alpha(r1, r2, 360, alpha) + with pytest.raises(ValueError): + openmc.model.CylinderSector.from_theta_alpha(r1, r2, -1, alpha) + + @pytest.mark.parametrize( "axis, plane_tb, plane_lr, axis_idx", [ ("x", "Z", "Y", 0), @@ -225,5 +314,3 @@ def test_isogonal_octagon(axis, plane_tb, plane_lr, axis_idx): # Make sure repr works repr(s) - -