mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Merge pull request #2027 from yardasol/cylinder-sector
`CylinderSector` composite surface
This commit is contained in:
Paul Romano
GitHub
commit
9395c656dc
3 changed files with 253 additions and 6 deletions
|
|
@ -25,6 +25,7 @@ Composite Surfaces
|
|||
:nosignatures:
|
||||
:template: myclass.rst
|
||||
|
||||
openmc.model.CylinderSector
|
||||
openmc.model.IsogonalOctagon
|
||||
openmc.model.RectangularParallelepiped
|
||||
openmc.model.RightCircularCylinder
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue