From 07be6306ce92bef7de529e592d0440ad226bddd1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 May 2016 16:37:11 -0500 Subject: [PATCH 1/9] Add Surface.evaluate() methods and Region.__contains__ methods --- openmc/region.py | 52 ++++++++++ openmc/surface.py | 246 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 298 insertions(+) diff --git a/openmc/region.py b/openmc/region.py index a2edbeedd..95f59546c 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -28,6 +28,10 @@ class Region(object): def __invert__(self): return Complement(self) + @abstractmethod + def __contains__(self, point): + return False + @abstractmethod def __str__(self): return '' @@ -229,6 +233,22 @@ class Intersection(Region): def __init__(self, *nodes): self.nodes = list(nodes) + def __contains__(self, point): + """Check whether a point is contained in the region. + + Parameters + ---------- + point : 3-tuple of float + Cartesian coordinates, :math:`(x',y',z')`, of the point + + Returns + ------- + bool + Whether the point is in the region + + """ + return all(point in n for n in self.nodes) + def __str__(self): return '(' + ' '.join(map(str, self.nodes)) + ')' @@ -281,6 +301,22 @@ class Union(Region): def __init__(self, *nodes): self.nodes = list(nodes) + def __contains__(self, point): + """Check whether a point is contained in the region. + + Parameters + ---------- + point : 3-tuple of float + Cartesian coordinates, :math:`(x',y',z')`, of the point + + Returns + ------- + bool + Whether the point is in the region + + """ + return any(point in n for n in self.nodes) + def __str__(self): return '(' + ' | '.join(map(str, self.nodes)) + ')' @@ -336,6 +372,22 @@ class Complement(Region): def __init__(self, node): self.node = node + def __contains__(self, point): + """Check whether a point is contained in the region. + + Parameters + ---------- + point : 3-tuple of float + Cartesian coordinates, :math:`(x',y',z')`, of the point + + Returns + ------- + bool + Whether the point is in the region + + """ + return point not in self.node + def __str__(self): return '~' + str(self.node) diff --git a/openmc/surface.py b/openmc/surface.py index 193780192..76f0d82e7 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -295,6 +295,25 @@ class Plane(Surface): self._periodic_surface = periodic_surface periodic_surface._periodic_surface = self + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`Ax' + By' + Cz' - d` + + """ + + x, y, z = point + return self.a*x + self.b*y + self.c*z - self.d + def create_xml_subelement(self): element = super(Plane, self).create_xml_subelement() @@ -392,6 +411,23 @@ class XPlane(Plane): return (np.array([self.x0, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`x' - x_0` + + """ + return point[0] - self.x0 + class YPlane(Plane): """A plane perpendicular to the y axis of the form :math:`y - y_0 = 0` @@ -481,6 +517,23 @@ class YPlane(Plane): return (np.array([-np.inf, self.y0, -np.inf]), np.array([np.inf, np.inf, np.inf])) + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`y' - y_0` + + """ + return point[1] - self.y0 + class ZPlane(Plane): """A plane perpendicular to the z axis of the form :math:`z - z_0 = 0` @@ -570,6 +623,23 @@ class ZPlane(Plane): return (np.array([-np.inf, -np.inf, self.z0]), np.array([np.inf, np.inf, np.inf])) + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`z' - z_0` + + """ + return point[2] - self.z0 + class Cylinder(Surface): """A cylinder whose length is parallel to the x-, y-, or z-axis. @@ -728,6 +798,25 @@ class XCylinder(Cylinder): return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`(y' - y_0)^2 + (z' - z_0)^2 - R^2` + + """ + y = point[1] - self.y0 + z = point[2] - self.z0 + return y**2 + z**2 - self.r**2 + class YCylinder(Cylinder): """An infinite cylinder whose length is parallel to the y-axis of the form @@ -831,6 +920,25 @@ class YCylinder(Cylinder): return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`(x' - x_0)^2 + (z' - z_0)^2 - R^2` + + """ + x = point[0] - self.x0 + z = point[2] - self.z0 + return x**2 + z**2 - self.r**2 + class ZCylinder(Cylinder): """An infinite cylinder whose length is parallel to the z-axis of the form @@ -934,6 +1042,25 @@ class ZCylinder(Cylinder): return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`(x' - x_0)^2 + (y' - y_0)^2 - R^2` + + """ + x = point[0] - self.x0 + y = point[1] - self.y0 + return x**2 + y**2 - self.r**2 + class Sphere(Surface): """A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = R^2`. @@ -1062,6 +1189,26 @@ class Sphere(Surface): return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, np.inf])) + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`(x' - x_0)^2 + (y' - y_0)^2 + (z' - z_0)^2 - R^2` + + """ + x = point[0] - self.x0 + y = point[1] - self.y0 + z = point[2] - self.z0 + return x**2 + y**2 + z**2 - self.r**2 + class Cone(Surface): """A conical surface parallel to the x-, y-, or z-axis. @@ -1214,6 +1361,26 @@ class XCone(Cone): self._type = 'x-cone' + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`(y' - y_0)^2 + (z' - z_0)^2 - R^2(x' - x_0)^2` + + """ + x = point[0] - self.x0 + y = point[1] - self.y0 + z = point[2] - self.z0 + return y**2 + z**2 - self.r2*x**2 + class YCone(Cone): """A cone parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 = @@ -1270,6 +1437,26 @@ class YCone(Cone): self._type = 'y-cone' + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`(x' - x_0)^2 + (z' - z_0)^2 - R^2(y' - y_0)^2` + + """ + x = point[0] - self.x0 + y = point[1] - self.y0 + z = point[2] - self.z0 + return x**2 + z**2 - self.r2*y**2 + class ZCone(Cone): """A cone parallel to the x-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 = @@ -1326,6 +1513,26 @@ class ZCone(Cone): self._type = 'z-cone' + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`(x' - x_0)^2 + (y' - y_0)^2 - R^2(z' - z_0)^2` + + """ + x = point[0] - self.x0 + y = point[1] - self.y0 + z = point[2] - self.z0 + return x**2 + y**2 - self.r2*z**2 + class Quadric(Surface): """A surface of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + @@ -1471,6 +1678,27 @@ class Quadric(Surface): check_type('k coefficient', k, Real) self._coefficients['k'] = k + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`Ax'^2 + By'^2 + Cz'^2 + Dx'y' + Ey'z' + Fx'z' + Gx' + Hy' + + Jz' + K = 0` + + """ + x, y, z = point + return x*(self.a*x + self.d*y + self.g) + \ + y*(self.b*y + self.e*z + self.h) + \ + z*(self.c*z + self.f*x + self.j) + self.k + class Halfspace(Region): """A positive or negative half-space region. @@ -1516,6 +1744,24 @@ class Halfspace(Region): def __invert__(self): return -self.surface if self.side == '+' else +self.surface + def __contains__(self, point): + """Check whether a point is contained in the half-space. + + Parameters + ---------- + point : 3-tuple of float + Cartesian coordinates, :math:`(x',y',z')`, of the point + + Returns + ------- + bool + Whether the point is in the half-space + + """ + + val = self.surface.evaluate(point) + return val >= 0. if self.side == '+' else val < 0. + @property def surface(self): return self._surface From e8d4dbd5f419ce52ff14355b9dbf59ef0a0e26b8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 May 2016 16:39:35 -0500 Subject: [PATCH 2/9] Add Intersection.__iter__ and Union.__iter__ methods --- openmc/region.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/openmc/region.py b/openmc/region.py index 95f59546c..9e1011271 100644 --- a/openmc/region.py +++ b/openmc/region.py @@ -223,7 +223,7 @@ class Intersection(Region): Attributes ---------- - nodes : tuple of openmc.Region + nodes : list of openmc.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 @@ -233,6 +233,10 @@ class Intersection(Region): def __init__(self, *nodes): self.nodes = list(nodes) + def __iter__(self): + for n in self.nodes: + yield n + def __contains__(self, point): """Check whether a point is contained in the region. @@ -301,6 +305,10 @@ class Union(Region): def __init__(self, *nodes): self.nodes = list(nodes) + def __iter__(self): + for n in self.nodes: + yield n + def __contains__(self, point): """Check whether a point is contained in the region. From 168f1269fdcd1ee73cc88090431c6686ab085b09 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 May 2016 10:25:27 -0500 Subject: [PATCH 3/9] Add Cell.__contains__ and Cell.rotation_matrix --- openmc/cell.py | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/openmc/cell.py b/openmc/cell.py index 8ddae6371..8f8ebd936 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -1,9 +1,12 @@ from collections import OrderedDict, Iterable +from math import cos, sin, pi from numbers import Real, Integral from xml.etree import ElementTree as ET import sys import warnings +import numpy as np + import openmc import openmc.checkvalue as cv from openmc.surface import Halfspace @@ -83,6 +86,7 @@ class Cell(object): self._type = None self._region = None self._rotation = None + self._rotation_matrix = None self._translation = None self._offsets = None self._distribcell_index = None @@ -92,6 +96,9 @@ class Cell(object): if region is not None: self.region = region + def __contains__(self, point): + return point in self.region + def __eq__(self, other): if not isinstance(other, Cell): return False @@ -124,6 +131,8 @@ class Cell(object): if isinstance(self._fill, openmc.Material): string += '{0: <16}{1}{2}\n'.format('\tMaterial', '=\t', self._fill._id) + elif isinstance(self._fill, basestring): + string += '{0: <16}=\tvoid\n'.format('\tMaterial') elif isinstance(self._fill, Iterable): string += '{0: <16}{1}'.format('\tMaterial', '=\t') string += '[' @@ -179,6 +188,10 @@ class Cell(object): def rotation(self): return self._rotation + @property + def rotation_matrix(self): + return self._rotation_matrix + @property def translation(self): return self._translation @@ -249,13 +262,23 @@ class Cell(object): cv.check_type('cell rotation', rotation, Iterable, Real) cv.check_length('cell rotation', rotation, 3) - self._rotation = rotation + self._rotation = np.asarray(rotation) + + # Save rotation matrix + phi, theta, psi = self.rotation*(-pi/180.) + c3, s3 = cos(phi), sin(phi) + c2, s2 = cos(theta), sin(theta) + c1, s1 = cos(psi), sin(psi) + self._rotation_matrix = np.array([ + [c1*c2, c1*s2*s3 - c3*s1, s1*s3 + c1*c3*s2], + [c2*s1, c1*c3 + s1*s2*s3, c3*s1*s2 - c1*s3], + [-s2, c2*s3, c2*c3]]) @translation.setter def translation(self, translation): cv.check_type('cell translation', translation, Iterable, Real) cv.check_length('cell translation', translation, 3) - self._translation = translation + self._translation = np.asarray(translation) @offsets.setter def offsets(self, offsets): From e40a369693e48d5f7b5a2e0d3edae5203f938e5a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 19 May 2016 20:14:44 -0500 Subject: [PATCH 4/9] Add geometry plotting capability and restructure lattice attributes --- .../pythonapi/examples/mgxs-part-iii.ipynb | 13 +- .../pythonapi/examples/mgxs-part-iv.ipynb | 1 - .../examples/pandas-dataframes.ipynb | 11 +- examples/python/lattice/nested/build-xml.py | 2 - examples/python/lattice/simple/build-xml.py | 1 - openmc/cell.py | 5 +- openmc/geometry.py | 17 + openmc/lattice.py | 444 +++++++++++++++--- openmc/opencg_compatible.py | 9 +- openmc/summary.py | 13 +- openmc/universe.py | 105 +++++ tests/input_set.py | 4 - .../test_asymmetric_lattice.py | 1 - tests/test_distribmat/test_distribmat.py | 1 - 14 files changed, 527 insertions(+), 100 deletions(-) diff --git a/docs/source/pythonapi/examples/mgxs-part-iii.ipynb b/docs/source/pythonapi/examples/mgxs-part-iii.ipynb index a38677945..5f0acde3f 100644 --- a/docs/source/pythonapi/examples/mgxs-part-iii.ipynb +++ b/docs/source/pythonapi/examples/mgxs-part-iii.ipynb @@ -260,7 +260,6 @@ "source": [ "# Create fuel assembly Lattice\n", "assembly = openmc.RectLattice(name='1.6% Fuel Assembly')\n", - "assembly.dimension = (17, 17)\n", "assembly.pitch = (1.26, 1.26)\n", "assembly.lower_left = [-1.26 * 17. / 2.0] * 2" ] @@ -1597,21 +1596,21 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 2", - "language": "python2", - "name": "python2" + "display_name": "Python 3", + "language": "python", + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.11" + "pygments_lexer": "ipython3", + "version": "3.5.1" } }, "nbformat": 4, diff --git a/docs/source/pythonapi/examples/mgxs-part-iv.ipynb b/docs/source/pythonapi/examples/mgxs-part-iv.ipynb index b330e7ace..65b3f44dc 100644 --- a/docs/source/pythonapi/examples/mgxs-part-iv.ipynb +++ b/docs/source/pythonapi/examples/mgxs-part-iv.ipynb @@ -244,7 +244,6 @@ "source": [ "# Create fuel assembly Lattice\n", "assembly = openmc.RectLattice(name='1.6% Fuel Assembly')\n", - "assembly.dimension = (17, 17)\n", "assembly.pitch = (1.26, 1.26)\n", "assembly.lower_left = [-1.26 * 17. / 2.0] * 2" ] diff --git a/docs/source/pythonapi/examples/pandas-dataframes.ipynb b/docs/source/pythonapi/examples/pandas-dataframes.ipynb index b88cf9949..2c222ad6e 100644 --- a/docs/source/pythonapi/examples/pandas-dataframes.ipynb +++ b/docs/source/pythonapi/examples/pandas-dataframes.ipynb @@ -199,7 +199,6 @@ "source": [ "# Create fuel assembly Lattice\n", "assembly = openmc.RectLattice(name='1.6% Fuel - 0BA')\n", - "assembly.dimension = (17, 17)\n", "assembly.pitch = (1.26, 1.26)\n", "assembly.lower_left = [-1.26 * 17. / 2.0] * 2\n", "assembly.universes = [[pin_cell_universe] * 17] * 17" @@ -2194,21 +2193,21 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.6" + "pygments_lexer": "ipython3", + "version": "3.5.1" } }, "nbformat": 4, diff --git a/examples/python/lattice/nested/build-xml.py b/examples/python/lattice/nested/build-xml.py index b2d611d34..a964d882c 100644 --- a/examples/python/lattice/nested/build-xml.py +++ b/examples/python/lattice/nested/build-xml.py @@ -98,14 +98,12 @@ univ4.add_cell(cell2) # Instantiate nested Lattices lattice1 = openmc.RectLattice(lattice_id=4, name='4x4 assembly') -lattice1.dimension = [2, 2] lattice1.lower_left = [-1., -1.] lattice1.pitch = [1., 1.] lattice1.universes = [[univ1, univ2], [univ2, univ3]] lattice2 = openmc.RectLattice(lattice_id=6, name='4x4 core') -lattice2.dimension = [2, 2] lattice2.lower_left = [-2., -2.] lattice2.pitch = [2., 2.] lattice2.universes = [[univ4, univ4], diff --git a/examples/python/lattice/simple/build-xml.py b/examples/python/lattice/simple/build-xml.py index 65c355479..4961b96b8 100644 --- a/examples/python/lattice/simple/build-xml.py +++ b/examples/python/lattice/simple/build-xml.py @@ -94,7 +94,6 @@ root.add_cell(cell1) # Instantiate a Lattice lattice = openmc.RectLattice(lattice_id=5) -lattice.dimension = [4, 4] lattice.lower_left = [-2., -2.] lattice.pitch = [1., 1.] lattice.universes = [[univ1, univ2, univ1, univ2], diff --git a/openmc/cell.py b/openmc/cell.py index 8f8ebd936..29f85754a 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -97,7 +97,10 @@ class Cell(object): self.region = region def __contains__(self, point): - return point in self.region + if self.region is None: + return True + else: + return point in self.region def __eq__(self, other): if not isinstance(other, Cell): diff --git a/openmc/geometry.py b/openmc/geometry.py index ed437f6e1..7eddadfc7 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -62,6 +62,23 @@ class Geometry(object): tree.write("geometry.xml", xml_declaration=True, encoding='utf-8', method="xml") + def find(self, point): + """Find cells/universes/lattices which contain a given point + + Parameters + ---------- + point : 3-tuple of float + Cartesian coordinatesof the point + + Returns + ------- + list + Sequence of universes, cells, and lattices which are traversed to + find the given point + + """ + return self.root_universe.find(point) + def get_cell_instance(self, path): """Return the instance number for the final cell in a geometry path. diff --git a/openmc/lattice.py b/openmc/lattice.py index af6c14a6a..d8deddc42 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -1,8 +1,12 @@ +from __future__ import division + import abc from collections import OrderedDict, Iterable +from math import sqrt, floor from numbers import Real, Integral from xml.etree import ElementTree as ET import sys +import warnings import numpy as np @@ -113,12 +117,6 @@ class Lattice(object): cv.check_type('outer universe', outer, openmc.Universe) self._outer = outer - @universes.setter - def universes(self, universes): - cv.check_iterable_type('lattice universes', universes, openmc.Universe, - min_depth=2, max_depth=3) - self._universes = np.asarray(universes) - def get_unique_universes(self): """Determine all unique universes in the lattice @@ -239,6 +237,11 @@ class Lattice(object): class RectLattice(Lattice): """A lattice consisting of rectangular prisms. + To completely define a rectangular lattice, the + :attr:`RectLattice.lower_left` :attr:`RectLattice.pitch`, + :attr:`RectLattice.outer`, and :attr:`RectLattice.universes` properties need + to be set. + Parameters ---------- lattice_id : int, optional @@ -253,12 +256,6 @@ class RectLattice(Lattice): Unique identifier for the lattice name : str Name of the lattice - dimension : Iterable of int - An array of two or three integers representing the number of lattice - cells in the x- and y- (and z-) directions, respectively. - lower_left : Iterable of float - The coordinates of the lower-left corner of the lattice. If the lattice - is two-dimensional, only the x- and y-coordinates are specified. pitch : Iterable of float Pitch of the lattice in the x, y, and (if applicable) z directions in cm. @@ -266,7 +263,25 @@ class RectLattice(Lattice): A universe to fill all space outside the lattice universes : Iterable of Iterable of openmc.Universe A two- or three-dimensional list/array of universes filling each element - of the lattice + of the lattice. The first dimension corresponds to the z-direction (if + applicable), the second dimension corresponds to the y-direction, and + the third dimension corresponds to the x-direction. Note that for the + y-direction, a higher index corresponds to a lower physical + y-value. Each z-slice in the array can be thought of as a top-down view + of the lattice. + lower_left : Iterable of float + The Cartesian coordinates of the lower-left corner of the lattice. If + the lattice is two-dimensional, only the x- and y-coordinates are + specified. + indices : list of tuple + A list of all possible (z,y,x) or (y,x) lattice element indices. These + indices correspond to indices in the :attr:`RectLattice.universes` + property. + ndim : int + The number of dimensions of the lattice + shape : Iterable of int + An array of two or three integers representing the number of lattice + cells in the x- and y- (and z-) directions, respectively. """ @@ -274,7 +289,6 @@ class RectLattice(Lattice): super(RectLattice, self).__init__(lattice_id, name) # Initialize Lattice class attributes - self._dimension = None self._lower_left = None self._offsets = None @@ -283,7 +297,7 @@ class RectLattice(Lattice): return False elif not super(RectLattice, self).__eq__(other): return False - elif self.dimension != other.dimension: + elif self.shape != other.shape: return False elif self.lower_left != other.lower_left: return False @@ -300,8 +314,8 @@ class RectLattice(Lattice): string = 'RectLattice\n' string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) - string += '{0: <16}{1}{2}\n'.format('\tDimension', '=\t', - self._dimension) + string += '{0: <16}{1}{2}\n'.format('\tShape', '=\t', + self.shape) string += '{0: <16}{1}{2}\n'.format('\tLower Left', '=\t', self._lower_left) string += '{0: <16}{1}{2}\n'.format('\tPitch', '=\t', self._pitch) @@ -320,7 +334,7 @@ class RectLattice(Lattice): string += '{0} '.format(universe._id) # Add a newline character every time we reach end of row of cells - if (i+1) % self._dimension[-1] == 0: + if (i+1) % self.shape[0] == 0: string += '\n' string = string.rstrip('\n') @@ -333,7 +347,7 @@ class RectLattice(Lattice): string += '{0} '.format(offset) # Add a newline character when we reach end of row of cells - if (i+1) % self._dimension[-1] == 0: + if (i+1) % self.shape[0] == 0: string += '\n' string = string.rstrip('\n') @@ -341,24 +355,29 @@ class RectLattice(Lattice): return string @property - def dimension(self): - return self._dimension + def indices(self): + if self.ndim == 2: + return list(np.broadcast(*np.ogrid[ + :self.shape[1], :self.shape[0]])) + else: + return list(np.broadcast(*np.ogrid[ + :self.shape[2], :self.shape[1], :self.shape[0]])) @property def lower_left(self): return self._lower_left + @property + def ndim(self): + return len(self.pitch) + @property def offsets(self): return self._offsets - @dimension.setter - def dimension(self, dimension): - cv.check_type('lattice dimension', dimension, Iterable, Integral) - cv.check_length('lattice dimension', dimension, 2, 3) - for dim in dimension: - cv.check_greater_than('lattice dimension', dim, 0) - self._dimension = dimension + @property + def shape(self): + return self._universes.shape[::-1] @lower_left.setter def lower_left(self, lower_left): @@ -379,8 +398,13 @@ class RectLattice(Lattice): cv.check_greater_than('lattice pitch', dim, 0.0) self._pitch = pitch - def get_cell_instance(self, path, distribcell_index): + @Lattice.universes.setter + def universes(self, universes): + cv.check_iterable_type('lattice universes', universes, openmc.Universe, + min_depth=2, max_depth=3) + self._universes = np.asarray(universes) + def get_cell_instance(self, path, distribcell_index): # Extract the lattice element from the path next_index = path.index('-') lat_id_indices = path[:next_index] @@ -395,7 +419,7 @@ class RectLattice(Lattice): lat_z = int(i.split(',')[2]) - 1 # For 2D Lattices - if len(self._dimension) == 2: + if self.ndim == 2: offset = self._offsets[lat_z, lat_y, lat_x, distribcell_index-1] offset += self._universes[lat_x][lat_y].get_cell_instance(path, distribcell_index) @@ -408,6 +432,128 @@ class RectLattice(Lattice): return offset + def find_element(self, point): + """Determine index of lattice element and local coordinates for a point + + Parameters + ---------- + point : Iterable of float + Cartesian coordinates of point + + Returns + ------- + 2- or 3-tuple of int + A tuple of the corresponding (x,y,z) lattice element indices + 3-tuple of float + Carestian coordinates of the point in the corresponding lattice + element coordinate system + + """ + ix = floor((point[0] - self._lower_left[0])/self._pitch[0]) + iy = floor((point[1] - self._lower_left[1])/self._pitch[1]) + if self.ndim == 2: + idx = (ix, iy) + else: + iz = floor((point[2] - self._lower_left[2])/self._pitch[2]) + idx = (ix, iy, iz) + return idx, self.get_local_coordinates(point, idx) + + def get_local_coordinates(self, point, idx): + """Determine local coordinates of a point within a lattice element + + Parameters + ---------- + point : Iterable of float + Cartesian coordinates of point + idx : Iterable of int + (x,y,z) indices of lattice element. If the lattice is 2D, the z + index can be omitted. + + Returns + ------- + 3-tuple of float + Cartesian coordinates of point in the lattice element coordinate + system + + """ + x = point[0] - (self._lower_left[0] + (idx[0] + 0.5)*self._pitch[0]) + y = point[1] - (self._lower_left[1] + (idx[1] + 0.5)*self._pitch[1]) + if self.ndim == 2: + z = point[2] + else: + z = point[2] - (self._lower_left[2] + (idx[2] + 0.5)*self._pitch[2]) + return (x, y, z) + + def get_universe_index(self, idx): + """Return index in the universes array corresponding to a lattice element index + + Parameters + ---------- + idx : Iterable of int + Lattice element indices in the :math:`(x,y,z)` coordinate system + + Returns + ------- + 2- or 3-tuple of int + Indices used when setting the :attr:`RectLattice.universes` property + + """ + max_y = self.shape[1] - 1 + if self.ndim == 2: + x, y = idx + return (max_y - y, x) + else: + x, y, z = idx + return (z, max_y - y, x) + + def is_valid_index(self, idx): + """Determine whether lattice element index is within defined range + + Parameters + ---------- + idx : Iterable of int + Lattice element indices in the :math:`(x,y,z)` coordinate system + + Returns + ------- + bool + Whether index is valid + + """ + if self.ndim == 2: + return (0 <= idx[0] < self.shape[0] and + 0 <= idx[1] < self.shape[1]) + else: + return (0 <= idx[0] < self.shape[0] and + 0 <= idx[1] < self.shape[1] and + 0 <= idx[2] < self.shape[2]) + + def find(self, point): + """Find cells/universes/lattices which contain a given point + + Parameters + ---------- + point : 3-tuple of float + Cartesian coordinatesof the point + + Returns + ------- + list + Sequence of universes, cells, and lattices which are traversed to + find the given point + + """ + idx, p = self.find_element(point) + if self.is_valid_index(idx): + idx_u = self.get_universe_index(idx) + u = self.universes[idx_u] + else: + if self.outer is not None: + u = self.outer + else: + return [] + return [(self, idx)] + u.find(p) + def create_xml_subelement(self, xml_element): # Determine if XML element already contains subelement for this Lattice @@ -436,7 +582,7 @@ class RectLattice(Lattice): # Export Lattice cell dimensions dimension = ET.SubElement(lattice_subelement, "dimension") - dimension.text = ' '.join(map(str, self._dimension)) + dimension.text = ' '.join(map(str, self.shape)) # Export Lattice lower left lower_left = ET.SubElement(lattice_subelement, "lower_left") @@ -446,10 +592,10 @@ class RectLattice(Lattice): universe_ids = '\n' # 3D Lattices - if len(self._dimension) == 3: - for z in range(self._dimension[2]): - for y in range(self._dimension[1]): - for x in range(self._dimension[0]): + if self.ndim == 3: + for z in range(self.shape[2]): + for y in range(self.shape[1]): + for x in range(self.shape[0]): universe = self._universes[z][y][x] # Append Universe ID to the Lattice XML subelement @@ -466,8 +612,8 @@ class RectLattice(Lattice): # 2D Lattices else: - for y in range(self._dimension[1]): - for x in range(self._dimension[0]): + for y in range(self.shape[1]): + for x in range(self.shape[0]): universe = self._universes[y][x] # Append Universe ID to Lattice XML subelement @@ -492,6 +638,10 @@ class RectLattice(Lattice): class HexLattice(Lattice): """A lattice consisting of hexagonal prisms. + To completely define a hexagonal lattice, the :attr:`HexLattice.center`, + :attr:`HexLattice.pitch`, :attr:`HexLattice.universes`, and + :attr:`HexLattice.outer` properties need to be set. + Parameters ---------- lattice_id : int, optional @@ -506,26 +656,31 @@ class HexLattice(Lattice): Unique identifier for the lattice name : str Name of the lattice - num_rings : int - Number of radial ring positions in the xy-plane - num_axial : int - Number of positions along the z-axis. - center : Iterable of float - Coordinates of the center of the lattice. If the lattice does not have - axial sections then only the x- and y-coordinates are specified pitch : Iterable of float Pitch of the lattice in cm. The first item in the iterable specifies the pitch in the radial direction and, if the lattice is 3D, the second item in the iterable specifies the pitch in the axial direction. outer : openmc.Universe A universe to fill all space outside the lattice - universes : Iterable of Iterable of openmc.Universe + universes : Nested Iterable of openmc.Universe A two- or three-dimensional list/array of universes filling each element of the lattice. Each sub-list corresponds to one ring of universes and should be ordered from outermost ring to innermost ring. The universes within each sub-list are ordered from the "top" and proceed in a clockwise fashion. The :meth:`HexLattice.show_indices` method can be used to help figure out indices for this property. + center : Iterable of float + Coordinates of the center of the lattice. If the lattice does not have + axial sections then only the x- and y-coordinates are specified + indices : list of tuple + A list of all possible (z,r,i) or (r,i) lattice element indices that are + possible, where z is the axial index, r is in the ring index (starting + from the outermost ring), and i is the index with a ring starting from + the top and proceeding clockwise. + num_rings : int + Number of radial ring positions in the xy-plane + num_axial : int + Number of positions along the z-axis. """ @@ -597,17 +752,15 @@ class HexLattice(Lattice): def center(self): return self._center - @num_rings.setter - def num_rings(self, num_rings): - cv.check_type('number of rings', num_rings, Integral) - cv.check_greater_than('number of rings', num_rings, 0) - self._num_rings = num_rings - - @num_axial.setter - def num_axial(self, num_axial): - cv.check_type('number of axial', num_axial, Integral) - cv.check_greater_than('number of axial', num_axial, 0) - self._num_axial = num_axial + @property + def indices(self): + if self.num_axial is None: + return [(r, i) for r in range(self._num_rings) + for i in range(max(6*(self._num_rings - 1 - r), 1))] + else: + return [(z, r, i) for z in range(self._num_axial) + for r in range(self._num_rings) + for i in range(max(6*(self._num_rings - 1 - r), 1))] @center.setter def center(self, center): @@ -625,8 +778,9 @@ class HexLattice(Lattice): @Lattice.universes.setter def universes(self, universes): - # Call Lattice.universes parent class setter property - Lattice.universes.fset(self, universes) + cv.check_iterable_type('lattice universes', universes, openmc.Universe, + min_depth=2, max_depth=3) + self._universes = universes # NOTE: This routine assumes that the user creates a "ragged" list of # lists, where each sub-list corresponds to one ring of Universes. @@ -649,14 +803,14 @@ class HexLattice(Lattice): # Set the number of axial positions. if n_dims == 3: - self.num_axial = len(self._universes) + self._num_axial = len(self._universes) else: self._num_axial = None # Set the number of rings and make sure this number is consistent for # all axial positions. if n_dims == 3: - self.num_rings = len(self._universes[0]) + self._num_rings = len(self._universes[0]) for rings in self._universes: if len(rings) != self._num_rings: msg = 'HexLattice ID={0:d} has an inconsistent number of ' \ @@ -664,7 +818,7 @@ class HexLattice(Lattice): raise ValueError(msg) else: - self.num_rings = len(self._universes) + self._num_rings = len(self._universes) # Make sure there are the correct number of elements in each ring. if n_dims == 3: @@ -705,6 +859,170 @@ class HexLattice(Lattice): 6*(self._num_rings - 1 - r)) raise ValueError(msg) + def find_element(self, point): + """Determine index of lattice element and local coordinates for a point + + Parameters + ---------- + point : Iterable of float + Cartesian coordinates of point + + Returns + ------- + 3-tuple of int + Indices of corresponding lattice element in (x,:math:`alpha`,z) + bases + numpy.ndarray + Carestian coordinates of the point in the corresponding lattice + element coordinate system + + """ + # Convert coordinates to skewed bases + x = point[0] - self._center[0] + y = point[1] - self._center[1] + if self._num_axial is None: + iz = 1 + else: + z = point[2] - self._center[2] + iz = floor(z/self._pitch[1] + 0.5*self._num_axial) + alpha = y - x/sqrt(3.) + ix = floor(x/(sqrt(0.75) * self._pitch[0])) + ia = floor(alpha/self._pitch[0]) + + # Check four lattice elements to see which one is closest based on local + # coordinates + d_min = np.inf + for idx in [(ix, ia, iz), (ix + 1, ia, iz), (ix, ia + 1, iz), + (ix + 1, ia + 1, iz)]: + p = self.get_local_coordinates(point, idx) + d = p[0]**2 + p[1]**2 + if d < d_min: + d_min = d + idx_min = idx + p_min = p + + return idx_min, p_min + + def get_local_coordinates(self, point, idx): + """Determine local coordinates of a point within a lattice element + + Parameters + ---------- + point : Iterable of float + Cartesian coordinates of point + idx : Iterable of int + Indices of lattice element in (x,:math:`alpha`,z) bases + + Returns + ------- + 3-tuple of float + Cartesian coordinates of point in the lattice element coordinate + system + + """ + x = point[0] - (self._center[0] + sqrt(0.75)*self._pitch[0]*idx[0]) + y = point[1] - (self._center[1] + (0.5*idx[0] + idx[1])*self._pitch[0]) + if self._num_axial is None: + z = point[2] + else: + z = point[2] - (self._center[2] + (idx[2] + 0.5 - 0.5*self._num_axial)* + self._pitch[1]) + return (x, y, z) + + def get_universe_index(self, idx): + """Return index in the universes array corresponding to a lattice element index + + Parameters + ---------- + idx : Iterable of int + Lattice element indices in the :math:`(x,\alpha,z)` coordinate + system + + Returns + ------- + 2- or 3-tuple of int + Indices used when setting the :attr:`HexLattice.universes` property + + """ + + # First we determine which ring the index corresponds to. + x = idx[0] + a = idx[1] + z = -a - x + g = max(abs(x), abs(a), abs(z)) + + # Next we use a clever method to figure out where along the ring we are. + i_ring = self._num_rings - 1 - g + if x >= 0: + if a >= 0: + i_within = x + else: + i_within = 2*g + z + else: + if a <= 0: + i_within = 3*g - x + else: + i_within = 5*g - z + + if self.num_axial is None: + return (i_ring, i_within) + else: + return (idx[2], i_ring, i_within) + + def is_valid_index(self, idx): + """Determine whether lattice element index is within defined range + + Parameters + ---------- + idx : Iterable of int + Lattice element indices in the :math:`(x,\alpha,z)` coordinate + system + + Returns + ------- + bool + Whether index is valid + + """ + x = idx[0] + y = idx[1] + z = 0 - y - x + g = max(abs(x), abs(y), abs(z)) + if self.num_axial is None: + return g < self.num_rings + else: + return g < self.num_rings and 0 <= idx[2] < self.num_axial + + def find(self, point): + """Find cells/universes/lattices which contain a given point + + Parameters + ---------- + point : 3-tuple of float + Cartesian coordinatesof the point + + Returns + ------- + list + Sequence of universes, cells, and lattices which are traversed to + find the given point + + """ + idx, p = self.find_element(point) + if self.is_valid_index(idx): + idx_u = self.get_universe_index(idx) + if self.num_axial is None: + u = self.universes[idx_u[0]][idx_u[1]] + else: + u = self.universes[idx_u[0]][idx_u[1]][idx_u[2]] + else: + if self.outer is not None: + u = self.outer + else: + return [] + + return [(self, idx)] + u.find(p) + def create_xml_subelement(self, xml_element): # Determine if XML element already contains subelement for this Lattice path = './hex_lattice[@id=\'{0}\']'.format(self._id) @@ -736,8 +1054,8 @@ class HexLattice(Lattice): lattice_subelement.set("n_axial", str(self._num_axial)) # Export Lattice cell center - dimension = ET.SubElement(lattice_subelement, "center") - dimension.text = ' '.join(map(str, self._center)) + center = ET.SubElement(lattice_subelement, "center") + center.text = ' '.join(map(str, self._center)) # Export the Lattice nested Universe IDs. diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index 562fe9cad..c6b364d7c 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -861,8 +861,8 @@ def get_opencg_lattice(openmc_lattice): universes = new_universes # Initialize an empty array for the OpenCG nested Universes in this Lattice - universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), - dtype=opencg.Universe) + universe_array = np.empty(tuple(np.array(dimension)[::-1]), + dtype=opencg.Universe) # Create OpenCG Universes for each unique nested Universe in this Lattice unique_universes = openmc_lattice.get_unique_universes() @@ -929,8 +929,8 @@ def get_openmc_lattice(opencg_lattice): outer = opencg_lattice.outside # Initialize an empty array for the OpenMC nested Universes in this Lattice - universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), - dtype=openmc.Universe) + universe_array = np.empty(tuple(np.array(dimension)[::-1]), + dtype=openmc.Universe) # Create OpenMC Universes for each unique nested Universe in this Lattice unique_universes = opencg_lattice.get_unique_universes() @@ -953,7 +953,6 @@ def get_openmc_lattice(opencg_lattice): np.array(dimension, dtype=np.float64))) / -2.0 openmc_lattice = openmc.RectLattice(lattice_id=lattice_id) - openmc_lattice.dimension = dimension openmc_lattice.pitch = width openmc_lattice.universes = universe_array openmc_lattice.lower_left = lower_left diff --git a/openmc/summary.py b/openmc/summary.py index 6fb6b7000..c3277809a 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -358,7 +358,6 @@ class Summary(object): # Create the Lattice lattice = openmc.RectLattice(lattice_id=lattice_id, name=name) - lattice.dimension = tuple(dimension) lattice.lower_left = lower_left lattice.pitch = pitch @@ -368,7 +367,7 @@ class Summary(object): # Build array of Universe pointers for the Lattice universes = \ - np.ndarray(tuple(universe_ids.shape), dtype=openmc.Universe) + np.empty(tuple(universe_ids.shape), dtype=openmc.Universe) for z in range(universe_ids.shape[0]): for y in range(universe_ids.shape[1]): @@ -403,8 +402,6 @@ class Summary(object): # Create the Lattice lattice = openmc.HexLattice(lattice_id=lattice_id, name=name) - lattice.num_rings = n_rings - lattice.num_axial = n_axial lattice.center = center lattice.pitch = pitch @@ -417,12 +414,12 @@ class Summary(object): # (x, alpha, z) to the Python API's format of a ragged nested # list of (z, ring, theta). universes = [] - for z in range(lattice.num_axial): + for z in range(n_axial): # Add a list for this axial level. universes.append([]) - x = lattice.num_rings - 1 - a = 2*lattice.num_rings - 2 - for r in range(lattice.num_rings - 1, 0, -1): + x = n_rings - 1 + a = 2*n_rings - 2 + for r in range(n_rings - 1, 0, -1): # Add a list for this ring. universes[-1].append([]) diff --git a/openmc/universe.py b/openmc/universe.py index 770e789da..d8e1c4dac 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -1,6 +1,7 @@ from collections import OrderedDict, Iterable from numbers import Integral from xml.etree import ElementTree as ET +import random import sys import warnings @@ -124,6 +125,110 @@ class Universe(object): else: self._name = '' + def find(self, point): + """Find cells/universes/lattices which contain a given point + + Parameters + ---------- + point : 3-tuple of float + Cartesian coordinatesof the point + + Returns + ------- + list + Sequence of universes, cells, and lattices which are traversed to + find the given point + + """ + p = np.asarray(point) + for cell in self._cells.values(): + if p in cell: + if cell._type in ('normal', 'void'): + return [self, cell] + elif cell._type == 'fill': + if cell.translation is not None: + p -= cell.translation + if cell.rotation is not None: + p[:] = cell.rotation_matrix.dot(p) + return [self, cell] + cell.fill.find(p) + else: + return [self, cell] + cell.fill.find(p) + return [] + + def plot(self, center=(0., 0., 0.), width=(1., 1.), pixels=(200, 200), + basis='xy', color_by='cell'): + """Display a slice plot of the universe. + + Parameters + ---------- + center : Iterable of float + Coordinates at the center of the plot + width : Iterable of float + Width of the plot in each basis direction + pixels : Iterable of int + Number of pixels to use in each basis direction + basis : {'xy', 'xz', 'yz'} + The basis directions for the plot + color_by : {'cell', 'material'} + Indicate whether the plot should be colored by cell or by material + + """ + import matplotlib.pyplot as plt + + if basis == 'xy': + x_min = center[0] - 0.5*width[0] + x_max = center[0] + 0.5*width[0] + y_min = center[1] - 0.5*width[1] + y_max = center[1] + 0.5*width[1] + elif basis == 'yz': + # The x-axis will correspond to physical y and the y-axis will correspond to physical z + x_min = center[1] - 0.5*width[0] + x_max = center[1] + 0.5*width[0] + y_min = center[2] - 0.5*width[1] + y_max = center[2] + 0.5*width[1] + elif basis == 'xz': + # The y-axis will correspond to physical z + x_min = center[0] - 0.5*width[0] + x_max = center[0] + 0.5*width[0] + y_min = center[2] - 0.5*width[1] + y_max = center[2] + 0.5*width[1] + + # Determine locations to determine cells at + x_coords = np.linspace(x_min, x_max, pixels[0], endpoint=False) + \ + 0.5*(x_max - x_min)/pixels[0] + y_coords = np.linspace(y_max, y_min, pixels[1], endpoint=False) - \ + 0.5*(y_max - y_min)/pixels[1] + + colors = {} + img = np.zeros(pixels + (4,)) # Use RGBA form + for i, x in enumerate(x_coords): + for j, y in enumerate(y_coords): + if basis == 'xy': + path = self.find((x, y, center[2])) + elif basis == 'yz': + path = self.find((center[0], x, y)) + elif basis == 'xz': + path = self.find((x, center[1], y)) + + if len(path) > 0: + try: + if color_by == 'cell': + uid = path[-1].id + elif color_by == 'material': + if path[-1].fill_type == 'material': + uid = path[-1].fill.id + else: + continue + except AttributeError: + continue + if uid not in colors: + colors[uid] = (random.random(), random.random(), + random.random(), 1.0) + img[j,i,:] = colors[uid] + + plt.imshow(img, extent=(x_min, x_max, y_min, y_max)) + plt.show() + def add_cell(self, cell): """Add a cell to the universe. diff --git a/tests/input_set.py b/tests/input_set.py index 2c6841e25..fefd2ca5f 100644 --- a/tests/input_set.py +++ b/tests/input_set.py @@ -350,7 +350,6 @@ class InputSet(object): # Define fuel lattices. l100 = openmc.RectLattice(name='Fuel assembly (lower half)', lattice_id=100) - l100.dimension = (17, 17) l100.lower_left = (-10.71, -10.71) l100.pitch = (1.26, 1.26) l100.universes = [ @@ -384,7 +383,6 @@ class InputSet(object): l101 = openmc.RectLattice(name='Fuel assembly (upper half)', lattice_id=101) - l101.dimension = (17, 17) l101.lower_left = (-10.71, -10.71) l101.pitch = (1.26, 1.26) l101.universes = [ @@ -444,7 +442,6 @@ class InputSet(object): # Define core lattices l200 = openmc.RectLattice(name='Core lattice (lower half)', lattice_id=200) - l200.dimension = (21, 21) l200.lower_left = (-224.91, -224.91) l200.pitch = (21.42, 21.42) l200.universes = [ @@ -472,7 +469,6 @@ class InputSet(object): l201 = openmc.RectLattice(name='Core lattice (lower half)', lattice_id=201) - l201.dimension = (21, 21) l201.lower_left = (-224.91, -224.91) l201.pitch = (21.42, 21.42) l201.universes = [ diff --git a/tests/test_asymmetric_lattice/test_asymmetric_lattice.py b/tests/test_asymmetric_lattice/test_asymmetric_lattice.py index 504cc4746..40def1f62 100644 --- a/tests/test_asymmetric_lattice/test_asymmetric_lattice.py +++ b/tests/test_asymmetric_lattice/test_asymmetric_lattice.py @@ -24,7 +24,6 @@ class AsymmetricLatticeTestHarness(PyAPITestHarness): # Construct a 3x3 lattice of fuel assemblies core_lat = openmc.RectLattice(name='3x3 Core Lattice', lattice_id=202) - core_lat.dimension = (3, 3) core_lat.lower_left = (-32.13, -32.13) core_lat.pitch = (21.42, 21.42) core_lat.universes = [[fuel, water, water], diff --git a/tests/test_distribmat/test_distribmat.py b/tests/test_distribmat/test_distribmat.py index d8f78c5cf..6700a96b6 100644 --- a/tests/test_distribmat/test_distribmat.py +++ b/tests/test_distribmat/test_distribmat.py @@ -53,7 +53,6 @@ class DistribmatTestHarness(PyAPITestHarness): fuel_univ.add_cells((c11, c12)) lat = openmc.RectLattice(lattice_id=101) - lat.dimension = [2, 2] lat.lower_left = [-2.0, -2.0] lat.pitch = [2.0, 2.0] lat.universes = [[fuel_univ]*2]*2 From 1f294c614e874ca92177aae8e8cd5f3de96834db Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 25 May 2016 09:57:42 -0500 Subject: [PATCH 5/9] Remove use of Cell.add_surface in OpenCG compatibility module --- openmc/opencg_compatible.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index c6b364d7c..b112ed327 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -1,4 +1,5 @@ import copy +import operator import numpy as np @@ -9,8 +10,6 @@ except ImportError: raise ImportError(msg) import openmc -from openmc.region import Intersection -from openmc.surface import Halfspace import openmc.checkvalue as cv @@ -467,13 +466,13 @@ def get_opencg_cell(openmc_cell): # half-spaces, i.e., no complex cells. region = openmc_cell.region if region is not None: - if isinstance(region, Halfspace): + if isinstance(region, openmc.Halfspace): surface = region.surface halfspace = -1 if region.side == '-' else 1 opencg_cell.add_surface(get_opencg_surface(surface), halfspace) - elif isinstance(region, Intersection): + elif isinstance(region, openmc.Intersection): for node in region.nodes: - if not isinstance(node, Halfspace): + if not isinstance(node, openmc.Halfspace): raise NotImplementedError("Complex cells not yet " "supported in OpenCG.") surface = node.surface @@ -697,12 +696,13 @@ def get_openmc_cell(opencg_cell): translation = np.asarray(opencg_cell.translation, dtype=np.float64) openmc_cell.translation = translation - surfaces = opencg_cell.surfaces - - for surface_id in surfaces: - surface = surfaces[surface_id][0] - halfspace = surfaces[surface_id][1] - openmc_cell.add_surface(get_openmc_surface(surface), halfspace) + surfaces = [] + operators = [] + for surface, halfspace in opencg_cell.surfaces.values(): + surfaces.append(get_openmc_surface(surface)) + operators.append(operator.neg if halfspace == -1 else operator.pos) + openmc_cell.region = openmc.Intersection( + *[op(s) for op, s in zip(operators, surfaces)]) # Add the OpenMC Cell to the global collection of all OpenMC Cells OPENMC_CELLS[cell_id] = openmc_cell From eb3f888c899b905429963edfc6fdbd5a75e1e071 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 25 May 2016 10:02:27 -0500 Subject: [PATCH 6/9] Allow root_universe to be passed to Geometry constructor --- examples/python/basic/build-xml.py | 3 +-- examples/python/boxes/build-xml.py | 3 +-- examples/python/lattice/hexagonal/build-xml.py | 3 +-- examples/python/lattice/nested/build-xml.py | 3 +-- examples/python/lattice/simple/build-xml.py | 3 +-- examples/python/pincell/build-xml.py | 3 +-- examples/python/pincell_multigroup/build-xml.py | 3 +-- examples/python/reflective/build-xml.py | 3 +-- openmc/geometry.py | 9 ++++++++- 9 files changed, 16 insertions(+), 17 deletions(-) diff --git a/examples/python/basic/build-xml.py b/examples/python/basic/build-xml.py index ffff03720..81aecc9f9 100644 --- a/examples/python/basic/build-xml.py +++ b/examples/python/basic/build-xml.py @@ -74,8 +74,7 @@ universe1.add_cells([cell2, cell3]) root.add_cells([cell1, cell4]) # Instantiate a Geometry, register the root Universe, and export to XML -geometry = openmc.Geometry() -geometry.root_universe = root +geometry = openmc.Geometry(root) geometry.export_to_xml() diff --git a/examples/python/boxes/build-xml.py b/examples/python/boxes/build-xml.py index 814f60beb..4be33dcf1 100644 --- a/examples/python/boxes/build-xml.py +++ b/examples/python/boxes/build-xml.py @@ -97,8 +97,7 @@ root = openmc.Universe(universe_id=0, name='root universe') root.add_cells([inner_box, middle_box, outer_box]) # Instantiate a Geometry, register the root Universe, and export to XML -geometry = openmc.Geometry() -geometry.root_universe = root +geometry = openmc.Geometry(root) geometry.export_to_xml() diff --git a/examples/python/lattice/hexagonal/build-xml.py b/examples/python/lattice/hexagonal/build-xml.py index ef3a12847..05cb2cb01 100644 --- a/examples/python/lattice/hexagonal/build-xml.py +++ b/examples/python/lattice/hexagonal/build-xml.py @@ -105,8 +105,7 @@ lattice.outer = univ2 cell1.fill = lattice # Instantiate a Geometry, register the root Universe, and export to XML -geometry = openmc.Geometry() -geometry.root_universe = root +geometry = openmc.Geometry(root) geometry.export_to_xml() diff --git a/examples/python/lattice/nested/build-xml.py b/examples/python/lattice/nested/build-xml.py index a964d882c..03cede9dc 100644 --- a/examples/python/lattice/nested/build-xml.py +++ b/examples/python/lattice/nested/build-xml.py @@ -114,8 +114,7 @@ cell1.fill = lattice2 cell2.fill = lattice1 # Instantiate a Geometry, register the root Universe, and export to XML -geometry = openmc.Geometry() -geometry.root_universe = root +geometry = openmc.Geometry(root) geometry.export_to_xml() diff --git a/examples/python/lattice/simple/build-xml.py b/examples/python/lattice/simple/build-xml.py index 4961b96b8..5a642d308 100644 --- a/examples/python/lattice/simple/build-xml.py +++ b/examples/python/lattice/simple/build-xml.py @@ -105,8 +105,7 @@ lattice.universes = [[univ1, univ2, univ1, univ2], cell1.fill = lattice # Instantiate a Geometry, register the root Universe, and export to XML -geometry = openmc.Geometry() -geometry.root_universe = root +geometry = openmc.Geometry(root) geometry.export_to_xml() diff --git a/examples/python/pincell/build-xml.py b/examples/python/pincell/build-xml.py index a3be3e97e..0afb2527f 100644 --- a/examples/python/pincell/build-xml.py +++ b/examples/python/pincell/build-xml.py @@ -149,8 +149,7 @@ root = openmc.Universe(universe_id=0, name='root universe') root.add_cells([fuel, gap, clad, water]) # Instantiate a Geometry, register the root Universe, and export to XML -geometry = openmc.Geometry() -geometry.root_universe = root +geometry = openmc.Geometry(root) geometry.export_to_xml() diff --git a/examples/python/pincell_multigroup/build-xml.py b/examples/python/pincell_multigroup/build-xml.py index 5ac5b376a..6dbfa336b 100644 --- a/examples/python/pincell_multigroup/build-xml.py +++ b/examples/python/pincell_multigroup/build-xml.py @@ -119,8 +119,7 @@ root = openmc.Universe(universe_id=0, name='root universe') root.add_cells([fuel, moderator]) # Instantiate a Geometry, register the root Universe, and export to XML -geometry = openmc.Geometry() -geometry.root_universe = root +geometry = openmc.Geometry(root) geometry.export_to_xml() diff --git a/examples/python/reflective/build-xml.py b/examples/python/reflective/build-xml.py index 4ecd0351f..949e57c8c 100644 --- a/examples/python/reflective/build-xml.py +++ b/examples/python/reflective/build-xml.py @@ -64,8 +64,7 @@ root = openmc.Universe(universe_id=0, name='root universe') root.add_cell(cell) # Instantiate a Geometry, register the root Universe, and export to XML -geometry = openmc.Geometry() -geometry.root_universe = root +geometry = openmc.Geometry(root) geometry.export_to_xml() diff --git a/openmc/geometry.py b/openmc/geometry.py index 7eddadfc7..4e85929fc 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -15,6 +15,11 @@ def reset_auto_ids(): class Geometry(object): """Geometry representing a collection of surfaces, cells, and universes. + Parameters + ---------- + root_universe : openmc.Universe, optional + Root universe which contains all others + Attributes ---------- root_universe : openmc.Universe @@ -22,9 +27,11 @@ class Geometry(object): """ - def __init__(self): + def __init__(self, root_universe=None): self._root_universe = None self._offsets = {} + if root_universe is not None: + self.root_universe = root_universe @property def root_universe(self): From 08b8082b7bc78aeb968c0ac49bf9636745471bdc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 29 May 2016 10:20:11 -0500 Subject: [PATCH 7/9] First round of fixes for @wbinventor comments on #656 --- docs/source/methods/geometry.rst | 2 + openmc/geometry.py | 2 +- openmc/lattice.py | 71 +++++++++++++++++++------------- openmc/universe.py | 11 ++++- 4 files changed, 55 insertions(+), 31 deletions(-) diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index f642cca10..bcd568ca3 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -437,6 +437,8 @@ where :math:`(x_0, y_0, z_0)` are the coordinates to the lower-left-bottom corner of the lattice, and :math:`p_0, p_1, p_2` are the pitches along the :math:`x`, :math:`y`, and :math:`z` axes, respectively. +.. _hexagonal_indexing: + Hexagonal Lattice Indexing -------------------------- diff --git a/openmc/geometry.py b/openmc/geometry.py index 4e85929fc..006151900 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -75,7 +75,7 @@ class Geometry(object): Parameters ---------- point : 3-tuple of float - Cartesian coordinatesof the point + Cartesian coordinates of the point Returns ------- diff --git a/openmc/lattice.py b/openmc/lattice.py index d8deddc42..0277d50cb 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -242,6 +242,14 @@ class RectLattice(Lattice): :attr:`RectLattice.outer`, and :attr:`RectLattice.universes` properties need to be set. + Most methods for this class use a natural indexing scheme wherein elements + are assigned an index corresponding to their position relative to the + (x,y,z) axes in a Cartesian coordinate system, i.e., an index of (0,0,0) in + the lattice gives the element whose x, y, and z coordinates are the + smallest. However, note that when universes are assigned to lattice elements + using the :attr:`RectLattice.universes` property, the array indices do not + correspond to natural indices. + Parameters ---------- lattice_id : int, optional @@ -449,12 +457,12 @@ class RectLattice(Lattice): element coordinate system """ - ix = floor((point[0] - self._lower_left[0])/self._pitch[0]) - iy = floor((point[1] - self._lower_left[1])/self._pitch[1]) + ix = floor((point[0] - self.lower_left[0])/self.pitch[0]) + iy = floor((point[1] - self.lower_left[1])/self.pitch[1]) if self.ndim == 2: idx = (ix, iy) else: - iz = floor((point[2] - self._lower_left[2])/self._pitch[2]) + iz = floor((point[2] - self.lower_left[2])/self.pitch[2]) idx = (ix, iy, iz) return idx, self.get_local_coordinates(point, idx) @@ -476,12 +484,12 @@ class RectLattice(Lattice): system """ - x = point[0] - (self._lower_left[0] + (idx[0] + 0.5)*self._pitch[0]) - y = point[1] - (self._lower_left[1] + (idx[1] + 0.5)*self._pitch[1]) + x = point[0] - (self.lower_left[0] + (idx[0] + 0.5)*self.pitch[0]) + y = point[1] - (self.lower_left[1] + (idx[1] + 0.5)*self.pitch[1]) if self.ndim == 2: z = point[2] else: - z = point[2] - (self._lower_left[2] + (idx[2] + 0.5)*self._pitch[2]) + z = point[2] - (self.lower_left[2] + (idx[2] + 0.5)*self.pitch[2]) return (x, y, z) def get_universe_index(self, idx): @@ -636,12 +644,19 @@ class RectLattice(Lattice): class HexLattice(Lattice): - """A lattice consisting of hexagonal prisms. + r"""A lattice consisting of hexagonal prisms. To completely define a hexagonal lattice, the :attr:`HexLattice.center`, :attr:`HexLattice.pitch`, :attr:`HexLattice.universes`, and :attr:`HexLattice.outer` properties need to be set. + Most methods for this class use a natural indexing scheme wherein elements + are assigned an index corresponding to their position relative to skewed + :math:`(x,\alpha,z)` axes as described fully in + :ref:`hexagonal_indexing`. However, note that when universes are assigned to + lattice elements using the :attr:`RectLattice.universes` property, the array + indices do not correspond to natural indices. + Parameters ---------- lattice_id : int, optional @@ -755,12 +770,12 @@ class HexLattice(Lattice): @property def indices(self): if self.num_axial is None: - return [(r, i) for r in range(self._num_rings) - for i in range(max(6*(self._num_rings - 1 - r), 1))] + return [(r, i) for r in range(self.num_rings) + for i in range(max(6*(self.num_rings - 1 - r), 1))] else: - return [(z, r, i) for z in range(self._num_axial) - for r in range(self._num_rings) - for i in range(max(6*(self._num_rings - 1 - r), 1))] + return [(z, r, i) for z in range(self.num_axial) + for r in range(self.num_rings) + for i in range(max(6*(self.num_rings - 1 - r), 1))] @center.setter def center(self, center): @@ -860,7 +875,7 @@ class HexLattice(Lattice): raise ValueError(msg) def find_element(self, point): - """Determine index of lattice element and local coordinates for a point + r"""Determine index of lattice element and local coordinates for a point Parameters ---------- @@ -870,7 +885,7 @@ class HexLattice(Lattice): Returns ------- 3-tuple of int - Indices of corresponding lattice element in (x,:math:`alpha`,z) + Indices of corresponding lattice element in :math:`(x,\alpha,z)` bases numpy.ndarray Carestian coordinates of the point in the corresponding lattice @@ -878,16 +893,16 @@ class HexLattice(Lattice): """ # Convert coordinates to skewed bases - x = point[0] - self._center[0] - y = point[1] - self._center[1] + x = point[0] - self.center[0] + y = point[1] - self.center[1] if self._num_axial is None: iz = 1 else: - z = point[2] - self._center[2] - iz = floor(z/self._pitch[1] + 0.5*self._num_axial) + z = point[2] - self.center[2] + iz = floor(z/self.pitch[1] + 0.5*self.num_axial) alpha = y - x/sqrt(3.) - ix = floor(x/(sqrt(0.75) * self._pitch[0])) - ia = floor(alpha/self._pitch[0]) + ix = floor(x/(sqrt(0.75) * self.pitch[0])) + ia = floor(alpha/self.pitch[0]) # Check four lattice elements to see which one is closest based on local # coordinates @@ -904,14 +919,14 @@ class HexLattice(Lattice): return idx_min, p_min def get_local_coordinates(self, point, idx): - """Determine local coordinates of a point within a lattice element + r"""Determine local coordinates of a point within a lattice element Parameters ---------- point : Iterable of float Cartesian coordinates of point idx : Iterable of int - Indices of lattice element in (x,:math:`alpha`,z) bases + Indices of lattice element in :math:`(x,\alpha,z)` bases Returns ------- @@ -920,17 +935,17 @@ class HexLattice(Lattice): system """ - x = point[0] - (self._center[0] + sqrt(0.75)*self._pitch[0]*idx[0]) - y = point[1] - (self._center[1] + (0.5*idx[0] + idx[1])*self._pitch[0]) + x = point[0] - (self.center[0] + sqrt(0.75)*self.pitch[0]*idx[0]) + y = point[1] - (self.center[1] + (0.5*idx[0] + idx[1])*self.pitch[0]) if self._num_axial is None: z = point[2] else: - z = point[2] - (self._center[2] + (idx[2] + 0.5 - 0.5*self._num_axial)* - self._pitch[1]) + z = point[2] - (self.center[2] + (idx[2] + 0.5 - 0.5*self.num_axial)* + self.pitch[1]) return (x, y, z) def get_universe_index(self, idx): - """Return index in the universes array corresponding to a lattice element index + r"""Return index in the universes array corresponding to a lattice element index Parameters ---------- @@ -970,7 +985,7 @@ class HexLattice(Lattice): return (idx[2], i_ring, i_within) def is_valid_index(self, idx): - """Determine whether lattice element index is within defined range + r"""Determine whether lattice element index is within defined range Parameters ---------- diff --git a/openmc/universe.py b/openmc/universe.py index d8e1c4dac..a28729d8e 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -13,7 +13,6 @@ import openmc.checkvalue as cv if sys.version_info[0] >= 3: basestring = str - # A dictionary for storing IDs of cell elements that have already been written, # used to optimize the writing process WRITTEN_IDS = {} @@ -156,7 +155,7 @@ class Universe(object): return [] def plot(self, center=(0., 0., 0.), width=(1., 1.), pixels=(200, 200), - basis='xy', color_by='cell'): + basis='xy', color_by='cell', seed=None): """Display a slice plot of the universe. Parameters @@ -171,10 +170,18 @@ class Universe(object): The basis directions for the plot color_by : {'cell', 'material'} Indicate whether the plot should be colored by cell or by material + seed : hashable object or None + Hashable object which is used to seed the random number generator + used to select colors. If None, the generator is seeded from the + current time. """ import matplotlib.pyplot as plt + # Seed the random number generator + if seed is not None: + random.seed(seed) + if basis == 'xy': x_min = center[0] - 0.5*width[0] x_max = center[0] + 0.5*width[0] From 0ee6a6169cb05ab3d79f57766533a9c0a69e57b1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 31 May 2016 06:47:27 -0500 Subject: [PATCH 8/9] Fix typo --- openmc/lattice.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/lattice.py b/openmc/lattice.py index 0277d50cb..7690104c7 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -654,7 +654,7 @@ class HexLattice(Lattice): are assigned an index corresponding to their position relative to skewed :math:`(x,\alpha,z)` axes as described fully in :ref:`hexagonal_indexing`. However, note that when universes are assigned to - lattice elements using the :attr:`RectLattice.universes` property, the array + lattice elements using the :attr:`HexLattice.universes` property, the array indices do not correspond to natural indices. Parameters From 2b219ae15f8f1a581e808f7d664f522353ba89cc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 31 May 2016 08:41:05 -0500 Subject: [PATCH 9/9] Refactor handling of cell fill types. Fix Geometry.get_* methods --- openmc/cell.py | 149 ++++++++++------------- openmc/geometry.py | 55 ++++----- openmc/summary.py | 4 +- tests/test_distribmat/results_true.dat | 2 +- tests/test_distribmat/test_distribmat.py | 2 +- 5 files changed, 92 insertions(+), 120 deletions(-) diff --git a/openmc/cell.py b/openmc/cell.py index 29f85754a..900c952e5 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -36,7 +36,7 @@ class Cell(object): automatically be assigned. name : str, optional Name of the cell. If not specified, the name is the empty string. - fill : openmc.Material or openmc.Universe or openmc.Lattice or 'void' or iterable of openmc.Material, optional + fill : openmc.Material or openmc.Universe or openmc.Lattice or None or iterable of openmc.Material, optional Indicates what the region of space is filled with region : openmc.Region, optional Region of space that is assigned to the cell. @@ -47,9 +47,13 @@ class Cell(object): Unique identifier for the cell name : str Name of the cell - fill : openmc.Material or openmc.Universe or openmc.Lattice or 'void' or iterable of openmc.Material - Indicates what the region of space is filled with - region : openmc.Region + fill : openmc.Material or openmc.Universe or openmc.Lattice or None or iterable of openmc.Material + Indicates what the region of space is filled with. If None, the cell is + treated as a void. An iterable of materials is used to fill repeated + instances of a cell with different materials. + fill_type : {'material', 'universe', 'lattice', 'distribmat', 'void'} + Indicates what the cell is filled with. + region : openmc.Region or None Region of space that is assigned to the cell. rotation : Iterable of float If the cell is filled with a universe, this array specifies the angles @@ -68,6 +72,9 @@ class Cell(object): \sin\theta \sin\psi \\ -\sin\theta & \sin\phi \cos\theta & \cos\phi \cos\theta \end{array} \right ] + rotation_matrix : numpy.ndarray + The rotation matrix defined by the angles specified in the + :attr:`Cell.rotation` property. translation : Iterable of float If the cell is filled with a universe, this array specifies a vector that is used to translate (shift) the universe. @@ -82,20 +89,14 @@ class Cell(object): # Initialize Cell class attributes self.id = cell_id self.name = name - self._fill = None - self._type = None - self._region = None + self.fill = fill + self.region = region self._rotation = None self._rotation_matrix = None self._translation = None self._offsets = None self._distribcell_index = None - if fill is not None: - self.fill = fill - if region is not None: - self.region = region - def __contains__(self, point): if self.region is None: return True @@ -128,35 +129,24 @@ class Cell(object): def __repr__(self): string = 'Cell\n' - string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) - string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) + string += '{: <16}=\t{}\n'.format('\tID', self.id) + string += '{: <16}=\t{}\n'.format('\tName', self.name) - if isinstance(self._fill, openmc.Material): - string += '{0: <16}{1}{2}\n'.format('\tMaterial', '=\t', - self._fill._id) - elif isinstance(self._fill, basestring): - string += '{0: <16}=\tvoid\n'.format('\tMaterial') - elif isinstance(self._fill, Iterable): - string += '{0: <16}{1}'.format('\tMaterial', '=\t') - string += '[' - string += ', '.join(['void' if m == 'void' else str(m.id) - for m in self.fill]) - string += ']\n' - elif isinstance(self._fill, (openmc.Universe, openmc.Lattice)): - string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t', - self._fill._id) + if self.fill_type == 'material': + string += '{: <16}=\tMaterial {}\n'.format('\tFill', self.fill.id) + elif self.fill_type == 'void': + string += '{: <16}=\tNone\n'.format('\tFill') + elif self.fill_type == 'distribmat': + string += '{: <16}=\t{}\n'.format('\tFill', list(map( + lambda m: m if m is None else m.id, self.fill))) else: - string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t', self._fill) + string += '{: <16}=\t{}\n'.format('\tFill', self.fill.id) - string += '{0: <16}{1}{2}\n'.format('\tRegion', '=\t', self._region) - - string += '{0: <16}{1}{2}\n'.format('\tRotation', '=\t', - self._rotation) - string += '{0: <16}{1}{2}\n'.format('\tTranslation', '=\t', - self._translation) - string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self._offsets) - string += '{0: <16}{1}{2}\n'.format('\tDistribcell index', '=\t', - self._distribcell_index) + string += '{: <16}=\t{}\n'.format('\tRegion', self.region) + string += '{: <16}=\t{}\n'.format('\tRotation', self.rotation) + string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation) + string += '{: <16}=\t{}\n'.format('\tOffset', self.offsets) + string += '{: <16}=\t{}\n'.format('\tDistribcell index', self.distribcell_index) return string @@ -180,8 +170,10 @@ class Cell(object): return 'universe' elif isinstance(self.fill, openmc.Lattice): return 'lattice' + elif isinstance(self.fill, Iterable): + return 'distribmat' else: - return None + return 'void' @property def region(self): @@ -228,33 +220,25 @@ class Cell(object): @fill.setter def fill(self, fill): - if isinstance(fill, basestring): - if fill.strip().lower() == 'void': - self._type = 'void' - else: + if fill is not None: + if isinstance(fill, basestring): + if fill.strip().lower() != 'void': + msg = 'Unable to set Cell ID="{0}" to use a non-Material ' \ + 'or Universe fill "{1}"'.format(self._id, fill) + raise ValueError(msg) + fill = None + + elif isinstance(fill, Iterable): + for i, f in enumerate(fill): + if f is not None: + cv.check_type('cell.fill[i]', f, openmc.Material) + + elif not isinstance(fill, (openmc.Material, openmc.Lattice, + openmc.Universe)): msg = 'Unable to set Cell ID="{0}" to use a non-Material or ' \ - 'Universe fill "{1}"'.format(self._id, fill) + 'Universe fill "{1}"'.format(self._id, fill) raise ValueError(msg) - elif isinstance(fill, openmc.Material): - self._type = 'normal' - - elif isinstance(fill, Iterable): - cv.check_type('cell.fill', fill, Iterable, - (openmc.Material, basestring)) - self._type = 'normal' - - elif isinstance(fill, openmc.Universe): - self._type = 'fill' - - elif isinstance(fill, openmc.Lattice): - self._type = 'lattice' - - else: - msg = 'Unable to set Cell ID="{0}" to use a non-Material or ' \ - 'Universe fill "{1}"'.format(self._id, fill) - raise ValueError(msg) - self._fill = fill @rotation.setter @@ -290,7 +274,8 @@ class Cell(object): @region.setter def region(self, region): - cv.check_type('cell region', region, Region) + if region is not None: + cv.check_type('cell region', region, Region) self._region = region @distribcell_index.setter @@ -345,11 +330,11 @@ class Cell(object): def get_cell_instance(self, path, distribcell_index): # If the Cell is filled by a Material - if self._type == 'normal' or self._type == 'void': + if self.fill_type in ('material', 'distribmat', 'void'): offset = 0 # If the Cell is filled by a Universe - elif self._type == 'fill': + elif self.fill_type == 'universe': offset = self.offsets[distribcell_index-1] offset += self.fill.get_cell_instance(path, distribcell_index) @@ -372,8 +357,8 @@ class Cell(object): nuclides = OrderedDict() - if self._type != 'void': - nuclides.update(self._fill.get_all_nuclides()) + if self.fill_type != 'void': + nuclides.update(self.fill.get_all_nuclides()) return nuclides @@ -391,8 +376,8 @@ class Cell(object): cells = OrderedDict() - if self._type == 'fill' or self._type == 'lattice': - cells.update(self._fill.get_all_cells()) + if self.fill_type in ('universe', 'lattice'): + cells.update(self.fill.get_all_cells()) return cells @@ -432,11 +417,11 @@ class Cell(object): universes = OrderedDict() - if self._type == 'fill': - universes[self._fill._id] = self._fill - universes.update(self._fill.get_all_universes()) - elif self._type == 'lattice': - universes.update(self._fill.get_all_universes()) + if self.fill_type == 'universe': + universes[self.fill.id] = self.fill + universes.update(self.fill.get_all_universes()) + elif self.fill_type == 'lattice': + universes.update(self.fill.get_all_universes()) return universes @@ -447,24 +432,20 @@ class Cell(object): if len(self._name) > 0: element.set("name", str(self.name)) - if isinstance(self.fill, basestring): + if self.fill_type == 'void': element.set("material", "void") - elif isinstance(self.fill, openmc.Material): + elif self.fill_type == 'material': element.set("material", str(self.fill.id)) - elif isinstance(self.fill, Iterable): - element.set("material", ' '.join([m if m == 'void' else str(m.id) + elif self.fill_type == 'distribmat': + element.set("material", ' '.join(['void' if m is None else str(m.id) for m in self.fill])) - elif isinstance(self.fill, (openmc.Universe, openmc.Lattice)): + elif self.fill_type in ('universe', 'lattice'): element.set("fill", str(self.fill.id)) self.fill.create_xml_subelement(xml_element) - else: - element.set("fill", str(self.fill)) - self.fill.create_xml_subelement(xml_element) - if self.region is not None: # Set the region attribute with the region specification element.set("region", str(self.region)) diff --git a/openmc/geometry.py b/openmc/geometry.py index 006151900..2625ed3fa 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -144,14 +144,8 @@ class Geometry(object): """ - all_cells = self._root_universe.get_all_cells() - cells = set() - - for cell in all_cells.values(): - if cell._type == 'normal': - cells.add(cell) - - cells = list(cells) + all_cells = self.root_universe.get_all_cells() + cells = list(set(all_cells.values())) cells.sort(key=lambda x: x.id) return cells @@ -166,12 +160,7 @@ class Geometry(object): """ all_universes = self._root_universe.get_all_universes() - universes = set() - - for universe in all_universes.values(): - universes.add(universe) - - universes = list(universes) + universes = list(set(all_universes.values())) universes.sort(key=lambda x: x.id) return universes @@ -204,15 +193,17 @@ class Geometry(object): """ material_cells = self.get_all_material_cells() - materials = set() + materials = [] for cell in material_cells: - if isinstance(cell.fill, Iterable): - for m in cell.fill: materials.add(m) - else: - materials.add(cell.fill) + if cell.fill_type == 'distribmat': + for m in cell.fill: + if m is not None and m not in materials: + materials.append(m) + elif cell.fill_type == 'material': + if cell.fill not in materials: + materials.append(cell.fill) - materials = list(materials) materials.sort(key=lambda x: x.id) return materials @@ -227,13 +218,13 @@ class Geometry(object): """ all_cells = self.get_all_cells() - material_cells = set() + material_cells = [] for cell in all_cells: - if cell._type == 'normal': - material_cells.add(cell) + if cell.fill_type in ('material', 'distribmat'): + if cell not in material_cells: + material_cells.append(cell) - material_cells = list(material_cells) material_cells.sort(key=lambda x: x.id) return material_cells @@ -248,15 +239,15 @@ class Geometry(object): """ all_universes = self.get_all_universes() - material_universes = set() + material_universes = [] for universe in all_universes: cells = universe.cells for cell in cells: - if cell._type == 'normal': - material_universes.add(universe) + if cell.fill_type in ('material', 'distribmat', 'void'): + if universe not in material_universes: + material_universes.append(universe) - material_universes = list(material_universes) material_universes.sort(key=lambda x: x.id) return material_universes @@ -271,13 +262,13 @@ class Geometry(object): """ cells = self.get_all_cells() - lattices = set() + lattices = [] for cell in cells: - if isinstance(cell.fill, openmc.Lattice): - lattices.add(cell.fill) + if cell.fill_type == 'lattice': + if cell.fill not in lattices: + lattices.append(cell.fill) - lattices = list(lattices) lattices.sort(key=lambda x: x.id) return lattices diff --git a/openmc/summary.py b/openmc/summary.py index c3277809a..2af9b6be6 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -493,13 +493,13 @@ class Summary(object): # Retrieve the object corresponding to the fill type and ID if fill_type == 'normal': if isinstance(fill_id, Iterable): - fill = [self.get_material_by_id(mat) if mat > 0 else 'void' + fill = [self.get_material_by_id(mat) if mat > 0 else None for mat in fill_id] else: if fill_id > 0: fill = self.get_material_by_id(fill_id) else: - fill = 'void' + fill = None elif fill_type == 'universe': fill = self.get_universe_by_id(fill_id) else: diff --git a/tests/test_distribmat/results_true.dat b/tests/test_distribmat/results_true.dat index 15a00ee7d..bf96784d7 100644 --- a/tests/test_distribmat/results_true.dat +++ b/tests/test_distribmat/results_true.dat @@ -3,7 +3,7 @@ k-combined: Cell ID = 11 Name = - Material = [2, 3, void, 2] + Fill = [2, 3, None, 2] Region = -10000 Rotation = None Translation = None diff --git a/tests/test_distribmat/test_distribmat.py b/tests/test_distribmat/test_distribmat.py index 6700a96b6..96d41c3fb 100644 --- a/tests/test_distribmat/test_distribmat.py +++ b/tests/test_distribmat/test_distribmat.py @@ -45,7 +45,7 @@ class DistribmatTestHarness(PyAPITestHarness): r0 = openmc.ZCylinder(R=0.3) c11 = openmc.Cell(cell_id=11) c11.region = -r0 - c11.fill = [dense_fuel, light_fuel, 'void', dense_fuel] + c11.fill = [dense_fuel, light_fuel, None, dense_fuel] c12 = openmc.Cell(cell_id=12) c12.region = +r0 c12.fill = moderator