Implemented Python API to auto-generate XML input files

src/utils/openmc/__init__.py

@@ -0,0 +1,14 @@
+#!/usr/bin/env python
+
+from openmc.element import *
+from openmc.geometry import *
+from openmc.nuclide import *
+from openmc.material import *
+from openmc.opencsg_compatible import *
+from openmc.plots import *
+from openmc.settings import *
+from openmc.surface import *
+from openmc.universe import *
+from openmc.tallies import *
+from openmc.executor import *
+#from statepoint import *

src/utils/openmc/checkvalue.py

@@ -0,0 +1,16 @@
+#!/usr/bin/env python
+
+import numpy as np
+
+
+def is_integer(val):
+  return isinstance(val, (int, np.int32, np.int64))
+
+
+def is_float(val):
+  return isinstance(val, (float, np.float32, np.float64))
+
+
+def is_string(val):
+  return isinstance(val, (str, np.str))
+

src/utils/openmc/clean_xml.py

@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+
+def sort_xml_elements(tree):
+
+  # Retrieve all children of the root XML node in the tree
+  elements = tree.getchildren()
+
+  # Initialize empty lists for the sorted and comment elements
+  sorted_elements = list()
+
+  # Initialize an empty set of tags (e.g., Surface, Cell, and Lattice)
+  tags = set()
+
+  # Find the unique tags in the tree
+  for element in elements:
+    tags.add(element.tag)
+
+  # Initialize an empty list for the comment elements
+  comment_elements = list()
+
+  # Find the comment elements and record their ordering within the
+  # tree using a precedence with respect to the subsequent nodes
+  for index, element in enumerate(elements):
+    next_element = None
+
+    if 'Comment' in str(element.tag):
+
+      if index < len(elements)-1:
+        next_element = elements[index+1]
+
+      comment_elements.append((element, next_element))
+
+  # Now iterate over all tags and order the elements within each tag
+  for tag in tags:
+
+    # Retrieve all of the elements for this tag
+    try:
+      tagged_elements = tree.findall(tag)
+    except:
+      continue
+
+    # Initialize an empty list of tuples to sort (id, element)
+    tagged_data = list()
+
+    # Retrieve the IDs for each of the elements
+    for element in tagged_elements:
+      key = element.get('id')
+
+      # If this element has an "ID" tag, append it to the list to sort
+      if not key is None:
+        tagged_data.append((key, element))
+
+    # Sort the elements according to the IDs for this tag
+    tagged_data.sort()
+    sorted_elements.extend(list(item[-1] for item in tagged_data))
+
+  # Add the comment elements while preserving the original precedence
+  for element, next_element in comment_elements:
+    index = sorted_elements.index(next_element)
+    sorted_elements.insert(index, element)
+
+  # Remove all of the sorted elements from the tree
+  for element in sorted_elements:
+    tree.remove(element)
+
+  # Add the sorted elements back to the tree in the proper order
+  tree.extend(sorted_elements)
+
+
+def clean_xml_indentation(element, level=0):
+  '''
+  copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint
+  it basically walks your tree and adds spaces and newlines so the tree is
+  printed in a nice way
+  '''
+
+  i = "\n" + level*"  "
+
+  if len(element):
+
+    if not element.text or not element.text.strip():
+      element.text = i + "  "
+
+    if not element.tail or not element.tail.strip():
+      element.tail = i
+
+    for element in element:
+      clean_xml_indentation(element, level+1)
+
+      if not element.tail or not element.tail.strip():
+        element.tail = i
+
+  else:
+    if level and (not element.tail or not element.tail.strip()):
+      element.tail = i

src/utils/openmc/constants.py

@@ -0,0 +1,125 @@
+#!/usr/bin/env python
+
+"""Dictionaries of integer-to-string mapppings from openmc/src/constants.F90"""
+
+SURFACE_TYPES = {1: 'x-plane',
+                 2: 'y-plane',
+                 3: 'z-plane',
+                 4: 'plane',
+                 5: 'x-cylinder',
+                 6: 'y-cylinder',
+                 7: 'z-cylinder',
+                 8: 'sphere',
+                 9: 'x-cone',
+                 10: 'y-cone',
+                 11: 'z-cone'}
+
+BC_TYPES = {0: 'transmission',
+            1: 'vacuum',
+            2: 'reflective',
+            3: 'periodic'}
+
+FILL_TYPES = {1: 'normal',
+              2: 'fill',
+              3: 'lattice'}
+
+LATTICE_TYPES = {1: 'rectangular',
+                 2: 'hexagonal'}
+
+ESTIMATOR_TYPES = {1: 'analog',
+                   2: 'tracklength'}
+
+FILTER_TYPES = {1: 'universe',
+                2: 'material',
+                3: 'cell',
+                4: 'cellborn',
+                5: 'surface',
+                6: 'mesh',
+                7: 'energy',
+                8: 'energyout',
+                9: 'distribcell'}
+
+SCORE_TYPES = {-1: 'flux',
+               -2: 'total',
+               -3: 'scatter',
+               -4: 'nu-scatter',
+               -5: 'scatter-n',
+               -6: 'scatter-pn',
+               -7: 'nu-scatter-n',
+               -8: 'nu-scatter-pn',
+               -9: 'transport',
+               -10: 'n1n',
+               -11: 'absorption',
+               -12: 'fission',
+               -13: 'nu-fission',
+               -14: 'kappa-fission',
+               -15: 'current',
+               -16: 'flux-yn',
+               -17: 'total-yn',
+               -18: 'scatter-yn',
+               -19: 'nu-scatter-yn',
+               -20: 'events',
+               1: '(n,total)',
+               2: '(n,elastic)',
+               4: '(n,level)',
+               11: '(n,2nd)',
+               16: '(n,2n)',
+               17: '(n,3n)',
+               18: '(n,fission)',
+               19: '(n,f)',
+               20: '(n,nf)',
+               21: '(n,2nf)',
+               22: '(n,na)',
+               23: '(n,n3a)',
+               24: '(n,2na)',
+               25: '(n,3na)',
+               28: '(n,np)',
+               29: '(n,n2a)',
+               30: '(n,2n2a)',
+               32: '(n,nd)',
+               33: '(n,nt)',
+               34: '(n,nHe-3)',
+               35: '(n,nd2a)',
+               36: '(n,nt2a)',
+               37: '(n,4n)',
+               38: '(n,3nf)',
+               41: '(n,2np)',
+               42: '(n,3np)',
+               44: '(n,n2p)',
+               45: '(n,npa)',
+               91: '(n,nc)',
+               101: '(n,disappear)',
+               102: '(n,gamma)',
+               103: '(n,p)',
+               104: '(n,d)',
+               105: '(n,t)',
+               106: '(n,3He)',
+               107: '(n,a)',
+               108: '(n,2a)',
+               109: '(n,3a)',
+               111: '(n,2p)',
+               112: '(n,pa)',
+               113: '(n,t2a)',
+               114: '(n,d2a)',
+               115: '(n,pd)',
+               116: '(n,pt)',
+               117: '(n,da)',
+               201: '(n,Xn)',
+               202: '(n,Xgamma)',
+               203: '(n,Xp)',
+               204: '(n,Xd)',
+               205: '(n,Xt)',
+               206: '(n,X3He)',
+               207: '(n,Xa)',
+               444: '(damage)',
+               649: '(n,pc)',
+               699: '(n,dc)',
+               749: '(n,tc)',
+               799: '(n,3Hec)',
+               849: '(n,tc)'}
+SCORE_TYPES.update({MT: '(n,n' + str(MT-50) + ')' for MT in range(51,91)})
+SCORE_TYPES.update({MT: '(n,p' + str(MT-600) + ')' for MT in range(600,649)})
+SCORE_TYPES.update({MT: '(n,d' + str(MT-650) + ')' for MT in range(650,699)})
+SCORE_TYPES.update({MT: '(n,t' + str(MT-700) + ')' for MT in range(700,749)})
+SCORE_TYPES.update({MT: '(n,3He' + str(MT-750) + ')' for MT in range(750,649)})
+SCORE_TYPES.update({MT: '(n,a' + str(MT-800) + ')' for MT in range(800,849)})

src/utils/openmc/element.py

@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+
+from openmc.checkvalue import *
+
+class Element(object):
+
+  def __init__(self, name='', xs=None):
+
+    # Initialize class attributes
+    self._name = ''
+    self._xs = None
+
+    # Set the Material class attributes
+    self.setName(name)
+
+    if not xs is None:
+      self.setXS(xs)
+
+
+  def __eq__(self, element2):
+
+    # check type
+    if not isinstance(element2, Element):
+      return False
+
+    # Check name
+    if self._name != element2._name:
+      return False
+
+    # Check xs
+    elif self._xs != element2._xs:
+      return False
+
+    else:
+      return True
+
+
+  def __hash__(self):
+    hashable = list()
+    hashable.append(self._name)
+    hashable.append(self._xs)
+    return hash(tuple(hashable))
+
+
+  def setName(self, name):
+
+    if not is_string(name):
+      msg = 'Unable to set name for Element with a non-string ' \
+            'value {0}'.format(name)
+      raise ValueError(msg)
+
+    self._name = name
+
+
+  def setXS(self, xs):
+
+    if not is_string(xs):
+      msg = 'Unable to set cross-section identifier xs for Element with ' \
+            'a non-string value {0}'.format(xs)
+      raise ValueError(msg)
+
+    self._xs = xs
+
+
+  def __repr__(self):
+
+    string = 'Element  -  {0}\n'.format(self._name)
+    string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs)
+    return string

src/utils/openmc/executor.py

@@ -0,0 +1,74 @@
+#!/usr/bin/env python
+
+from openmc.checkvalue import *
+import subprocess
+import os
+
+
+FNULL = open(os.devnull, 'w')
+
+
+class Executor(object):
+
+
+  def __init__(self):
+    self._working_directory = '.'
+
+
+  def setWorkingDirectory(self, working_directory):
+
+    if not is_string(working_directory):
+      msg = 'Unable to set Executor\'s working directory to {0} ' \
+            'since it is not a string'.format(working_directory)
+      raise ValueError(msg)
+
+    elif not os.path.isdir(working_directory):
+      msg = 'Unable to set Executor\'s working directory to {0} ' \
+            'which does not exist'.format(working_directory)
+      raise ValueError(msg)
+
+    self._working_directory = working_directory
+
+
+  def plotGeometry(self, output=True):
+
+    if output:
+      subprocess.check_call('openmc -p', shell=True,
+                            cwd=self._working_directory)
+    else:
+      subprocess.check_call('openmc -p', shell=True, stdout=FNULL,
+                            cwd=self._working_directory)
+
+
+  def runSimulation(self, particles=None, threads=None, geometry_debug=False,
+                    restart_file=None, tracks=False, mpi_procs=1, output=True):
+
+    post_args = ' '
+    pre_args = ''
+
+    if is_integer(particles) and particles > 0:
+      post_args += '-n {0} '.format(particles)
+
+    if is_integer(threads) and threads > 0:
+      post_args += '-s {0} '.format(threads)
+
+    if geometry_debug:
+      post_args += '-g '
+
+    if is_string(restart_file):
+      post_args += '-r {0} '.format(restart_file)
+
+    if tracks:
+      post_args += '-t'
+
+    if is_integer(mpi_procs) and mpi_procs > 1:
+      pre_args += 'mpirun -n {0} '.format(mpi_procs)
+
+    command = pre_args + 'openmc ' + post_args
+
+    if output:
+      subprocess.check_call(command, shell=True,
+                            cwd=self._working_directory)
+    else:
+      subprocess.check_call(command, shell=True, stdout=FNULL,
+                            cwd=self._working_directory)

src/utils/openmc/geometry.py

@@ -0,0 +1,160 @@
+#!/usr/bin/env python
+
+import openmc
+from openmc.clean_xml import *
+from xml.etree import ElementTree as ET
+
+def reset_auto_ids():
+  openmc.reset_auto_material_id()
+  openmc.reset_auto_surface_id()
+  openmc.reset_auto_cell_id()
+  openmc.reset_auto_universe_id()
+
+
+class Geometry(object):
+
+  def __init__(self):
+
+    # Initialize Geometry class attributes
+    self._root_universe = None
+    self._offsets = dict()
+
+
+  def getOffset(self, path, filter_offset):
+    """
+    Returns the corresponding location in the results array for a given
+    path and filter number.
+
+    Parameters
+    ----------
+    path : list
+        A list of IDs that form the path to the target. It should begin with 0
+        for the base universe, and should cover every universe, cell, and
+        lattice passed through. For the case of the lattice, a tuple should be
+        provided to indicate which coordinates in the lattice should be entered.
+        This should be in the form: (lat_id, i_x, i_y, i_z)
+
+    filter_offset : int
+        An integer that specifies which offset map the filter is using
+
+    """
+
+    # Return memoize'd offset if possible
+    if self._offsets.has_key((path, filter_offset)):
+      offset = self._offsets[(path, filter_offset)]
+
+    # Begin recursive call to compute the offset starting with the base Universe
+    else:
+      offset = self._root_universe.getOffset(path, filter_offset)
+      self._offsets[(path, filter_offset)] = offset
+
+    # Return the final offset
+    return offset
+
+
+  def getAllCells(self):
+    return self._root_universe.getAllCells()
+
+
+  def getAllUniverses(self):
+    return self._root_universe.getAllUniverses()
+
+
+  def getAllNuclides(self):
+
+    nuclides = dict()
+    materials = self.getAllMaterials()
+
+    for material in materials:
+      nuclides.update(material.getAllNuclides())
+
+    return nuclides
+
+
+  def getAllMaterials(self):
+
+    material_cells = self.getAllMaterialCells()
+    materials = set()
+
+    for cell in material_cells:
+      materials.add(cell._fill)
+
+    return list(materials)
+
+
+  def getAllMaterialCells(self):
+
+    all_cells = self.getAllCells()
+    material_cells = set()
+
+    for cell_id, cell in all_cells.items():
+      if cell._type == 'normal':
+        material_cells.add(cell)
+
+    return list(material_cells)
+
+
+  def getAllMaterialUniverses(self):
+
+    all_universes = self.getAllUniverses()
+    material_universes = set()
+
+    for universe_id, universe in all_universes.items():
+
+      cells = universe._cells
+
+      for cell_id, cell in cells.items():
+        if cell._type == 'normal':
+          material_universes.add(universe)
+
+    return list(material_universes)
+
+
+  def setRootUniverse(self, root_universe):
+
+    if not isinstance(root_universe, openmc.Universe):
+      msg = 'Unable to add root Universe {0} to Geometry since ' \
+            'it is not a Universe'.format(root_universe)
+      raise ValueError(msg)
+
+    elif root_universe._id != 0:
+      msg = 'Unable to add root Universe {0} to Geometry since it has ' \
+            'ID={1} instead of ID=0'.format(root_universe, root_universe._id)
+      raise ValueError(msg)
+
+    self._root_universe = root_universe
+
+
+
+class GeometryFile(object):
+
+  def __init__(self):
+
+    # Initialize GeometryFile class attributes
+    self._geometry = None
+    self._geometry_file = ET.Element("geometry")
+
+
+  def setGeometry(self, geometry):
+
+    if not isinstance(geometry, Geometry):
+      msg = 'Unable to set the Geometry to {0} for the GeometryFile ' \
+            'since it is not a Geometry object'.format(geometry)
+      raise ValueError(msg)
+
+    self._geometry = geometry
+
+
+  def exportToXML(self):
+
+    root_universe = self._geometry._root_universe
+    root_universe.createXMLSubElement(self._geometry_file)
+
+    # Clean the indentation in the file to be user-readable
+    sort_xml_elements(self._geometry_file)
+    clean_xml_indentation(self._geometry_file)
+
+    # Write the XML Tree to the materials.xml file
+    tree = ET.ElementTree(self._geometry_file)
+    tree.write("geometry.xml", xml_declaration=True,
+               encoding='utf-8', method="xml")

src/utils/openmc/material.py

@@ -0,0 +1,413 @@
+#!/usr/bin/env python
+
+import openmc
+from openmc.checkvalue import *
+from openmc.clean_xml import *
+from xml.etree import ElementTree as ET
+from collections import MappingView
+from copy import deepcopy
+import numpy as np
+
+
+# A list of all IDs for all Materials created
+MATERIAL_IDS = list()
+
+# A static variable for auto-generated Material IDs
+AUTO_MATERIAL_ID = 10000
+
+def reset_auto_material_id():
+  global AUTO_MATERIAL_ID, MATERIAL_IDS
+  AUTO_MATERIAL_ID = 10000
+  MATERIAL_IDS = list()
+
+
+# Units for density supported by OpenMC
+DENSITY_UNITS = ['g/cm3', 'kg/cm3', 'at/b-cm', 'at/cm3', 'sum']
+
+# ENDF temperatures
+ENDF_TEMPS = np.array([300, 600, 700, 900, 1200, 1500])
+
+# ENDF ZAIDs
+ENDF_ZAIDS = np.array(['70c', '71c', '72c', '73c', '74c'])
+
+# Constant for density when not needed
+NO_DENSITY = 99999.
+
+
+class Material(object):
+
+  def __init__(self, material_id=None, name=''):
+
+    # Initialize class attributes
+    self._id = None
+    self._name = ''
+    self._density = None
+    self._density_units = ''
+
+    # A dictionary of Nuclides
+    # Keys     - Nuclide names
+    # Values   - tuple (nuclide, percent, percent type)
+    self._nuclides = dict()
+
+    # A dictionary of Elements
+    # Keys     - Element names
+    # Values   - tuple (element, percent, percent type)
+    self._elements = dict()
+
+    # If specified, a list of tuples of (table name, xs identifier)
+    self._sab = list()
+
+    # Set the Material class attributes
+    self.setId(material_id)
+    self.setName(name)
+
+
+  def setId(self, material_id=None):
+
+    global MATERIAL_IDS
+
+    # If the Material already has an ID, remove it from global list
+    if not self._id is None:
+      MATERIAL_IDS.remove(self._id)
+
+    if material_id is None:
+      global AUTO_MATERIAL_ID
+      self._id = AUTO_MATERIAL_ID
+      MATERIAL_IDS.append(AUTO_MATERIAL_ID)
+      AUTO_MATERIAL_ID += 1
+
+    # Check that the ID is an integer and wasn't already used
+    elif not is_integer(material_id):
+      msg = 'Unable to set a non-integer Material ID {0}'.format(material_id)
+      raise ValueError(msg)
+
+    elif material_id in MATERIAL_IDS:
+      msg = 'Unable to set Material ID to {0} since a Material with ' \
+            'this ID was already initialized'.format(material_id)
+      raise ValueError(msg)
+
+    elif material_id < 0:
+      msg = 'Unable to set Material ID to {0} since it must be a ' \
+            'non-negative integer'.format(material_id)
+      raise ValueError(msg)
+
+    else:
+      self._id = material_id
+      MATERIAL_IDS.append(material_id)
+
+
+  def setName(self, name):
+
+    if not is_string(name):
+      msg = 'Unable to set name for Material ID={0} with a non-string ' \
+            'value {1}'.format(self._id, name)
+      raise ValueError(msg)
+
+    else:
+      self._name = name
+
+
+  def setDensity(self, units, density=NO_DENSITY):
+
+    if not is_float(density):
+      msg = 'Unable to set the density for Material ID={0} to a ' \
+            'non-floating point value {1}'.format(self._id, density)
+      raise ValueError(msg)
+
+    elif not units in DENSITY_UNITS:
+      msg = 'Unable to set the density for Material ID={0} with ' \
+            'units {1}'.format(self._id, units)
+      raise ValueError(msg)
+
+    if density == NO_DENSITY and units is not 'sum':
+      msg = 'Unable to set the density Material ID={0} because '\
+            'a density must be set when not using sum unit'.format(self._id)
+      raise ValueError(msg)
+
+    self._density = density
+    self._density_units = units
+
+
+  def addNuclide(self, nuclide, percent, percent_type='ao'):
+
+    if not isinstance(nuclide, openmc.Nuclide):
+      msg = 'Unable to add an Nuclide to Material ID={0} with a ' \
+            'non-Nuclide value {1}'.format(self._id, nuclide)
+      raise ValueError(msg)
+
+    elif not is_float(percent):
+      msg = 'Unable to add an Nuclide to Material ID={0} with a non-floating ' \
+            'non-floating point percent value {1}'.format(self._id, percent)
+      raise ValueError(msg)
+
+    elif not percent_type in ['ao', 'wo', 'at/g-cm']:
+      msg = 'Unable to add an Nuclide to Material ID={0} with a ' \
+            'percent type {1}'.format(self._id, percent_type)
+      raise ValueError(msg)
+
+    # Copy this Nuclide to separate it from the same Nuclide in other Materials
+    nuclide = deepcopy(nuclide)
+
+    self._nuclides[nuclide._name] = (nuclide, percent, percent_type)
+
+
+  def removeNuclide(self, nuclide):
+
+    if not isinstance(nuclide, openmc.Nuclide):
+      msg = 'Unable to remove a Nuclide {0} in Material ID={1} ' \
+            'since it is not a Nuclide'.format(self._id, nuclide)
+      raise ValueError(msg)
+
+    # If the Material contains the Nuclide, delete it
+    if self._nuclides.has_key(nuclide._name):
+      del self._nuclides[nuclide._name]
+
+
+  def addElement(self, element, percent, percent_type='ao'):
+
+    if not isinstance(element, openmc.Element):
+      msg = 'Unable to add an Element to Material ID={0} with a ' \
+            'non-Element value {1}'.format(self._id, element)
+      raise ValueError(msg)
+
+    if not is_float(percent):
+      msg = 'Unable to add an Element to Material ID={0} with a non-floating ' \
+            'point percent value {1}'.format(self._id, percent)
+      raise ValueError(msg)
+
+    if not percent_type in ['ao', 'wo']:
+      msg = 'Unable to add an Element to Material ID={0} with a percent ' \
+            'percent type {1}'.format(self._id, percent_type)
+      raise ValueError(msg)
+
+    # Copy this Element to separate it from the same Element in other Materials
+    element = deepcopy(element)
+
+    self._elements[element._name] = (element, percent, percent_type)
+
+
+  def removeElement(self, element):
+
+    # If the Material contains the Element, delete it
+    if self._elements.has_key(element._name):
+      del self._elements[element._name]
+
+
+  def addSAlphaBeta(self, name, xs):
+
+    if not is_string(name):
+      msg = 'Unable to add an S(a,b) table to Material ID={0} with a ' \
+            'non-string table name {1}'.format(self._id, name)
+      raise ValueError(msg)
+
+    if not is_string(xs):
+      msg = 'Unable to add an S(a,b) table to Material ID={0} with a ' \
+            'non-string cross-section identifier {1}'.format(self._id, xs)
+      raise ValueError(msg)
+
+    self._sab.append((name, xs))
+
+
+  def getAllNuclides(self):
+
+    nuclides = dict()
+
+    for nuclide_name, nuclide_tuple in self._nuclides.items():
+      nuclide = nuclide_tuple[0]
+      density = nuclide_tuple[1]
+      nuclides[nuclide._name] = (nuclide, density)
+
+    return nuclides
+
+
+  def _repr__(self):
+
+    string = 'Material\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}'.format('\tDensity', '=\t', self._density)
+    string += ' [{0}]\n'.format(self._density_units)
+
+    string += '{0: <16}'.format('\tS(a,b) Tables') + '\n'
+
+    for sab in self._sab:
+      string += '{0: <16}{1}[{2}{3}]\n'.format('\tS(a,b)', '=\t', sab[0], sab[1])
+
+    string += '{0: <16}'.format('\tNuclides') + '\n'
+
+    for nuclide in self._nuclides:
+      percent = self._nuclides[nuclide][1]
+      percent_type = self._nuclides[nuclide][2]
+      string += '{0: <16}'.format('\t{0}'.format(nuclide))
+      string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type)
+
+    string += '{0: <16}\n'.format('\tElements')
+
+    for element in self._elements:
+      percent = self._nuclides[element][1]
+      percent_type = self._nuclides[element][2]
+      string += '{0: >16}'.format('\t{0}'.format(element))
+      string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type)
+
+    return string
+
+
+  def getNuclideXML(self, nuclide):
+
+    xml_element = ET.Element("nuclide")
+    xml_element.set("name", nuclide[0]._name)
+
+    if nuclide[2] is 'ao':
+      xml_element.set("ao", str(nuclide[1]))
+    else:
+      xml_element.set("wo", str(nuclide[1]))
+
+    if not nuclide[0]._xs is None:
+      xml_element.set("xs", nuclide[0]._xs)
+
+    return xml_element
+
+
+  def getElementXML(self, element):
+
+    xml_element = ET.Element("element")
+    xml_element.set("name", str(element[0]._name))
+
+    if element[2] is 'ao':
+      xml_element.set("ao", str(element[1]))
+    else:
+      xml_element.set("wo", str(element[1]))
+
+    return xml_element
+
+
+  def getNuclidesXML(self, nuclides):
+
+    xml_elements = list()
+
+    for nuclide in nuclides.values():
+      xml_elements.append(self.getNuclideXML(nuclide))
+
+    return xml_elements
+
+
+  def getElementsXML(self, elements):
+
+    xml_elements = list()
+
+    for element in elements.values():
+      xml_elements.append(self.getElementXML(element))
+
+    return xml_elements
+
+
+  def getMaterialXML(self):
+
+    # Create Material XML element
+    element = ET.Element("material")
+    element.set("id", str(self._id))
+
+    # Create density XML subelement
+    subelement = ET.SubElement(element, "density")
+    if self._density_units is not 'sum':
+      subelement.set("value", str(self._density))
+    subelement.set("units", self._density_units)
+
+    # Create nuclide XML subelements
+    subelements = self.getNuclidesXML(self._nuclides)
+    for subelement in subelements:
+      element.append(subelement)
+
+    # Create element XML subelements
+    subelements = self.getElementsXML(self._elements)
+    for subelement in subelements:
+      element.append(subelement)
+
+    if len(self._sab) > 0:
+      for sab in self._sab:
+        subelement = ET.SubElement(element, "sab")
+        subelement.set("name", sab[0])
+        subelement.set("xs", sab[1])
+
+    return element
+
+
+class MaterialsFile(object):
+
+  def __init__(self):
+
+    # Initialize MaterialsFile class attributes
+    self._materials = list()
+    self._default_xs = None
+    self._materials_file = ET.Element("materials")
+
+
+  def addMaterial(self, material):
+
+    if not isinstance(material, Material):
+      msg = 'Unable to add a non-Material {0} to the ' \
+            'MaterialsFile'.format(material)
+      raise ValueError(msg)
+
+    self._materials.append(material)
+
+
+  def addMaterials(self, materials):
+
+    if not isinstance(materials, (tuple, list, MappingView)):
+      msg = 'Unable to create OpenMC materials.xml file from {0} which ' \
+            'is not a Python tuple/list'.format(materials)
+      raise ValueError(msg)
+
+    for material in materials:
+      self.addMaterial(material)
+
+
+  def removeMaterial(self, material):
+
+    if not isinstance(material, Material):
+      msg = 'Unable to remove a non-Material {0} from the ' \
+            'MaterialsFile'.format(material)
+      raise ValueError(msg)
+
+    self._materials.remove(material)
+
+
+  def setDefaultXS(self, xs):
+
+    if not is_string(xs):
+      msg = 'Unable to set default xs to a non-string value'.format(xs)
+      raise ValueError(msg)
+
+    self._default_xs = xs
+
+
+  def createMaterialSubelements(self):
+
+    subelement = ET.SubElement(self._materials_file, "default_xs")
+
+    if not self._default_xs is None:
+      subelement.text = self._default_xs
+
+    for material in self._materials:
+      xml_element = material.getMaterialXML()
+
+      if len(material._name) > 0:
+        self._materials_file.append(ET.Comment(material._name))
+
+      self._materials_file.append(xml_element)
+
+
+  def exportToXML(self):
+
+    self.createMaterialSubelements()
+
+    # Clean the indentation in the file to be user-readable
+    sort_xml_elements(self._materials_file)
+    clean_xml_indentation(self._materials_file)
+
+    # Write the XML Tree to the materials.xml file
+    tree = ET.ElementTree(self._materials_file)
+    tree.write("materials.xml", xml_declaration=True,
+               encoding='utf-8', method="xml")

src/utils/openmc/nuclide.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+
+from openmc.checkvalue import *
+
+
+class Nuclide(object):
+
+  def __init__(self, name='', xs=None):
+
+    # Initialize class attributes
+    self._name = ''
+    self._xs = None
+
+    # Set the Material class attributes
+    self.setName(name)
+
+    if not xs is None:
+      self.setXS(xs)
+
+
+  def __eq__(self, nuclide2):
+
+    # Check type
+    if not isinstance(nuclide2, Nuclide):
+      return False
+
+    # Check name
+    elif self._name != nuclide2._name:
+      return False
+
+    # Check xs
+    elif self._xs != nuclide2._xs:
+      return False
+
+    else:
+      return True
+
+
+  def __hash__(self):
+    hashable = list()
+    hashable.append(self._name)
+    hashable.append(self._xs)
+    return hash(tuple(hashable))
+
+
+  def setName(self, name):
+
+    if not is_string(name):
+      msg = 'Unable to set name for Nuclide with a non-string ' \
+            'value {0}'.format(name)
+      raise ValueError(msg)
+
+    self._name = name
+
+
+  def setXS(self, xs):
+
+    if not is_string(xs):
+      msg = 'Unable to set cross-section identifier xs for Nuclide with ' \
+            'a non-string value {0}'.format(xs)
+      raise ValueError(msg)
+
+    self._xs = xs
+
+
+  def __repr__(self):
+
+    string = 'Nuclide  -  {0}\n'.format(self._name)
+    string += '{0: <16}{1}{2}\n'.format('\tXS', '=\t', self._xs)
+    return string

src/utils/openmc/opencsg_compatible.py

@@ -0,0 +1,820 @@
+#!/usr/bin/env python
+
+import openmc
+import opencsg
+import copy
+import numpy as np
+
+
+# A dictionary of all OpenMC Materials created
+# Keys    - Material IDs
+# Values  - Materials
+OPENMC_MATERIALS = dict()
+
+# A dictionary of all OpenCSG Materials created
+# Keys    - Material IDs
+# Values  - Materials
+OPENCSG_MATERIALS = dict()
+
+# A dictionary of all OpenMC Surfaces created
+# Keys    - Surface IDs
+# Values  - Surfaces
+OPENMC_SURFACES = dict()
+
+# A dictionary of all OpenCSG Surfaces created
+# Keys    - Surface IDs
+# Values  - Surfaces
+OPENCSG_SURFACES = dict()
+
+# A dictionary of all OpenMC Cells created
+# Keys    - Cell IDs
+# Values  - Cells
+OPENMC_CELLS = dict()
+
+# A dictionary of all OpenCSG Cells created
+# Keys    - Cell IDs
+# Values  - Cells
+OPENCSG_CELLS = dict()
+
+# A dictionary of all OpenMC Universes created
+# Keys    - Universes IDs
+# Values  - Universes
+OPENMC_UNIVERSES = dict()
+
+# A dictionary of all OpenCSG Universes created
+# Keys    - Universes IDs
+# Values  - Universes
+OPENCSG_UNIVERSES = dict()
+
+# A dictionary of all OpenMC Lattices created
+# Keys    - Lattice IDs
+# Values  - Lattices
+OPENMC_LATTICES = dict()
+
+# A dictionary of all OpenCSG Lattices created
+# Keys    - Lattice IDs
+# Values  - Lattices
+OPENCSG_LATTICES = dict()
+
+
+
+def get_opencsg_material(openmc_material):
+
+  if not isinstance(openmc_material, openmc.Material):
+    msg = 'Unable to create an OpenCSG Material from {0} ' \
+          'which is not an OpenMC Material'.format(openmc_material)
+    raise ValueError(msg)
+
+  global OPENCSG_MATERIALS
+  material_id = openmc_material._id
+
+  # If this Material was already created, use it
+  if material_id in OPENCSG_MATERIALS.keys():
+    return OPENCSG_MATERIALS[material_id]
+
+  # Create an OpenCSG Material to represent this OpenMC Material
+  name = openmc_material._name
+  opencsg_material = opencsg.Material(material_id=material_id, name=name)
+
+  # Add the OpenMC Material to the global collection of all OpenMC Materials
+  OPENMC_MATERIALS[material_id] = openmc_material
+
+  # Add the OpenCSG Material to the global collection of all OpenCSG Materials
+  OPENCSG_MATERIALS[material_id] = opencsg_material
+
+  return opencsg_material
+
+
+def get_openmc_material(opencsg_material):
+
+  if not isinstance(opencsg_material, opencsg.Material):
+    msg = 'Unable to create an OpenMC Material from {0} ' \
+          'which is not an OpenCSG Material'.format(opencsg_material)
+    raise ValueError(msg)
+
+  global OPENMC_MATERIALS
+  material_id = opencsg_material._id
+
+  # If this Material was already created, use it
+  if material_id in OPENMC_MATERIALS.keys():
+    return OPENMC_MATERIALS[material_id]
+
+  # Create an OpenMC Material to represent this OpenCSG Material
+  name = opencsg_material._name
+  openmc_material = openmc.Material(material_id=material_id, name=name)
+
+  # Add the OpenMC Material to the global collection of all OpenMC Materials
+  OPENMC_MATERIALS[material_id] = openmc_material
+
+  # Add the OpenCSG Material to the global collection of all OpenCSG Materials
+  OPENCSG_MATERIALS[material_id] = opencsg_material
+
+  return openmc_material
+
+
+def is_opencsg_surface_compatible(opencsg_surface):
+
+  if not isinstance(opencsg_surface, opencsg.Surface):
+    msg = 'Unable to check if OpenCSG Surface is compatible' \
+          'since {0} is not a Surface'.format(opencsg_surface)
+    raise ValueError(msg)
+
+  if opencsg_surface._type in ['x-squareprism', 'y-squareprism', 'z-squareprism']:
+    return False
+  else:
+    return True
+
+
+def get_opencsg_surface(openmc_surface):
+
+  if not isinstance(openmc_surface, openmc.Surface):
+    msg = 'Unable to create an OpenCSG Surface from {0} ' \
+          'which is not an OpenMC Surface'.format(openmc_surface)
+    raise ValueError(msg)
+
+  global OPENCSG_SURFACES
+  surface_id = openmc_surface._id
+
+  # If this Material was already created, use it
+  if surface_id in OPENCSG_SURFACES.keys():
+    return OPENCSG_SURFACES[surface_id]
+
+  # Create an OpenCSG Surface to represent this OpenMC Surface
+  name = openmc_surface._name
+
+  # Correct for OpenMC's syntax for Surfaces dividing Cells
+  boundary = openmc_surface._bc_type
+  if boundary == 'transmission':
+    boundary = 'interface'
+
+  opencsg_surface = None
+
+  if openmc_surface._type == 'plane':
+    A = openmc_surface._coeffs['A']
+    B = openmc_surface._coeffs['B']
+    C = openmc_surface._coeffs['C']
+    D = openmc_surface._coeffs['D']
+    opencsg_surface = opencsg.Plane(surface_id, name, boundary, A, B, C, D)
+
+  elif openmc_surface._type == 'x-plane':
+    x0 = openmc_surface._coeffs['x0']
+    opencsg_surface = opencsg.XPlane(surface_id, name, boundary, x0)
+
+  elif openmc_surface._type == 'y-plane':
+    y0 = openmc_surface._coeffs['y0']
+    opencsg_surface = opencsg.YPlane(surface_id, name, boundary, y0)
+
+  elif openmc_surface._type == 'z-plane':
+    z0 = openmc_surface._coeffs['z0']
+    opencsg_surface = opencsg.ZPlane(surface_id, name, boundary, z0)
+
+  elif openmc_surface._type == 'x-cylinder':
+    y0 = openmc_surface._coeffs['y0']
+    z0 = openmc_surface._coeffs['z0']
+    R = openmc_surface._coeffs['R']
+    opencsg_surface = opencsg.XCylinder(surface_id, name, boundary, y0, z0, R)
+
+  elif openmc_surface._type == 'y-cylinder':
+    x0 = openmc_surface._coeffs['x0']
+    z0 = openmc_surface._coeffs['z0']
+    R = openmc_surface._coeffs['R']
+    opencsg_surface = opencsg.YCylinder(surface_id, name, boundary, x0, z0, R)
+
+  elif openmc_surface._type == 'z-cylinder':
+    x0 = openmc_surface._coeffs['x0']
+    y0 = openmc_surface._coeffs['y0']
+    R = openmc_surface._coeffs['R']
+    opencsg_surface = opencsg.ZCylinder(surface_id, name, boundary, x0, y0, R)
+
+  # Add the OpenMC Surface to the global collection of all OpenMC Surfaces
+  OPENMC_SURFACES[surface_id] = openmc_surface
+
+  # Add the OpenCSG Surface to the global collection of all OpenCSG Surfaces
+  OPENCSG_SURFACES[surface_id] = opencsg_surface
+
+  return opencsg_surface
+
+
+def get_openmc_surface(opencsg_surface):
+
+  if not isinstance(opencsg_surface, opencsg.Surface):
+    msg = 'Unable to create an OpenMC Surface from {0} which ' \
+          'is not an OpenCSG Surface'.format(opencsg_surface)
+    raise ValueError(msg)
+
+  global openmc_surface
+  surface_id = opencsg_surface._id
+
+  # If this Surface was already created, use it
+  if surface_id in OPENMC_SURFACES.keys():
+    return OPENMC_SURFACES[surface_id]
+
+  # Create an OpenMC Surface to represent this OpenCSG Surface
+  name = opencsg_surface._name
+
+  # Correct for OpenMC's syntax for Surfaces dividing Cells
+  boundary = opencsg_surface._boundary_type
+  if boundary == 'interface':
+    boundary = 'transmission'
+
+  if opencsg_surface._type == 'plane':
+    A = opencsg_surface._coeffs['A']
+    B = opencsg_surface._coeffs['B']
+    C = opencsg_surface._coeffs['C']
+    D = opencsg_surface._coeffs['D']
+    openmc_surface = openmc.Plane(surface_id, boundary, A, B, C, D, name)
+
+  elif opencsg_surface._type == 'x-plane':
+    x0 = opencsg_surface._coeffs['x0']
+    openmc_surface = openmc.XPlane(surface_id, boundary, x0, name)
+
+  elif opencsg_surface._type == 'y-plane':
+    y0 = opencsg_surface._coeffs['y0']
+    openmc_surface = openmc.YPlane(surface_id, boundary, y0, name)
+
+  elif opencsg_surface._type == 'z-plane':
+    z0 = opencsg_surface._coeffs['z0']
+    openmc_surface = openmc.ZPlane(surface_id, boundary, z0, name)
+
+  elif opencsg_surface._type == 'x-cylinder':
+    y0 = opencsg_surface._coeffs['y0']
+    z0 = opencsg_surface._coeffs['z0']
+    R = opencsg_surface._coeffs['R']
+    openmc_surface = openmc.XCylinder(surface_id, boundary, y0, z0, R, name)
+
+  elif opencsg_surface._type == 'y-cylinder':
+    x0 = opencsg_surface._coeffs['x0']
+    z0 = opencsg_surface._coeffs['z0']
+    R = opencsg_surface._coeffs['R']
+    openmc_surface = openmc.YCylinder(surface_id, boundary, x0, z0, R, name)
+
+  elif opencsg_surface._type == 'z-cylinder':
+    x0 = opencsg_surface._coeffs['x0']
+    y0 = opencsg_surface._coeffs['y0']
+    R = opencsg_surface._coeffs['R']
+    openmc_surface = openmc.ZCylinder(surface_id, boundary, x0, y0, R, name)
+
+  else:
+    msg = 'Unable to create an OpenMC Surface from an OpenCSG ' \
+          'Surface of type {0} since it is not a compatible ' \
+          'Surface type in OpenMC'.format(opencsg_surface._type)
+    raise ValueError(msg)
+
+
+  # Add the OpenMC Surface to the global collection of all OpenMC Surfaces
+  OPENMC_SURFACES[surface_id] = openmc_surface
+
+  # Add the OpenCSG Surface to the global collection of all OpenCSG Surfaces
+  OPENCSG_SURFACES[surface_id] = opencsg_surface
+
+  return openmc_surface
+
+
+def get_compatible_opencsg_surfaces(opencsg_surface):
+
+  if not isinstance(opencsg_surface, opencsg.Surface):
+    msg = 'Unable to create an OpenMC Surface from {0} which ' \
+          'is not an OpenCSG Surface'.format(opencsg_surface)
+    raise ValueError(msg)
+
+  global OPENMC_SURFACES
+  surface_id = opencsg_surface._id
+
+  # If this Surface was already created, use it
+  if surface_id in OPENMC_SURFACES.keys():
+    return OPENMC_SURFACES[surface_id]
+
+  # Create an OpenMC Surface to represent this OpenCSG Surface
+  name = opencsg_surface._name
+  boundary = opencsg_surface._boundary_type
+
+  if opencsg_surface._type == 'x-squareprism':
+    y0 = opencsg_surface._coeffs['y0']
+    z0 = opencsg_surface._coeffs['z0']
+    R = opencsg_surface._coeffs['R']
+
+    # Create a list of the four planes we need
+    left = opencsg.YPlane(name=name, boundary=boundary, y0=y0-R)
+    right = opencsg.YPlane(name=name, boundary=boundary, y0=y0+R)
+    bottom = opencsg.ZPlane(name=name, boundary=boundary, z0=z0-R)
+    top = opencsg.ZPlane(name=name, boundary=boundary, z0=z0+R)
+    surfaces = [left, right, bottom, top]
+
+  elif opencsg_surface._type == 'y-squareprism':
+    x0 = opencsg_surface._coeffs['x0']
+    z0 = opencsg_surface._coeffs['z0']
+    R = opencsg_surface._coeffs['R']
+
+    # Create a list of the four planes we need
+    left = opencsg.XPlane(name=name, boundary=boundary, x0=x0-R)
+    right = opencsg.XPlane(name=name, boundary=boundary, x0=x0+R)
+    bottom = opencsg.ZPlane(name=name, boundary=boundary, z0=z0-R)
+    top = opencsg.ZPlane(name=name, boundary=boundary, z0=z0+R)
+    surfaces = [left, right, bottom, top]
+
+  elif opencsg_surface._type == 'z-squareprism':
+    x0 = opencsg_surface._coeffs['x0']
+    y0 = opencsg_surface._coeffs['y0']
+    R = opencsg_surface._coeffs['R']
+
+    # Create a list of the four planes we need
+    left = opencsg.XPlane(name=name, boundary=boundary, x0=x0-R)
+    right = opencsg.XPlane(name=name, boundary=boundary, x0=x0+R)
+    bottom = opencsg.YPlane(name=name, boundary=boundary, y0=y0-R)
+    top = opencsg.YPlane(name=name, boundary=boundary, y0=y0+R)
+    surfaces = [left, right, bottom, top]
+
+  else:
+    msg = 'Unable to create a compatible OpenMC Surface an OpenCSG ' \
+          'Surface of type {0} since it already a compatible ' \
+          'Surface type in OpenMC'.format(opencsg_surface._type)
+    raise ValueError(msg)
+
+  # Add the OpenMC Surface(s) to the global collection of all OpenMC Surfaces
+  OPENMC_SURFACES[surface_id] = surfaces
+
+  # Add the OpenCSG Surface to the global collection of all OpenCSG Surfaces
+  OPENCSG_SURFACES[surface_id] = opencsg_surface
+
+  return surfaces
+
+
+def get_opencsg_cell(openmc_cell):
+
+  if not isinstance(openmc_cell, openmc.Cell):
+    msg = 'Unable to create an OpenCSG Cell from {0} which ' \
+          'is not an OpenMC Cell'.format(openmc_cell)
+    raise ValueError(msg)
+
+  global OPENCSG_CELLS
+  cell_id = openmc_cell._id
+
+  # If this Cell was already created, use it
+  if cell_id in OPENCSG_CELLS.keys():
+    return OPENCSG_CELLS[cell_id]
+
+  # Create an OpenCSG Cell to represent this OpenMC Cell
+  name = openmc_cell._name
+  opencsg_cell = opencsg.Cell(cell_id, name)
+
+  fill = openmc_cell._fill
+
+  if (openmc_cell._type == 'normal'):
+    opencsg_cell.setFill(get_opencsg_material(fill))
+  elif (openmc_cell._type == 'fill'):
+    opencsg_cell.setFill(get_opencsg_universe(fill))
+  else:
+    opencsg_cell.setFill(get_opencsg_lattice(fill))
+
+  surfaces = openmc_cell._surfaces
+
+  for surface_id in surfaces.keys():
+    surface = surfaces[surface_id][0]
+    halfspace = surfaces[surface_id][1]
+    opencsg_cell.addSurface(get_opencsg_surface(surface), halfspace)
+
+  # Add the OpenMC Cell to the global collection of all OpenMC Cells
+  OPENMC_CELLS[cell_id] = openmc_cell
+
+  # Add the OpenCSG Cell to the global collection of all OpenCSG Cells
+  OPENCSG_CELLS[cell_id] = opencsg_cell
+
+  return opencsg_cell
+
+
+def get_compatible_opencsg_cells(opencsg_cell, opencsg_surface, halfspace):
+
+  if not isinstance(opencsg_cell, opencsg.Cell):
+    msg = 'Unable to create compatible OpenMC Cell from {0} which ' \
+          'is not an OpenCSG Cell'.format(opencsg_cell)
+    raise ValueError(msg)
+
+  elif not isinstance(opencsg_surface, opencsg.Surface):
+    msg = 'Unable to create compatible OpenMC Cell since {0} is ' \
+          'not an OpenCSG Surface'.format(opencsg_surface)
+    raise ValueError(msg)
+
+  elif not halfspace in [-1, +1]:
+    msg = 'Unable to create compatible Cell since {0}' \
+          'is not a +/-1 halfspace'.format(halfspace)
+    raise ValueError(msg)
+
+  # Initialize an empty list for the new compatible cells
+  compatible_cells = list()
+
+  # SquarePrism Surfaces
+  if opencsg_surface._type in ['x-squareprism', 'y-squareprism', 'z-squareprism']:
+
+    # Get the compatible Surfaces (XPlanes and YPlanes)
+    compatible_surfaces = get_compatible_opencsg_surfaces(opencsg_surface)
+
+    opencsg_cell.removeSurface(opencsg_surface)
+
+    # If Cell is inside SquarePrism, add "inside" of Surface halfspaces
+    if halfspace == -1:
+      opencsg_cell.addSurface(compatible_surfaces[0], +1)
+      opencsg_cell.addSurface(compatible_surfaces[1], -1)
+      opencsg_cell.addSurface(compatible_surfaces[2], +1)
+      opencsg_cell.addSurface(compatible_surfaces[3], -1)
+      compatible_cells.append(opencsg_cell)
+
+    # If Cell is outside SquarePrism, add "outside" of Surface halfspaces
+    else:
+
+      # Create 8 Cell clones to represent each of the disjoint planar
+      # Surface halfspace intersections
+      num_clones = 8
+
+      for clone_id in range(num_clones):
+
+        # Create a cloned OpenCSG Cell with Surfaces compatible with OpenMC
+        clone = opencsg_cell.clone()
+        compatible_cells.append(clone)
+
+        # Top left subcell - add left XPlane, top YPlane
+        if clone_id == 0:
+          clone.addSurface(compatible_surfaces[0], -1)
+          clone.addSurface(compatible_surfaces[3], +1)
+
+        # Top center subcell - add top YPlane, left/right XPlanes
+        elif clone_id == 1:
+          clone.addSurface(compatible_surfaces[0], +1)
+          clone.addSurface(compatible_surfaces[1], -1)
+          clone.addSurface(compatible_surfaces[3], +1)
+
+        # Top right subcell - add top YPlane, right XPlane
+        elif clone_id == 2:
+          clone.addSurface(compatible_surfaces[1], +1)
+          clone.addSurface(compatible_surfaces[3], +1)
+
+        # Right center subcell - add right XPlane, top/bottom YPlanes
+        elif clone_id == 3:
+          clone.addSurface(compatible_surfaces[1], +1)
+          clone.addSurface(compatible_surfaces[3], -1)
+          clone.addSurface(compatible_surfaces[2], +1)
+
+        # Bottom right subcell - add right XPlane, bottom YPlane
+        elif clone_id == 4:
+          clone.addSurface(compatible_surfaces[1], +1)
+          clone.addSurface(compatible_surfaces[2], -1)
+
+        # Bottom center subcell - add bottom YPlane, left/right XPlanes
+        elif clone_id == 5:
+          clone.addSurface(compatible_surfaces[0], +1)
+          clone.addSurface(compatible_surfaces[1], -1)
+          clone.addSurface(compatible_surfaces[2], -1)
+
+        # Bottom left subcell - add bottom YPlane, left XPlane
+        elif clone_id == 6:
+          clone.addSurface(compatible_surfaces[0], -1)
+          clone.addSurface(compatible_surfaces[2], -1)
+
+        # Left center subcell - add left XPlane, top/bottom YPlanes
+        elif clone_id == 7:
+          clone.addSurface(compatible_surfaces[0], -1)
+          clone.addSurface(compatible_surfaces[3], -1)
+          clone.addSurface(compatible_surfaces[2], +1)
+
+  # Remove redundant Surfaces from the Cells
+  for cell in compatible_cells:
+    cell.removeRedundantSurfaces()
+
+  # Return the list of OpenMC compatible OpenCSG Cells
+  return compatible_cells
+
+
+def make_opencsg_cells_compatible(opencsg_universe):
+
+  if not isinstance(opencsg_universe, opencsg.Universe):
+    msg = 'Unable to make compatible OpenCSG Cells for {0} which ' \
+          'is not an OpenCSG Universe'.format(opencsg_universe)
+    raise ValueError(msg)
+
+  # Check all OpenCSG Cells in this Universe for compatibility with OpenMC
+  opencsg_cells = opencsg_universe._cells
+
+  for cell_id, opencsg_cell in opencsg_cells.items():
+
+    # Check each of the OpenCSG Surfaces for OpenMC compatibility
+    surfaces = opencsg_cell._surfaces
+
+    for surface_id in surfaces.keys():
+      surface = surfaces[surface_id][0]
+      halfspace = surfaces[surface_id][1]
+
+      # If this Surface is not compatible with OpenMC, create compatible
+      # OpenCSG cells with a compatible version of this OpenCSG Surface
+      if not is_opencsg_surface_compatible(surface):
+
+        # Get the one or more OpenCSG Cells that are compatible with OpenMC
+        # NOTE: This does not necessarily make the OpenCSG fully compatible.
+        #       It only removes the incompatible Surface and replaces it with
+        #       compatible OpenCSG Surface(s). The recursive call at the end
+        #       of this block is necessary in the event that there are more
+        #       incompatible Surfaces in this Cell that are not accounted for.
+        cells = get_compatible_opencsg_cells(opencsg_cell, surface, halfspace)
+
+        # Remove the non-compatible OpenCSG Cell from the Universe
+        opencsg_universe.removeCell(opencsg_cell)
+
+        # Add the compatible OpenCSG Cells to the Universe
+        opencsg_universe.addCells(cells)
+
+        # Make recursive call to look at the updated state of the
+        # OpenCSG Universe and return
+        return make_opencsg_cells_compatible(opencsg_universe)
+
+  # If all OpenCSG Cells in the OpenCSG Universe are compatible, return
+  return
+
+
+
+def get_openmc_cell(opencsg_cell):
+
+  if not isinstance(opencsg_cell, opencsg.Cell):
+    msg = 'Unable to create an OpenMC Cell from {0} which ' \
+          'is not an OpenCSG Cell'.format(opencsg_cell)
+    raise ValueError(msg)
+
+  global OPENMC_CELLS
+  cell_id = opencsg_cell._id
+
+  # If this Cell was already created, use it
+  if cell_id in OPENMC_CELLS.keys():
+    return OPENMC_CELLS[cell_id]
+
+  # Create an OpenCSG Cell to represent this OpenMC Cell
+  name = opencsg_cell._name
+  openmc_cell = openmc.Cell(cell_id, name)
+
+  fill = opencsg_cell._fill
+
+  if (opencsg_cell._type == 'universe'):
+    openmc_cell.setFill(get_openmc_universe(fill))
+  elif (opencsg_cell._type == 'lattice'):
+    openmc_cell.setFill(get_openmc_lattice(fill))
+  else:
+    openmc_cell.setFill(get_openmc_material(fill))
+
+  surfaces = opencsg_cell._surfaces
+
+  for surface_id in surfaces.keys():
+    surface = surfaces[surface_id][0]
+    halfspace = surfaces[surface_id][1]
+    openmc_cell.addSurface(get_openmc_surface(surface), halfspace)
+
+  # Add the OpenMC Cell to the global collection of all OpenMC Cells
+  OPENMC_CELLS[cell_id] = openmc_cell
+
+  # Add the OpenCSG Cell to the global collection of all OpenCSG Cells
+  OPENCSG_CELLS[cell_id] = opencsg_cell
+
+  return openmc_cell
+
+
+
+def get_opencsg_universe(openmc_universe):
+
+  if not isinstance(openmc_universe, openmc.Universe):
+    msg = 'Unable to create an OpenCSG Universe from {0} which ' \
+          'is not an OpenMC Universe'.format(openmc_universe)
+    raise ValueError(msg)
+
+  global OPENCSG_UNIVERSES
+  universe_id = openmc_universe._id
+
+  # If this Universe was already created, use it
+  if universe_id in OPENCSG_UNIVERSES.keys():
+    return OPENCSG_UNIVERSES[universe_id]
+
+  # Create an OpenCSG Universe to represent this OpenMC Universe
+  name = openmc_universe._name
+  opencsg_universe = opencsg.Universe(universe_id, name)
+
+  # Convert all OpenMC Cells in this Universe to OpenCSG Cells
+  openmc_cells = openmc_universe._cells
+
+  for cell_id, openmc_cell in openmc_cells.items():
+    opencsg_cell = get_opencsg_cell(openmc_cell)
+    opencsg_universe.addCell(opencsg_cell)
+
+  # Add the OpenMC Universe to the global collection of all OpenMC Universes
+  OPENMC_UNIVERSES[universe_id] = openmc_universe
+
+  # Add the OpenCSG Universe to the global collection of all OpenCSG Universes
+  OPENCSG_UNIVERSES[universe_id] = opencsg_universe
+
+  return opencsg_universe
+
+
+def get_openmc_universe(opencsg_universe):
+
+  if not isinstance(opencsg_universe, opencsg.Universe):
+    msg = 'Unable to create an OpenMC Universe from {0} which ' \
+          'is not an OpenCSG Universe'.format(opencsg_universe)
+    raise ValueError(msg)
+
+  global OPENMC_UNIVERSES
+  universe_id = opencsg_universe._id
+
+  # If this Universe was already created, use it
+  if universe_id in OPENMC_UNIVERSES.keys():
+    return OPENMC_UNIVERSES[universe_id]
+
+  # Make all OpenCSG Cells and Surfaces in this Universe compatible with OpenMC
+  make_opencsg_cells_compatible(opencsg_universe)
+
+  # Create an OpenMC Universe to represent this OpenCSg Universe
+  name = opencsg_universe._name
+  openmc_universe = openmc.Universe(universe_id, name)
+
+  # Convert all OpenCSG Cells in this Universe to OpenMC Cells
+  opencsg_cells = opencsg_universe._cells
+
+  for cell_id, opencsg_cell in opencsg_cells.items():
+    openmc_cell = get_openmc_cell(opencsg_cell)
+    openmc_universe.addCell(openmc_cell)
+
+  # Add the OpenMC Universe to the global collection of all OpenMC Universes
+  OPENMC_UNIVERSES[universe_id] = openmc_universe
+
+  # Add the OpenCSG Universe to the global collection of all OpenCSG Universes
+  OPENCSG_UNIVERSES[universe_id] = opencsg_universe
+
+  return openmc_universe
+
+
+def get_opencsg_lattice(openmc_lattice):
+
+  if not isinstance(openmc_lattice, openmc.Lattice):
+    msg = 'Unable to create an OpenCSG Lattice from {0} which ' \
+          'is not an OpenMC Lattice'.format(openmc_lattice)
+    raise ValueError(msg)
+
+  global OPENCSG_LATTICES
+  lattice_id = openmc_lattice._id
+
+  # If this Lattice was already created, use it
+  if lattice_id in OPENCSG_LATTICES.keys():
+    return OPENCSG_LATTICES[lattice_id]
+
+  # Create an OpenCSG Lattice to represent this OpenMC Lattice
+  name = openmc_lattice._name
+  dimension = openmc_lattice._dimension
+  width = openmc_lattice._width
+  lower_left = openmc_lattice._lower_left
+  universes = openmc_lattice._universes
+
+  # Initialize an empty array for the OpenCSG nested Universes in this Lattice
+  universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), \
+                              dtype=opencsg.Universe)
+
+  # Create OpenCSG Universes for each unique nested Universe in this Lattice
+  unique_universes = openmc_lattice.getUniqueUniverses()
+
+  for universe_id, universe in unique_universes.items():
+    unique_universes[universe_id] = get_opencsg_universe(universe)
+
+  # Build the nested Universe array
+  for z in range(dimension[2]):
+    for y in range(dimension[1]):
+      for x in range(dimension[0]):
+        universe_id = universes[x][y][z]._id
+        universe_array[z][y][x] = unique_universes[universe_id]
+
+  opencsg_lattice = opencsg.Lattice(lattice_id, name)
+  opencsg_lattice.setDimension(dimension)
+  opencsg_lattice.setWidth(width)
+  opencsg_lattice.setUniverses(universe_array)
+
+  offset = np.array(lower_left, dtype=np.float64) - \
+           ((np.array(width, dtype=np.float64) * \
+             np.array(dimension, dtype=np.float64))) / -2.0
+  opencsg_lattice.setOffset(offset)
+
+  # Add the OpenMC Lattice to the global collection of all OpenMC Lattices
+  OPENMC_LATTICES[lattice_id] = openmc_lattice
+
+  # Add the OpenCSG Lattice to the global collection of all OpenCSG Lattices
+  OPENCSG_LATTICES[lattice_id] = opencsg_lattice
+
+  return opencsg_lattice
+
+
+def get_openmc_lattice(opencsg_lattice):
+
+  if not isinstance(opencsg_lattice, opencsg.Lattice):
+    msg = 'Unable to create an OpenMC Lattice from {0} which ' \
+          'is not an OpenCSG Lattice'.format(opencsg_lattice)
+    raise ValueError(msg)
+
+  global OPENMC_LATTICES
+  lattice_id = opencsg_lattice._id
+
+  # If this Lattice was already created, use it
+  if lattice_id in OPENMC_LATTICES.keys():
+    return OPENMC_LATTICES[lattice_id]
+
+  dimension = opencsg_lattice._dimension
+  width = opencsg_lattice._width
+  offset = opencsg_lattice._offset
+  universes = opencsg_lattice._universes
+
+  # Initialize an empty array for the OpenMC nested Universes in this Lattice
+  universe_array = np.ndarray(tuple(np.array(dimension)), \
+                              dtype=openmc.Universe)
+
+  # Create OpenMC Universes for each unique nested Universe in this Lattice
+  unique_universes = opencsg_lattice.getUniqueUniverses()
+
+  for universe_id, universe in unique_universes.items():
+    unique_universes[universe_id] = get_openmc_universe(universe)
+
+  # Build the nested Universe array
+  for z in range(dimension[2]):
+    for y in range(dimension[1]):
+      for x in range(dimension[0]):
+        universe_id = universes[z][y][x]._id
+        universe_array[x][y][z] = unique_universes[universe_id]
+
+  openmc_lattice = openmc.Lattice(lattice_id=lattice_id)
+  openmc_lattice.setDimension(dimension)
+  openmc_lattice.setWidth(width)
+  openmc_lattice.setUniverses(universe_array)
+
+  lower_left = np.array(offset, dtype=np.float64) + \
+               ((np.array(width, dtype=np.float64) * \
+                 np.array(dimension, dtype=np.float64))) / -2.0
+  openmc_lattice.setLowerLeft(lower_left)
+
+  # Add the OpenMC Lattice to the global collection of all OpenMC Lattices
+  OPENMC_LATTICES[lattice_id] = openmc_lattice
+
+  # Add the OpenCSG Lattice to the global collection of all OpenCSG Lattices
+  OPENCSG_LATTICES[lattice_id] = opencsg_lattice
+
+  return openmc_lattice
+
+
+def get_opencsg_geometry(openmc_geometry):
+
+  if not isinstance(openmc_geometry, openmc.Geometry):
+    msg = 'Unable to get OpenCSG geometry from {0} which is ' \
+          'not an OpenMC Geometry object'.format(openmc_geometry)
+    raise ValueError(msg)
+
+  # Clear dictionaries and auto-generated IDs
+  OPENMC_SURFACES.clear()
+  OPENCSG_SURFACES.clear()
+  OPENMC_CELLS.clear()
+  OPENCSG_CELLS.clear()
+  OPENMC_UNIVERSES.clear()
+  OPENCSG_UNIVERSES.clear()
+  OPENMC_LATTICES.clear()
+  OPENCSG_LATTICES.clear()
+
+  opencsg.geometry.reset_auto_ids()
+
+  openmc_root_universe = openmc_geometry._root_universe
+  opencsg_root_universe = get_opencsg_universe(openmc_root_universe)
+
+  opencsg_geometry = opencsg.Geometry()
+  opencsg_geometry.setRootUniverse(opencsg_root_universe)
+  opencsg_geometry.initializeCellOffsets()
+
+  return opencsg_geometry
+
+
+def get_openmc_geometry(opencsg_geometry):
+
+  if not isinstance(opencsg_geometry, opencsg.Geometry):
+    msg = 'Unable to get OpenMC geometry from {0} which is ' \
+          'not an OpenCSG Geometry object'.format(opencsg_geometry)
+    raise ValueError(msg)
+
+  # Deep copy the goemetry since it may be modified to make all Surfaces
+  # compatible with OpenMC's specifications
+  opencsg_geometry = copy.deepcopy(opencsg_geometry)
+
+  # Update Cell bounding boxes in Geometry
+  opencsg_geometry.updateBoundingBoxes()
+
+  # Clear dictionaries and auto-generated ID
+  OPENMC_SURFACES.clear()
+  OPENCSG_SURFACES.clear()
+  OPENMC_CELLS.clear()
+  OPENCSG_CELLS.clear()
+  OPENMC_UNIVERSES.clear()
+  OPENCSG_UNIVERSES.clear()
+  OPENMC_LATTICES.clear()
+  OPENCSG_LATTICES.clear()
+
+  openmc.reset_auto_ids()
+
+  opencsg_root_universe = opencsg_geometry._root_universe
+  openmc_root_universe = get_openmc_universe(opencsg_root_universe)
+
+  openmc_geometry = openmc.Geometry()
+  openmc_geometry.setRootUniverse(openmc_root_universe)
+
+  return openmc_geometry

src/utils/openmc/plots.py

@@ -0,0 +1,434 @@
+#!/usr/bin/env python
+
+from openmc.checkvalue import *
+from openmc.clean_xml import *
+from xml.etree import ElementTree as ET
+import numpy as np
+
+
+# A static variable for auto-generated Plot IDs
+AUTO_PLOT_ID = 10000
+
+def reset_auto_plot_id():
+  global AUTO_PLOT_ID
+  AUTO_PLOT_ID = 10000
+
+
+BASES = ['xy', 'xz', 'yz']
+
+
+class Plot(object):
+
+  def __init__(self, plot_id=None, name=''):
+
+    # Initialize Plot class attributes
+    self._id = None
+    self._name = ''
+    self._width = [4.0, 4.0]
+    self._pixels = [1000, 1000]
+    self._origin = [0., 0., 0.]
+    self._filename = 'plot'
+    self._color = 'cell'
+    self._type = 'slice'
+    self._basis = 'xy'
+    self._background = [0, 0, 0]
+    self._mask_components = None
+    self._mask_background = None
+    self._col_spec = None
+
+    self.setId(plot_id)
+    self.setName(name)
+
+
+  def setId(self, plot_id=None):
+
+    if plot_id is None:
+      global AUTO_PLOT_ID
+      self._id = AUTO_PLOT_ID
+      AUTO_PLOT_ID += 1
+
+    # Check that the ID is an integer and wasn't already used
+    elif not is_integer(plot_id):
+      msg = 'Unable to set a non-integer Plot ID {0}'.format(plot_id)
+      raise ValueError(msg)
+
+    elif plot_id < 0:
+      msg = 'Unable to set Plot ID to {0} since it must be a ' \
+            'non-negative integer'.format(plot_id)
+      raise ValueError(msg)
+
+    else:
+      self._id = plot_id
+
+
+  def setName(self, name):
+
+    if not is_string(name):
+      msg = 'Unable to set name for Plot ID={0} with a non-string ' \
+            'value {1}'.format(self._id, name)
+      raise ValueError(msg)
+
+    else:
+      self._name = name
+
+
+  def setWidth(self, width):
+
+    if not isinstance(width, (tuple, list, np.ndarray)):
+      msg = 'Unable to create Plot ID={0} with width {1} which is not ' \
+            'a Python tuple/list or NumPy array'.format(self._id, width)
+      raise ValueError(msg)
+
+    elif len(width) != 2 and len(width) != 3:
+      msg = 'Unable to create Plot ID={0} with width {1} since only 2D ' \
+            'and 3D plots are supported'.format(self._id, width)
+      raise ValueError(msg)
+
+    for dim in width:
+      if not is_integer(dim) and not is_float(dim):
+        msg = 'Unable to create Plot ID={0} with width {1} since each ' \
+              'element must be a floating point value or ' \
+              'integer'.format(self._id, width)
+        raise ValueError(msg)
+
+    self._width = width
+
+
+  def setOrigin(self, origin):
+
+    if not isinstance(origin, (tuple, list, np.ndarray)):
+      msg = 'Unable to create Plot ID={0} with origin {1} which is not a ' \
+            'Python tuple/list or NumPy array'.format(self._id, origin)
+      raise ValueError(msg)
+
+    elif len(origin) != 3:
+      msg = 'Unable to create Plot ID={0} with origin {1} since only a 3D ' \
+            'coordinate must be input'.format(self._id, origin)
+      raise ValueError(msg)
+
+
+    for dim in origin:
+      if not is_integer(dim) and not is_float(dim):
+        msg = 'Unable to create Plot ID={0} with origin {1} since each ' \
+              'element must be a floating point value or ' \
+              'integer'.format(self._id, origin)
+        raise ValueError(msg)
+
+    self._origin = origin
+
+
+  def setPixels(self, pixels):
+
+    if not isinstance(pixels, (tuple, list, np.ndarray)):
+      msg = 'Unable to create Plot ID={0} with pixels {1} which is not a ' \
+            'Python tuple/list or NumPy array'.format(self._id, pixels)
+      raise ValueError(msg)
+
+    elif len(pixels) != 2 and len(pixels) != 3:
+      msg = 'Unable to create Plot ID={0} with pixels {1} since only 2D ' \
+            'and 3D plots are supported'.format(self._id, pixels)
+      raise ValueError(msg)
+
+    for dim in pixels:
+
+      if not is_integer(dim):
+        msg  = 'Unable to create Plot ID={0} with pixel value {1} which is ' \
+               'not an integer'.format(self._id, dim)
+        raise ValueError(msg)
+
+      elif dim < 0:
+        msg = 'Unable to create Plot ID={0} with pixel value {1} which ' \
+              'is less than 0'.format(self._id, dim)
+        raise ValueError(msg)
+
+    self._pixels = pixels
+
+
+  def setFilename(self, filename):
+
+    if not is_string(filename):
+      msg = 'Unable to create Plot ID={0} with filename {1} which is ' \
+            'not a string'.format(self._id, filename)
+      raise ValueError(msg)
+
+    self._filename = filename
+
+
+  def setColor(self, color):
+
+    if not is_string(color):
+      msg = 'Unable to create Plot ID={0} with color {1} which is not ' \
+            'a string'.format(self._id, color)
+      raise ValueError(msg)
+
+    elif not color in ['cell', 'mat']:
+      msg = 'Unable to create Plot ID={0} with color {1} which is not ' \
+            'a cell or mat'.format(self._id, color)
+      raise ValueError(msg)
+
+    self._color = color
+
+
+  def setType(self, type):
+
+    if not is_string(type):
+      msg = 'Unable to create Plot ID={0} with type {1} which is not ' \
+            'a string'.format(self._id, type)
+      raise ValueError(msg)
+
+    elif not type in ['slice', 'voxel']:
+      msg = 'Unable to create Plot ID={0} with type {1} which is not ' \
+            'slice or voxel'.format(self._id, type)
+      raise ValueError(msg)
+
+    self._type = type
+
+
+  def setBasis(self, basis):
+
+    if not is_string(basis):
+      msg = 'Unable to create Plot ID={0} with basis {1} which is not ' \
+            'a string'.format(self._id, basis)
+      raise ValueError(msg)
+
+    elif not basis in ['xy', 'xz', 'yz']:
+      msg = 'Unable to create Plot ID={0} with basis {1} which is not ' \
+            'xy, xz, or yz'.format(self._id, basis)
+      raise ValueError(msg)
+
+    self._basis = basis
+
+
+  def setBackground(self, background):
+
+    if not isinstance(background, (tuple, list, np.ndarray)):
+      msg = 'Unable to create Plot ID={0} with background {1} which is not ' \
+            'a Python tuple/list or NumPy array'.format(self._id, background)
+      raise ValueError(msg)
+
+    elif len(background) != 3:
+      msg = 'Unable to create Plot ID={0} with background {1} which is ' \
+            'not 3 integer RGB values'.format(self._id, background)
+      raise ValueError(msg)
+
+    for rgb in background:
+
+      if not is_integer(rgb):
+        msg = 'Unable to create Plot ID={0} with background RGB value {1} ' \
+              'which is not an integer'.format(self._id, rgb)
+        raise ValueError(msg)
+
+      elif rgb < 0 or rgb > 255:
+        msg = 'Unable to create Plot ID={0} with background RGB value {1} ' \
+              'which is not between 0 and 255'.format(self._id, rgb)
+        raise ValueError(msg)
+
+    self._background = background
+
+
+  def setColSpec(self, col_spec):
+
+    if not isinstance(col_spec, dict):
+      msg= 'Unable to create Plot ID={0} with col_spec parameter {1} ' \
+           'which is not a Python dictionary of IDs to ' \
+           'pixels'.format(self._id, col_spec)
+      raise ValueError(msg)
+
+    for key in col_spec.keys():
+
+      if not is_integer(key):
+        msg = 'Unable to create Plot ID={0} with col_spec ID {1} ' \
+              'which is not an integer'.format(self._id, key)
+        raise ValueError(msg)
+
+      elif key < 0:
+        msg = 'Unable to create Plot ID={0} with col_spec ID {1} which ' \
+              'is less than 0'.format(self._id, key)
+        raise ValueError(msg)
+
+      elif not isinstance(col_spec[key], (tuple, list, np.ndarray)):
+        msg = 'Unable to create Plot ID={0} with col_spec RGB values {1} ' \
+              'which is not a Python tuple/list or NumPy ' \
+              'array'.format(self._id, col_spec[key])
+        raise ValueError(msg)
+
+      elif len(col_spec[key]) != 3:
+        msg = 'Unable to create Plot ID={0} with col_spec RGB values of ' \
+              'length {1} since 3 values must be ' \
+              'input'.format(self._id, len(col_spec[key]))
+        raise ValueError(msg)
+
+    self._col_spec = col_spec
+
+
+  def setMaskComponents(self, mask_components):
+
+    if not isinstance(mask_components, (list, tuple, np.ndarray)):
+      msg = 'Unable to create Plot ID={0} with mask components {1} ' \
+            'which is not a Python tuple/list or NumPy ' \
+            'array'.format(self._id, mask_components)
+      raise ValueError(msg)
+
+    for component in mask_components:
+      if not is_integer(component):
+        msg = 'Unable to create Plot ID={0} with mask component {1} ' \
+              'which is not an integer'.format(self._id, component)
+        raise ValueError(msg)
+
+      elif component < 0:
+        msg = 'Unable to create Plot ID={0} with mask component {1} ' \
+              'which is less than 0'.format(self._id, component)
+        raise ValueError(msg)
+
+    self._mask_components = mask_components
+
+
+  def setMaskBackground(self, mask_background):
+
+    if not isinstance(mask_background, (list, tuple, np.ndarray)):
+      msg = 'Unable to create Plot ID={0} with mask background {1} ' \
+            'which is not a Python tuple/list or NumPy ' \
+            'array'.format(self._id, mask_background)
+      raise ValueError(msg)
+
+    elif len(mask_background) != 3 and len(mask_background) != 0:
+      msg = 'Unable to create Plot ID={0} with mask background {1} since ' \
+            '3 RGB values must be input'.format(self._id, mask_background)
+      raise ValueError(msg)
+
+    for rgb in mask_background:
+
+      if not is_integer(rgb):
+        msg = 'Unable to create Plot ID={0} with mask background RGB ' \
+              'value {1} which is not an integer'.format(self._id, rgb)
+        raise ValueError(msg)
+
+      elif rgb < 0 or rgb > 255:
+        msg = 'Unable to create Plot ID={0} with mask bacground RGB ' \
+              'value {1} which is not between 0 and 255'.format(self._id, rgb)
+        raise ValueError(msg)
+
+    self._mask_background = mask_background
+
+
+  def __repr__(self):
+
+    string = 'Plot\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('\tFilename', '=\t', self._filename)
+    string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type)
+    string += '{0: <16}{1}{2}\n'.format('\tBasis', '=\t', self._basis)
+    string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self._width)
+    string += '{0: <16}{1}{2}\n'.format('\tOrigin', '=\t', self._origin)
+    string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._origin)
+    string += '{0: <16}{1}{2}\n'.format('\tColor', '=\t', self._color)
+    string += '{0: <16}{1}{2}\n'.format('\tMask', '=\t', self._mask_components)
+    string += '{0: <16}{1}{2}\n'.format('\tMask',  '=\t', self._mask_background)
+    string += '{0: <16}{1}{2}\n'.format('\tCol Spec', '=\t', self._col_spec)
+    return string
+
+
+  def getPlotXML(self):
+
+    element = ET.Element("plot")
+    element.set("id", str(self._id))
+    element.set("filename", self._filename)
+    element.set("color", self._color)
+    element.set("type", self._type)
+
+    if self._type is 'slice':
+      element.set("basis", self._basis)
+
+    subelement = ET.SubElement(element, "origin")
+    text = ''
+    for coord in self._origin:
+      text += str(coord) + ' '
+    subelement.text = text.rstrip(' ')
+
+    subelement = ET.SubElement(element, "width")
+    text = ''
+    for dim in self._width:
+      text += str(dim) + ' '
+    subelement.text = text.rstrip(' ')
+
+    subelement = ET.SubElement(element, "pixels")
+    text = ''
+    for dim in self._pixels:
+      text += str(dim) + ' '
+    subelement.text = text.rstrip(' ')
+
+    if not self._mask_background is None:
+      subelement = ET.SubElement(element, "background")
+      text = ''
+      for rgb in self._background:
+        text += str(rgb) + ' '
+      subelement.text = text.rstrip(' ')
+
+    if not self._col_spec is None:
+
+      for key in self._col_spec.keys():
+        subelement = ET.SubElement(element, "col_spec")
+        subelement.set("id", '{0}'.format(key))
+        value = self._col_spec[key]
+        subelement.set("rgb", '{0} {1} {2}'.format(value[0], value[1], value[2]))
+
+    if not self._mask_components is None:
+      subelement = ET.SubElement(element, "mask")
+
+      text = ''
+      for id in self._mask_components:
+        text += str(id) + ' '
+      subelement.set("components", text.rstrip(' '))
+
+      rgb = self._mask_background
+      subelement.set("background", '{0} {1} {2}'.format(rgb[0], rgb[1], rgb[2]))
+
+    return element
+
+
+class PlotsFile(object):
+
+  def __init__(self):
+
+    # Initialize PlotsFile class attributes
+    self._plots = list()
+    self._plots_file = ET.Element("plots")
+
+
+  def addPlot(self, plot):
+
+    if not isinstance(plot, Plot):
+      msg = 'Unable to add a non-Plot {0} to the PlotsFile'.format(plot)
+      raise ValueError(msg)
+
+    self._plots.append(plot)
+
+
+  def removePlot(self, plot):
+    self._plots.remove(plot)
+
+
+  def createPlotSubelements(self):
+
+    for plot in self._plots:
+
+      xml_element = plot.getPlotXML()
+
+      if len(plot._name) > 0:
+        self._plots_file.append(ET.Comment(plot._name))
+
+      self._plots_file.append(xml_element)
+
+
+  def exportToXML(self):
+
+    self.createPlotSubelements()
+
+    # Clean the indentation in the file to be user-readable
+    clean_xml_indentation(self._plots_file)
+
+    # Write the XML Tree to the plots.xml file
+    tree = ET.ElementTree(self._plots_file)
+    tree.write("plots.xml", xml_declaration=True,
+               encoding='utf-8', method="xml")

src/utils/openmc/surface.py

@@ -0,0 +1,593 @@
+#!/usr/bin/env python
+
+from openmc.checkvalue import *
+from openmc.constants import BC_TYPES
+from xml.etree import ElementTree as ET
+
+
+# A static variable for auto-generated Surface IDs
+AUTO_SURFACE_ID = 10000
+
+def reset_auto_surface_id():
+  global AUTO_SURFACE_ID
+  AUTO_SURFACE_ID = 10000
+
+
+
+class Surface(object):
+
+  def __init__(self, surface_id=None, bc_type='transmission', name=''):
+
+    # Initialize class attributes
+    self._id = None
+    self._name = ''
+    self._type = ''
+    self._bc_type = ''
+
+    # A dictionary of the quadratic surface coefficients
+    # Key    - coefficeint name
+    # Value  - coefficient value
+    self._coeffs = dict()
+
+    # An ordered list of the coefficient names to export to XML in the
+    # proper order
+    self._coeff_keys = list()
+
+    self.setId(surface_id)
+    self.setBoundaryType(bc_type)
+    self.setName(name)
+
+
+  def setId(self, surface_id=None):
+
+    if surface_id is None:
+      global AUTO_SURFACE_ID
+      self._id = AUTO_SURFACE_ID
+      AUTO_SURFACE_ID += 1
+
+    # Check that the ID is an integer and wasn't already used
+    elif not is_integer(surface_id):
+      msg = 'Unable to set a non-integer Surface ID {0}'.format(surface_id)
+      raise ValueError(msg)
+
+    elif surface_id < 0:
+      msg = 'Unable to set Surface ID to {0} since it must be a ' \
+            'non-negative integer'.format(surface_id)
+      raise ValueError(msg)
+
+    else:
+      self._id = surface_id
+
+
+  def setName(self, name):
+
+    if not is_string(name):
+      msg = 'Unable to set name for Surface ID={0} with a non-string ' \
+            'value {1}'.format(self._id, name)
+      raise ValueError(msg)
+
+    else:
+      self._name = name
+
+
+  def setBoundaryType(self, bc_type):
+
+    if not is_string(bc_type):
+      msg = 'Unable to set boundary type for Surface ID={0} with a ' \
+            'non-string value {1}'.format(self._id, bc_type)
+      raise ValueError(msg)
+
+    elif not bc_type in BC_TYPES.values():
+      msg = 'Unable to set boundary type for Surface ID={0} to {1} which ' \
+            'is not trasmission, vacuum or reflective'.format(bc_type)
+      raise ValueError(msg)
+
+    else:
+      self._bc_type = bc_type
+
+
+  def __repr__(self):
+
+    string = 'Surface\n'
+    string += '{0: <16}{1}{2}\n'.format('\tID', '=\t{0}', self._id)
+    string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
+    string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._type)
+    string += '{0: <16}{1}{2}\n'.format('\tBoundary', '=\t', self._bc_type)
+
+    coeffs = '{0: <16}'.format('\tCoefficients') + '\n'
+
+    for coeff in self._coeffs.keys():
+      coeffs += '{0: <16}{1}{2}\n'.format(coeff, '=\t', self._coeffs[coeff])
+
+    string += coeffs
+
+    return string
+
+
+  def createXMLSubElement(self):
+
+    element = ET.Element("surface")
+    element.set("id", str(self._id))
+    element.set("type", self._type)
+    element.set("boundary", self._bc_type)
+
+    coeffs = ''
+
+    for coeff in self._coeff_keys:
+      coeffs += '{0} '.format(self._coeffs[coeff])
+
+    element.set("coeffs", coeffs.rstrip(' '))
+
+    return element
+
+
+
+class Plane(Surface):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               A=None, B=None, C=None, D=None, name='',):
+
+    # Initialize Plane class attributes
+    super(Plane, self).__init__(surface_id, bc_type, name=name)
+
+    self._A = None
+    self._B = None
+    self._C = None
+    self._D = None
+    self._type = 'plane'
+    self._coeff_keys = ['A', 'B', 'C', 'D']
+
+    if not A is None:
+      self.setA(A)
+
+    if not B is None:
+      self.setB(B)
+
+    if not C is None:
+      self.setC(C)
+
+    if not D is None:
+      self.setD(D)
+
+
+  def setA(self, A):
+
+    if not is_integer(A) and not is_float(A):
+      msg = 'Unable to set A coefficient for Plane ID={0} to a non-integer ' \
+            'value {1}'.format(self._id, A)
+      raise ValueError(msg)
+
+    self._coeffs['A'] = A
+
+
+  def setB(self, B):
+
+    if not is_integer(B) and not is_float(B):
+      msg = 'Unable to set B coefficient for Plane ID={0} to a non-integer ' \
+            'value {1}'.format(self._id, B)
+      raise ValueError(msg)
+
+    self._coeffs['B'] = B
+
+
+  def setC(self, C):
+
+    if not is_integer(C) and not is_float(C):
+      msg = 'Unable to set C coefficient for Plane ID={0} to a non-integer ' \
+            'value {1}'.format(self._id, C)
+      raise ValueError(msg)
+
+    self._coeffs['C'] = C
+
+
+  def setD(self, D):
+
+    if not is_integer(D) and not is_float(D):
+      msg = 'Unable to set D coefficient for Plane ID={0} to a non-integer ' \
+            'value {1}'.format(self._id, D)
+      raise ValueError(msg)
+
+    self._coeffs['D'] = D
+
+
+
+class XPlane(Plane):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               x0=None, name=''):
+
+    # Initialize XPlane class attributes
+    super(XPlane, self).__init__(surface_id, bc_type, name=name)
+
+    self._x0 = None
+    self._type = 'x-plane'
+    self._coeff_keys = ['x0']
+
+    if not x0 is None:
+      self.setX0(x0)
+
+
+  def setX0(self, x0):
+
+    if not is_integer(x0) and not is_float(x0):
+      msg = 'Unable to set x0 coefficient for XPlane ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, x0)
+      raise ValueError(msg)
+
+    self._coeffs['x0'] = x0
+
+
+
+class YPlane(Plane):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               y0=None, name=''):
+
+    # Initialize YPlane class attributes
+    super(YPlane, self).__init__(surface_id, bc_type, name=name)
+
+    self._y0 = None
+    self._type = 'y-plane'
+    self._coeff_keys = ['y0']
+
+    if not y0 is None:
+      self.setY0(y0)
+
+
+  def setY0(self, y0):
+
+    if not is_integer(y0) and not is_float(y0):
+      msg = 'Unable to set y0 coefficient for XPlane ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, y0)
+      raise ValueError(msg)
+
+    self._coeffs['y0'] = y0
+
+
+
+class ZPlane(Plane):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               z0=None, name=''):
+
+    # Initialize ZPlane class attributes
+    super(ZPlane, self).__init__(surface_id, bc_type, name=name)
+
+    self._z0 = None
+    self._type = 'z-plane'
+    self._coeff_keys = ['z0']
+
+    if not z0 is None:
+      self.setZ0(z0)
+
+
+  def setZ0(self, z0):
+
+    if not is_integer(z0) and not is_float(z0):
+      msg = 'Unable to set z0 coefficient for ZPlane ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, z0)
+      raise ValueError(msg)
+
+    self._coeffs['z0'] = z0
+
+
+
+class Cylinder(Surface):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               R=None, name=''):
+
+    # Initialize Cylinder class attributes
+    super(Cylinder, self).__init__(surface_id, bc_type, name=name)
+
+    self._R = None
+    self._coeff_keys = ['R']
+
+    if not R is None:
+      self.setR(R)
+
+
+  def setR(self, R):
+
+    if not is_integer(R) and not is_float(R):
+      msg = 'Unable to set R coefficient for Cylinder ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, R)
+      raise ValueError(msg)
+
+    self._coeffs['R'] = R
+
+
+
+class XCylinder(Cylinder):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               y0=None, z0=None, R=None, name=''):
+
+    # Initialize XCylinder class attributes
+    super(XCylinder, self).__init__(surface_id, bc_type, R, name=name)
+
+    self._type = 'x-cylinder'
+    self._y0 = None
+    self._z0 = None
+    self._coeff_keys = ['y0', 'z0', 'R']
+
+    if not y0 is None:
+      self.setY0(y0)
+
+    if not z0 is None:
+      self.setZ0(z0)
+
+
+  def setY0(self, y0):
+
+    if not is_integer(y0) and not is_float(y0):
+      msg = 'Unable to set y0 coefficient for XCylinder ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, y0)
+      raise ValueError(msg)
+
+    self._coeffs['y0'] = y0
+
+
+  def setZ0(self, z0):
+
+    if not is_integer(z0) and not is_float(z0):
+      msg = 'Unable to set z0 coefficient for XCylinder ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, z0)
+      raise ValueError(msg)
+
+    self._coeffs['z0'] = z0
+
+
+
+class YCylinder(Cylinder):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               x0=None, z0=None, R=None, name=''):
+
+    # Initialize YCylinder class attributes
+    super(YCylinder, self).__init__(surface_id, bc_type, R, name=name)
+
+    self._type = 'y-cylinder'
+    self._x0 = None
+    self._z0 = None
+    self._coeff_keys = ['x0', 'z0', 'R']
+
+    if not x0 is None:
+      self.setX0(x0)
+
+    if not z0 is None:
+      self.setZ0(z0)
+
+
+  def setX0(self, x0):
+
+    if not is_integer(x0) and not is_float(x0):
+      msg = 'Unable to set x0 coefficient for YCylinder ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, x0)
+      raise ValueError(msg)
+
+    self._coeffs['x0'] = x0
+
+
+  def setZ0(self, z0):
+
+    if not is_integer(z0) and not is_float(z0):
+      msg = 'Unable to set z0 coefficient for YCylinder ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, z0)
+      raise ValueError(msg)
+
+    self._coeffs['z0'] = z0
+
+
+class ZCylinder(Cylinder):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               x0=None, y0=None, R=None, name=''):
+
+    # Initialize ZCylinder class attributes
+    # Initialize YPlane class attributes
+    super(ZCylinder, self).__init__(surface_id, bc_type, R, name=name)
+
+    self._type = 'z-cylinder'
+    self._x0 = None
+    self._y0 = None
+    self._coeff_keys = ['x0', 'y0', 'R']
+
+    if not x0 is None:
+      self.setX0(x0)
+
+    if not y0 is None:
+      self.setY0(y0)
+
+
+  def setX0(self, x0):
+
+    if not is_integer(x0) and not is_float(x0):
+      msg = 'Unable to set x0 coefficient for ZCylinder ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, x0)
+      raise ValueError(msg)
+
+    self._coeffs['x0'] = x0
+
+
+  def setY0(self, y0):
+
+    if not is_integer(y0) and not is_float(y0):
+      msg = 'Unable to set y0 coefficient for ZCylinder ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, y0)
+      raise ValueError(msg)
+
+    self._coeffs['y0'] = y0
+
+
+
+class Sphere(Surface):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               x0=None, y0=None, z0=None, R=None, name=''):
+
+    # Initialize Sphere class attributes
+    super(Sphere, self).__init__(surface_id, bc_type, name=name)
+
+    self._type = 'sphere'
+    self._x0 = None
+    self._y0 = None
+    self._z0 = None
+    self._R = None
+    self._coeff_keys = ['x0', 'y0', 'z0', 'R']
+
+    if not x0 is None:
+      self.setX0(x0)
+
+    if not y0 is None:
+      self.setY0(y0)
+
+    if not z0 is None:
+      self.setZ0(z0)
+
+    if not R is None:
+      self.setZ0(R)
+
+
+  def setX0(self, x0):
+
+    if not is_integer(x0) and not is_float(x0):
+      msg = 'Unable to set x0 coefficient for Sphere ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, x0)
+      raise ValueError(msg)
+
+    self._coeffs['x0'] = x0
+
+
+  def setY0(self, y0):
+
+    if not is_integer(y0) and not is_float(y0):
+      msg = 'Unable to set y0 coefficient for Sphere ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, y0)
+      raise ValueError(msg)
+
+    self._coeffs['y0'] = y0
+
+
+  def setZ0(self, z0):
+
+    if not is_integer(z0) and not is_float(z0):
+      msg = 'Unable to set z0 coefficient for Sphere ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, z0)
+      raise ValueError(msg)
+
+    self._coeffs['z0'] = z0
+
+
+  def setR(self, R):
+
+    if not is_integer(R) and not is_float(R):
+      msg = 'Unable to set R coefficient for Sphere ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, R)
+      raise ValueError(msg)
+
+    self._coeffs['R'] = R
+
+
+
+class Cone(Surface):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               x0=None, y0=None, z0=None, R2=None, name=''):
+
+    # Initialize Cone class attributes
+    super(Cone, self).__init__(surface_id, bc_type, name=name)
+
+    self._x0 = None
+    self._y0 = None
+    self._z0 = None
+    self._R2 = None
+    self._coeff_keys = ['x0', 'y0', 'z0', 'R2']
+
+    if not x0 is None:
+      self.setX0(x0)
+
+    if not y0 is None:
+      self.setY0(y0)
+
+    if not z0 is None:
+      self.setZ0(z0)
+
+    if not R2 is None:
+      self.setZ0(R2)
+
+
+  def setX0(self, x0):
+
+    if not is_integer(x0) and not is_float(x0):
+      msg = 'Unable to set x0 coefficient for Cone ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, x0)
+      raise ValueError(msg)
+
+    self._coeffs['x0'] = x0
+
+
+  def setY0(self, y0):
+
+    if not is_integer(y0) and not is_float(y0):
+      msg = 'Unable to set y0 coefficient for Cone ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, y0)
+      raise ValueError(msg)
+
+    self._coeffs['y0'] = y0
+
+
+  def setZ0(self, z0):
+
+    if not is_integer(z0) and not is_float(z0):
+      msg = 'Unable to set z0 coefficient for Cone ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, z0)
+      raise ValueError(msg)
+
+    self._coeffs['z0'] = z0
+
+
+  def setR2(self, R2):
+
+    if not is_integer(R2) and not is_float(R2):
+      msg = 'Unable to set R^2 coefficient for Cone ID={0} to a ' \
+            'non-integer value {1}'.format(self._id, R2)
+      raise ValueError(msg)
+
+    self._coeffs['R2'] = R2
+
+
+
+class XCone(Cone):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               x0=None, y0=None, z0=None, R2=None, name=''):
+
+    # Initialize XCone class attributes
+    super(XCone, self).__init__(surface_id, bc_type, x0, y0, z0, R2, name=name)
+
+    self._type = 'x-cone'
+
+
+
+class YCone(Cone):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               x0=None, y0=None, z0=None, R2=None, name=''):
+
+    # Initialize YCone class attributes
+    super(YCone, self).__init__(surface_id, bc_type, x0, y0, z0, R2, name=name)
+
+    self._type = 'y-cone'
+
+
+
+class ZCone(Cone):
+
+  def __init__(self, surface_id=None, bc_type='transmission',
+               x0=None, y0=None, z0=None, R2=None, name=''):
+
+    # Initialize ZCone class attributes
+    super(ZCone, self).__init__(surface_id, bc_type, x0, y0, z0, R2, name=name)
+
+    self._type = 'z-cone'

src/utils/openmc/universe.py

@@ -0,0 +1,932 @@
+#!/usr/bin/env python
+
+import openmc
+from openmc.checkvalue import *
+from xml.etree import ElementTree as ET
+from collections import OrderedDict
+import operator
+import numpy as np
+
+
+################################################################################
+####################################  Cell  ####################################
+################################################################################
+
+# A static variable for auto-generated Cell IDs
+AUTO_CELL_ID = 10000
+
+def reset_auto_cell_id():
+  global AUTO_CELL_ID
+  AUTO_CELL_ID = 10000
+
+
+class Cell(object):
+
+  def __init__(self, cell_id=None, name=''):
+
+    # Initialize Cell class attributes
+    self._id = None
+    self._name = None
+    self._fill = None
+    self._type = None
+    self._surfaces = dict()
+    self._rotation = None
+    self._translation = None
+    self._offset = None
+
+    self.setId(cell_id)
+    self.setName(name)
+
+
+  def getOffset(self, path, filter_offset):
+
+    # Get the current element and remove it from the list
+    cell_id = path[0]
+    path = path[1:]
+
+    # If the Cell is filled by a Material
+    if self._type == 'normal':
+      offset = 0
+
+    # If the Cell is filled by a Universe
+    elif self._type == 'fill':
+      offset = self._offset[filter_offset-1]
+      offset += self._fill.getOffset(path, filter_offset)
+
+    # If the Cell is filled by a Lattice
+    else:
+      offset = self._fill.getOffset(path, filter_offset)
+
+    return offset
+
+    # Make a recursive call to the Universe filling this Cell
+    offset = self._cells[cell_id].getOffset(path, filter_offset)
+
+    # Return the offset computed at all nested Universe levels below this one
+    return offset
+
+
+  def setId(self, cell_id=None):
+
+    if cell_id is None:
+      global AUTO_CELL_ID
+      self._id = AUTO_CELL_ID
+      AUTO_CELL_ID += 1
+
+    # Check that the ID is an integer and wasn't already used
+    elif not is_integer(cell_id):
+       msg = 'Unable to set a non-integer Cell ID {0}'.format(cell_id)
+       raise ValueError(msg)
+
+    elif cell_id < 0:
+      msg = 'Unable to set Cell ID to {0} since it must be a ' \
+            'non-negative integer'.format(cell_id)
+      raise ValueError(msg)
+
+    else:
+      self._id = cell_id
+
+
+  def setName(self, name):
+
+    if not isinstance(name, str):
+      msg = 'Unable to set name for Cell ID={0} with a non-string ' \
+            'value {1}'.format(self._id, name)
+      raise ValueError(msg)
+
+    else:
+      self._name = name
+
+
+  def setFill(self, fill):
+
+    if not isinstance(fill, (openmc.Material, Universe, Lattice)) \
+      and fill != 'void':
+      msg = 'Unable to set Cell ID={0} to use a a non-Material or ' \
+            'Universe fill {1}'.format(self._id, fill)
+      raise ValueError(msg)
+
+    self._fill = fill
+
+    if isinstance(fill, Lattice):
+      self._type = 'lattice'
+    elif isinstance(fill, Universe):
+      self._type = 'fill'
+    elif fill == 'void':
+      self._type = 'normal'
+    else:
+      self._type = 'normal'
+
+
+  def addSurface(self, surface, halfspace):
+
+    if not isinstance(surface, openmc.Surface):
+      msg = 'Unable to add Surface {0} to Cell ID={1} since it is ' \
+            'not a Surface object'.format(surface, self._id)
+      raise ValueError(msg)
+
+    if not halfspace in [-1, +1]:
+      msg = 'Unable to add Surface {0} to Cell ID={1} with halfspace {2} ' \
+            'since it is not +/-1'.format(surface, self._id, halfspace)
+      raise ValueError(msg)
+
+    # If the Cell does not already contain the Surface, add it
+    if not surface._id in self._surfaces.keys():
+      self._surfaces[surface._id] = (surface, halfspace)
+
+
+  def removeSurface(self, surface):
+
+    if not isinstance(surface, openmc.Surface):
+      msg = 'Unable to remove Surface {0} from Cell ID={1} since it is ' \
+            'not a Surface object'.format(surface, self._id)
+      raise ValueError(msg)
+
+    # If the Cell contains the Surface, delete it
+    if surface._id in self._surfaces:
+      del self._surfaces[surface._id]
+
+
+  def setSurfaces(self, surfaces, halfspaces):
+
+    if not isinstance(surfaces, (tuple, list, np.ndarray)):
+      msg = 'Unable to set Cell ID={0} with Surfaces {1} since it is not a ' \
+            'Python tuple/list or NumPy array'.format(self._id, surfaces)
+      raise ValueError(msg)
+
+    if not isinstance(halfspaces, (tuple, list, np.ndarray)):
+      msg = 'Unable to set Cell ID={0} with Surface halfspaces {1} ' \
+            'since it is not a Python tuple/list or NumPy ' \
+            'array'.format(self._id, halfspaces)
+      raise ValueError(msg)
+
+    for surface in surfaces:
+      self.addSurface(surface)
+
+
+  def setRotation(self, rotation):
+
+    if not isinstance(rotation, (tuple, list, np.ndarray)):
+      msg = 'Unable to add rotation {0} to Cell ID={1} since it is not a ' \
+            'Python tuple/list or NumPy array'.format(rotation, self._id)
+      raise ValueError(msg)
+
+    elif len(rotation) != 3:
+      msg = 'Unable to add rotation {0} to Cell ID={1} since it does not ' \
+            'contain 3 values'.format(rotation, self._id)
+      raise ValueError(msg)
+
+    for axis in rotation:
+
+      if not is_integer(axis) and not is_float(axis):
+        msg = 'Unable to add rotation {0} to Cell ID={1} since it is not ' \
+              'an integer or floating point value'.format(axis, self._id)
+        raise ValueError(msg)
+
+      self._rotation = rotation
+
+
+  def setTranslation(self, translation):
+
+    if not isinstance(translation, (tuple, list, np.ndarray)):
+      msg = 'Unable to add translation {0} to Cell ID={1} since it is not a ' \
+            'Python tuple/list or NumPy array'.format(translation, self._id)
+      raise ValueError(msg)
+
+    elif len(translation) != 3:
+      msg = 'Unable to add translation {0} to Cell ID={1} since it does not ' \
+            'contain 3 values'.format(translation, self._id)
+      raise ValueError(msg)
+
+    for axis in translation:
+      if not is_integer(axis) and not is_float(axis):
+        msg = 'Unable to add translation {0} to Cell ID={1} since it is ' \
+              'not an integer or floating point value'.format(axis, self._id)
+        raise ValueError(msg)
+
+    self._translation = translation
+
+  def setOffset(self, offset):
+
+    if not isinstance(offset, (tuple, list, np.ndarray)):
+      msg = 'Unable to set offset {0} to Cell ID={1} since it is not a ' \
+            'Python tuple/list or NumPy array'.format(offset, self._id)
+      raise ValueError(msg)
+
+    self._offset = offset
+
+
+  def getAllNuclides(self):
+
+    nuclides = dict()
+
+    if self._type != 'void':
+      nuclides.update(self._fill.getAllNuclides())
+
+    return nuclides
+
+
+  def getAllCells(self):
+
+    cells = dict()
+
+    if self._type == 'fill' or self._type == 'lattice':
+      cells.update(self._fill.getAllCells())
+
+    return cells
+
+
+  def getAllUniverses(self):
+
+    universes = dict()
+
+    if self._type == 'fill':
+      universes[self._fill._id] = self._fill
+      universes.update(self._fill.getAllUniverses())
+    elif self._type == 'lattice':
+      universes.update(self._fill.getAllUniverses())
+
+    return universes
+
+
+  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)
+
+    if isinstance(self._fill, openmc.Material):
+      string += '{0: <16}{1}{2}\n'.format('\tMaterial', '=\t', self._fill._id)
+    elif isinstance(self._fill, (Universe, Lattice)):
+      string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t', self._fill._id)
+    else:
+      string += '{0: <16}{1}{2}\n'.format('\tFill', '=\t', self._fill)
+
+    string += '{0: <16}{1}\n'.format('\tSurfaces', '=\t')
+
+    for surface_id in self._surfaces:
+      halfspace = self._surfaces[surface_id][1]
+      string += '{0} '.format(halfspace * surface_id)
+
+    string = string.rstrip(' ') + '\n'
+
+    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._offset)
+
+    return string
+
+
+  def createXMLSubElement(self, xml_element):
+
+    element = ET.Element("cell")
+    element.set("id", str(self._id))
+
+    if isinstance(self._fill, openmc.Material):
+      element.set("material", str(self._fill._id))
+
+    elif isinstance(self._fill, (Universe, Lattice)):
+      element.set("fill", str(self._fill._id))
+      self._fill.createXMLSubElement(xml_element)
+
+    else:
+      element.set("fill", str(self._fill))
+      self._fill.createXMLSubElement(xml_element)
+
+
+    if not self._surfaces is None:
+
+      surfaces = ''
+
+      for surface_id in self._surfaces.keys():
+
+        # Determine if XML element already includes subelement for this Surface
+        path = './surface[@id=\'{0}\']'.format(surface_id)
+        test = xml_element.find(path)
+
+        # If the element does not contain the Surface subelement, then add it
+        if test is None:
+
+          # Create the XML subelement for this Surface and add it to the element
+          surface = self._surfaces[surface_id][0]
+          surface_subelement = surface.createXMLSubElement()
+
+          if len(surface._name) > 0:
+            xml_element.append(ET.Comment(surface._name))
+
+          xml_element.append(surface_subelement)
+
+        # Append to the Surfaces's XML attribute the halfspace and Surface ID
+        halfspace = self._surfaces[surface_id][1]
+        surfaces += '{0} '.format(halfspace * surface_id)
+
+      element.set("surfaces", surfaces.rstrip(' '))
+
+
+    if not self._translation is None:
+
+      translation  = ''
+
+      for axis in self._translation:
+        translation += '{0} '.format(axis)
+
+      element.set("translation", translation.rstrip(' '))
+
+
+    if not self._rotation is None:
+
+      rotation  = ''
+
+      for axis in self._rotation:
+        rotation += '{0} '.format(axis)
+
+      element.set("rotation", rotation.rstrip(' '))
+
+    return element
+
+
+
+################################################################################
+###################################  Universe  #################################
+################################################################################
+
+# A static variable for auto-generated Lattice (Universe) IDs
+AUTO_UNIVERSE_ID = 10000
+
+def reset_auto_universe_id():
+  global AUTO_UNIVERSE_ID
+  AUTO_UNIVERSE_ID = 10000
+
+
+class Universe(object):
+
+  def __init__(self, universe_id=None, name=''):
+
+    # Initialize Cell class attributes
+    self._id = None
+    self._name = None
+
+    # Keys   - Cell IDs
+    # Values - Cells
+    self._cells = dict()
+
+    # Keys   - Cell IDs
+    # Values - Offsets
+    self._cell_offsets = OrderedDict()
+    self._num_regions = 0
+
+    self.setId(universe_id)
+    self.setName(name)
+
+
+  def setId(self, universe_id=None):
+
+    if universe_id is None:
+      global AUTO_UNIVERSE_ID
+      self._id = AUTO_UNIVERSE_ID
+      AUTO_UNIVERSE_ID += 1
+
+    # Check that the ID is an integer and wasn't already used
+    elif not is_integer(universe_id):
+      msg = 'Unable to set Universe ID to a non-integer {0}'.format(universe_id)
+      raise ValueError(msg)
+
+    elif universe_id < 0:
+      msg = 'Unable to set Universe ID to {0} since it must be a ' \
+            'non-negative integer'.format(universe_id)
+      raise ValueError(msg)
+
+    else:
+      self._id = universe_id
+
+
+  def setName(self, name):
+
+    if not is_string(name):
+      msg = 'Unable to set name for Universe ID={0} with a non-string ' \
+            'value {1}'.format(self._id, name)
+      raise ValueError(msg)
+
+    else:
+      self._name = name
+
+
+  def addCell(self, cell):
+
+    if not isinstance(cell, Cell):
+      msg = 'Unable to add a Cell to Universe ID={0} since {1} is not ' \
+            'a Cell'.format(self._id, cell)
+      raise ValueError(msg)
+
+    cell_id = cell._id
+
+    if not cell_id in self._cells.keys():
+      self._cells[cell_id] = cell
+
+
+  def addCells(self, cells):
+
+    if not isinstance(cells, (list, tuple, np.ndarray)):
+      msg = 'Unable to add Cells to Universe ID={0} since {1} is not a ' \
+            'Python tuple/list or NumPy array'.format(self._id, cells)
+      raise ValueError(msg)
+
+    for i in range(len(cells)):
+      self.addCell(cells[i])
+
+
+  def removeCell(self, cell):
+
+    if not isinstance(cell, Cell):
+      msg = 'Unable to remove a Cell from Universe ID={0} since {1} is ' \
+            'not a Cell'.format(self._id, cell)
+      raise ValueError(msg)
+
+    cell_id = cell.getId()
+
+    # If the Cell is in the Universe's list of Cells, delete it
+    if cell_id in self._cells.keys():
+      del self._cells[cell_id]
+
+
+  def clearCells(self):
+    self._cells.clear()
+
+
+  def getOffset(self, path, filter_offset):
+
+    # Get the current element and remove it from the list
+    path = path[1:]
+
+    # Get the Cell ID
+    cell_id = path[0]
+
+    # Make a recursive call to the Cell within this Universe
+    offset = self._cells[cell_id].getOffset(path, filter_offset)
+
+    # Return the offset computed at all nested Universe levels below this one
+    return offset
+
+
+  def getAllNuclides(self):
+
+    nuclides = dict()
+
+    # Append all Nuclides in each Cell in the Universe to the dictionary
+    for cell_id, cell in self._cells.items():
+      nuclides.update(cell.getAllNuclides())
+
+    return nuclides
+
+
+  def getAllCells(self):
+
+    cells = dict()
+
+    # Add this Universe's cells to the dictionary
+    cells.update(self._cells)
+
+    # Append all Cells in each Cell in the Universe to the dictionary
+    for cell_id, cell in self._cells.items():
+      cells.update(cell.getAllCells())
+
+    return cells
+
+
+  def getAllUniverses(self):
+
+    # Get all Cells in this Universe
+    cells = self.getAllCells()
+
+    universes = dict()
+
+    # Append all Universes containing each Cell to the dictionary
+    for cell_id, cell in cells.items():
+      universes.update(cell.getAllUniverses())
+
+    return universes
+
+
+  def __repr__(self):
+
+    string = 'Universe\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('\tCells', '=\t', self._cells.keys())
+    string += '{0: <16}{1}{2}\n'.format('\t# Regions', '=\t', self._num_regions)
+    return string
+
+
+  def createXMLSubElement(self, xml_element):
+
+    # Iterate over all Cells
+    for cell_id, cell in self._cells.items():
+
+      # Determine if XML element already contains subelement for this Cell
+      path = './cell[@id=\'{0}\']'.format(cell_id)
+      test = xml_element.find(path)
+
+      # If the element does not contain the Cell subelement, then add it
+      if test is None:
+
+        # Create the XML subelement for this Cell and everything beneath it
+        cell_subelement = cell.createXMLSubElement(xml_element)
+
+        # Append the Universe ID to the subelement and add it to the element
+        cell_subelement.set("universe", str(self._id))
+
+        if len(cell._name) > 0:
+          xml_element.append(ET.Comment(cell._name))
+
+        xml_element.append(cell_subelement)
+
+
+
+################################################################################
+###################################  Lattice  ##################################
+################################################################################
+
+
+class Lattice(object):
+
+  def __init__(self, lattice_id=None, name='', type='rectangular'):
+
+    # Initialize Lattice class attributes
+    self._id = None
+    self._name = None
+    self._type = ''
+    self._dimension = None
+    self._lower_left = None
+    self._width = None
+    self._outside = None
+    self._universes = None
+    self._offsets = None
+
+
+    self.setId(lattice_id)
+    self.setName(name)
+    self.setType(type)
+
+
+  def setId(self, lattice_id=None):
+
+    if lattice_id is None:
+      global AUTO_UNIVERSE_ID
+      self._id = AUTO_UNIVERSE_ID
+      AUTO_UNIVERSE_ID += 1
+
+    # Check that the ID is an integer and wasn't already used
+    elif not is_integer(lattice_id):
+      msg = 'Unable to set a non-integer Lattice ID {0}'.format(lattice_id)
+      raise ValueError(msg)
+
+    elif lattice_id < 0:
+      msg = 'Unable to set Lattice ID to {0} since it must be a ' \
+           'non-negative integer'.format(lattice_id)
+      raise ValueError(msg)
+
+    else:
+      self._id = lattice_id
+
+
+  def setName(self, name):
+
+    if not is_string(name):
+      msg = 'Unable to set name for Lattice ID={0} with a non-string ' \
+            'value {1}'.format(self._id, name)
+      raise ValueError(msg)
+
+    else:
+      self._name = name
+
+
+  def setType(self, type):
+
+    if not is_string(type):
+      msg = 'Unable to set type for Lattice ID={0} with a non-string ' \
+           'value {1}'.format(self._id, type)
+      raise ValueError(msg)
+
+    elif not type in ['rectangular', 'hexagonal']:
+      msg = 'Unable to set type for Lattice ID={0} as {1} since it is not ' \
+           'rectangular or hexagonal'.format(self._id, type)
+      raise ValueError(msg)
+
+    self._type = type
+
+
+  def setDimension(self, dimension):
+
+    if not isinstance(dimension, (tuple, list, np.ndarray)):
+      msg = 'Unable to set Lattice ID={0} dimension to {1} since it is not ' \
+            'a Python tuple/list or NumPy array'.format(self._id, dimension)
+      raise ValueError(msg)
+
+    elif len(dimension) != 2 and len(dimension) != 3:
+      msg = 'Unable to set Lattice ID={0} dimension to {1} since it does ' \
+            'not contain 2 or 3 coordinates'.format(self._id, dimension)
+      raise ValueError(msg)
+
+    for dim in dimension:
+
+      if not is_integer(dim) and not is_float(dim):
+        msg = 'Unable to set Lattice ID={0} dimension to {1} since it is ' \
+             'not an integer or floating point value'.format(self._id, dim)
+        raise ValueError(msg)
+
+
+      elif dim < 0:
+        msg = 'Unable to set Lattice ID={0} dimension to {1} since it ' \
+              'is a negative value'.format(self._id, dim)
+        raise ValueError(msg)
+
+    self._dimension = dimension
+
+
+  def setLowerLeft(self, lower_left):
+
+    if not isinstance(lower_left, (tuple, list, np.ndarray)):
+      msg = 'Unable to set Lattice ID={0} lower_left to {1} since it is not ' \
+            'a Python tuple/list or NumPy array'.format(self._id, lower_left)
+      raise ValueError(msg)
+
+    elif len(lower_left) != 2 and len(lower_left) != 3:
+      msg = 'Unable to set Lattice ID={0} lower_left to {1} since it does ' \
+            'not contain 2 or 3 coordinates'.format(self._id, lower_left)
+      raise ValueError(msg)
+
+
+    for dim in lower_left:
+
+      if not is_integer(dim) and not is_float(dim):
+        msg = 'Unable to set Lattice ID={0} lower_left to {1} since it is ' \
+              'is not an integer or floating point value'.format(self._id, dim)
+        raise ValueError(msg)
+
+    self._lower_left = lower_left
+
+
+  def setWidth(self, width):
+
+    if not isinstance(width, (tuple, list, np.ndarray)):
+      msg = 'Unable to set Lattice ID={0} width to {1} since it is not a ' \
+             'Python tuple/list or NumPy array'.format(self._id, width)
+      raise ValueError(msg)
+
+
+    elif len(width) != 2 and len(width) != 3:
+      msg = 'Unable to set Lattice ID={0} width to {1} since it does ' \
+            'not contain 2 or 3 coordinates'.format(self._id, width)
+      raise ValueError(msg)
+
+    for dim in width:
+
+      if not is_integer(dim) and not is_float(dim):
+        msg = 'Unable to set Lattice ID={0} width to {1} since it is not an ' \
+              'an integer or floating point value'.format(self._id, dim)
+        raise ValueError(msg)
+
+      elif dim < 0:
+        msg = 'Unable to set Lattice ID={0} width to {1} since it ' \
+              'is a negative value'.format(self._id, dim)
+        raise ValueError(msg)
+
+    self._width = width
+
+
+  def setOutside(self, outside):
+
+    if not isinstance(outside, Universe):
+      msg = 'Unable to set Lattice ID={0} outside universe to {1} ' \
+            'since it is not a Universe object'.format(self._id, outside)
+      raise ValueError(msg)
+
+    self._outside = outside
+
+
+  def setUniverses(self, universes):
+
+    if not isinstance(universes, (tuple, list, np.ndarray)):
+      msg = 'Unable to set Lattice ID={0} universes to {1} since it is not ' \
+            'a Python tuple/list or NumPy array'.format(self._id, universes)
+      raise ValueError(msg)
+
+    self._universes = np.asarray(universes, dtype=Universe)
+
+
+  def setOffsets(self, offsets):
+
+    if not isinstance(offsets, (tuple, list, np.ndarray)):
+      msg = 'Unable to set Lattice ID={0} offsets to {1} since it is not ' \
+            'a Python tuple/list or NumPy array'.format(self._id, offsets)
+      raise ValueError(msg)
+
+    self._offsets = offsets
+
+
+  def getOffset(self, path, filter_offset):
+
+    # Get the current element and remove it from the list
+    i = path[0]
+    path = path[1:]
+
+    # For 2D Lattices
+    if len(self._dimension) == 2:
+      offset = self._offsets[i[1]-1, i[2]-1, 0, filter_offset-1]
+      offset += self._universes[i[1]][i[2]].getOffset(path, filter_offset)
+
+    # For 3D Lattices
+    else:
+      offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, filter_offset-1]
+      offset += self._universes[i[1]-1][i[2]-1][i[3]-1].getOffset(path, \
+                                                                  filter_offset)
+
+    return offset
+
+
+  def getUniqueUniverses(self):
+
+    unique_universes = np.unique(self._universes.ravel())
+    universes = dict()
+
+    for universe in unique_universes:
+      universes[universe._id] = universe
+
+    return universes
+
+
+  def getAllNuclides(self):
+
+    nuclides = dict()
+
+    # Get all unique Universes contained in each of the lattice cells
+    unique_universes = self.getUniqueUniverses()
+
+    # Append all Universes containing each cell to the dictionary
+    for universe_id, universe in unique_universes.items():
+      nuclides.update(universe.getAllNuclides())
+
+    return nuclides
+
+
+  def getAllCells(self):
+
+    cells = dict()
+    unique_universes = self.getUniqueUniverses()
+
+    for universe_id, universe in unique_universes.items():
+      cells.update(universe.getAllCells())
+
+    return cells
+
+
+  def getAllUniverses(self):
+
+    # Initialize a dictionary of all Universes contained by the Lattice
+    # in each nested Universe level
+    all_universes = dict()
+
+    # Get all unique Universes contained in each of the lattice cells
+    unique_universes = self.getUniqueUniverses()
+
+    # Add the unique Universes filling each Lattice cell
+    all_universes.update(unique_universes)
+
+    # Append all Universes containing each cell to the dictionary
+    for universe_id, universe in unique_universes.items():
+      all_universes.update(universe.getAllUniverses())
+
+    return all_universes
+
+
+  def __repr__(self):
+
+    string = 'Lattice\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('\tType', '=\t', self._type)
+    string += '{0: <16}{1}{2}\n'.format('\tDimension', '=\t', self._dimension)
+    string += '{0: <16}{1}{2}\n'.format('\tLower Left', '=\t', self._lower_left)
+    string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self._width)
+
+    if self._outside is not None:
+      string += '{0: <16}{1}{2}\n'.format('\tOutside', '=\t', self._outside._id)
+    else:
+      string += '{0: <16}{1}{2}\n'.format('\tOutside', '=\t', self._outside)
+
+    string += '{0: <16}\n'.format('\tUniverses')
+
+    # Lattice nested Universe IDs - column major for Fortran
+    for i, universe in enumerate(np.ravel(self._universes)):
+      string += '{0} '.format(universe._id)
+
+      # Add a newline character every time we reach the end of a row of cells
+      if (i+1) % self._dimension[-1] == 0:
+        string += '\n'
+
+    string = string.rstrip('\n')
+
+    if self._offsets is not None:
+      string += '{0: <16}\n'.format('\tOffsets')
+
+      # Lattice cell offsets
+      for i, offset in enumerate(np.ravel(self._offsets)):
+
+        string += '{0} '.format(offset)
+
+        # Add a newline character every time we reach the end of a row of cells
+        if (i+1) % self._dimension[-1] == 0:
+          string += '\n'
+
+      string = string.rstrip('\n')
+
+    return string
+
+
+  def createXMLSubElement(self, xml_element):
+
+    # Determine if XML element already contains subelement for this Lattice
+    path = './lattice[@id=\'{0}\']'.format(self._id)
+    test = xml_element.find(path)
+
+    # If the element does contain the Lattice subelement, then return
+    if not test is None:
+      return
+
+    lattice_subelement = ET.Element("lattice")
+    lattice_subelement.set("id", str(self._id))
+    lattice_subelement.set("type", self._type)
+
+    # Export Lattice cell dimensions
+    dimension = ET.SubElement(lattice_subelement, "dimension")
+    dimension.text = '{0} {1} {2}'.format(self._dimension[0], \
+                                          self._dimension[1], \
+                                          self._dimension[2])
+
+    # Export Lattice lower left
+    lower_left = ET.SubElement(lattice_subelement, "lower_left")
+    lower_left.text = '{0} {1} {2}'.format(self._lower_left[0], \
+                                           self._lower_left[1], \
+                                           self._lower_left[2])
+
+    # Export the Lattice cell width/height
+    width = ET.SubElement(lattice_subelement, "width")
+    width.text = '{0} {1} {2}'.format(self._width[0], \
+                                      self._width[1], \
+                                      self._width[2])
+
+    # Export the Lattice outside Universe (if specified)
+    if self._outside is not None:
+      outside = ET.SubElement(lattice_subelement, "outside")
+      outside.text = '{0}'.format(self._outside._id)
+
+    # Export the Lattice nested Universe IDs - column major for Fortran
+    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]):
+
+            universe = self._universes[x][y][z]
+
+            # Append Universe ID to the string for the Lattice XML subelement
+            universe_ids += '{0} '.format(universe._id)
+
+            # Create XML subelement for this Universe and everything beneath it
+            universe.createXMLSubElement(xml_element)
+
+          # Add newline character every time we reach the end of row of cells
+          universe_ids += '\n'
+
+        # Add newline character every time we reach the end of row of cells
+        universe_ids += '\n'
+
+    # 2D Lattices
+    else:
+      for y in range(self._dimension[1]):
+        for x in range(self._dimension[0]):
+
+          universe = self._universes[x][y]
+
+          # Append Universe ID to the string for the Lattice XML subelement
+          universe_ids += '{0} '.format(universe._id)
+
+          # Create XML subelement for this Universe and everything beneath it
+          universe.createXMLSubElement(xml_element)
+
+        # Add newline character every time we reach the end of row of cells
+        universe_ids += '\n'
+
+    # Remove trailing newline character from Universe IDs string
+    universe_ids = universe_ids.rstrip('\n')
+
+    universes = ET.SubElement(lattice_subelement, "universes")
+    universes.text = universe_ids
+
+    if len(self._name) > 0:
+      xml_element.append(ET.Comment(self._name))
+
+    # Append the XML subelement for this Lattice to the XML element
+    xml_element.append(lattice_subelement)

src/utils/setup.py

@@ -9,3 +9,11 @@ setup(name='statepoint',
       author_email='paul.k.romano@gmail.com',
       url='https://github.com/mit-crpg/openmc',
       py_modules=['statepoint'])
+
+setup(name='openmc',
+      version='0.6.0',
+      description='OpenMC Python API',
+      author='Will Boyd',
+      author_email='wbinventor@gmail.com',
+      url='https://github.com/mit-crpg/openmc',
+      packages=['openmc'])