diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index e6cf1e1cc..f1633f3f9 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -91,18 +91,6 @@ Many of the above classes are derived from several abstract classes: openmc.Region openmc.Lattice -Two helper function are also available to create rectangular and hexagonal -prisms defined by the intersection of four and six surface half-spaces, -respectively. - -.. autosummary:: - :toctree: generated - :nosignatures: - :template: myfunction.rst - - openmc.get_hexagonal_prism - openmc.get_rectangular_prism - .. _pythonapi_tallies: Constructing Tallies diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 55a5d17b1..4ba247468 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -5,11 +5,18 @@ Convenience Functions --------------------- +Several helper functions are available here. Ther first two create rectangular +and hexagonal prisms defined by the intersection of four and six surface +half-spaces, respectively. The last function takes a sequence of surfaces and +returns the regions that separate them. + .. autosummary:: :toctree: generated :nosignatures: :template: myfunction.rst + openmc.model.get_hexagonal_prism + openmc.model.get_rectangular_prism openmc.model.subdivide TRISO Fuel Modeling diff --git a/openmc/__init__.py b/openmc/__init__.py index d692ebbae..8fb9bcf37 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -28,4 +28,7 @@ from openmc.mixin import * from openmc.plotter import * from openmc.search import * +# Import a few convencience functions that used to be here +from openmc.model import get_rectangular_prism, get_hexagonal_prism + __version__ = '0.9.0' diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py index d9755b0e8..589765d2c 100644 --- a/openmc/model/funcs.py +++ b/openmc/model/funcs.py @@ -1,3 +1,257 @@ +from __future__ import division +from collections import Iterable, OrderedDict +from math import sqrt +from numbers import Real + +from openmc import XPlane, YPlane, Plane, ZCylinder +from openmc.checkvalue import check_type, check_value + + +def get_rectangular_prism(width, height, axis='z', origin=(0., 0.), + boundary_type='transmission', corner_radius=0.): + """Get an infinite rectangular prism from four planar surfaces. + + Parameters + ---------- + width: float + Prism width in units of cm. The width is aligned with the y, x, + or x axes for prisms parallel to the x, y, or z axis, respectively. + 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. + axis : {'x', 'y', 'z'} + Axis with which the infinite length of the prism should be aligned. + Defaults to 'z'. + origin: Iterable of two floats + Origin of the prism. The two floats correspond to (y,z), (x,z) or + (x,y) for prisms 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 rectangular prism (default is 'transmission'). + corner_radius: float + Prism corner radius in units of cm. Defaults to 0. + + Returns + ------- + openmc.Region + The inside of a rectangular prism + + """ + + 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) + + # 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': + x1, x2 = 'x', 'z' + else: + x1, x2 = 'x', 'y' + + # Get cylinder class corresponding to given axis + cyl = globals()['{}Cylinder'.format(axis.upper())] + + # 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]) + prism = +min_x1 & -max_x1 & +min_x2 & -max_x2 + + # Handle rounded corners if given + if corner_radius > 0.: + args = {'R': corner_radius, 'boundary_type': boundary_type} + + args[x1 + '0'] = origin[0] - width/2 + corner_radius + args[x2 + '0'] = origin[1] - height/2 + corner_radius + x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args) + + args[x1 + '0'] = origin[0] - width/2 + corner_radius + args[x2 + '0'] = origin[1] - height/2 + corner_radius + x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args) + + args[x1 + '0'] = origin[0] - width/2 + corner_radius + args[x2 + '0'] = origin[1] + height/2 - corner_radius + x1_min_x2_max = cyl(name='{} min {} max'.format(x1, x2), **args) + + args[x1 + '0'] = origin[0] + width/2 - corner_radius + args[x2 + '0'] = origin[1] - height/2 + corner_radius + x1_max_x2_min = cyl(name='{} max {} min'.format(x1, x2), **args) + + args[x1 + '0'] = origin[0] + width/2 - corner_radius + 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) + + corners = (+x1_min_x2_min & -x1_min & -x2_min) | \ + (+x1_min_x2_max & -x1_min & +x2_max) | \ + (+x1_max_x2_min & +x1_max & -x2_min) | \ + (+x1_max_x2_max & +x1_max & +x2_max) + + prism = prism & ~corners + + return prism + + +def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.), + boundary_type='transmission', corner_radius=0.): + """Create a hexagon region from six surface planes. + + Parameters + ---------- + edge_length : float + Length of a side of the hexagon in cm + 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 + parallel to the y-axis. + origin: Iterable of two floats + Origin of the prism. Defaults to (0., 0.). + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + Boundary condition that defines the behavior for particles hitting the + surfaces comprising the hexagonal prism (default is 'transmission'). + corner_radius: float + Prism corner radius in units of cm. Defaults to 0. + + Returns + ------- + openmc.Region + The inside of a hexagonal prism + + """ + + l = edge_length + x, y = origin + + if orientation == 'y': + right = XPlane(x0=x + sqrt(3.)/2*l, boundary_type=boundary_type) + left = XPlane(x0=x - sqrt(3.)/2*l, boundary_type=boundary_type) + c = sqrt(3.)/3. + + # y = -x/sqrt(3) + a + upper_right = Plane(A=c, B=1., D=l+x*c+y, boundary_type=boundary_type) + + # y = x/sqrt(3) + a + upper_left = Plane(A=-c, B=1., D=l-x*c+y, boundary_type=boundary_type) + + # y = x/sqrt(3) - a + lower_right = Plane(A=-c, B=1., D=-l-x*c+y, boundary_type=boundary_type) + + # y = -x/sqrt(3) - a + lower_left = Plane(A=c, B=1., D=-l+x*c+y, boundary_type=boundary_type) + + prism = -right & +left & -upper_right & -upper_left & \ + +lower_right & +lower_left + + if boundary_type == 'periodic': + right.periodic_surface = left + upper_right.periodic_surface = lower_left + 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) + c = sqrt(3.) + + # y = -sqrt(3)*(x - a) + upper_right = Plane(A=c, B=1., D=c*l+x*c+y, boundary_type=boundary_type) + + # y = sqrt(3)*(x + a) + lower_right = Plane(A=-c, B=1., D=-c*l-x*c+y, + boundary_type=boundary_type) + + # y = -sqrt(3)*(x + a) + lower_left = Plane(A=c, B=1., D=-c*l+x*c+y, boundary_type=boundary_type) + + # y = sqrt(3)*(x + a) + upper_left = Plane(A=-c, B=1., D=c*l-x*c+y, boundary_type=boundary_type) + + prism = -top & +bottom & -upper_right & +lower_right & \ + +lower_left & -upper_left + + if boundary_type == 'periodic': + top.periodic_surface = bottom + upper_right.periodic_surface = lower_left + lower_right.periodic_surface = upper_left + + # Handle rounded corners if given + if corner_radius > 0.: + if boundary_type == 'periodic': + raise ValueError('Periodic boundary conditions not permitted when ' + 'rounded corners are used.') + + c = sqrt(3.)/2 + 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), + boundary_type=boundary_type) + + if orientation == 'x': + x_min_y_min_in = cyl1(name='x min y min in', x0=x-t/2, y0=y-c*t) + x_min_y_max_in = cyl1(name='x min y max in', x0=x+t/2, y0=y-c*t) + x_max_y_min_in = cyl1(name='x max y min in', x0=x-t/2, y0=y+c*t) + x_max_y_max_in = cyl1(name='x max y max in', x0=x+t/2, y0=y+c*t) + x_min_in = cyl1(name='x min in', x0=x-t, y0=y) + x_max_in = cyl1(name='x max in', x0=x+t, y0=y) + + x_min_y_min_out = cyl2(name='x min y min out', x0=x-l/2, y0=y-c*l) + x_min_y_max_out = cyl2(name='x min y max out', x0=x+l/2, y0=y-c*l) + x_max_y_min_out = cyl2(name='x max y min out', x0=x-l/2, y0=y+c*l) + x_max_y_max_out = cyl2(name='x max y max out', x0=x+l/2, y0=y+c*l) + x_min_out = cyl2(name='x min out', x0=x-l, y0=y) + x_max_out = cyl2(name='x max out', x0=x+l, y0=y) + + 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) + + elif orientation == 'y': + x_min_y_min_in = cyl1(name='x min y min in', x0=x-c*t, y0=y-t/2) + x_min_y_max_in = cyl1(name='x min y max in', x0=x-c*t, y0=y+t/2) + x_max_y_min_in = cyl1(name='x max y min in', x0=x+c*t, y0=y-t/2) + x_max_y_max_in = cyl1(name='x max y max in', x0=x+c*t, y0=y+t/2) + y_min_in = cyl1(name='y min in', x0=x, y0=y-t) + y_max_in = cyl1(name='y max in', x0=x, y0=y+t) + + x_min_y_min_out = cyl2(name='x min y min out', x0=x-c*l, y0=y-l/2) + x_min_y_max_out = cyl2(name='x min y max out', x0=x-c*l, y0=y+l/2) + x_max_y_min_out = cyl2(name='x max y min out', x0=x+c*l, y0=y-l/2) + x_max_y_max_out = cyl2(name='x max y max out', x0=x+c*l, y0=y+l/2) + y_min_out = cyl2(name='y min out', x0=x, y0=y-l) + y_max_out = cyl2(name='y max out', x0=x, y0=y+l) + + 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) + + prism = prism & ~corners + + return prism + + def subdivide(surfaces): """Create regions separated by a series of surfaces. diff --git a/openmc/surface.py b/openmc/surface.py index b415ecb93..961ac4d45 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1,23 +1,19 @@ from __future__ import division from abc import ABCMeta -from collections import Iterable, OrderedDict +from collections import OrderedDict from copy import deepcopy from functools import partial from numbers import Real, Integral from xml.etree import ElementTree as ET -from math import sqrt from six import add_metaclass, string_types import numpy as np -from openmc.checkvalue import check_type, check_value, check_greater_than +from openmc.checkvalue import check_type, check_value from openmc.region import Region, Intersection, Union from openmc.mixin import IDManagerMixin -# A static variable for auto-generated Surface IDs -AUTO_SURFACE_ID = 10000 - _BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic'] @@ -1944,248 +1940,3 @@ class Halfspace(Region): clone = deepcopy(self) clone.surface = self.surface.clone(memo) return clone - - -def get_rectangular_prism(width, height, axis='z', origin=(0., 0.), - boundary_type='transmission', corner_radius=0.): - """Get an infinite rectangular prism from four planar surfaces. - - Parameters - ---------- - width: float - Prism width in units of cm. The width is aligned with the y, x, - or x axes for prisms parallel to the x, y, or z axis, respectively. - 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. - axis : {'x', 'y', 'z'} - Axis with which the infinite length of the prism should be aligned. - Defaults to 'z'. - origin: Iterable of two floats - Origin of the prism. The two floats correspond to (y,z), (x,z) or - (x,y) for prisms 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 rectangular prism (default is 'transmission'). - corner_radius: float - Prism corner radius in units of cm. Defaults to 0. - - Returns - ------- - openmc.Region - The inside of a rectangular prism - - """ - - 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) - - # 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': - x1, x2 = 'x', 'z' - else: - x1, x2 = 'x', 'y' - - # Get cylinder class corresponding to given axis - cyl = globals()['{}Cylinder'.format(axis.upper())] - - # 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]) - prism = +min_x1 & -max_x1 & +min_x2 & -max_x2 - - # Handle rounded corners if given - if corner_radius > 0.: - args = {'R': corner_radius, 'boundary_type': boundary_type} - - args[x1 + '0'] = origin[0] - width/2 + corner_radius - args[x2 + '0'] = origin[1] - height/2 + corner_radius - x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args) - - args[x1 + '0'] = origin[0] - width/2 + corner_radius - args[x2 + '0'] = origin[1] - height/2 + corner_radius - x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args) - - args[x1 + '0'] = origin[0] - width/2 + corner_radius - args[x2 + '0'] = origin[1] + height/2 - corner_radius - x1_min_x2_max = cyl(name='{} min {} max'.format(x1, x2), **args) - - args[x1 + '0'] = origin[0] + width/2 - corner_radius - args[x2 + '0'] = origin[1] - height/2 + corner_radius - x1_max_x2_min = cyl(name='{} max {} min'.format(x1, x2), **args) - - args[x1 + '0'] = origin[0] + width/2 - corner_radius - 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) - - corners = (+x1_min_x2_min & -x1_min & -x2_min) | \ - (+x1_min_x2_max & -x1_min & +x2_max) | \ - (+x1_max_x2_min & +x1_max & -x2_min) | \ - (+x1_max_x2_max & +x1_max & +x2_max) - - prism = prism & ~corners - - return prism - - -def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.), - boundary_type='transmission', corner_radius=0.): - """Create a hexagon region from six surface planes. - - Parameters - ---------- - edge_length : float - Length of a side of the hexagon in cm - 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 - parallel to the y-axis. - origin: Iterable of two floats - Origin of the prism. Defaults to (0., 0.). - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} - Boundary condition that defines the behavior for particles hitting the - surfaces comprising the hexagonal prism (default is 'transmission'). - corner_radius: float - Prism corner radius in units of cm. Defaults to 0. - - Returns - ------- - openmc.Region - The inside of a hexagonal prism - - """ - - l = edge_length - x, y = origin - - if orientation == 'y': - right = XPlane(x0=x + sqrt(3.)/2*l, boundary_type=boundary_type) - left = XPlane(x0=x - sqrt(3.)/2*l, boundary_type=boundary_type) - c = sqrt(3.)/3. - - # y = -x/sqrt(3) + a - upper_right = Plane(A=c, B=1., D=l+x*c+y, boundary_type=boundary_type) - - # y = x/sqrt(3) + a - upper_left = Plane(A=-c, B=1., D=l-x*c+y, boundary_type=boundary_type) - - # y = x/sqrt(3) - a - lower_right = Plane(A=-c, B=1., D=-l-x*c+y, boundary_type=boundary_type) - - # y = -x/sqrt(3) - a - lower_left = Plane(A=c, B=1., D=-l+x*c+y, boundary_type=boundary_type) - - prism = -right & +left & -upper_right & -upper_left & \ - +lower_right & +lower_left - - if boundary_type == 'periodic': - right.periodic_surface = left - upper_right.periodic_surface = lower_left - 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) - c = sqrt(3.) - - # y = -sqrt(3)*(x - a) - upper_right = Plane(A=c, B=1., D=c*l+x*c+y, boundary_type=boundary_type) - - # y = sqrt(3)*(x + a) - lower_right = Plane(A=-c, B=1., D=-c*l-x*c+y, - boundary_type=boundary_type) - - # y = -sqrt(3)*(x + a) - lower_left = Plane(A=c, B=1., D=-c*l+x*c+y, boundary_type=boundary_type) - - # y = sqrt(3)*(x + a) - upper_left = Plane(A=-c, B=1., D=c*l-x*c+y, boundary_type=boundary_type) - - prism = -top & +bottom & -upper_right & +lower_right & \ - +lower_left & -upper_left - - if boundary_type == 'periodic': - top.periodic_surface = bottom - upper_right.periodic_surface = lower_left - lower_right.periodic_surface = upper_left - - # Handle rounded corners if given - if corner_radius > 0.: - if boundary_type == 'periodic': - raise ValueError('Periodic boundary conditions not permitted when ' - 'rounded corners are used.') - - c = sqrt(3.)/2 - 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), - boundary_type=boundary_type) - - if orientation == 'x': - x_min_y_min_in = cyl1(name='x min y min in', x0=x-t/2, y0=y-c*t) - x_min_y_max_in = cyl1(name='x min y max in', x0=x+t/2, y0=y-c*t) - x_max_y_min_in = cyl1(name='x max y min in', x0=x-t/2, y0=y+c*t) - x_max_y_max_in = cyl1(name='x max y max in', x0=x+t/2, y0=y+c*t) - x_min_in = cyl1(name='x min in', x0=x-t, y0=y) - x_max_in = cyl1(name='x max in', x0=x+t, y0=y) - - x_min_y_min_out = cyl2(name='x min y min out', x0=x-l/2, y0=y-c*l) - x_min_y_max_out = cyl2(name='x min y max out', x0=x+l/2, y0=y-c*l) - x_max_y_min_out = cyl2(name='x max y min out', x0=x-l/2, y0=y+c*l) - x_max_y_max_out = cyl2(name='x max y max out', x0=x+l/2, y0=y+c*l) - x_min_out = cyl2(name='x min out', x0=x-l, y0=y) - x_max_out = cyl2(name='x max out', x0=x+l, y0=y) - - 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) - - elif orientation == 'y': - x_min_y_min_in = cyl1(name='x min y min in', x0=x-c*t, y0=y-t/2) - x_min_y_max_in = cyl1(name='x min y max in', x0=x-c*t, y0=y+t/2) - x_max_y_min_in = cyl1(name='x max y min in', x0=x+c*t, y0=y-t/2) - x_max_y_max_in = cyl1(name='x max y max in', x0=x+c*t, y0=y+t/2) - y_min_in = cyl1(name='y min in', x0=x, y0=y-t) - y_max_in = cyl1(name='y max in', x0=x, y0=y+t) - - x_min_y_min_out = cyl2(name='x min y min out', x0=x-c*l, y0=y-l/2) - x_min_y_max_out = cyl2(name='x min y max out', x0=x-c*l, y0=y+l/2) - x_max_y_min_out = cyl2(name='x max y min out', x0=x+c*l, y0=y-l/2) - x_max_y_max_out = cyl2(name='x max y max out', x0=x+c*l, y0=y+l/2) - y_min_out = cyl2(name='y min out', x0=x, y0=y-l) - y_max_out = cyl2(name='y max out', x0=x, y0=y+l) - - 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) - - prism = prism & ~corners - - return prism