mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-21 22:45:33 -04:00
1962 lines
76 KiB
Python
1962 lines
76 KiB
Python
from __future__ import annotations
|
|
from abc import ABC, abstractmethod
|
|
from collections.abc import Iterable, Sequence
|
|
from copy import copy
|
|
from functools import partial
|
|
from math import sqrt, pi, sin, cos, isclose
|
|
from numbers import Real
|
|
import warnings
|
|
import operator
|
|
|
|
import numpy as np
|
|
from scipy.spatial import ConvexHull, Delaunay
|
|
|
|
import openmc
|
|
from openmc.checkvalue import (check_greater_than, check_value, check_less_than,
|
|
check_iterable_type, check_length, check_type)
|
|
|
|
|
|
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 = self if inplace else 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 component_surfaces(self):
|
|
return [getattr(self, name) for name in self._surface_names]
|
|
|
|
@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
|
|
classes = (XConeOneSided, YConeOneSided, ZConeOneSided, Vessel)
|
|
for name in self._surface_names:
|
|
if isinstance(self, classes) and name.startswith('plane'):
|
|
continue
|
|
getattr(self, name).boundary_type = boundary_type
|
|
|
|
def __repr__(self):
|
|
return f"<{type(self).__name__} at 0x{id(self):x}>"
|
|
|
|
@property
|
|
@abstractmethod
|
|
def _surface_names(self):
|
|
"""Iterable of attribute names corresponding to underlying surfaces."""
|
|
|
|
@abstractmethod
|
|
def __neg__(self) -> openmc.Region:
|
|
"""Return the negative half-space of the composite surface."""
|
|
|
|
def __pos__(self) -> openmc.Region:
|
|
"""Return the positive half-space of the composite surface."""
|
|
return ~(-self)
|
|
|
|
|
|
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.
|
|
|
|
.. versionadded:: 0.13.1
|
|
|
|
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), (x, z), 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('r2 must be greater than r1.')
|
|
|
|
if theta2 <= theta1:
|
|
raise ValueError('theta2 must be greater than theta1.')
|
|
|
|
# Determine whether the angle between theta1 and theta2 is a reflex
|
|
# angle, in which case we need to use a union between the planar
|
|
# half-spaces
|
|
self._reflex = (theta2 - theta1 > 180.0)
|
|
|
|
phi1 = pi / 180 * theta1
|
|
phi2 = pi / 180 * theta2
|
|
|
|
# Coords for axis-perpendicular planes
|
|
p1 = np.array([center[0], center[1], 1.])
|
|
|
|
p2_plane1 = np.array([r1 * cos(phi1) + center[0], r1 * sin(phi1) + center[1], 0.])
|
|
p3_plane1 = np.array([r2 * cos(phi1) + center[0], r2 * sin(phi1) + center[1], 0.])
|
|
|
|
p2_plane2 = np.array([r1 * cos(phi2) + center[0], r1 * sin(phi2)+ center[1], 0.])
|
|
p3_plane2 = np.array([r2 * cos(phi2) + center[0], r2 * sin(phi2)+ center[1], 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 = [0, 2, 1]
|
|
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)
|
|
|
|
# Reorder the points to correspond to the correct central axis
|
|
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)
|
|
if axis == 'y':
|
|
self.plane1.flip_normal()
|
|
self.plane2.flip_normal()
|
|
|
|
@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), (x, z), 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('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):
|
|
if self._reflex:
|
|
return -self.outer_cyl & +self.inner_cyl & (-self.plane1 | +self.plane2)
|
|
else:
|
|
return -self.outer_cyl & +self.inner_cyl & -self.plane1 & +self.plane2
|
|
|
|
|
|
class IsogonalOctagon(CompositeSurface):
|
|
r"""Infinite isogonal octagon composite surface
|
|
|
|
An isogonal octagon is composed of eight planar surfaces. The prism is
|
|
parallel to the x, y, or z axis. The remaining two axes (y and z, x and z,
|
|
or x and y) serve as a basis for constructing the surfaces. Two surfaces
|
|
are parallel to the first basis axis, two surfaces are parallel
|
|
to the second basis axis, and the remaining four surfaces intersect both
|
|
basis axes at 45 degree angles.
|
|
|
|
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 octagonal prism.
|
|
|
|
.. versionadded:: 0.13.1
|
|
|
|
Parameters
|
|
----------
|
|
center : iterable of float
|
|
Coordinate for the central axis of the octagon in the
|
|
(y, z), (x, z), or (x, y) basis depending on the axis parameter.
|
|
r1 : float
|
|
Half-width of octagon across its basis axis-parallel sides in units
|
|
of cm. Must be less than :math:`r_2\sqrt{2}`.
|
|
r2 : float
|
|
Half-width of octagon across its basis axis intersecting sides in
|
|
units of cm. Must be less than than :math:`r_1\sqrt{2}`.
|
|
axis : {'x', 'y', 'z'}
|
|
Central axis of octagon. Defaults to 'z'
|
|
**kwargs
|
|
Keyword arguments passed to underlying plane classes
|
|
|
|
Attributes
|
|
----------
|
|
top : openmc.ZPlane, openmc.XPlane, or openmc.YPlane
|
|
Top planar surface of octagon
|
|
bottom : openmc.ZPlane, openmc.XPlane, or openmc.YPlane
|
|
Bottom planar surface of octagon
|
|
right : openmc.YPlane, openmc.ZPlane, or openmc.XPlane
|
|
Right planar surface of octagon
|
|
left : openmc.YPlane, openmc.ZPlane, or openmc.XPlane
|
|
Left planar surface of octagon
|
|
upper_right : openmc.Plane
|
|
Upper right planar surface of octagon
|
|
lower_right : openmc.Plane
|
|
Lower right planar surface of octagon
|
|
lower_left : openmc.Plane
|
|
Lower left planar surface of octagon
|
|
upper_left : openmc.Plane
|
|
Upper left planar surface of octagon
|
|
|
|
"""
|
|
|
|
_surface_names = ('top', 'bottom',
|
|
'upper_right', 'lower_left',
|
|
'right', 'left',
|
|
'lower_right', 'upper_left')
|
|
|
|
def __init__(self, center, r1, r2, axis='z', **kwargs):
|
|
c1, c2 = center
|
|
|
|
# Coordinates for axis-perpendicular planes
|
|
cright = c1 + r1
|
|
cleft = c1 - r1
|
|
|
|
ctop = c2 + r1
|
|
cbottom = c2 - r1
|
|
|
|
# Side lengths
|
|
if r2 > r1 * sqrt(2):
|
|
raise ValueError('r2 is greater than sqrt(2) * r1. Octagon' +
|
|
' may be erroneous.')
|
|
if r1 > r2 * sqrt(2):
|
|
raise ValueError('r1 is greater than sqrt(2) * r2. Octagon' +
|
|
' may be erroneous.')
|
|
|
|
L_basis_ax = (r2 * sqrt(2) - r1)
|
|
|
|
# Coordinates for quadrant planes
|
|
p1_ur = np.array([L_basis_ax, r1, 0.])
|
|
p2_ur = np.array([r1, L_basis_ax, 0.])
|
|
p3_ur = np.array([r1, L_basis_ax, 1.])
|
|
|
|
p1_lr = np.array([r1, -L_basis_ax, 0.])
|
|
p2_lr = np.array([L_basis_ax, -r1, 0.])
|
|
p3_lr = np.array([L_basis_ax, -r1, 1.])
|
|
|
|
p1_ll = -p1_ur
|
|
p2_ll = -p2_ur
|
|
p3_ll = -p3_ur
|
|
|
|
p1_ul = -p1_lr
|
|
p2_ul = -p2_lr
|
|
p3_ul = -p3_lr
|
|
|
|
points = [p1_ur, p2_ur, p3_ur, p1_lr, p2_lr, p3_lr,
|
|
p1_ll, p2_ll, p3_ll, p1_ul, p2_ul, p3_ul]
|
|
|
|
# Orientation specific variables
|
|
if axis == 'z':
|
|
coord_map = [0, 1, 2]
|
|
self.top = openmc.YPlane(ctop, **kwargs)
|
|
self.bottom = openmc.YPlane(cbottom, **kwargs)
|
|
self.right = openmc.XPlane(cright, **kwargs)
|
|
self.left = openmc.XPlane(cleft, **kwargs)
|
|
elif axis == 'y':
|
|
coord_map = [0, 2, 1]
|
|
self.top = openmc.ZPlane(ctop, **kwargs)
|
|
self.bottom = openmc.ZPlane(cbottom, **kwargs)
|
|
self.right = openmc.XPlane(cright, **kwargs)
|
|
self.left = openmc.XPlane(cleft, **kwargs)
|
|
elif axis == 'x':
|
|
coord_map = [2, 0, 1]
|
|
self.top = openmc.ZPlane(ctop, **kwargs)
|
|
self.bottom = openmc.ZPlane(cbottom, **kwargs)
|
|
self.right = openmc.YPlane(cright, **kwargs)
|
|
self.left = openmc.YPlane(cleft, **kwargs)
|
|
self.axis = axis
|
|
|
|
# Put our coordinates in (x,y,z) order and add the offset
|
|
for p in points:
|
|
p[0] += c1
|
|
p[1] += c2
|
|
p[:] = p[coord_map]
|
|
|
|
self.upper_right = openmc.Plane.from_points(p1_ur, p2_ur, p3_ur,
|
|
**kwargs)
|
|
self.lower_right = openmc.Plane.from_points(p1_lr, p2_lr, p3_lr,
|
|
**kwargs)
|
|
self.lower_left = openmc.Plane.from_points(p1_ll, p2_ll, p3_ll,
|
|
**kwargs)
|
|
self.upper_left = openmc.Plane.from_points(p1_ul, p2_ul, p3_ul,
|
|
**kwargs)
|
|
|
|
def __neg__(self):
|
|
if self.axis == 'y':
|
|
region = -self.top & +self.bottom & -self.right & +self.left & \
|
|
-self.upper_right & -self.lower_right & +self.lower_left & \
|
|
+self.upper_left
|
|
else:
|
|
region = -self.top & +self.bottom & -self.right & +self.left & \
|
|
+self.upper_right & +self.lower_right & -self.lower_left & \
|
|
-self.upper_left
|
|
|
|
return region
|
|
|
|
def __pos__(self):
|
|
if self.axis == 'y':
|
|
region = +self.top | -self.bottom | +self.right | -self.left | \
|
|
+self.upper_right | +self.lower_right | -self.lower_left | \
|
|
-self.upper_left
|
|
else:
|
|
region = +self.top | -self.bottom | +self.right | -self.left | \
|
|
-self.upper_right | -self.lower_right | +self.lower_left | \
|
|
+self.upper_left
|
|
return region
|
|
|
|
|
|
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
|
|
upper_fillet_radius : float
|
|
Upper edge fillet radius in [cm].
|
|
lower_fillet_radius : float
|
|
Lower edge fillet radius in [cm].
|
|
**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
|
|
upper_fillet_torus : openmc.Torus
|
|
Surface that creates the filleted edge for the upper end of the
|
|
cylinder. Only present if :attr:`upper_fillet_radius` is set.
|
|
upper_fillet_cylinder : openmc.Cylinder
|
|
Surface that bounds :attr:`upper_fillet_torus` radially. Only present
|
|
if :attr:`upper_fillet_radius` is set.
|
|
upper_fillet_plane : openmc.Plane
|
|
Surface that bounds :attr:`upper_fillet_torus` axially. Only present if
|
|
:attr:`upper_fillet_radius` is set.
|
|
lower_fillet_torus : openmc.Torus
|
|
Surface that creates the filleted edge for the lower end of the
|
|
cylinder. Only present if :attr:`lower_fillet_radius` is set.
|
|
lower_fillet_cylinder : openmc.Cylinder
|
|
Surface that bounds :attr:`lower_fillet_torus` radially. Only present
|
|
if :attr:`lower_fillet_radius` is set.
|
|
lower_fillet_plane : openmc.Plane
|
|
Surface that bounds :attr:`lower_fillet_torus` axially. Only present if
|
|
:attr:`lower_fillet_radius` is set.
|
|
|
|
"""
|
|
_surface_names = ('cyl', 'bottom', 'top')
|
|
|
|
def __init__(self, center_base, height, radius, axis='z',
|
|
upper_fillet_radius=0., lower_fillet_radius=0., **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'))
|
|
check_type('upper_fillet_radius', upper_fillet_radius, float)
|
|
check_less_than('upper_fillet_radius', upper_fillet_radius,
|
|
radius, equality=True)
|
|
check_type('lower_fillet_radius', lower_fillet_radius, float)
|
|
check_less_than('lower_fillet_radius', lower_fillet_radius,
|
|
radius, equality=True)
|
|
|
|
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)
|
|
x1, x2 = 'y', 'z'
|
|
axcoord, axcoord1, axcoord2 = 0, 1, 2
|
|
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)
|
|
x1, x2 = 'x', 'z'
|
|
axcoord, axcoord1, axcoord2 = 1, 0, 2
|
|
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)
|
|
x1, x2 = 'x', 'y'
|
|
axcoord, axcoord1, axcoord2 = 2, 0, 1
|
|
|
|
def _create_fillet_objects(axis_args, height, center_base, radius, fillet_radius, pos='upper'):
|
|
axis, x1, x2, axcoord, axcoord1, axcoord2 = axis_args
|
|
fillet_ext = height / 2 - fillet_radius
|
|
sign = 1
|
|
if pos == 'lower':
|
|
sign = -1
|
|
coord = center_base[axcoord] + (height / 2) + sign * fillet_ext
|
|
|
|
# cylinder
|
|
cyl_name = f'{pos}_min'
|
|
cylinder_args = {
|
|
x1 + '0': center_base[axcoord1],
|
|
x2 + '0': center_base[axcoord2],
|
|
'r': radius - fillet_radius
|
|
}
|
|
cls = getattr(openmc, f'{axis.upper()}Cylinder')
|
|
cyl = cls(name=f'{cyl_name} {axis}', **cylinder_args)
|
|
|
|
#torus
|
|
tor_name = f'{axis} {pos}'
|
|
tor_args = {
|
|
'a': radius - fillet_radius,
|
|
'b': fillet_radius,
|
|
'c': fillet_radius,
|
|
x1 + '0': center_base[axcoord1],
|
|
x2 + '0': center_base[axcoord2],
|
|
axis + '0': coord
|
|
}
|
|
cls = getattr(openmc, f'{axis.upper()}Torus')
|
|
torus = cls(name=tor_name, **tor_args)
|
|
|
|
# plane
|
|
p_name = f'{pos} ext'
|
|
p_args = {axis + '0': coord}
|
|
cls = getattr(openmc, f'{axis.upper()}Plane')
|
|
plane = cls(name=p_name, **p_args)
|
|
|
|
return cyl, torus, plane
|
|
|
|
if upper_fillet_radius > 0. or lower_fillet_radius > 0.:
|
|
if 'boundary_type' in kwargs:
|
|
if kwargs['boundary_type'] == 'periodic':
|
|
raise ValueError('Periodic boundary conditions not permitted when '
|
|
'rounded corners are used.')
|
|
|
|
axis_args = (axis, x1, x2, axcoord, axcoord1, axcoord2)
|
|
if upper_fillet_radius > 0.:
|
|
cylinder, torus, plane = _create_fillet_objects(
|
|
axis_args, height, center_base, radius, upper_fillet_radius)
|
|
self.upper_fillet_cylinder = cylinder
|
|
self.upper_fillet_torus = torus
|
|
self.upper_fillet_plane = plane
|
|
self._surface_names += ('upper_fillet_cylinder',
|
|
'upper_fillet_torus',
|
|
'upper_fillet_plane')
|
|
|
|
if lower_fillet_radius > 0.:
|
|
cylinder, torus, plane = _create_fillet_objects(
|
|
axis_args, height, center_base, radius, lower_fillet_radius,
|
|
pos='lower'
|
|
)
|
|
self.lower_fillet_cylinder = cylinder
|
|
self.lower_fillet_torus = torus
|
|
self.lower_fillet_plane = plane
|
|
|
|
self._surface_names += ('lower_fillet_cylinder',
|
|
'lower_fillet_torus',
|
|
'lower_fillet_plane')
|
|
|
|
def _get_fillet(self):
|
|
upper_fillet = self._get_upper_fillet()
|
|
lower_fillet = self._get_lower_fillet()
|
|
has_upper_fillet = upper_fillet is not None
|
|
has_lower_fillet = lower_fillet is not None
|
|
if has_lower_fillet and has_upper_fillet:
|
|
fillet = lower_fillet | upper_fillet
|
|
elif has_upper_fillet and not has_lower_fillet:
|
|
fillet = upper_fillet
|
|
elif not has_upper_fillet and has_lower_fillet:
|
|
fillet = lower_fillet
|
|
else:
|
|
fillet = None
|
|
return fillet
|
|
|
|
def _get_upper_fillet(self):
|
|
has_upper_fillet = hasattr(self, 'upper_fillet_plane')
|
|
if has_upper_fillet:
|
|
upper_fillet = +self.upper_fillet_cylinder & +self.upper_fillet_torus & +self.upper_fillet_plane
|
|
else:
|
|
upper_fillet = None
|
|
return upper_fillet
|
|
|
|
def _get_lower_fillet(self):
|
|
has_lower_fillet = hasattr(self, 'lower_fillet_plane')
|
|
if has_lower_fillet:
|
|
lower_fillet = +self.lower_fillet_cylinder & +self.lower_fillet_torus & -self.lower_fillet_plane
|
|
else:
|
|
lower_fillet = None
|
|
return lower_fillet
|
|
|
|
def __neg__(self):
|
|
prism = -self.cyl & +self.bottom & -self.top
|
|
fillet = self._get_fillet()
|
|
if fillet is not None:
|
|
prism = prism & ~fillet
|
|
return prism
|
|
|
|
def __pos__(self):
|
|
prism = +self.cyl | -self.bottom | +self.top
|
|
fillet = self._get_fillet()
|
|
if fillet is not None:
|
|
prism = prism | fillet
|
|
return prism
|
|
|
|
|
|
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.xmax & +self.xmin & -self.ymax & +self.ymin & -self.zmax & +self.zmin
|
|
|
|
def __pos__(self):
|
|
return +self.xmax | -self.xmin | +self.ymax | -self.ymin | +self.zmax | -self.zmin
|
|
|
|
|
|
class OrthogonalBox(CompositeSurface):
|
|
"""Arbitrarily oriented orthogonal box
|
|
|
|
This composite surface is composed of four or six planar surfaces that form
|
|
an arbitrarily oriented orthogonal box when combined.
|
|
|
|
Parameters
|
|
----------
|
|
v : iterable of float
|
|
(x,y,z) coordinates of a corner of the box
|
|
a1 : iterable of float
|
|
Vector of first side starting from ``v``
|
|
a2 : iterable of float
|
|
Vector of second side starting from ``v``
|
|
a3 : iterable of float, optional
|
|
Vector of third side starting from ``v``. When not specified, it is
|
|
assumed that the box will be infinite along the vector normal to the
|
|
plane specified by ``a1`` and ``a2``.
|
|
**kwargs
|
|
Keyword arguments passed to underlying plane classes
|
|
|
|
Attributes
|
|
----------
|
|
ax1_min, ax1_max : openmc.Plane
|
|
Planes representing minimum and maximum along first axis
|
|
ax2_min, ax2_max : openmc.Plane
|
|
Planes representing minimum and maximum along second axis
|
|
ax3_min, ax3_max : openmc.Plane
|
|
Planes representing minimum and maximum along third axis
|
|
|
|
"""
|
|
_surface_names = ('ax1_min', 'ax1_max', 'ax2_min', 'ax2_max', 'ax3_min', 'ax3_max')
|
|
|
|
def __init__(self, v, a1, a2, a3=None, **kwargs):
|
|
v = np.array(v)
|
|
a1 = np.array(a1)
|
|
a2 = np.array(a2)
|
|
if has_a3 := a3 is not None:
|
|
a3 = np.array(a3)
|
|
else:
|
|
a3 = np.cross(a1, a2) # normal to plane specified by a1 and a2
|
|
|
|
# Generate corners of box
|
|
p1 = v
|
|
p2 = v + a1
|
|
p3 = v + a2
|
|
p4 = v + a3
|
|
p5 = v + a1 + a2
|
|
p6 = v + a2 + a3
|
|
p7 = v + a1 + a3
|
|
|
|
# Generate 6 planes of box
|
|
self.ax1_min = openmc.Plane.from_points(p1, p3, p4, **kwargs)
|
|
self.ax1_max = openmc.Plane.from_points(p2, p5, p7, **kwargs)
|
|
self.ax2_min = openmc.Plane.from_points(p1, p4, p2, **kwargs)
|
|
self.ax2_max = openmc.Plane.from_points(p3, p6, p5, **kwargs)
|
|
if has_a3:
|
|
self.ax3_min = openmc.Plane.from_points(p1, p2, p3, **kwargs)
|
|
self.ax3_max = openmc.Plane.from_points(p4, p7, p6, **kwargs)
|
|
|
|
# Make sure a point inside the box produces the correct senses. If not,
|
|
# flip the plane coefficients so it does.
|
|
mid_point = v + (a1 + a2 + a3)/2
|
|
nums = (1, 2, 3) if has_a3 else (1, 2)
|
|
for num in nums:
|
|
min_surf = getattr(self, f'ax{num}_min')
|
|
max_surf = getattr(self, f'ax{num}_max')
|
|
if mid_point in -min_surf:
|
|
min_surf.flip_normal()
|
|
if mid_point in +max_surf:
|
|
max_surf.flip_normal()
|
|
|
|
def __neg__(self):
|
|
region = (+self.ax1_min & -self.ax1_max &
|
|
+self.ax2_min & -self.ax2_max)
|
|
if hasattr(self, 'ax3_min'):
|
|
region &= (+self.ax3_min & -self.ax3_max)
|
|
return region
|
|
|
|
|
|
class XConeOneSided(CompositeSurface):
|
|
r"""One-sided cone parallel the x-axis
|
|
|
|
A one-sided cone is composed of a normal cone surface and a "disambiguation"
|
|
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 in [cm].
|
|
y0 : float, optional
|
|
y-coordinate of the apex in [cm].
|
|
z0 : float, optional
|
|
z-coordinate of the apex in [cm].
|
|
r2 : float, optional
|
|
The square of the slope of the cone. It is defined as
|
|
:math:`\left(\frac{r}{h}\right)^2` for a radius, :math:`r` and an axial
|
|
distance :math:`h` from the apex. An easy way to define this quantity is
|
|
to take the square of the radius of the cone (in cm) 1 cm from the apex.
|
|
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
|
|
Disambiguation 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)
|
|
|
|
|
|
class YConeOneSided(CompositeSurface):
|
|
r"""One-sided cone parallel the y-axis
|
|
|
|
A one-sided cone is composed of a normal cone surface and a "disambiguation"
|
|
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 in [cm].
|
|
y0 : float, optional
|
|
y-coordinate of the apex in [cm].
|
|
z0 : float, optional
|
|
z-coordinate of the apex in [cm].
|
|
r2 : float, optional
|
|
The square of the slope of the cone. It is defined as
|
|
:math:`\left(\frac{r}{h}\right)^2` for a radius, :math:`r` and an axial
|
|
distance :math:`h` from the apex. An easy way to define this quantity is
|
|
to take the square of the radius of the cone (in cm) 1 cm from the apex.
|
|
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
|
|
Disambiguation 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__
|
|
|
|
|
|
class ZConeOneSided(CompositeSurface):
|
|
r"""One-sided cone parallel the z-axis
|
|
|
|
A one-sided cone is composed of a normal cone surface and a "disambiguation"
|
|
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 in [cm].
|
|
y0 : float, optional
|
|
y-coordinate of the apex in [cm].
|
|
z0 : float, optional
|
|
z-coordinate of the apex in [cm].
|
|
r2 : float, optional
|
|
The square of the slope of the cone. It is defined as
|
|
:math:`\left(\frac{r}{h}\right)^2` for a radius, :math:`r` and an axial
|
|
distance :math:`h` from the apex. An easy way to define this quantity is
|
|
to take the square of the radius of the cone (in cm) 1 cm from the apex.
|
|
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
|
|
Disambiguation 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__
|
|
|
|
|
|
class Polygon(CompositeSurface):
|
|
"""Polygon formed from a path of closed points.
|
|
|
|
.. versionadded:: 0.13.3
|
|
|
|
Parameters
|
|
----------
|
|
points : np.ndarray
|
|
An Nx2 array of points defining the vertices of the polygon.
|
|
basis : {'rz', 'xy', 'yz', 'xz'}, optional
|
|
2D basis set for the polygon. The polygon is two dimensional and has
|
|
infinite extent in the third (unspecified) dimension. For example, the
|
|
'xy' basis produces a polygon with infinite extent in the +/- z
|
|
direction. For the 'rz' basis the phi extent is infinite, thus forming
|
|
an axisymmetric surface.
|
|
|
|
Attributes
|
|
----------
|
|
points : np.ndarray
|
|
An Nx2 array of points defining the vertices of the polygon.
|
|
basis : {'rz', 'xy', 'yz', 'xz'}
|
|
2D basis set for the polygon.
|
|
regions : list of openmc.Region
|
|
A list of :class:`openmc.Region` objects, one for each of the convex polygons
|
|
formed during the decomposition of the input polygon.
|
|
region : openmc.Union
|
|
The union of all the regions comprising the polygon.
|
|
"""
|
|
|
|
def __init__(self, points, basis='rz'):
|
|
check_value('basis', basis, ('xy', 'yz', 'xz', 'rz'))
|
|
self._basis = basis
|
|
|
|
# Create a constrained triangulation of the validated points.
|
|
# The constrained triangulation is set to the _tri attribute
|
|
self._constrain_triangulation(self._validate_points(points))
|
|
|
|
# Decompose the polygon into groups of simplices forming convex subsets
|
|
# and get the sets of (surface, operator) pairs defining the polygon
|
|
self._surfsets = self._decompose_polygon_into_convex_sets()
|
|
|
|
# Set surface names as required by CompositeSurface protocol
|
|
surfnames = []
|
|
i = 0
|
|
for surfset in self._surfsets:
|
|
for surf, op, on_boundary in surfset:
|
|
if on_boundary:
|
|
setattr(self, f'surface_{i}', surf)
|
|
surfnames.append(f'surface_{i}')
|
|
i += 1
|
|
self._surfnames = tuple(surfnames)
|
|
|
|
# Generate a list of regions whose union represents the polygon.
|
|
regions = []
|
|
for surfs_ops in self._surfsets:
|
|
regions.append([op(surf) for surf, op, _ in surfs_ops])
|
|
self._regions = [openmc.Intersection(regs) for regs in regions]
|
|
|
|
# Create the union of all the convex subsets
|
|
self._region = openmc.Union(self._regions)
|
|
|
|
def __neg__(self):
|
|
return self._region
|
|
|
|
@property
|
|
def _surface_names(self):
|
|
return self._surfnames
|
|
|
|
@CompositeSurface.boundary_type.setter
|
|
def boundary_type(self, boundary_type):
|
|
if boundary_type != 'transmission':
|
|
warnings.warn("Setting boundary_type to a value other than "
|
|
"'transmission' on Polygon composite surfaces can "
|
|
"result in unintended behavior. Please use the "
|
|
"regions property of the Polygon to generate "
|
|
"individual openmc.Cell objects to avoid unwanted "
|
|
"behavior.")
|
|
for name in self._surface_names:
|
|
getattr(self, name).boundary_type = boundary_type
|
|
|
|
@property
|
|
def points(self):
|
|
return self._tri.points
|
|
|
|
@property
|
|
def basis(self):
|
|
return self._basis
|
|
|
|
@property
|
|
def _normals(self):
|
|
"""Generate the outward normal unit vectors for the polygon."""
|
|
# Rotation matrix for 90 degree clockwise rotation (-90 degrees about z
|
|
# axis for an 'xy' basis).
|
|
rotation = np.array([[0., 1.], [-1., 0.]])
|
|
# Get the unit vectors that point from one point in the polygon to the
|
|
# next given that they are ordered counterclockwise and that the final
|
|
# point is connected to the first point
|
|
tangents = np.diff(self.points, axis=0, append=[self.points[0, :]])
|
|
tangents /= np.linalg.norm(tangents, axis=-1, keepdims=True)
|
|
# Rotate the tangent vectors clockwise by 90 degrees, which for a
|
|
# counter-clockwise ordered polygon will produce the outward normal
|
|
# vectors.
|
|
return rotation.dot(tangents.T).T
|
|
|
|
@property
|
|
def _equations(self):
|
|
normals = self._normals
|
|
equations = np.empty((normals.shape[0], 3))
|
|
equations[:, :2] = normals
|
|
equations[:, 2] = -np.sum(normals*self.points, axis=-1)
|
|
return equations
|
|
|
|
@property
|
|
def regions(self):
|
|
return self._regions
|
|
|
|
@property
|
|
def region(self):
|
|
return self._region
|
|
|
|
def _validate_points(self, points):
|
|
"""Ensure the closed path defined by points does not intersect and is
|
|
oriented counter-clockwise.
|
|
|
|
Parameters
|
|
----------
|
|
points : np.ndarray (Nx2)
|
|
An Nx2 array of coordinate pairs describing the vertices.
|
|
|
|
Returns
|
|
-------
|
|
ordered_points : the input points ordered counter-clockwise
|
|
"""
|
|
points = np.asarray(points, dtype=float)
|
|
check_iterable_type('points', points, float, min_depth=2, max_depth=2)
|
|
check_length('points', points[0, :], 2, 2)
|
|
|
|
# If the last point is the same as the first, remove it and make sure
|
|
# there are still at least 3 points for a valid polygon.
|
|
if np.allclose(points[0, :], points[-1, :]):
|
|
points = points[:-1, :]
|
|
check_length('points', points, 3)
|
|
|
|
if len(points) != len(np.unique(points, axis=0)):
|
|
raise ValueError('Duplicate points were detected in the Polygon input')
|
|
|
|
# Order the points counter-clockwise (necessary for offset method)
|
|
# Calculates twice the signed area of the polygon using the "Shoelace
|
|
# Formula" https://en.wikipedia.org/wiki/Shoelace_formula
|
|
# If signed area is positive the curve is oriented counter-clockwise.
|
|
# If the signed area is negative the curve is oriented clockwise.
|
|
xpts, ypts = points.T
|
|
if np.sum(ypts*(np.roll(xpts, 1) - np.roll(xpts, -1))) < 0:
|
|
points = points[::-1, :]
|
|
|
|
# Check if polygon is self-intersecting by comparing edges pairwise
|
|
n = len(points)
|
|
for i in range(n):
|
|
p0 = np.append(points[i, :], 0)
|
|
p1 = np.append(points[(i + 1) % n, :], 0)
|
|
for j in range(i + 1, n):
|
|
p2 = np.append(points[j, :], 0)
|
|
p3 = np.append(points[(j + 1) % n, :], 0)
|
|
# Compute orientation of p0 wrt p2->p3 line segment
|
|
cp0 = np.cross(p3-p0, p2-p0)[-1]
|
|
# Compute orientation of p1 wrt p2->p3 line segment
|
|
cp1 = np.cross(p3-p1, p2-p1)[-1]
|
|
# Compute orientation of p2 wrt p0->p1 line segment
|
|
cp2 = np.cross(p1-p2, p0-p2)[-1]
|
|
# Compute orientation of p3 wrt p0->p1 line segment
|
|
cp3 = np.cross(p1-p3, p0-p3)[-1]
|
|
|
|
# Group cross products in an array and find out how many are 0
|
|
cross_products = np.array([[cp0, cp1], [cp2, cp3]])
|
|
cps_near_zero = np.isclose(cross_products, 0).astype(int)
|
|
num_zeros = np.sum(cps_near_zero)
|
|
|
|
# Topologies of 2 finite line segments categorized by the number
|
|
# of zero-valued cross products:
|
|
#
|
|
# 0: No 3 points lie on the same line
|
|
# 1: 1 point lies on the same line defined by the other line
|
|
# segment, but is not coincident with either of the points
|
|
# 2: 2 points are coincident, but the line segments are not
|
|
# collinear which guarantees no intersection
|
|
# 3: not possible, except maybe floating point issues?
|
|
# 4: Both line segments are collinear, simply need to check if
|
|
# they overlap or not
|
|
# adapted from algorithm linked below and modified to only
|
|
# consider intersections on the interior of line segments as
|
|
# proper intersections: i.e. segments sharing end points do not
|
|
# count as intersections.
|
|
# https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/
|
|
|
|
if num_zeros == 0:
|
|
# If the orientations of p0 and p1 have opposite signs
|
|
# and the orientations of p2 and p3 have opposite signs
|
|
# then there is an intersection.
|
|
if all(np.prod(cross_products, axis=-1) < 0):
|
|
raise ValueError('Polygon cannot be self-intersecting')
|
|
continue
|
|
|
|
elif num_zeros == 1:
|
|
# determine which line segment has 2 out of the 3 collinear
|
|
# points
|
|
idx = np.argwhere(np.sum(cps_near_zero, axis=-1) == 0)
|
|
if np.prod(cross_products[idx, :]) < 0:
|
|
raise ValueError('Polygon cannot be self-intersecting')
|
|
continue
|
|
|
|
elif num_zeros == 2:
|
|
continue
|
|
|
|
elif num_zeros == 3:
|
|
warnings.warn('Unclear if Polygon is self-intersecting')
|
|
continue
|
|
|
|
else:
|
|
# All 4 cross products are zero
|
|
# Determine number of unique points, x span and y span for
|
|
# both line segments
|
|
xmin1, xmax1 = min(p0[0], p1[0]), max(p0[0], p1[0])
|
|
ymin1, ymax1 = min(p0[1], p1[1]), max(p0[1], p1[1])
|
|
xmin2, xmax2 = min(p2[0], p3[0]), max(p2[0], p3[0])
|
|
ymin2, ymax2 = min(p2[1], p3[1]), max(p2[1], p3[1])
|
|
xlap = xmin1 < xmax2 and xmin2 < xmax1
|
|
ylap = ymin1 < ymax2 and ymin2 < ymax1
|
|
if xlap or ylap:
|
|
raise ValueError('Polygon cannot be self-intersecting')
|
|
continue
|
|
|
|
return points
|
|
|
|
def _constrain_triangulation(self, points, depth=0):
|
|
"""Generate a constrained triangulation by ensuring all edges of the
|
|
Polygon are contained within the simplices.
|
|
|
|
Parameters
|
|
----------
|
|
points : np.ndarray (Nx2)
|
|
An Nx2 array of coordinate pairs describing the vertices. These
|
|
points represent a planar straight line graph.
|
|
|
|
Returns
|
|
-------
|
|
None
|
|
"""
|
|
# Only attempt the triangulation up to 5 times.
|
|
if depth > 4:
|
|
raise RuntimeError('Could not create a valid triangulation after 5'
|
|
' attempts')
|
|
|
|
tri = Delaunay(points, qhull_options='QJ')
|
|
# Loop through the boundary edges of the polygon. If an edge is not
|
|
# included in the triangulation, break it into two line segments.
|
|
n = len(points)
|
|
new_pts = []
|
|
for i, j in zip(range(n), range(1, n + 1)):
|
|
# If both vertices of any edge are not found in any simplex, insert
|
|
# a new point between them.
|
|
if not any([i in s and j % n in s for s in tri.simplices]):
|
|
newpt = (points[i, :] + points[j % n, :]) / 2
|
|
new_pts.append((j, newpt))
|
|
|
|
# If all the edges are included in the triangulation set it, otherwise
|
|
# try again with additional points inserted on offending edges.
|
|
if not new_pts:
|
|
self._tri = tri
|
|
else:
|
|
for i, pt in new_pts[::-1]:
|
|
points = np.insert(points, i, pt, axis=0)
|
|
self._constrain_triangulation(points, depth=depth + 1)
|
|
|
|
def _group_simplices(self, neighbor_map, group=None):
|
|
"""Generate a convex grouping of simplices.
|
|
|
|
Parameters
|
|
----------
|
|
neighbor_map : dict
|
|
A map whose keys are simplex indices for simplices inside the polygon
|
|
and whose values are a list of simplex indices that neighbor this
|
|
simplex and are also inside the polygon.
|
|
group : list
|
|
A list of simplex indices that comprise the current convex group.
|
|
|
|
Returns
|
|
-------
|
|
group : list
|
|
The list of simplex indices that comprise the complete convex group.
|
|
"""
|
|
# If neighbor_map is empty there's nothing left to do
|
|
if not neighbor_map:
|
|
return group
|
|
# If group is empty, grab the next simplex in the dictionary and recurse
|
|
if group is None:
|
|
# Start with smallest neighbor lists
|
|
sidx = sorted(neighbor_map.items(), key=lambda item: len(item[1]))[0][0]
|
|
return self._group_simplices(neighbor_map, group=[sidx])
|
|
# Otherwise use the last simplex in the group
|
|
else:
|
|
sidx = group[-1]
|
|
# Remove current simplex from dictionary since it is in a group
|
|
neighbors = neighbor_map.pop(sidx, [])
|
|
# For each neighbor check if it is part of the same convex
|
|
# hull as the rest of the group. If yes, recurse. If no, continue on.
|
|
for n in neighbors:
|
|
if n in group or neighbor_map.get(n) is None:
|
|
continue
|
|
test_group = group + [n]
|
|
test_point_idx = np.unique(self._tri.simplices[test_group, :])
|
|
test_points = self.points[test_point_idx]
|
|
test_hull = ConvexHull(test_points, qhull_options='Qc')
|
|
pts_on_hull = len(test_hull.vertices) + len(test_hull.coplanar)
|
|
# If test_points are convex (including coplanar) keep adding to
|
|
# this group
|
|
if len(test_points) == pts_on_hull:
|
|
group = self._group_simplices(neighbor_map, group=test_group)
|
|
return group
|
|
|
|
def _get_convex_hull_surfs(self, qhull):
|
|
"""Generate a list of surfaces given by a set of linear equations
|
|
|
|
Parameters
|
|
----------
|
|
qhull : scipy.spatial.ConvexHull
|
|
A ConvexHull object representing the sub-region of the polygon.
|
|
|
|
Returns
|
|
-------
|
|
surfs_ops : list of (surface, operator) tuples
|
|
|
|
"""
|
|
basis = self.basis
|
|
boundary_eqns = self._equations
|
|
# Collect surface/operator pairs such that the intersection of the
|
|
# regions defined by these pairs is the inside of the polygon.
|
|
surfs_ops = []
|
|
# hull facet equation: dx*x + dy*y + c = 0
|
|
for dx, dy, c in qhull.equations:
|
|
# check if this facet is on the boundary of the polygon
|
|
facet_eq = np.array([dx, dy, c])
|
|
on_boundary = any([np.allclose(facet_eq, eq) for eq in boundary_eqns])
|
|
# Check if the facet is horizontal
|
|
if isclose(dx, 0, abs_tol=1e-8):
|
|
if basis in ('xz', 'yz', 'rz'):
|
|
surf = openmc.ZPlane(z0=-c/dy)
|
|
else:
|
|
surf = openmc.YPlane(y0=-c/dy)
|
|
# if (0, 1).(dx, dy) < 0 we want positive halfspace instead
|
|
op = operator.pos if dy < 0 else operator.neg
|
|
# Check if the facet is vertical
|
|
elif isclose(dy, 0, abs_tol=1e-8):
|
|
if basis in ('xy', 'xz'):
|
|
surf = openmc.XPlane(x0=-c/dx)
|
|
elif basis == 'yz':
|
|
surf = openmc.YPlane(y0=-c/dx)
|
|
else:
|
|
surf = openmc.ZCylinder(r=-c/dx)
|
|
# if (1, 0).(dx, dy) < 0 we want positive halfspace instead
|
|
op = operator.pos if dx < 0 else operator.neg
|
|
# Otherwise the facet is at an angle
|
|
else:
|
|
op = operator.neg
|
|
if basis == 'xy':
|
|
surf = openmc.Plane(a=dx, b=dy, c=0.0, d=-c)
|
|
elif basis == 'yz':
|
|
surf = openmc.Plane(a=0.0, b=dx, c=dy, d=-c)
|
|
elif basis == 'xz':
|
|
surf = openmc.Plane(a=dx, b=0.0, c=dy, d=-c)
|
|
else:
|
|
y0 = -c/dy
|
|
r2 = dy**2 / dx**2
|
|
# Check if the *slope* of the facet is positive. If dy/dx < 0
|
|
# then we want up to be True for the one-sided cones.
|
|
up = dy / dx < 0
|
|
surf = openmc.model.ZConeOneSided(z0=y0, r2=r2, up=up)
|
|
# if (1, -1).(dx, dy) < 0 for up cones we want positive halfspace
|
|
# if (1, 1).(dx, dy) < 0 for down cones we want positive halfspace
|
|
# otherwise we keep the negative halfspace operator
|
|
if (up and dx - dy < 0) or (not up and dx + dy < 0):
|
|
op = operator.pos
|
|
|
|
surfs_ops.append((surf, op, on_boundary))
|
|
|
|
return surfs_ops
|
|
|
|
def _decompose_polygon_into_convex_sets(self):
|
|
"""Decompose the Polygon into a set of convex polygons.
|
|
|
|
Returns
|
|
-------
|
|
surfsets : a list of lists of surface, operator pairs
|
|
"""
|
|
from matplotlib.path import Path
|
|
|
|
# Get centroids of all the simplices and determine if they are inside
|
|
# the polygon defined by input vertices or not.
|
|
centroids = np.mean(self.points[self._tri.simplices], axis=1)
|
|
in_polygon = Path(self.points).contains_points(centroids)
|
|
self._in_polygon = in_polygon
|
|
|
|
# Build a map with keys of simplex indices inside the polygon whose
|
|
# values are lists of that simplex's neighbors also inside the
|
|
# polygon
|
|
neighbor_map = {}
|
|
for i, nlist in enumerate(self._tri.neighbors):
|
|
if not in_polygon[i]:
|
|
continue
|
|
neighbor_map[i] = [n for n in nlist if in_polygon[n] and n >=0]
|
|
|
|
# Get the groups of simplices forming convex polygons whose union
|
|
# comprises the full input polygon. While there are still simplices
|
|
# left in the neighbor map, group them together into convex sets.
|
|
groups = []
|
|
while neighbor_map:
|
|
groups.append(self._group_simplices(neighbor_map))
|
|
self._groups = groups
|
|
|
|
# Generate lists of (surface, operator) pairs for each convex
|
|
# sub-region.
|
|
surfsets = []
|
|
for group in groups:
|
|
# Find all the unique points in the convex group of simplices,
|
|
# generate the convex hull and find the resulting surfaces and
|
|
# unary operators that represent this convex subset of the polygon.
|
|
idx = np.unique(self._tri.simplices[group, :])
|
|
qhull = ConvexHull(self.points[idx, :])
|
|
surf_ops = self._get_convex_hull_surfs(qhull)
|
|
surfsets.append(surf_ops)
|
|
return surfsets
|
|
|
|
def offset(self, distance: float | Sequence[float] | np.ndarray) -> Polygon:
|
|
"""Offset this polygon by a set distance
|
|
|
|
Parameters
|
|
----------
|
|
distance : float or sequence of float or np.ndarray
|
|
The distance to offset the polygon by. Positive is outward
|
|
(expanding) and negative is inward (shrinking). If a float is
|
|
provided, the same offset is applied to all vertices. If a list or
|
|
tuple is provided, each vertex gets a different offset. If an
|
|
iterable or numpy array is provided, each vertex gets a different
|
|
offset.
|
|
|
|
Returns
|
|
-------
|
|
offset_polygon : openmc.model.Polygon
|
|
"""
|
|
|
|
if isinstance(distance, float):
|
|
distance = np.full(len(self.points), distance)
|
|
elif isinstance(distance, Sequence):
|
|
distance = np.array(distance)
|
|
elif not isinstance(distance, np.ndarray):
|
|
raise TypeError("Distance must be a float or sequence of float.")
|
|
|
|
if len(distance) != len(self.points):
|
|
raise ValueError(
|
|
f"Length of distance {len(distance)} array must "
|
|
f"match number of polygon points {len(self.points)}"
|
|
)
|
|
|
|
normals = np.insert(self._normals, 0, self._normals[-1, :], axis=0)
|
|
cos2theta = np.sum(normals[1:, :]*normals[:-1, :], axis=-1, keepdims=True)
|
|
costheta = np.cos(np.arccos(cos2theta) / 2)
|
|
nvec = (normals[1:, :] + normals[:-1, :])
|
|
unit_nvec = nvec / np.linalg.norm(nvec, axis=-1, keepdims=True)
|
|
disp_vec = distance[:, np.newaxis] / costheta * unit_nvec
|
|
|
|
return type(self)(self.points + disp_vec, basis=self.basis)
|
|
|
|
|
|
class CruciformPrism(CompositeSurface):
|
|
"""Generalized cruciform prism
|
|
|
|
This surface represents a prism parallel to an axis formed by planes at
|
|
multiple distances from the center. Equivalent to the 'gcross' derived
|
|
surface in Serpent.
|
|
|
|
.. versionadded:: 0.14.0
|
|
|
|
Parameters
|
|
----------
|
|
distances : iterable of float
|
|
A monotonically increasing (or decreasing) iterable of distances in [cm]
|
|
that form the planes of the generalized cruciform.
|
|
center : iterable of float
|
|
The center of the prism in the two non-parallel axes (e.g., (x, y) when
|
|
axis is 'z') in [cm]
|
|
axis : {'x', 'y', 'z'}
|
|
Axis to which the prism is parallel
|
|
**kwargs
|
|
Keyword arguments passed to underlying plane classes
|
|
|
|
"""
|
|
|
|
def __init__(self, distances, center=(0., 0.), axis='z', **kwargs):
|
|
x0, y0 = center
|
|
self.distances = distances
|
|
|
|
if axis == 'x':
|
|
cls_horizontal = openmc.YPlane
|
|
cls_vertical = openmc.ZPlane
|
|
elif axis == 'y':
|
|
cls_horizontal = openmc.XPlane
|
|
cls_vertical = openmc.ZPlane
|
|
elif axis == 'z':
|
|
cls_horizontal = openmc.XPlane
|
|
cls_vertical = openmc.YPlane
|
|
else:
|
|
raise ValueError("axis must be 'x', 'y', or 'z'")
|
|
|
|
# Create each planar surface
|
|
surfnames = []
|
|
for i, d in enumerate(distances):
|
|
setattr(self, f'hmin{i}', cls_horizontal(x0 - d, **kwargs))
|
|
setattr(self, f'hmax{i}', cls_horizontal(x0 + d, **kwargs))
|
|
setattr(self, f'vmin{i}', cls_vertical(y0 - d, **kwargs))
|
|
setattr(self, f'vmax{i}', cls_vertical(y0 + d, **kwargs))
|
|
surfnames.extend([f'hmin{i}', f'hmax{i}', f'vmin{i}', f'vmax{i}'])
|
|
|
|
# Set _surfnames to satisfy CompositeSurface protocol
|
|
self._surfnames = tuple(surfnames)
|
|
|
|
@property
|
|
def _surface_names(self):
|
|
return self._surfnames
|
|
|
|
@property
|
|
def distances(self):
|
|
return self._distances
|
|
|
|
@distances.setter
|
|
def distances(self, values):
|
|
values = np.array(values, dtype=float)
|
|
# check for positive values
|
|
if not (values > 0).all():
|
|
raise ValueError("distances must be positive")
|
|
# Check for monotonicity
|
|
if (values[1:] > values[:-1]).all() or (values[1:] < values[:-1]).all():
|
|
self._distances = values
|
|
else:
|
|
raise ValueError("distances must be monotonic")
|
|
|
|
def __neg__(self):
|
|
n = len(self.distances)
|
|
regions = []
|
|
for i in range(n):
|
|
regions.append(
|
|
+getattr(self, f'hmin{i}') &
|
|
-getattr(self, f'hmax{i}') &
|
|
+getattr(self, f'vmin{n-1-i}') &
|
|
-getattr(self, f'vmax{n-1-i}')
|
|
)
|
|
return openmc.Union(regions)
|
|
|
|
|
|
# Define function to create a plane on given axis
|
|
def _plane(axis, name, value, boundary_type='transmission', albedo=1.0):
|
|
cls = getattr(openmc, f'{axis.upper()}Plane')
|
|
return cls(value, name=f'{name} {axis}',
|
|
boundary_type=boundary_type, albedo=albedo)
|
|
|
|
|
|
class RectangularPrism(CompositeSurface):
|
|
"""Infinite rectangular prism bounded by four planar surfaces.
|
|
|
|
.. versionadded:: 0.14.0
|
|
|
|
Parameters
|
|
----------
|
|
width : float
|
|
Prism width in units of [cm]. The width is aligned with the x, x, or z
|
|
axes for prisms parallel to the x, y, or z axis, respectively.
|
|
height : float
|
|
Prism height in units of [cm]. The height is aligned with the x, y, or z
|
|
axes for prisms parallel to the x, y, or z axis, respectively.
|
|
axis : {'x', 'y', 'z'}
|
|
Axis with which the infinite length of the prism should be aligned.
|
|
origin : Iterable of two floats
|
|
Origin of the prism. The two floats correspond to (y,z), (x,z) or (x,y)
|
|
for prisms parallel to the x, y or z axis, respectively.
|
|
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}
|
|
Boundary condition that defines the behavior for particles hitting the
|
|
surfaces comprising the rectangular prism.
|
|
albedo : float, optional
|
|
Albedo of the prism's surfaces as a ratio of particle weight after
|
|
interaction with the surface to the initial weight. Values must be
|
|
positive. Only applicable if the boundary type is 'reflective',
|
|
'periodic', or 'white'.
|
|
corner_radius : float
|
|
Prism corner radius in units of [cm].
|
|
|
|
"""
|
|
_surface_names = ('min_x1', 'max_x1', 'min_x2', 'max_x2')
|
|
|
|
def __init__(
|
|
self,
|
|
width: float,
|
|
height: float,
|
|
axis: str = 'z',
|
|
origin: Sequence[float] = (0., 0.),
|
|
boundary_type: str = 'transmission',
|
|
albedo: float = 1.,
|
|
corner_radius: float = 0.
|
|
):
|
|
check_type('width', width, Real)
|
|
check_type('height', height, Real)
|
|
check_type('albedo', albedo, Real)
|
|
check_type('corner_radius', corner_radius, Real)
|
|
check_value('axis', axis, ('x', 'y', 'z'))
|
|
check_type('origin', origin, Iterable, Real)
|
|
|
|
if axis == 'x':
|
|
x1, x2 = 'y', 'z'
|
|
elif axis == 'y':
|
|
x1, x2 = 'x', 'z'
|
|
else:
|
|
x1, x2 = 'x', 'y'
|
|
|
|
# Get cylinder class corresponding to given axis
|
|
cyl = getattr(openmc, f'{axis.upper()}Cylinder')
|
|
|
|
# Create container for boundary arguments
|
|
bc_args = {'boundary_type': boundary_type, 'albedo': albedo}
|
|
|
|
# Create rectangular region
|
|
self.min_x1 = _plane(x1, 'minimum', -width/2 + origin[0], **bc_args)
|
|
self.max_x1 = _plane(x1, 'maximum', width/2 + origin[0], **bc_args)
|
|
self.min_x2 = _plane(x2, 'minimum', -height/2 + origin[1], **bc_args)
|
|
self.max_x2 = _plane(x2, 'maximum', height/2 + origin[1], **bc_args)
|
|
if boundary_type == 'periodic':
|
|
self.min_x1.periodic_surface = self.max_x1
|
|
self.min_x2.periodic_surface = self.max_x2
|
|
|
|
# Handle rounded corners if given
|
|
if corner_radius > 0.:
|
|
if boundary_type == 'periodic':
|
|
raise ValueError('Periodic boundary conditions not permitted when '
|
|
'rounded corners are used.')
|
|
|
|
args = {'r': corner_radius, 'boundary_type': boundary_type, 'albedo': albedo}
|
|
|
|
args[x1 + '0'] = origin[0] - width/2 + corner_radius
|
|
args[x2 + '0'] = origin[1] - height/2 + corner_radius
|
|
self.x1_min_x2_min = cyl(name=f'{x1} min {x2} min', **args)
|
|
|
|
args[x1 + '0'] = origin[0] - width/2 + corner_radius
|
|
args[x2 + '0'] = origin[1] + height/2 - corner_radius
|
|
self.x1_min_x2_max = cyl(name=f'{x1} min {x2} max', **args)
|
|
|
|
args[x1 + '0'] = origin[0] + width/2 - corner_radius
|
|
args[x2 + '0'] = origin[1] - height/2 + corner_radius
|
|
self.x1_max_x2_min = cyl(name=f'{x1} max {x2} min', **args)
|
|
|
|
args[x1 + '0'] = origin[0] + width/2 - corner_radius
|
|
args[x2 + '0'] = origin[1] + height/2 - corner_radius
|
|
self.x1_max_x2_max = cyl(name=f'{x1} max {x2} max', **args)
|
|
|
|
self.x1_min = _plane(x1, 'min', -width/2 + origin[0] + corner_radius,
|
|
**bc_args)
|
|
self.x1_max = _plane(x1, 'max', width/2 + origin[0] - corner_radius,
|
|
**bc_args)
|
|
self.x2_min = _plane(x2, 'min', -height/2 + origin[1] + corner_radius,
|
|
**bc_args)
|
|
self.x2_max = _plane(x2, 'max', height/2 + origin[1] - corner_radius,
|
|
**bc_args)
|
|
self._surface_names += (
|
|
'x1_min_x2_min', 'x1_min_x2_max', 'x1_max_x2_min',
|
|
'x1_max_x2_max', 'x1_min', 'x1_max', 'x2_min', 'x2_max'
|
|
)
|
|
|
|
def __neg__(self):
|
|
prism = +self.min_x1 & -self.max_x1 & +self.min_x2 & -self.max_x2
|
|
|
|
# Cut out corners if a corner radius was given
|
|
if hasattr(self, 'x1_min'):
|
|
corners = (
|
|
(+self.x1_min_x2_min & -self.x1_min & -self.x2_min) |
|
|
(+self.x1_min_x2_max & -self.x1_min & +self.x2_max) |
|
|
(+self.x1_max_x2_min & +self.x1_max & -self.x2_min) |
|
|
(+self.x1_max_x2_max & +self.x1_max & +self.x2_max)
|
|
)
|
|
prism &= ~corners
|
|
|
|
return prism
|
|
|
|
|
|
class HexagonalPrism(CompositeSurface):
|
|
"""Hexagonal prism comoposed of six planar surfaces
|
|
|
|
.. versionadded:: 0.14.0
|
|
|
|
Parameters
|
|
----------
|
|
edge_length : float
|
|
Length of a side of the hexagon in [cm]
|
|
orientation : {'x', 'y'}
|
|
An 'x' orientation means that two sides of the hexagon are parallel to
|
|
the x-axis and a 'y' orientation means that two sides of the hexagon are
|
|
parallel to the y-axis.
|
|
origin : Iterable of two floats
|
|
Origin of the prism.
|
|
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}
|
|
Boundary condition that defines the behavior for particles hitting the
|
|
surfaces comprising the hexagonal prism.
|
|
albedo : float, optional
|
|
Albedo of the prism's surfaces as a ratio of particle weight after
|
|
interaction with the surface to the initial weight. Values must be
|
|
positive. Only applicable if the boundary type is 'reflective',
|
|
'periodic', or 'white'.
|
|
corner_radius : float
|
|
Prism corner radius in units of [cm].
|
|
|
|
"""
|
|
_surface_names = ('plane_max', 'plane_min', 'upper_right', 'upper_left',
|
|
'lower_right', 'lower_left')
|
|
|
|
def __init__(
|
|
self,
|
|
edge_length: float = 1.,
|
|
orientation: str = 'y',
|
|
origin: Sequence[float] = (0., 0.),
|
|
boundary_type: str = 'transmission',
|
|
albedo: float = 1.,
|
|
corner_radius: float = 0.
|
|
):
|
|
check_type('edge_length', edge_length, Real)
|
|
check_type('albedo', albedo, Real)
|
|
check_type('corner_radius', corner_radius, Real)
|
|
check_value('orientation', orientation, ('x', 'y'))
|
|
check_type('origin', origin, Iterable, Real)
|
|
|
|
l = edge_length
|
|
x, y = origin
|
|
|
|
# Create container for boundary arguments
|
|
bc_args = {'boundary_type': boundary_type, 'albedo': albedo}
|
|
|
|
if orientation == 'y':
|
|
# Left and right planes
|
|
self.plane_max = openmc.XPlane(x + sqrt(3.)/2*l, **bc_args)
|
|
self.plane_min = openmc.XPlane(x - sqrt(3.)/2*l, **bc_args)
|
|
c = sqrt(3.)/3.
|
|
|
|
# y = -x/sqrt(3) + a
|
|
self.upper_right = openmc.Plane(a=c, b=1., d=l+x*c+y, **bc_args)
|
|
|
|
# y = x/sqrt(3) + a
|
|
self.upper_left = openmc.Plane(a=-c, b=1., d=l-x*c+y, **bc_args)
|
|
|
|
# y = x/sqrt(3) - a
|
|
self.lower_right = openmc.Plane(a=-c, b=1., d=-l-x*c+y, **bc_args)
|
|
|
|
# y = -x/sqrt(3) - a
|
|
self.lower_left = openmc.Plane(a=c, b=1., d=-l+x*c+y, **bc_args)
|
|
|
|
elif orientation == 'x':
|
|
self.plane_max = openmc.YPlane(y + sqrt(3.)/2*l, **bc_args)
|
|
self.plane_min = openmc.YPlane(y - sqrt(3.)/2*l, **bc_args)
|
|
c = sqrt(3.)
|
|
|
|
# Upper-right surface: y = -sqrt(3)*(x - a)
|
|
self.upper_right = openmc.Plane(a=c, b=1., d=c*l+x*c+y, **bc_args)
|
|
|
|
# Lower-right surface: y = sqrt(3)*(x + a)
|
|
self.lower_right = openmc.Plane(a=-c, b=1., d=-c*l-x*c+y, **bc_args)
|
|
|
|
# Lower-left surface: y = -sqrt(3)*(x + a)
|
|
self.lower_left = openmc.Plane(a=c, b=1., d=-c*l+x*c+y, **bc_args)
|
|
|
|
# Upper-left surface: y = sqrt(3)*(x + a)
|
|
self.upper_left = openmc.Plane(a=-c, b=1., d=c*l-x*c+y, **bc_args)
|
|
|
|
# Handle periodic boundary conditions
|
|
if boundary_type == 'periodic':
|
|
self.plane_min.periodic_surface = self.plane_max
|
|
self.upper_right.periodic_surface = self.lower_left
|
|
self.lower_right.periodic_surface = self.upper_left
|
|
|
|
# Handle rounded corners if given
|
|
if corner_radius > 0.:
|
|
if boundary_type == 'periodic':
|
|
raise ValueError('Periodic boundary conditions not permitted '
|
|
'when rounded corners are used.')
|
|
|
|
c = sqrt(3.)/2
|
|
t = l - corner_radius/c
|
|
|
|
# Cylinder with corner radius and boundary type pre-applied
|
|
cyl1 = partial(openmc.ZCylinder, r=corner_radius, **bc_args)
|
|
cyl2 = partial(openmc.ZCylinder, r=corner_radius/(2*c), **bc_args)
|
|
|
|
if orientation == 'x':
|
|
self.x_min_y_min_in = cyl1(name='x min y min in', x0=x-t/2, y0=y-c*t)
|
|
self.x_min_y_max_in = cyl1(name='x min y max in', x0=x+t/2, y0=y-c*t)
|
|
self.x_max_y_min_in = cyl1(name='x max y min in', x0=x-t/2, y0=y+c*t)
|
|
self.x_max_y_max_in = cyl1(name='x max y max in', x0=x+t/2, y0=y+c*t)
|
|
self.min_in = cyl1(name='x min in', x0=x-t, y0=y)
|
|
self.max_in = cyl1(name='x max in', x0=x+t, y0=y)
|
|
|
|
self.x_min_y_min_out = cyl2(name='x min y min out', x0=x-l/2, y0=y-c*l)
|
|
self.x_min_y_max_out = cyl2(name='x min y max out', x0=x+l/2, y0=y-c*l)
|
|
self.x_max_y_min_out = cyl2(name='x max y min out', x0=x-l/2, y0=y+c*l)
|
|
self.x_max_y_max_out = cyl2(name='x max y max out', x0=x+l/2, y0=y+c*l)
|
|
self.min_out = cyl2(name='x min out', x0=x-l, y0=y)
|
|
self.max_out = cyl2(name='x max out', x0=x+l, y0=y)
|
|
|
|
elif orientation == 'y':
|
|
self.x_min_y_min_in = cyl1(name='x min y min in', x0=x-c*t, y0=y-t/2)
|
|
self.x_min_y_max_in = cyl1(name='x min y max in', x0=x-c*t, y0=y+t/2)
|
|
self.x_max_y_min_in = cyl1(name='x max y min in', x0=x+c*t, y0=y-t/2)
|
|
self.x_max_y_max_in = cyl1(name='x max y max in', x0=x+c*t, y0=y+t/2)
|
|
self.min_in = cyl1(name='y min in', x0=x, y0=y-t)
|
|
self.max_in = cyl1(name='y max in', x0=x, y0=y+t)
|
|
|
|
self.x_min_y_min_out = cyl2(name='x min y min out', x0=x-c*l, y0=y-l/2)
|
|
self.x_min_y_max_out = cyl2(name='x min y max out', x0=x-c*l, y0=y+l/2)
|
|
self.x_max_y_min_out = cyl2(name='x max y min out', x0=x+c*l, y0=y-l/2)
|
|
self.x_max_y_max_out = cyl2(name='x max y max out', x0=x+c*l, y0=y+l/2)
|
|
self.min_out = cyl2(name='y min out', x0=x, y0=y-l)
|
|
self.max_out = cyl2(name='y max out', x0=x, y0=y+l)
|
|
|
|
# Add to tuple of surface names
|
|
for s in ('in', 'out'):
|
|
self._surface_names += (
|
|
f'x_min_y_min_{s}', f'x_min_y_max_{s}',
|
|
f'x_max_y_min_{s}', f'x_max_y_max_{s}',
|
|
f'min_{s}', f'max_{s}')
|
|
|
|
def __neg__(self) -> openmc.Region:
|
|
prism = (
|
|
-self.plane_max & +self.plane_min &
|
|
-self.upper_right & -self.upper_left &
|
|
+self.lower_right & +self.lower_left
|
|
)
|
|
|
|
# Cut out corners if a corner radius was given
|
|
if hasattr(self, 'min_in'):
|
|
corners = (
|
|
+self.x_min_y_min_in & -self.x_min_y_min_out |
|
|
+self.x_min_y_max_in & -self.x_min_y_max_out |
|
|
+self.x_max_y_min_in & -self.x_max_y_min_out |
|
|
+self.x_max_y_max_in & -self.x_max_y_max_out |
|
|
+self.min_in & -self.min_out |
|
|
+self.max_in & -self.max_out
|
|
)
|
|
prism &= ~corners
|
|
|
|
return prism
|
|
|
|
|
|
def _rotation_matrix(v1, v2):
|
|
"""Compute rotation matrix that would rotate v1 into v2.
|
|
|
|
Parameters
|
|
----------
|
|
v1 : numpy.ndarray
|
|
Unrotated vector
|
|
v2 : numpy.ndarray
|
|
Rotated vector
|
|
|
|
Returns
|
|
-------
|
|
3x3 rotation matrix
|
|
|
|
"""
|
|
# Normalize vectors and compute cosine
|
|
u1 = v1 / np.linalg.norm(v1)
|
|
u2 = v2 / np.linalg.norm(v2)
|
|
cos_angle = np.dot(u1, u2)
|
|
|
|
I = np.identity(3)
|
|
|
|
# Handle special case where vectors are parallel or anti-parallel
|
|
if isclose(abs(cos_angle), 1.0, rel_tol=1e-8):
|
|
return np.sign(cos_angle)*I
|
|
else:
|
|
# Calculate rotation angle
|
|
sin_angle = np.sqrt(1 - cos_angle*cos_angle)
|
|
|
|
# Calculate axis of rotation
|
|
axis = np.cross(u1, u2)
|
|
axis /= np.linalg.norm(axis)
|
|
|
|
# Create cross-product matrix K
|
|
kx, ky, kz = axis
|
|
K = np.array([
|
|
[0.0, -kz, ky],
|
|
[kz, 0.0, -kx],
|
|
[-ky, kx, 0.0]
|
|
])
|
|
|
|
# Create rotation matrix using Rodrigues' rotation formula
|
|
return I + K * sin_angle + (K @ K) * (1 - cos_angle)
|
|
|
|
|
|
class ConicalFrustum(CompositeSurface):
|
|
"""Conical frustum.
|
|
|
|
A conical frustum, also known as a right truncated cone, is a cone that is
|
|
truncated by two parallel planes that are perpendicular to the axis of the
|
|
cone. The lower and upper base of the conical frustum are circular faces.
|
|
This surface is equivalent to the TRC macrobody in MCNP.
|
|
|
|
.. versionadded:: 0.15.1
|
|
|
|
Parameters
|
|
----------
|
|
center_base : iterable of float
|
|
Cartesian coordinates of the center of the bottom planar face.
|
|
axis : iterable of float
|
|
Vector from the center of the bottom planar face to the center of the
|
|
top planar face that defines the axis of the cone. The length of this
|
|
vector is the height of the conical frustum.
|
|
r1 : float
|
|
Radius of the lower cone base
|
|
r2 : float
|
|
Radius of the upper cone base
|
|
**kwargs
|
|
Keyword arguments passed to underlying plane classes
|
|
|
|
Attributes
|
|
----------
|
|
cone : openmc.Cone
|
|
Cone surface
|
|
plane_bottom : openmc.Plane
|
|
Plane surface defining the bottom of the frustum
|
|
plane_top : openmc.Plane
|
|
Plane surface defining the top of the frustum
|
|
|
|
"""
|
|
_surface_names = ('cone', 'plane_bottom', 'plane_top')
|
|
|
|
def __init__(self, center_base: Sequence[float], axis: Sequence[float],
|
|
r1: float, r2: float, **kwargs):
|
|
center_base = np.array(center_base)
|
|
axis = np.array(axis)
|
|
|
|
# Determine length of axis height vector
|
|
h = np.linalg.norm(axis)
|
|
|
|
# To create the frustum oriented with the correct axis, first we will
|
|
# create a cone along the z axis and then rotate it according to the
|
|
# given axis. Thus, we first need to determine the apex using the z axis
|
|
# as a reference.
|
|
x0, y0, z0 = center_base
|
|
if r1 != r2:
|
|
apex = z0 + r1*h/(r1 - r2)
|
|
r_sq = ((r1 - r2)/h)**2
|
|
cone = openmc.ZCone(x0, y0, apex, r2=r_sq, **kwargs)
|
|
else:
|
|
# In the degenerate case r1 == r2, the cone becomes a cylinder
|
|
cone = openmc.ZCylinder(x0, y0, r1, **kwargs)
|
|
|
|
# Create the parallel planes
|
|
plane_bottom = openmc.ZPlane(z0, **kwargs)
|
|
plane_top = openmc.ZPlane(z0 + h, **kwargs)
|
|
|
|
# Determine rotation matrix corresponding to specified axis
|
|
u = np.array([0., 0., 1.])
|
|
rotation = _rotation_matrix(u, axis)
|
|
|
|
# Rotate the surfaces
|
|
self.cone = cone.rotate(rotation, pivot=center_base)
|
|
self.plane_bottom = plane_bottom.rotate(rotation, pivot=center_base)
|
|
self.plane_top = plane_top.rotate(rotation, pivot=center_base)
|
|
|
|
def __neg__(self) -> openmc.Region:
|
|
return +self.plane_bottom & -self.plane_top & -self.cone
|
|
|
|
|
|
class Vessel(CompositeSurface):
|
|
"""Vessel composed of cylinder with semi-ellipsoid top and bottom.
|
|
|
|
This composite surface is represented by a finite cylinder with ellipsoidal
|
|
top and bottom surfaces. This surface is equivalent to the 'vesesl' surface
|
|
in Serpent.
|
|
|
|
.. versionadded:: 0.15.1
|
|
|
|
Parameters
|
|
----------
|
|
r : float
|
|
Radius of vessel.
|
|
p1 : float
|
|
Minimum coordinate for cylindrical part of vessel.
|
|
p2 : float
|
|
Maximum coordinate for cylindrical part of vessel.
|
|
h1 : float
|
|
Height of bottom ellipsoidal part of vessel.
|
|
h2 : float
|
|
Height of top ellipsoidal part of vessel.
|
|
center : 2-tuple of float
|
|
Coordinate for central axis of the cylinder in the (y, z), (x, z), or
|
|
(x, y) basis. Defaults to (0,0).
|
|
axis : {'x', 'y', 'z'}
|
|
Central axis of the cylinder.
|
|
|
|
"""
|
|
|
|
_surface_names = ('cyl', 'plane_bottom', 'plane_top', 'bottom', 'top')
|
|
|
|
def __init__(self, r: float, p1: float, p2: float, h1: float, h2: float,
|
|
center: Sequence[float] = (0., 0.), axis: str = 'z', **kwargs):
|
|
if p1 >= p2:
|
|
raise ValueError('p1 must be less than p2')
|
|
check_value('axis', axis, {'x', 'y', 'z'})
|
|
|
|
c1, c2 = center
|
|
cyl_class = getattr(openmc, f'{axis.upper()}Cylinder')
|
|
plane_class = getattr(openmc, f'{axis.upper()}Plane')
|
|
self.cyl = cyl_class(c1, c2, r, **kwargs)
|
|
self.plane_bottom = plane_class(p1)
|
|
self.plane_top = plane_class(p2)
|
|
|
|
# General equation for an ellipsoid:
|
|
# (x-x₀)²/r² + (y-y₀)²/r² + (z-z₀)²/h² = 1
|
|
# (x-x₀)² + (y-y₀)² + (z-z₀)²s² = r²
|
|
# Let s = r/h:
|
|
# (x² - 2x₀x + x₀²) + (y² - 2y₀y + y₀²) + (z² - 2z₀z + z₀²)s² = r²
|
|
# x² + y² + s²z² - 2x₀x - 2y₀y - 2s²z₀z + (x₀² + y₀² + z₀²s² - r²) = 0
|
|
|
|
sb = (r/h1)
|
|
st = (r/h2)
|
|
kwargs['a'] = kwargs['b'] = kwargs['c'] = 1.0
|
|
kwargs_bottom = kwargs
|
|
kwargs_top = kwargs.copy()
|
|
|
|
sb2 = sb*sb
|
|
st2 = st*st
|
|
kwargs_bottom['k'] = c1*c1 + c2*c2 + p1*p1*sb2 - r*r
|
|
kwargs_top['k'] = c1*c1 + c2*c2 + p2*p2*st2 - r*r
|
|
|
|
if axis == 'x':
|
|
kwargs_bottom['a'] *= sb2
|
|
kwargs_top['a'] *= st2
|
|
kwargs_bottom['g'] = -2*p1*sb2
|
|
kwargs_top['g'] = -2*p2*st2
|
|
kwargs_top['h'] = kwargs_bottom['h'] = -2*c1
|
|
kwargs_top['j'] = kwargs_bottom['j'] = -2*c2
|
|
elif axis == 'y':
|
|
kwargs_bottom['b'] *= sb2
|
|
kwargs_top['b'] *= st2
|
|
kwargs_top['g'] = kwargs_bottom['g'] = -2*c1
|
|
kwargs_bottom['h'] = -2*p1*sb2
|
|
kwargs_top['h'] = -2*p2*st2
|
|
kwargs_top['j'] = kwargs_bottom['j'] = -2*c2
|
|
elif axis == 'z':
|
|
kwargs_bottom['c'] *= sb2
|
|
kwargs_top['c'] *= st2
|
|
kwargs_top['g'] = kwargs_bottom['g'] = -2*c1
|
|
kwargs_top['h'] = kwargs_bottom['h'] = -2*c2
|
|
kwargs_bottom['j'] = -2*p1*sb2
|
|
kwargs_top['j'] = -2*p2*st2
|
|
|
|
self.bottom = openmc.Quadric(**kwargs_bottom)
|
|
self.top = openmc.Quadric(**kwargs_top)
|
|
|
|
def __neg__(self):
|
|
return ((-self.cyl & +self.plane_bottom & -self.plane_top) |
|
|
(-self.bottom & -self.plane_bottom) |
|
|
(-self.top & +self.plane_top))
|