diff --git a/examples/python/boxes/build-xml.py b/examples/python/boxes/build-xml.py index 4bac9fff4d..9c28d37bb3 100644 --- a/examples/python/boxes/build-xml.py +++ b/examples/python/boxes/build-xml.py @@ -1,3 +1,5 @@ +import numpy as np + import openmc ############################################################################### @@ -115,7 +117,7 @@ settings_file = openmc.SettingsFile() settings_file.batches = batches settings_file.inactive = inactive settings_file.particles = particles -settings_file.set_source_space('point', [0., 0., 0.]) +settings_file.set_source_space('box', np.concatenate(outer_cube.bounding_box)) settings_file.export_to_xml() ############################################################################### diff --git a/examples/python/reflective/build-xml.py b/examples/python/reflective/build-xml.py index 9eab4aec32..44b544d20d 100644 --- a/examples/python/reflective/build-xml.py +++ b/examples/python/reflective/build-xml.py @@ -1,3 +1,5 @@ +import numpy as np + import openmc ############################################################################### @@ -82,5 +84,5 @@ settings_file = openmc.SettingsFile() settings_file.batches = batches settings_file.inactive = inactive settings_file.particles = particles -settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1]) +settings_file.set_source_space('box', np.concatenate(cell.region.bounding_box)) settings_file.export_to_xml() diff --git a/examples/xml/boxes/settings.xml b/examples/xml/boxes/settings.xml index 0ac26ec4d4..eff7c1c105 100644 --- a/examples/xml/boxes/settings.xml +++ b/examples/xml/boxes/settings.xml @@ -10,7 +10,7 @@ - + diff --git a/openmc/region.py b/openmc/region.py index 936e6e5121..b7cfbca7b8 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -1,6 +1,8 @@ from abc import ABCMeta, abstractmethod from collections import Iterable +import numpy as np + from openmc.checkvalue import check_type @@ -218,6 +220,8 @@ class Intersection(Region): ---------- nodes : tuple of Region Regions to take the intersection of + bounding_box : tuple of numpy.array + Lower-left and upper-right coordinates of an axis-aligned bounding box """ @@ -231,6 +235,16 @@ class Intersection(Region): def nodes(self): return self._nodes + @property + def bounding_box(self): + lower_left = np.array([-np.inf, -np.inf, -np.inf]) + upper_right = np.array([np.inf, np.inf, np.inf]) + for n in self.nodes: + lower_left_n, upper_right_n = n.bounding_box + lower_left[:] = np.maximum(lower_left, lower_left_n) + upper_right[:] = np.minimum(upper_right, upper_right_n) + return lower_left, upper_right + @nodes.setter def nodes(self, nodes): check_type('nodes', nodes, Iterable, Region) @@ -257,6 +271,8 @@ class Union(Region): ---------- nodes : tuple of Region Regions to take the union of + bounding_box : tuple of numpy.array + Lower-left and upper-right coordinates of an axis-aligned bounding box """ @@ -270,6 +286,16 @@ class Union(Region): def nodes(self): return self._nodes + @property + def bounding_box(self): + lower_left = np.array([np.inf, np.inf, np.inf]) + upper_right = np.array([-np.inf, -np.inf, -np.inf]) + for n in self.nodes: + lower_left_n, upper_right_n = n.bounding_box + lower_left[:] = np.minimum(lower_left, lower_left_n) + upper_right[:] = np.maximum(upper_right, upper_right_n) + return lower_left, upper_right + @nodes.setter def nodes(self, nodes): check_type('nodes', nodes, Iterable, Region) @@ -300,6 +326,8 @@ class Complement(Region): ---------- node : Region Regions to take the complement of + bounding_box : tuple of numpy.array + Lower-left and upper-right coordinates of an axis-aligned bounding box """ @@ -317,3 +345,16 @@ class Complement(Region): def node(self, node): check_type('node', node, Region) self._node = node + + @property + def bounding_box(self): + # Use De Morgan's laws to distribute the complement operator so that it + # only applies to surface half-spaces, thus allowing us to calculate the + # bounding box in the usual recursive manner. + if isinstance(self.node, Union): + temp_region = Intersection(*[~n for n in self.node.nodes]) + elif isinstance(self.node, Intersection): + temp_region = Union(*[~n for n in self.node.nodes]) + else: + temp_region = ~n + return temp_region.bounding_box diff --git a/openmc/surface.py b/openmc/surface.py index 279246b029..8dc45209be 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -3,6 +3,8 @@ from numbers import Real, Integral from xml.etree import ElementTree as ET import sys +import numpy as np + from openmc.checkvalue import check_type, check_value, check_greater_than from openmc.region import Region @@ -136,6 +138,33 @@ class Surface(object): check_value('boundary type', boundary_type, _BC_TYPES) self._boundary_type = boundary_type + def bounding_box(self, side): + """Determine an axis-aligned bounding box. + + An axis-aligned bounding box for surface half-spaces is represented by + its lower-left and upper-right coordinates. If the half-space is + unbounded in a particular direction, numpy.inf is used to represent + infinity. + + Parameters + ---------- + side : {'+', '-'} + Indicates the negative or positive half-space + + Returns + ------- + numpy.array + Lower-left coordinates of the axis-aligned bounding box for the + desired half-space + numpy.array + Upper-right coordinates of the axis-aligned bounding box for the + desired half-space + + """ + + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + def create_xml_subelement(self): element = ET.Element("surface") element.set("id", str(self._id)) @@ -194,6 +223,10 @@ class Plane(Surface): self._type = 'plane' self._coeff_keys = ['A', 'B', 'C', 'D'] + self._coeffs['A'] = 1. + self._coeffs['B'] = 0. + self._coeffs['C'] = 0. + self._coeffs['D'] = 0. if A is not None: self.a = A @@ -276,6 +309,7 @@ class XPlane(Plane): self._type = 'x-plane' self._coeff_keys = ['x0'] + self._coeffs['x0'] = 0. if x0 is not None: self.x0 = x0 @@ -289,6 +323,37 @@ class XPlane(Plane): check_type('x0 coefficient', x0, Real) self._coeffs['x0'] = x0 + def bounding_box(self, side): + """Determine an axis-aligned bounding box. + + An axis-aligned bounding box for surface half-spaces is represented by + its lower-left and upper-right coordinates. For the x-plane surface, the + half-spaces are unbounded in their y- and z- directions. To represent + infinity, numpy.inf is used. + + Parameters + ---------- + side : {'+', '-'} + Indicates the negative or positive half-space + + Returns + ------- + numpy.array + Lower-left coordinates of the axis-aligned bounding box for the + desired half-space + numpy.array + Upper-right coordinates of the axis-aligned bounding box for the + desired half-space + + """ + + if side == '-': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([self.x0, np.inf, np.inf])) + elif side == '+': + return (np.array([self.x0, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + class YPlane(Plane): """A plane perpendicular to the y axis, i.e. a surface of the form :math:`y - @@ -322,6 +387,7 @@ class YPlane(Plane): self._type = 'y-plane' self._coeff_keys = ['y0'] + self._coeffs['y0'] = 0. if y0 is not None: self.y0 = y0 @@ -335,6 +401,37 @@ class YPlane(Plane): check_type('y0 coefficient', y0, Real) self._coeffs['y0'] = y0 + def bounding_box(self, side): + """Determine an axis-aligned bounding box. + + An axis-aligned bounding box for surface half-spaces is represented by + its lower-left and upper-right coordinates. For the y-plane surface, the + half-spaces are unbounded in their x- and z- directions. To represent + infinity, numpy.inf is used. + + Parameters + ---------- + side : {'+', '-'} + Indicates the negative or positive half-space + + Returns + ------- + numpy.array + Lower-left coordinates of the axis-aligned bounding box for the + desired half-space + numpy.array + Upper-right coordinates of the axis-aligned bounding box for the + desired half-space + + """ + + if side == '-': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, self.y0, np.inf])) + elif side == '+': + return (np.array([-np.inf, self.y0, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + class ZPlane(Plane): """A plane perpendicular to the z axis, i.e. a surface of the form :math:`z - @@ -368,6 +465,7 @@ class ZPlane(Plane): self._type = 'z-plane' self._coeff_keys = ['z0'] + self._coeffs['z0'] = 0. if z0 is not None: self.z0 = z0 @@ -381,6 +479,37 @@ class ZPlane(Plane): check_type('z0 coefficient', z0, Real) self._coeffs['z0'] = z0 + def bounding_box(self, side): + """Determine an axis-aligned bounding box. + + An axis-aligned bounding box for surface half-spaces is represented by + its lower-left and upper-right coordinates. For the z-plane surface, the + half-spaces are unbounded in their x- and y- directions. To represent + infinity, numpy.inf is used. + + Parameters + ---------- + side : {'+', '-'} + Indicates the negative or positive half-space + + Returns + ------- + numpy.array + Lower-left coordinates of the axis-aligned bounding box for the + desired half-space + numpy.array + Upper-right coordinates of the axis-aligned bounding box for the + desired half-space + + """ + + if side == '-': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, self.z0])) + elif side == '+': + return (np.array([-np.inf, -np.inf, self.z0]), + np.array([np.inf, np.inf, np.inf])) + class Cylinder(Surface): """A cylinder whose length is parallel to the x-, y-, or z-axis. @@ -415,6 +544,7 @@ class Cylinder(Surface): super(Cylinder, self).__init__(surface_id, boundary_type, name=name) self._coeff_keys = ['R'] + self._coeffs['R'] = 1. if R is not None: self.r = R @@ -468,6 +598,8 @@ class XCylinder(Cylinder): self._type = 'x-cylinder' self._coeff_keys = ['y0', 'z0', 'R'] + self._coeffs['y0'] = 0. + self._coeffs['z0'] = 0. if y0 is not None: self.y0 = y0 @@ -493,6 +625,38 @@ class XCylinder(Cylinder): check_type('z0 coefficient', z0, Real) self._coeffs['z0'] = z0 + def bounding_box(self, side): + """Determine an axis-aligned bounding box. + + An axis-aligned bounding box for surface half-spaces is represented by + its lower-left and upper-right coordinates. For the x-cylinder surface, + the negative half-space is unbounded in the x- direction and the + positive half-space is unbounded in all directions. To represent + infinity, numpy.inf is used. + + Parameters + ---------- + side : {'+', '-'} + Indicates the negative or positive half-space + + Returns + ------- + numpy.array + Lower-left coordinates of the axis-aligned bounding box for the + desired half-space + numpy.array + Upper-right coordinates of the axis-aligned bounding box for the + desired half-space + + """ + + if side == '-': + return (np.array([-np.inf, self.y0 - self.r, self.z0 - self.r]), + np.array([np.inf, self.y0 + self.r, self.z0 + self.r])) + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + class YCylinder(Cylinder): """An infinite cylinder whose length is parallel to the y-axis. This is a @@ -533,6 +697,8 @@ class YCylinder(Cylinder): self._type = 'y-cylinder' self._coeff_keys = ['x0', 'z0', 'R'] + self._coeffs['x0'] = 0. + self._coeffs['z0'] = 0. if x0 is not None: self.x0 = x0 @@ -558,6 +724,38 @@ class YCylinder(Cylinder): check_type('z0 coefficient', z0, Real) self._coeffs['z0'] = z0 + def bounding_box(self, side): + """Determine an axis-aligned bounding box. + + An axis-aligned bounding box for surface half-spaces is represented by + its lower-left and upper-right coordinates. For the y-cylinder surface, + the negative half-space is unbounded in the y- direction and the + positive half-space is unbounded in all directions. To represent + infinity, numpy.inf is used. + + Parameters + ---------- + side : {'+', '-'} + Indicates the negative or positive half-space + + Returns + ------- + numpy.array + Lower-left coordinates of the axis-aligned bounding box for the + desired half-space + numpy.array + Upper-right coordinates of the axis-aligned bounding box for the + desired half-space + + """ + + if side == '-': + return (np.array([self.x0 - self.r, -np.inf, self.z0 - self.r]), + np.array([self.x0 + self.r, np.inf, self.z0 + self.r])) + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + class ZCylinder(Cylinder): """An infinite cylinder whose length is parallel to the z-axis. This is a @@ -598,6 +796,8 @@ class ZCylinder(Cylinder): self._type = 'z-cylinder' self._coeff_keys = ['x0', 'y0', 'R'] + self._coeffs['x0'] = 0. + self._coeffs['y0'] = 0. if x0 is not None: self.x0 = x0 @@ -623,6 +823,38 @@ class ZCylinder(Cylinder): check_type('y0 coefficient', y0, Real) self._coeffs['y0'] = y0 + def bounding_box(self, side): + """Determine an axis-aligned bounding box. + + An axis-aligned bounding box for surface half-spaces is represented by + its lower-left and upper-right coordinates. For the z-cylinder surface, + the negative half-space is unbounded in the z- direction and the + positive half-space is unbounded in all directions. To represent + infinity, numpy.inf is used. + + Parameters + ---------- + side : {'+', '-'} + Indicates the negative or positive half-space + + Returns + ------- + numpy.array + Lower-left coordinates of the axis-aligned bounding box for the + desired half-space + numpy.array + Upper-right coordinates of the axis-aligned bounding box for the + desired half-space + + """ + + if side == '-': + return (np.array([self.x0 - self.r, self.y0 - self.r, -np.inf]), + np.array([self.x0 + self.r, self.y0 + self.r, np.inf])) + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + class Sphere(Surface): """A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = R^2`. @@ -667,6 +899,10 @@ class Sphere(Surface): self._type = 'sphere' self._coeff_keys = ['x0', 'y0', 'z0', 'R'] + self._coeffs['x0'] = 0. + self._coeffs['y0'] = 0. + self._coeffs['z0'] = 0. + self._coeffs['R'] = 1. if x0 is not None: self.x0 = x0 @@ -716,6 +952,39 @@ class Sphere(Surface): check_type('R coefficient', R, Real) self._coeffs['R'] = R + def bounding_box(self, side): + """Determine an axis-aligned bounding box. + + An axis-aligned bounding box for surface half-spaces is represented by + its lower-left and upper-right coordinates. The positive half-space of a + sphere is unbounded in all directions. To represent infinity, numpy.inf + is used. + + Parameters + ---------- + side : {'+', '-'} + Indicates the negative or positive half-space + + Returns + ------- + numpy.array + Lower-left coordinates of the axis-aligned bounding box for the + desired half-space + numpy.array + Upper-right coordinates of the axis-aligned bounding box for the + desired half-space + + """ + + if side == '-': + return (np.array([self.x0 - self.r, self.y0 - self.r, + self.z0 - self.r]), + np.array([self.x0 + self.r, self.y0 + self.r, + self.z0 + self.r])) + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + class Cone(Surface): """A conical surface parallel to the x-, y-, or z-axis. @@ -761,6 +1030,10 @@ class Cone(Surface): super(Cone, self).__init__(surface_id, boundary_type, name=name) self._coeff_keys = ['x0', 'y0', 'z0', 'R2'] + self._coeffs['x0'] = 0. + self._coeffs['y0'] = 0. + self._coeffs['z0'] = 0. + self._coeffs['R2'] = 1. if x0 is not None: self.x0 = x0 @@ -982,6 +1255,8 @@ class Quadric(Surface): self._type = 'quadric' self._coeff_keys = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k'] + for key in self._coeff_keys: + self._coeffs[key] = 0. if a is not None: self.a = a @@ -1127,6 +1402,8 @@ class Halfspace(Region): Surface which divides Euclidean space. side : {'+', '-'} Indicates whether the positive or negative half-space is used. + bounding_box : tuple of numpy.array + Lower-left and upper-right coordinates of an axis-aligned bounding box """ @@ -1155,6 +1432,10 @@ class Halfspace(Region): check_value('side', side, ('+', '-')) self._side = side + @property + def bounding_box(self): + return self.surface.bounding_box(self.side) + def __str__(self): return '-' + str(self.surface.id) if self.side == '-' \ else str(self.surface.id)