crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #646 from paulromano/lattice-improvements

Browse source 
Improve clarity of constructing hexagonal lattices
This commit is contained in:
Sterling Harper 2016-05-14 21:22:38 -04:00
commit cebc7a3494
8 changed files with 214 additions and 18 deletions
Download patch file Download diff file Expand all files Collapse all files

10
docs/source/pythonapi/index.rst
View file

@ -122,6 +122,16 @@ Many of the above classes are derived from several abstract classes:
openmc.Region
openmc.Lattice
One function is also available to create a hexagonal region defined by the
intersection of six surface half-spaces.
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.make_hexagon_region
Constructing Tallies
--------------------

14
docs/source/usersguide/input.rst
View file

@ -1033,6 +1033,20 @@ Each ``<cell>`` element can have the following attributes or sub-elements:
<cell fill="..." rotation="0 0 90" />
The rotation applied is an intrinsic rotation whose Tait-Bryan angles are
given as those specified about the x, y, and z axes respectively. That is to
say, if the angles are :math:`(\phi, \theta, \psi)`, then the rotation
matrix applied is :math:`R_z(\psi) R_y(\theta) R_x(\phi)` or
.. math::
\left [ \begin{array}{ccc} \cos\theta \cos\psi & -\cos\theta \sin\psi +
\sin\phi \sin\theta \cos\psi & \sin\phi \sin\psi + \cos\phi \sin\theta
\cos\psi \\ \cos\theta \sin\psi & \cos\phi \cos\psi + \sin\phi \sin\theta
\sin\psi & -\sin\phi \cos\psi + \cos\phi \sin\theta \sin\psi \\
-\sin\theta & \sin\phi \cos\theta & \cos\phi \cos\theta \end{array}
\right ]
*Default*: None
:translation:

25
openmc/cell.py
View file

@ -23,7 +23,7 @@ def reset_auto_cell_id():
class Cell(object):
"""A region of space defined as the intersection of half-space created by
r"""A region of space defined as the intersection of half-space created by
quadric surfaces.
Parameters
@ -48,11 +48,24 @@ class Cell(object):
Indicates what the region of space is filled with
region : openmc.Region
Region of space that is assigned to the cell.
rotation : numpy.ndarray
rotation : Iterable of float
If the cell is filled with a universe, this array specifies the angles
in degrees about the x, y, and z axes that the filled universe should be
rotated.
translation : numpy.ndarray
rotated. The rotation applied is an intrinsic rotation with specified
Tait-Bryan angles. That is to say, if the angles are :math:`(\phi,
\theta, \psi)`, then the rotation matrix applied is :math:`R_z(\psi)
R_y(\theta) R_x(\phi)` or
.. math::
\left [ \begin{array}{ccc} \cos\theta \cos\psi & -\cos\theta \sin\psi
+ \sin\phi \sin\theta \cos\psi & \sin\phi \sin\psi + \cos\phi
\sin\theta \cos\psi \\ \cos\theta \sin\psi & \cos\phi \cos\psi +
\sin\phi \sin\theta \sin\psi & -\sin\phi \cos\psi + \cos\phi
\sin\theta \sin\psi \\ -\sin\theta & \sin\phi \cos\theta & \cos\phi
\cos\theta \end{array} \right ]
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.
offsets : ndarray
@ -230,6 +243,10 @@ class Cell(object):
@rotation.setter
def rotation(self, rotation):
if not isinstance(self.fill, openmc.Universe):
raise RuntimeError('Cell rotation can only be applied if the cell '
'is filled with a Universe')
cv.check_type('cell rotation', rotation, Iterable, Real)
cv.check_length('cell rotation', rotation, 3)
self._rotation = rotation

5
openmc/element.py
View file

@ -1,7 +1,7 @@
import sys
import openmc
from openmc.checkvalue import check_type
from openmc.checkvalue import check_type, check_length
from openmc.data import natural_abundance
if sys.version_info[0] >= 3:
@ -99,7 +99,8 @@ class Element(object):
@name.setter
def name(self, name):
check_type('name', name, basestring)
check_type('element name', name, basestring)
check_length('element name', name, 1, 2)
self._name = name
@scattering.setter

125
openmc/lattice.py
View file

@ -30,8 +30,8 @@ class Lattice(object):
Unique identifier for the lattice
name : str
Name of the lattice
pitch : float
Pitch of the lattice in cm
pitch : Iterable of float
Pitch of the lattice in each direction in cm
outer : openmc.Universe
A universe to fill all space outside the lattice
universes : Iterable of Iterable of openmc.Universe
@ -259,8 +259,9 @@ class RectLattice(Lattice):
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 : float
Pitch of the lattice in cm
pitch : Iterable of float
Pitch of the lattice in the x, y, and (if applicable) z directions in
cm.
outer : openmc.Universe
A universe to fill all space outside the lattice
universes : Iterable of Iterable of openmc.Universe
@ -512,13 +513,19 @@ class HexLattice(Lattice):
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 : float
Pitch of the lattice in cm
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
A two- or three-dimensional list/array of universes filling each element
of the lattice
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.
"""
@ -882,3 +889,107 @@ class HexLattice(Lattice):
# Join the rows together and return the string.
universe_ids = '\n'.join(rows)
return universe_ids
@staticmethod
def show_indices(num_rings):
"""Return a diagram of the hexagonal lattice layout with indices.
This method can be used to show the proper indices to be used when
setting the :attr:`HexLattice.universes` property. For example, running
this method with num_rings=3 will return the following diagram::
(0, 0)
(0,11) (0, 1)
(0,10) (1, 0) (0, 2)
(1, 5) (1, 1)
(0, 9) (2, 0) (0, 3)
(1, 4) (1, 2)
(0, 8) (1, 3) (0, 4)
(0, 7) (0, 5)
(0, 6)
Parameters
----------
num_rings : int
Number of rings in the hexagonal lattice
Returns
-------
str
Diagram of the hexagonal lattice showing indices
"""
# Find the largest string and count the number of digits so we can
# properly pad the output string later
largest_index = 6*(num_rings - 1)
n_digits_index = len(str(largest_index))
n_digits_ring = len(str(num_rings - 1))
str_form = '({{:{}}},{{:{}}})'.format(n_digits_ring, n_digits_index)
pad = ' '*(n_digits_index + n_digits_ring + 3)
# Initialize the list for each row.
rows = [[] for i in range(1 + 4 * (num_rings-1))]
middle = 2 * (num_rings - 1)
# Start with the degenerate first ring.
rows[middle] = [str_form.format(num_rings - 1, 0)]
# Add universes one ring at a time.
for r in range(1, num_rings):
# r_prime increments down while r increments up.
r_prime = num_rings - 1 - r
theta = 0
y = middle + 2*r
for i in range(r):
# Climb down the top-right.
rows[y].append(str_form.format(r_prime, theta))
y -= 1
theta += 1
for i in range(r):
# Climb down the right.
rows[y].append(str_form.format(r_prime, theta))
y -= 2
theta += 1
for i in range(r):
# Climb down the bottom-right.
rows[y].append(str_form.format(r_prime, theta))
y -= 1
theta += 1
for i in range(r):
# Climb up the bottom-left.
rows[y].insert(0, str_form.format(r_prime, theta))
y += 1
theta += 1
for i in range(r):
# Climb up the left.
rows[y].insert(0, str_form.format(r_prime, theta))
y += 2
theta += 1
for i in range(r):
# Climb up the top-left.
rows[y].insert(0, str_form.format(r_prime, theta))
y += 1
theta += 1
# Flip the rows and join each row into a single string.
rows = [pad.join(x) for x in rows[::-1]]
# Pad the beginning of the rows so they line up properly.
for y in range(num_rings - 1):
rows[y] = (num_rings - 1 - y)*pad + rows[y]
rows[-1 - y] = (num_rings - 1 - y)*pad + rows[-1 - y]
for y in range(num_rings % 2, num_rings, 2):
rows[middle + y] = pad + rows[middle + y]
if y != 0:
rows[middle - y] = pad + rows[middle - y]
# Join the rows together and return the string.
return '\n'.join(rows)

6
openmc/stats/multivariate.py
View file

@ -328,8 +328,8 @@ class Point(Spatial):
Parameters
----------
xyz : Iterable of float
Cartesian coordinates of location
xyz : Iterable of float, optional
Cartesian coordinates of location. Defaults to (0., 0., 0.).
Attributes
----------
@ -338,7 +338,7 @@ class Point(Spatial):
"""
def __init__(self, xyz):
def __init__(self, xyz=(0., 0., 0.)):
super(Point, self).__init__()
self.xyz = xyz

2
openmc/summary.py
View file

@ -267,7 +267,7 @@ class Summary(object):
rotation = \
self._f['geometry/cells'][key]['rotation'][...]
rotation = np.asarray(rotation, dtype=np.int)
cell.rotation = rotation
cell._rotation = rotation
# Store Cell fill information for after Universe/Lattice creation
self._cell_fills[index] = (fill_type, fill)

45
openmc/surface.py
View file

@ -2,11 +2,12 @@ from abc import ABCMeta
from numbers import Real, Integral
from xml.etree import ElementTree as ET
import sys
from math import sqrt
import numpy as np
from openmc.checkvalue import check_type, check_value, check_greater_than
from openmc.region import Region
from openmc.region import Region, Intersection
if sys.version_info[0] >= 3:
basestring = str
@ -1503,3 +1504,45 @@ class Halfspace(Region):
def __str__(self):
return '-' + str(self.surface.id) if self.side == '-' \
else str(self.surface.id)
def make_hexagon_region(edge_length=1., orientation='y'):
"""Create a hexagon region from six surface planes.
Parameters
----------
edge_length : float
Length of a side of the hexagon in cm
orientation : {'x', 'y'}
An 'x' orientation means that two sides of the hexagon are parallel to
the x-axis and a 'y' orientation means that two sides of the hexagon are
parallel to the y-axis.
Returns
-------
openmc.Region
The inside of a hexagonal prism
"""
l = edge_length
if orientation == 'x':
right = XPlane(x0=sqrt(3.)/2.*l)
left = XPlane(x0=-sqrt(3.)/2.*l)
c = sqrt(3.)/3.
ur = Plane(A=c, B=1., D=l) # y = -x/sqrt(3) + a
ul = Plane(A=-c, B=1., D=l) # y = x/sqrt(3) + a
lr = Plane(A=-c, B=1., D=-l) # y = x/sqrt(3) - a
ll = Plane(A=c, B=1., D=-l) # y = -x/sqrt(3) - a
return Intersection(-right, +left, -ur, -ul, +lr, +ll)
elif orientation == 'y':
top = YPlane(y0=sqrt(3.)/2.*l)
bottom = YPlane(y0=-sqrt(3.)/2.*l)
c = sqrt(3.)
ur = Plane(A=c, B=1., D=c*l) # y = -sqrt(3)*(x - a)
lr = Plane(A=-c, B=1., D=-c*l) # y = sqrt(3)*(x + a)
ll = Plane(A=c, B=1., D=-c*l) # y = -sqrt(3)*(x + a)
ul = Plane(A=-c, B=1., D=c*l) # y = sqrt(3)*(x + a)
return Intersection(-top, +bottom, -ur, +lr, +ll, -ul)