move fillet functionality to RightCircularCylinder

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