From f4501527f74166b8f71e4c6966863cacace6a827 Mon Sep 17 00:00:00 2001 From: yardasol Date: Fri, 24 Feb 2023 16:31:30 -0600 Subject: [PATCH] move fillet functionality to RightCircularCylinder --- openmc/model/surface_composite.py | 103 +++++++++++++++++++++++++++++- 1 file changed, 100 insertions(+), 3 deletions(-) diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py index 002fad9d9e8..4ab7fffd7d5 100644 --- a/openmc/model/surface_composite.py +++ b/openmc/model/surface_composite.py @@ -373,6 +373,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 @@ -393,24 +397,117 @@ def __init__(self, center_base, height, radius, axis='z', **kwargs): 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, Real) + check_less_than('upper_fillet_radius', upper_fillet_radius, + radius, equality=True) + check_type('lower_fillet_radius', lower_fillet_radius, Real) + 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 + + 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.: + upper_fillet_cyl, upper_fillet_tor, upper_fillet_plane = \ + self._create_fillet_objects(axis_args, + center_base, + radius, + upper_fillet_radius) + upper_fillet = +upper_fillet_cyl & +upper_fillet_torus & +upper_fillet_plane + else: + upper_fillet = None + + if lower_fillet_radius > 0.: + lower_fillet_cyl, lower_fillet_tor, lower_fillet_plane = \ + self._create_fillet_objects(axis_args, + center_base, + radius, + lower_fillet_radius, + pos='lower') + lower_fillet = +lower_fillet_cyl & +lower_fillet_torus & -lower_fillet_plane + else: + lower_fillet = None + + if lower_fillet is not None and upper_fillet is not None: + self._fillet = lower_fillet | upper_fillet + elif lower_fillet is None: + self._fillet = upper_fillet + else: + self._fillet = lower_fillet + + def _create_fillet_objects(self, axis_args, 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 + (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'{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 = tor(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 def __neg__(self): - return -self.cyl & +self.bottom & -self.top + prism = -self.cyl & +self.bottom & -self.top + if self._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 + if self._fillet is not None: + prism = prism | fillet + return prism class RectangularParallelepiped(CompositeSurface): @@ -807,7 +904,7 @@ def _group_simplices(self, neighbor_map, group=None): if group is None: sidx = next(iter(neighbor_map)) return self._group_simplices(neighbor_map, group=[sidx]) - # Otherwise use the last simplex in the group + # Otherwise use the last simplex in the group else: sidx = group[-1] # Remove current simplex from dictionary since it is in a group