mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Merge pull request #656 from paulromano/pyapi-geometry-improvements
Python API geometry enhancements
This commit is contained in:
Will Boyd
commit
a787ea8818
24 changed files with 1000 additions and 249 deletions
|
|
@ -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
|
||||
--------------------------
|
||||
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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"
|
||||
]
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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],
|
||||
|
|
@ -116,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()
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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],
|
||||
|
|
@ -106,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()
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
||||
|
||||
|
|
|
|||
173
openmc/cell.py
173
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
|
||||
|
|
@ -33,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.
|
||||
|
|
@ -44,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
|
||||
|
|
@ -65,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.
|
||||
|
|
@ -79,18 +89,19 @@ 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
|
||||
else:
|
||||
return point in self.region
|
||||
|
||||
def __eq__(self, other):
|
||||
if not isinstance(other, Cell):
|
||||
|
|
@ -118,33 +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, 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
|
||||
|
||||
|
|
@ -168,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):
|
||||
|
|
@ -179,6 +183,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
|
||||
|
|
@ -212,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
|
||||
|
|
@ -249,13 +249,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):
|
||||
|
|
@ -264,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
|
||||
|
|
@ -319,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)
|
||||
|
||||
|
|
@ -346,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
|
||||
|
||||
|
|
@ -365,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
|
||||
|
||||
|
|
@ -406,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
|
||||
|
||||
|
|
@ -421,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))
|
||||
|
|
|
|||
|
|
@ -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):
|
||||
|
|
@ -62,6 +69,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 coordinates of 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.
|
||||
|
||||
|
|
@ -120,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
|
||||
|
||||
|
|
@ -142,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
|
||||
|
||||
|
|
@ -180,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
|
||||
|
||||
|
|
@ -203,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
|
||||
|
||||
|
|
@ -224,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
|
||||
|
||||
|
|
@ -247,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
|
||||
|
||||
|
|
|
|||
|
|
@ -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,19 @@ 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.
|
||||
|
||||
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
|
||||
|
|
@ -253,12 +264,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 +271,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 +297,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 +305,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 +322,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 +342,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 +355,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 +363,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 +406,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 +427,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 +440,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 +590,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 +600,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 +620,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
|
||||
|
|
@ -490,7 +644,18 @@ 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:`HexLattice.universes` property, the array
|
||||
indices do not correspond to natural indices.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
|
|
@ -506,26 +671,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 +767,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 +793,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 +818,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 +833,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 +874,170 @@ class HexLattice(Lattice):
|
|||
6*(self._num_rings - 1 - r))
|
||||
raise ValueError(msg)
|
||||
|
||||
def find_element(self, point):
|
||||
r"""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 :math:`(x,\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):
|
||||
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 :math:`(x,\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):
|
||||
r"""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):
|
||||
r"""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 +1069,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.
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 ''
|
||||
|
|
@ -219,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
|
||||
|
|
@ -229,6 +233,26 @@ 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.
|
||||
|
||||
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 +305,26 @@ 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.
|
||||
|
||||
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 +380,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)
|
||||
|
||||
|
|
|
|||
|
|
@ -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([])
|
||||
|
||||
|
|
@ -496,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:
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
@ -12,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 = {}
|
||||
|
|
@ -124,6 +124,118 @@ 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', seed=None):
|
||||
"""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
|
||||
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]
|
||||
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.
|
||||
|
||||
|
|
|
|||
|
|
@ -354,7 +354,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 = [
|
||||
|
|
@ -388,7 +387,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 = [
|
||||
|
|
@ -448,7 +446,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 = [
|
||||
|
|
@ -476,7 +473,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 = [
|
||||
|
|
|
|||
|
|
@ -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],
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue