Skip to content

Commit

Permalink
Cylindrical prism with optional filleted edges (#2309)
Browse files Browse the repository at this point in the history 
Co-authored-by: Paul Romano <[email protected]>
Co-authored-by: Patrick Shriwise <[email protected]>
Co-authored-by: Paul Wilson <[email protected]>
  • Loading branch information
@yardasol @paulromano
@pshriwise @gonuke
4 people committed Jul 18, 2023
1 parent 3f9cd0d commit 922b3f9
Show file tree
Hide file tree
Showing 3 changed files with 226 additions and 72 deletions.
68 changes: 37 additions & 31 deletions openmc/model/funcs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,9 @@
from operator import attrgetter
from warnings import warn

from openmc import (
XPlane, YPlane, Plane, ZCylinder, Cylinder, XCylinder,
YCylinder, Universe, Cell)
from ..checkvalue import (
check_type, check_value, check_length, check_less_than,
check_iterable_type)
from openmc import Plane, Cylinder, Universe, Cell
from ..checkvalue import (check_type, check_value, check_length,
check_less_than, check_iterable_type)
import openmc.data


Expand Down Expand Up @@ -111,6 +108,13 @@ def borated_water(boron_ppm, temperature=293., pressure=0.1013, temp_unit='K',
return out


# Define function to create a plane on given axis
def _plane(axis, name, value, boundary_type='transmission'):
cls = getattr(openmc, f'{axis.upper()}Plane')
return cls(value, name=f'{name} {axis}',
boundary_type=boundary_type)


def rectangular_prism(width, height, axis='z', origin=(0., 0.),
boundary_type='transmission', corner_radius=0.):
"""Get an infinite rectangular prism from four planar surfaces.
Expand Down Expand Up @@ -153,13 +157,6 @@ def rectangular_prism(width, height, axis='z', origin=(0., 0.),
check_value('axis', axis, ['x', 'y', 'z'])
check_type('origin', origin, Iterable, Real)

# Define function to create a plane on given axis
def plane(axis, name, value):
cls = globals()['{}Plane'.format(axis.upper())]
return cls(name='{} {}'.format(name, axis),
boundary_type=boundary_type,
**{axis + '0': value})

if axis == 'x':
x1, x2 = 'y', 'z'
elif axis == 'y':
Expand All @@ -168,13 +165,17 @@ def plane(axis, name, value):
x1, x2 = 'x', 'y'

# Get cylinder class corresponding to given axis
cyl = globals()['{}Cylinder'.format(axis.upper())]
cyl = getattr(openmc, f'{axis.upper()}Cylinder')

# Create rectangular region
min_x1 = plane(x1, 'minimum', -width/2 + origin[0])
max_x1 = plane(x1, 'maximum', width/2 + origin[0])
min_x2 = plane(x2, 'minimum', -height/2 + origin[1])
max_x2 = plane(x2, 'maximum', height/2 + origin[1])
min_x1 = _plane(x1, 'minimum', -width/2 + origin[0],
boundary_type=boundary_type)
max_x1 = _plane(x1, 'maximum', width/2 + origin[0],
boundary_type=boundary_type)
min_x2 = _plane(x2, 'minimum', -height/2 + origin[1],
boundary_type=boundary_type)
max_x2 = _plane(x2, 'maximum', height/2 + origin[1],
boundary_type=boundary_type)
if boundary_type == 'periodic':
min_x1.periodic_surface = max_x1
min_x2.periodic_surface = max_x2
Expand Down Expand Up @@ -208,10 +209,14 @@ def plane(axis, name, value):
args[x2 + '0'] = origin[1] + height/2 - corner_radius
x1_max_x2_max = cyl(name='{} max {} max'.format(x1, x2), **args)

x1_min = plane(x1, 'min', -width/2 + origin[0] + corner_radius)
x1_max = plane(x1, 'max', width/2 + origin[0] - corner_radius)
x2_min = plane(x2, 'min', -height/2 + origin[1] + corner_radius)
x2_max = plane(x2, 'max', height/2 + origin[1] - corner_radius)
x1_min = _plane(x1, 'min', -width/2 + origin[0] + corner_radius,
boundary_type=boundary_type)
x1_max = _plane(x1, 'max', width/2 + origin[0] - corner_radius,
boundary_type=boundary_type)
x2_min = _plane(x2, 'min', -height/2 + origin[1] + corner_radius,
boundary_type=boundary_type)
x2_max = _plane(x2, 'max', height/2 + origin[1] - corner_radius,
boundary_type=boundary_type)

corners = (+x1_min_x2_min & -x1_min & -x2_min) | \
(+x1_min_x2_max & -x1_min & +x2_max) | \
Expand Down Expand Up @@ -265,8 +270,8 @@ def hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
x, y = origin

if orientation == 'y':
right = XPlane(x + sqrt(3.)/2*l, boundary_type=boundary_type)
left = XPlane(x - sqrt(3.)/2*l, boundary_type=boundary_type)
right = openmc.XPlane(x + sqrt(3.)/2*l, boundary_type=boundary_type)
left = openmc.XPlane(x - sqrt(3.)/2*l, boundary_type=boundary_type)
c = sqrt(3.)/3.

# y = -x/sqrt(3) + a
Expand All @@ -290,8 +295,8 @@ def hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
lower_right.periodic_surface = upper_left

elif orientation == 'x':
top = YPlane(y0=y + sqrt(3.)/2*l, boundary_type=boundary_type)
bottom = YPlane(y0=y - sqrt(3.)/2*l, boundary_type=boundary_type)
top = openmc.YPlane(y0=y + sqrt(3.)/2*l, boundary_type=boundary_type)
bottom = openmc.YPlane(y0=y - sqrt(3.)/2*l, boundary_type=boundary_type)
c = sqrt(3.)

# y = -sqrt(3)*(x - a)
Expand Down Expand Up @@ -325,8 +330,9 @@ def hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
t = l - corner_radius/c

# Cylinder with corner radius and boundary type pre-applied
cyl1 = partial(ZCylinder, r=corner_radius, boundary_type=boundary_type)
cyl2 = partial(ZCylinder, r=corner_radius/(2*c),
cyl1 = partial(openmc.ZCylinder, r=corner_radius,
boundary_type=boundary_type)
cyl2 = partial(openmc.ZCylinder, r=corner_radius/(2*c),
boundary_type=boundary_type)

if orientation == 'x':
Expand Down Expand Up @@ -462,11 +468,11 @@ def pin(surfaces, items, subdivisions=None, divide_vols=True,
check_iterable_type("surfaces", surfaces[1:], surf_type)

# Check for increasing radii and equal centers
if surf_type is ZCylinder:
if surf_type is openmc.ZCylinder:
center_getter = attrgetter("x0", "y0")
elif surf_type is YCylinder:
elif surf_type is openmc.YCylinder:
center_getter = attrgetter("x0", "z0")
elif surf_type is XCylinder:
elif surf_type is openmc.XCylinder:
center_getter = attrgetter("z0", "y0")
else:
raise TypeError(
Expand Down
154 changes: 149 additions & 5 deletions openmc/model/surface_composite.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,8 @@
from scipy.spatial import ConvexHull, Delaunay

import openmc
from openmc.checkvalue import (check_greater_than, check_value,
check_iterable_type, check_length)
from openmc.checkvalue import (check_greater_than, check_value, check_less_than,
check_iterable_type, check_length, check_type)


class CompositeSurface(ABC):
Expand Down Expand Up @@ -372,6 +372,10 @@ class RightCircularCylinder(CompositeSurface):
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
Expand All @@ -383,33 +387,173 @@ class RightCircularCylinder(CompositeSurface):
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', **kwargs):
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):
return -self.cyl & +self.bottom & -self.top
prism = -self.cyl & +self.bottom & -self.top
fillet = self._get_fillet()
if fillet is not None:
prism = prism & ~fillet
return prism

def __pos__(self):
return +self.cyl | -self.bottom | +self.top
prism = +self.cyl | -self.bottom | +self.top
fillet = self._get_fillet()
if fillet is not None:
prism = prism | fillet
return prism


class RectangularParallelepiped(CompositeSurface):
Expand Down
Loading

0 comments on commit 922b3f9

Please sign in to comment.