Cylindrical prism with optional filleted edges

docs/source/pythonapi/model.rst

@@ -13,6 +13,7 @@ Convenience Functions
    openmc.model.borated_water
    openmc.model.hexagonal_prism
    openmc.model.rectangular_prism
+   openmc.model.cylindrical_prism
    openmc.model.subdivide
    openmc.model.pin
 

openmc/model/funcs.py

@@ -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
 
 
@@ -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.
@@ -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':
@@ -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
@@ -208,10 +209,10 @@ 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)
+        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)
 
         corners = (+x1_min_x2_min & -x1_min & -x2_min) | \
                   (+x1_min_x2_max & -x1_min & +x2_max) | \
@@ -265,8 +266,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
@@ -290,8 +291,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)
@@ -325,8 +326,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':
@@ -462,11 +464,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(

openmc/model/surface_composite.py

@@ -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):
@@ -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
 
@@ -383,33 +387,183 @@ 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
+
+        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_cylinder,
+                upper_fillet_torus,
+                upper_fillet_plane) = \
+                    self._create_fillet_objects(axis_args,
+                                                height,
+                                                center_base,
+                                                radius,
+                                                upper_fillet_radius)
+                self.upper_fillet_cylinder = upper_fillet_cylinder
+                self.upper_fillet_torus = upper_fillet_torus
+                self.upper_fillet_plane = upper_fillet_plane
+                self._surface_names += ('upper_fillet_cylinder',
+                                        'upper_fillet_torus',
+                                        'upper_fillet_plane')
+
+            if lower_fillet_radius > 0.:
+                (lower_fillet_cylinder,
+                lower_fillet_torus,
+                lower_fillet_plane) = \
+                    self._create_fillet_objects(axis_args,
+                                                height,
+                                                center_base,
+                                                radius,
+                                                lower_fillet_radius,
+                                                pos='lower')
+                self.lower_fillet_cylinder = lower_fillet_cylinder
+                self.lower_fillet_torus = lower_fillet_torus
+                self.lower_fillet_plane = lower_fillet_plane
+
+                self._surface_names += ('lower_fillet_cylinder',
+                                        'lower_fillet_torus',
+                                        'lower_fillet_plane')
+
+    def _create_fillet_objects(self, 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
+
+    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):

tests/unit_tests/test_geometry.py

@@ -3,6 +3,7 @@
 
 import numpy as np
 import openmc
+import openmc.lib
 import pytest