From 455c655c41bba608ed7c6cfd517b9f0e38caca15 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Wed, 17 May 2017 10:59:55 -0400 Subject: [PATCH] moved corner radius into get_rectangular_prism and get_hexagonal_prism functions --- openmc/surface.py | 362 +++++++++++++++++++++++++++++++++------------- 1 file changed, 260 insertions(+), 102 deletions(-) diff --git a/openmc/surface.py b/openmc/surface.py index 3a445c54c1..eb786515d2 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1960,7 +1960,7 @@ class Halfspace(Region): return clone -def get_rectangular_prism(width, height, axis='z', origin=(0., 0.), +def get_rectangular_prism(width, height, corner_radius=0., axis='z', origin=(0., 0.), boundary_type='transmission'): """Get an infinite rectangular prism from four planar surfaces. @@ -1972,6 +1972,8 @@ def get_rectangular_prism(width, height, axis='z', origin=(0., 0.), height: float Prism height in units of cm. The height is aligned with the z, z, or y axes for prisms parallel to the x, y, or z axis, respectively. + corner_radius: float + Prism corner radius in units of cm. Defaults to 0. axis : {'x', 'y', 'z'} Axis with which the infinite length of the prism should be aligned. Defaults to 'z'. @@ -1992,6 +1994,7 @@ def get_rectangular_prism(width, height, axis='z', origin=(0., 0.), check_type('width', width, Real) check_type('height', height, Real) + check_type('corner_radius', corner_radius, Real) check_value('axis', axis, ['x','y','z']) check_type('origin', origin, Iterable, Real) @@ -2026,10 +2029,120 @@ def get_rectangular_prism(width, height, axis='z', origin=(0., 0.), boundary_type=boundary_type) prism = +min_x & -max_x & +min_y & -max_y + if corner_radius > 0.: + if axis == 'x': + y_min_z_min = XCylinder(name='y min z min', + y0=origin[0] - width /2. + corner_radius, + z0=origin[1] - height/2. + corner_radius, + R=corner_radius, + boundary_type=boundary_type) + y_min_z_max = XCylinder(name='y min z max', + y0=origin[0] - width /2. + corner_radius, + z0=origin[1] + height/2. - corner_radius, + R=corner_radius, + boundary_type=boundary_type) + y_max_z_min = XCylinder(name='y max z min', + y0=origin[0] + width /2. - corner_radius, + z0=origin[1] - height/2. + corner_radius, + R=corner_radius, + boundary_type=boundary_type) + y_max_z_max = XCylinder(name='y max z max', + y0=origin[0] + width /2. - corner_radius, + z0=origin[1] + height/2. - corner_radius, + R=corner_radius, + boundary_type=boundary_type) + y_min = YPlane(name='min y', y0=-width/2.+origin[0]+corner_radius, + boundary_type=boundary_type) + y_max = YPlane(name='max y', y0=+width/2.+origin[0]-corner_radius, + boundary_type=boundary_type) + z_min = ZPlane(name='min z', z0=-height/2.+origin[1]+corner_radius, + boundary_type=boundary_type) + z_max = ZPlane(name='max z', z0=+height/2.+origin[1]-corner_radius, + boundary_type=boundary_type) + + corners = (+y_min_z_min & -y_min & -z_min) | \ + (+y_min_z_max & -y_min & +z_max) | \ + (+y_max_z_min & +y_max & -z_min) | \ + (+y_max_z_max & +y_max & +z_max) + + elif axis == 'y': + x_min_z_min = YCylinder(name='x min z min', + x0=origin[0] - width /2. + corner_radius, + z0=origin[1] - height/2. + corner_radius, + R=corner_radius, + boundary_type=boundary_type) + x_min_z_max = YCylinder(name='x min z max', + x0=origin[0] - width /2. + corner_radius, + z0=origin[1] + height/2. - corner_radius, + R=corner_radius, + boundary_type=boundary_type) + x_max_z_min = YCylinder(name='x max z min', + x0=origin[0] + width /2. - corner_radius, + z0=origin[1] - height/2. + corner_radius, + R=corner_radius, + boundary_type=boundary_type) + x_max_z_max = YCylinder(name='x max z max', + x0=origin[0] + width /2. - corner_radius, + z0=origin[1] + height/2. - corner_radius, + R=corner_radius, + boundary_type=boundary_type) + + x_min = XPlane(name='min x', x0=-width/2.+origin[0]+corner_radius, + boundary_type=boundary_type) + x_max = XPlane(name='max x', x0=+width/2.+origin[0]-corner_radius, + boundary_type=boundary_type) + z_min = ZPlane(name='min z', z0=-height/2.+origin[1]+corner_radius, + boundary_type=boundary_type) + z_max = ZPlane(name='max z', z0=+height/2.+origin[1]-corner_radius, + boundary_type=boundary_type) + + corners = (+x_min_z_min & -x_min & -z_min) | \ + (+x_min_z_max & -x_min & +z_max) | \ + (+x_max_z_min & +x_max & -z_min) | \ + (+x_max_z_max & +x_max & +z_max) + + else: + x_min_y_min = ZCylinder(name='x min y min', + x0=origin[0] - width /2. + corner_radius, + y0=origin[1] - height/2. + corner_radius, + R=corner_radius, + boundary_type=boundary_type) + x_min_y_max = ZCylinder(name='x min y max', + x0=origin[0] - width /2. + corner_radius, + y0=origin[1] + height/2. - corner_radius, + R=corner_radius, + boundary_type=boundary_type) + x_max_y_min = ZCylinder(name='x max y min', + x0=origin[0] + width /2. - corner_radius, + y0=origin[1] - height/2. + corner_radius, + R=corner_radius, + boundary_type=boundary_type) + x_max_y_max = ZCylinder(name='x max y max', + x0=origin[0] + width /2. - corner_radius, + y0=origin[1] + height/2. - corner_radius, + R=corner_radius, + boundary_type=boundary_type) + + x_min = XPlane(name='min x', x0=-width/2.+origin[0]+corner_radius, + boundary_type=boundary_type) + x_max = XPlane(name='max x', x0=+width/2.+origin[0]-corner_radius, + boundary_type=boundary_type) + y_min = YPlane(name='min y', y0=-height/2.+origin[1]+corner_radius, + boundary_type=boundary_type) + y_max = YPlane(name='max y', y0=+height/2.+origin[1]-corner_radius, + boundary_type=boundary_type) + + corners = (+x_min_y_min & -x_min & -y_min) | \ + (+x_min_y_max & -x_min & +y_max) | \ + (+x_max_y_min & +x_max & -y_min) | \ + (+x_max_y_max & +x_max & +y_max) + + prism = prism & ~corners + return prism -def get_hexagonal_prism(edge_length=1., orientation='y', +def get_hexagonal_prism(edge_length=1., corner_radius=0., orientation='y', boundary_type='transmission'): """Create a hexagon region from six surface planes. @@ -2037,6 +2150,8 @@ def get_hexagonal_prism(edge_length=1., orientation='y', ---------- edge_length : float Length of a side of the hexagon in cm + corner_radius: float + Prism corner radius in units of cm. Defaults to 0. 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 @@ -2070,8 +2185,8 @@ def get_hexagonal_prism(edge_length=1., orientation='y', # y = -x/sqrt(3) - a lower_left = Plane(A=c, B=1., D=-l, boundary_type=boundary_type) - return Intersection(-right, +left, -upper_right, -upper_left, - +lower_right, +lower_left) + prism = -right & +left & -upper_right & -upper_left & \ + +lower_right & +lower_left elif orientation == 'x': top = YPlane(y0=sqrt(3.)/2.*l) @@ -2089,109 +2204,152 @@ def get_hexagonal_prism(edge_length=1., orientation='y', # y = sqrt(3)*(x + a) upper_left = Plane(A=-c, B=1., D=c*l, boundary_type=boundary_type) - return Intersection(-top, +bottom, -upper_right, +lower_right, - +lower_left, -upper_left) + prism = -top & +bottom & -upper_right & +lower_right & \ + +lower_left & -upper_left + if corner_radius > 0.: + if orientation == 'x': + c = sqrt(3.) + x_min_y_min_in = ZCylinder(name='x min y min in', + x0=-0.5 * (l - 2./c*corner_radius), + y0=-c/2. * (l - 2./c*corner_radius), + R=corner_radius, + boundary_type=boundary_type) + x_min_y_max_in = ZCylinder(name='x min y max in', + x0= 0.5 * (l - 2./c*corner_radius), + y0=-c/2. * (l - 2./c*corner_radius), + R=corner_radius, + boundary_type=boundary_type) + x_max_y_min_in = ZCylinder(name='x max y min in', + x0=-0.5 * (l - 2./c*corner_radius), + y0= c/2. * (l - 2./c*corner_radius), + R=corner_radius, + boundary_type=boundary_type) + x_max_y_max_in = ZCylinder(name='x max y max in', + x0= 0.5 * (l - 2./c*corner_radius), + y0= c/2. * (l - 2./c*corner_radius), + R=corner_radius, + boundary_type=boundary_type) + x_min_in = ZCylinder(name='x min in', + x0=-(l - 2./c*corner_radius), + y0=0., + R=corner_radius, + boundary_type=boundary_type) + x_max_in = ZCylinder(name='x max in', + x0= (l - 2./c*corner_radius), + y0=0., + R=corner_radius, + boundary_type=boundary_type) -def get_square_cylinder(r_cylinder, r_corner, axis='z', origin=(0., 0.), - boundary_type='transmission'): - """Get an infinite square cylinder using a collection of planes and cylinders. + x_min_y_min_out = ZCylinder(name='x min y min out', + x0=-0.5 * l, + y0=-c/2. * l, + R=corner_radius/c, + boundary_type=boundary_type) + x_min_y_max_out = ZCylinder(name='x min y max out', + x0=+0.5 * l, + y0=-c/2. * l, + R=corner_radius/c, + boundary_type=boundary_type) + x_max_y_min_out = ZCylinder(name='x max y min out', + x0=-0.5 * l, + y0= c/2. * l, + R=corner_radius/c, + boundary_type=boundary_type) + x_max_y_max_out = ZCylinder(name='x max y max out', + x0= 0.5 * l, + y0= c/2. * l, + R=corner_radius/c, + boundary_type=boundary_type) + x_min_out = ZCylinder(name='x min out', + x0=-l, + y0=0., + R=corner_radius/c, + boundary_type=boundary_type) + x_max_out = ZCylinder(name='x max out', + x0= l, + y0=0., + R=corner_radius/c, + boundary_type=boundary_type) - Parameters - ---------- - r_cylinder: float - Square cylinder radius in units of cm. The square cylinder is aligned - with the x, y, or z axis. - r_corner: float - Square cylinder corner radius in units of cm. - axis : {'x', 'y', 'z'} - Axis with which the infinite length of the square cylinder should be - aligned. - Defaults to 'z'. - origin: Iterable of two floats - Origin of the square cylinder. The two floats correspond to (y,z), (x,z) - or (x,y) for square cylinders parallel to the x, y or z axis, - respectively. - Defaults to (0., 0.). - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} - Boundary condition that defines the behavior for particles hitting the - surfaces comprising the square cylinder (default is 'transmission'). + corners = (+x_min_y_min_in & -x_min_y_min_out) | \ + (+x_min_y_max_in & -x_min_y_max_out) | \ + (+x_max_y_min_in & -x_max_y_min_out) | \ + (+x_max_y_max_in & -x_max_y_max_out) | \ + (+x_min_in & -x_min_out ) | \ + (+x_max_in & -x_max_out ) - Returns - ------- - openmc.Region - The inside of a square cylinder + elif orientation == 'y': + c = sqrt(3.) + x_min_y_min_in = ZCylinder(name='x min y min in', + x0=-c/2. * (l - 2./c*corner_radius), + y0=-0.5 * (l - 2./c*corner_radius), + R=corner_radius, + boundary_type=boundary_type) + x_min_y_max_in = ZCylinder(name='x min y max in', + x0=-c/2. * (l - 2./c*corner_radius), + y0= 0.5 * (l - 2./c*corner_radius), + R=corner_radius, + boundary_type=boundary_type) + x_max_y_min_in = ZCylinder(name='x max y min in', + x0= c/2. * (l - 2./c*corner_radius), + y0=-0.5 * (l - 2./c*corner_radius), + R=corner_radius, + boundary_type=boundary_type) + x_max_y_max_in = ZCylinder(name='x max y max in', + x0= c/2. * (l - 2./c*corner_radius), + y0= 0.5 * (l - 2./c*corner_radius), + R=corner_radius, + boundary_type=boundary_type) + y_min_in = ZCylinder(name='y min in', + x0=0., + y0=-(l - 2./c*corner_radius), + R=corner_radius, + boundary_type=boundary_type) + y_max_in = ZCylinder(name='y max in', + x0=0., + y0= (l - 2./c*corner_radius), + R=corner_radius, + boundary_type=boundary_type) - """ + x_min_y_min_out = ZCylinder(name='x min y min out', + x0=-c/2. * l, + y0=-0.5 * l, + R=corner_radius/c, + boundary_type=boundary_type) + x_min_y_max_out = ZCylinder(name='x min y max out', + x0=-c/2. * l, + y0=+0.5 * l, + R=corner_radius/c, + boundary_type=boundary_type) + x_max_y_min_out = ZCylinder(name='x max y min out', + x0= c/2. * l, + y0=-0.5 * l, + R=corner_radius/c, + boundary_type=boundary_type) + x_max_y_max_out = ZCylinder(name='x max y max out', + x0= c/2. * l, + y0= 0.5 * l, + R=corner_radius/c, + boundary_type=boundary_type) + y_min_out = ZCylinder(name='y min out', + x0=0., + y0=-l, + R=corner_radius/c, + boundary_type=boundary_type) + y_max_out = ZCylinder(name='y max out', + x0=0., + y0= l, + R=corner_radius/c, + boundary_type=boundary_type) - check_type('r_cylinder', r_cylinder, Real) - check_type('r_corner', r_corner, Real) - check_value('axis', axis, ['x','y','z']) - check_type('origin', origin, Iterable, Real) + corners = (+x_min_y_min_in & -x_min_y_min_out) | \ + (+x_min_y_max_in & -x_min_y_max_out) | \ + (+x_max_y_min_in & -x_max_y_min_out) | \ + (+x_max_y_max_in & -x_max_y_max_out) | \ + (+y_min_in & -y_min_out ) | \ + (+y_max_in & -y_max_out ) - if axis == 'x': - y_min_z_min = XCylinder(name='y min z min', y0=-r_cylinder+r_corner, - z0=-r_cylinder+r_corner, R=r_corner, - boundary_type=boundary_type) - y_min_z_max = XCylinder(name='y min z max', y0=-r_cylinder+r_corner, - z0= r_cylinder-r_corner, R=r_corner, - boundary_type=boundary_type) - y_max_z_min = XCylinder(name='y max z min', y0= r_cylinder-r_corner, - z0=-r_cylinder+r_corner, R=r_corner, - boundary_type=boundary_type) - y_max_z_max = XCylinder(name='y max z max', y0= r_cylinder-r_corner, - z0= r_cylinder-r_corner, R=r_corner, - boundary_type=boundary_type) - corners = -y_min_z_min | -y_min_z_max | -y_max_z_min | -y_max_z_max - elif axis == 'y': - x_min_z_min = YCylinder(name='x min z min', x0=-r_cylinder+r_corner, - z0=-r_cylinder+r_corner, R=r_corner, - boundary_type=boundary_type) - x_min_z_max = YCylinder(name='x min z max', x0=-r_cylinder+r_corner, - z0= r_cylinder-r_corner, R=r_corner, - boundary_type=boundary_type) - x_max_z_min = YCylinder(name='x max z min', x0= r_cylinder-r_corner, - z0=-r_cylinder+r_corner, R=r_corner, - boundary_type=boundary_type) - x_max_z_max = YCylinder(name='x max z max', x0= r_cylinder-r_corner, - z0= r_cylinder-r_corner, R=r_corner, - boundary_type=boundary_type) - corners = -x_min_z_min | -x_min_z_max | -x_max_z_min | -x_max_z_max - elif axis == 'z': - x_min_y_min = ZCylinder(name='x min y min', x0=-r_cylinder+r_corner, - y0=-r_cylinder+r_corner, R=r_corner, - boundary_type=boundary_type) - x_min_y_max = ZCylinder(name='x min y max', x0=-r_cylinder+r_corner, - y0= r_cylinder-r_corner, R=r_corner, - boundary_type=boundary_type) - x_max_y_min = ZCylinder(name='x max y min', x0= r_cylinder-r_corner, - y0=-r_cylinder+r_corner, R=r_corner, - boundary_type=boundary_type) - x_max_y_max = ZCylinder(name='x max y max', x0= r_cylinder-r_corner, - y0= r_cylinder-r_corner, R=r_corner, - boundary_type=boundary_type) - corners = -x_min_y_min | -x_min_y_max | -x_max_y_min | -x_max_y_max + prism = prism & ~corners - prism_1 = get_rectangular_prism(width=2*r_cylinder, - height=2*(r_cylinder - r_corner), - axis=axis, - origin=origin, - boundary_type=boundary_type) - origin_2 = list(origin) - origin_2[1] += r_cylinder - r_corner/2. - prism_2 = get_rectangular_prism(width=2*(r_cylinder - r_corner), - height=r_corner, - axis=axis, - origin=origin_2, - boundary_type=boundary_type) - origin_3 = list(origin) - origin_3[1] -= r_cylinder - r_corner/2. - prism_3 = get_rectangular_prism(width=2*(r_cylinder - r_corner), - height=r_corner, - axis=axis, - origin=origin_3, - boundary_type=boundary_type) - - sqc = corners | prism_1 | prism_2 | prism_3 - - return sqc + return prism