diff --git a/openmc/region.py b/openmc/region.py index 936e6e5121..7b640a07f4 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): + ll = np.array([-np.inf, -np.inf, -np.inf]) + ur = np.array([np.inf, np.inf, np.inf]) + for n in self.nodes: + ll_n, ur_n = n.bounding_box + ll[:] = np.maximum(ll, ll_n) + ur[:] = np.minimum(ur, ur_n) + return ll, ur + @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): + ll = np.array([np.inf, np.inf, np.inf]) + ur = np.array([-np.inf, -np.inf, -np.inf]) + for n in self.nodes: + ll_n, ur_n = n.bounding_box + ll[:] = np.minimum(ll, ll_n) + ur[:] = np.maximum(ur, ur_n) + return ll, ur + @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..83465ac6d0 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. 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 + + """ + + 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. 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 + + """ + + 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. 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 + + """ + + 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,37 @@ 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. 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 + + """ + + if side == '-': + return (np.array([-np.inf, self.y0 - self.r, self.z0 - self.r]), + np.array([np.inf, self.y0 + self.r, self.y0 + 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 +696,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 +723,37 @@ 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. 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 + + """ + + if side == '-': + return (np.array([self.x0 - self.r, -np.inf, self.z0 - self.r]), + np.array([self.x0 + self.r, np.inf, self.y0 + 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 +794,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 +821,37 @@ 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. 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 + + """ + + 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 +896,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 +949,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. 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 + + """ + + 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 +1027,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 +1252,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 +1399,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 +1429,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)