diff --git a/examples/python/lattice/hexagonal/build-xml.py b/examples/python/lattice/hexagonal/build-xml.py
new file mode 100644
index 0000000000..6e0a58bc1b
--- /dev/null
+++ b/examples/python/lattice/hexagonal/build-xml.py
@@ -0,0 +1,158 @@
+import openmc
+
+
+###############################################################################
+#                      Simulation Input File Parameters
+###############################################################################
+
+# OpenMC simulation parameters
+batches = 20
+inactive = 10
+particles = 10000
+
+
+###############################################################################
+#                 Exporting to OpenMC materials.xml File
+###############################################################################
+
+# Instantiate some Nuclides
+h1 = openmc.Nuclide('H-1')
+o16 = openmc.Nuclide('O-16')
+u235 = openmc.Nuclide('U-235')
+fe56 = openmc.Nuclide('Fe-56')
+
+# Instantiate some Materials and register the appropriate Nuclides
+fuel = openmc.Material(material_id=1, name='fuel')
+fuel.set_density('g/cc', 4.5)
+fuel.add_nuclide(u235, 1.)
+
+moderator = openmc.Material(material_id=2, name='moderator')
+moderator.set_density('g/cc', 1.0)
+moderator.add_nuclide(h1, 2.)
+moderator.add_nuclide(o16, 1.)
+moderator.add_s_alpha_beta('HH2O', '71t')
+
+iron = openmc.Material(material_id=3, name='iron')
+iron.set_density('g/cc', 7.9)
+iron.add_nuclide(fe56, 1.)
+
+# Instantiate a MaterialsFile, register all Materials, and export to XML
+materials_file = openmc.MaterialsFile()
+materials_file.set_default_xs('71c')
+materials_file.add_materials([moderator, fuel, iron])
+materials_file.export_to_xml()
+
+
+###############################################################################
+#                 Exporting to OpenMC geometry.xml File
+###############################################################################
+
+# Instantiate Surfaces
+left = openmc.XPlane(surface_id=1, x0=-3, name='left')
+right = openmc.XPlane(surface_id=2, x0=3, name='right')
+bottom = openmc.YPlane(surface_id=3, y0=-4, name='bottom')
+top = openmc.YPlane(surface_id=4, y0=4, name='top')
+fuel_surf = openmc.ZCylinder(surface_id=5, x0=0, y0=0, R=0.4)
+
+left.set_boundary_type('vacuum')
+right.set_boundary_type('vacuum')
+top.set_boundary_type('vacuum')
+bottom.set_boundary_type('vacuum')
+
+# Instantiate Cells
+cell1 = openmc.Cell(cell_id=1, name='Cell 1')
+cell2 = openmc.Cell(cell_id=101, name='cell 2')
+cell3 = openmc.Cell(cell_id=102, name='cell 3')
+cell4 = openmc.Cell(cell_id=500, name='cell 4')
+cell5 = openmc.Cell(cell_id=600, name='cell 5')
+cell6 = openmc.Cell(cell_id=601, name='cell 6')
+
+# Register Surfaces with Cells
+cell1.add_surface(left, halfspace=+1)
+cell1.add_surface(right, halfspace=-1)
+cell1.add_surface(bottom, halfspace=+1)
+cell1.add_surface(top, halfspace=-1)
+cell2.add_surface(fuel_surf, halfspace=-1)
+cell3.add_surface(fuel_surf, halfspace=+1)
+cell5.add_surface(fuel_surf, halfspace=-1)
+cell6.add_surface(fuel_surf, halfspace=+1)
+
+# Register Materials with Cells
+cell2.set_fill(fuel)
+cell3.set_fill(moderator)
+cell4.set_fill(moderator)
+cell5.set_fill(iron)
+cell6.set_fill(moderator)
+
+# Instantiate Universe
+univ1 = openmc.Universe(universe_id=1)
+univ2 = openmc.Universe(universe_id=3)
+univ3 = openmc.Universe(universe_id=4)
+root = openmc.Universe(universe_id=0, name='root universe')
+
+# Register Cells with Universe
+univ1.add_cells([cell2, cell3])
+univ2.add_cells([cell4])
+univ3.add_cells([cell5, cell6])
+root.add_cell(cell1)
+
+# Instantiate a Lattice
+lattice = openmc.HexLattice(lattice_id=5)
+lattice.set_center([0., 0., 0.])
+lattice.set_pitch([1., 2.])
+lattice.set_universes([
+     [ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ],
+     [ [univ2] + [univ1]*11, [univ2] + [univ1]*5, [univ1] ],
+     [ [univ2] + [univ3]*11, [univ2] + [univ3]*5, [univ3] ]])
+lattice.set_outer(univ2)
+
+# Fill Cell with the Lattice
+cell1.set_fill(lattice)
+
+# Instantiate a Geometry and register the root Universe
+geometry = openmc.Geometry()
+geometry.set_root_universe(root)
+
+# Instantiate a GeometryFile, register Geometry, and export to XML
+geometry_file = openmc.GeometryFile()
+geometry_file.set_geometry(geometry)
+geometry_file.export_to_xml()
+
+
+###############################################################################
+#                   Exporting to OpenMC settings.xml File
+###############################################################################
+
+# Instantiate a SettingsFile, set all runtime parameters, and export to XML
+settings_file = openmc.SettingsFile()
+settings_file.set_batches(batches)
+settings_file.set_inactive(inactive)
+settings_file.set_particles(particles)
+settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
+settings_file.export_to_xml()
+
+
+###############################################################################
+#                   Exporting to OpenMC plots.xml File
+###############################################################################
+
+plot_xy = openmc.Plot(plot_id=1)
+plot_xy.set_filename('plot_xy')
+plot_xy.set_origin([0, 0, 0])
+plot_xy.set_width([6, 6])
+plot_xy.set_pixels([400, 400])
+plot_xy.set_color('mat')
+
+plot_yz = openmc.Plot(plot_id=2)
+plot_yz.set_filename('plot_yz')
+plot_yz.set_basis('yz')
+plot_yz.set_origin([0, 0, 0])
+plot_yz.set_width([8, 8])
+plot_yz.set_pixels([400, 400])
+plot_yz.set_color('mat')
+
+# Instantiate a PlotsFile, add Plot, and export to XML
+plot_file = openmc.PlotsFile()
+plot_file.add_plot(plot_xy)
+plot_file.add_plot(plot_yz)
+plot_file.export_to_xml()
diff --git a/src/utils/openmc/universe.py b/src/utils/openmc/universe.py
index 57956a4991..e70b6effd5 100644
--- a/src/utils/openmc/universe.py
+++ b/src/utils/openmc/universe.py
@@ -610,36 +610,6 @@ class Lattice(object):
             self._name = name
 
 
-    def set_pitch(self, pitch):
-
-        if not isinstance(pitch, (tuple, list, np.ndarray)):
-            msg = 'Unable to set Lattice ID={0} pitch to {1} since ' \
-                  'it is not a Python tuple/list or NumPy ' \
-                  'array'.format(self._id, pitch)
-            raise ValueError(msg)
-
-
-        elif len(pitch) != 2 and len(pitch) != 3:
-            msg = 'Unable to set Lattice ID={0} pitch to {1} since it does ' \
-                  'not contain 2 or 3 coordinates'.format(self._id, pitch)
-            raise ValueError(msg)
-
-        for dim in pitch:
-
-            if not is_integer(dim) and not is_float(dim):
-                msg = 'Unable to set Lattice ID={0} pitch 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} pitch to {1} since it ' \
-                      'is a negative value'.format(self._id, dim)
-                raise ValueError(msg)
-
-        self._pitch = pitch
-
-
     def set_outer(self, outer):
 
         if not isinstance(outer, Universe):
@@ -716,37 +686,6 @@ class Lattice(object):
         return all_universes
 
 
-    def create_xml_subelement(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))
-
-        # Export the Lattice cell pitch
-        if len(self._pitch) == 3:
-            pitch = ET.SubElement(lattice_subelement, "pitch")
-            pitch.text = '{0} {1} {2}'.format(self._pitch[0], \
-                                              self._pitch[1], \
-                                              self._pitch[2])
-        else:
-            pitch = ET.SubElement(lattice_subelement, "pitch")
-            pitch.text = '{0} {1}'.format(self._pitch[0], \
-                                          self._pitch[1])
-
-        # Export the Lattice outer Universe (if specified)
-        if self._outer is not None:
-            outer = ET.SubElement(lattice_subelement, "outer")
-            outer.text = '{0}'.format(self._outer._id)
-
-        return lattice_subelement
-
 
 
 ################################################################################
@@ -834,6 +773,36 @@ class RectLattice(Lattice):
         self._offsets = offsets
 
 
+    def set_pitch(self, pitch):
+
+        if not isinstance(pitch, (tuple, list, np.ndarray)):
+            msg = 'Unable to set Lattice ID={0} pitch to {1} since ' \
+                  'it is not a Python tuple/list or NumPy ' \
+                  'array'.format(self._id, pitch)
+            raise ValueError(msg)
+
+
+        elif len(pitch) != 2 and len(pitch) != 3:
+            msg = 'Unable to set Lattice ID={0} pitch to {1} since it does ' \
+                  'not contain 2 or 3 coordinates'.format(self._id, pitch)
+            raise ValueError(msg)
+
+        for dim in pitch:
+
+            if not is_integer(dim) and not is_float(dim):
+                msg = 'Unable to set Lattice ID={0} pitch 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} pitch to {1} since it ' \
+                      'is a negative value'.format(self._id, dim)
+                raise ValueError(msg)
+
+        self._pitch = pitch
+
+
     def get_offset(self, path, filter_offset):
 
         # Get the current element and remove it from the list
@@ -911,8 +880,24 @@ class RectLattice(Lattice):
         if not test is None:
             return
 
-        lattice_subelement = \
-             super(RectLattice, self).create_xml_subelement(xml_element)
+        lattice_subelement = ET.Element("lattice")
+        lattice_subelement.set("id", str(self._id))
+
+        # Export the Lattice cell pitch
+        if len(self._pitch) == 3:
+            pitch = ET.SubElement(lattice_subelement, "pitch")
+            pitch.text = '{0} {1} {2}'.format(self._pitch[0], \
+                                              self._pitch[1], \
+                                              self._pitch[2])
+        else:
+            pitch = ET.SubElement(lattice_subelement, "pitch")
+            pitch.text = '{0} {1}'.format(self._pitch[0], \
+                                          self._pitch[1])
+
+        # Export the Lattice outer Universe (if specified)
+        if self._outer is not None:
+            outer = ET.SubElement(lattice_subelement, "outer")
+            outer.text = '{0}'.format(self._outer._id)
 
         # Export Lattice cell dimensions
         if len(self._dimension) == 3:
@@ -1050,6 +1035,36 @@ class HexLattice(Lattice):
         self._center = center
 
 
+    def set_pitch(self, pitch):
+
+        if not isinstance(pitch, (tuple, list, np.ndarray)):
+            msg = 'Unable to set Lattice ID={0} pitch to {1} since ' \
+                  'it is not a Python tuple/list or NumPy ' \
+                  'array'.format(self._id, pitch)
+            raise ValueError(msg)
+
+
+        elif len(pitch) != 1 and len(pitch) != 2:
+            msg = 'Unable to set Lattice ID={0} pitch to {1} since it does ' \
+                  'not contain 2 or 3 coordinates'.format(self._id, pitch)
+            raise ValueError(msg)
+
+        for dim in pitch:
+
+            if not is_integer(dim) and not is_float(dim):
+                msg = 'Unable to set Lattice ID={0} pitch 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} pitch to {1} since it ' \
+                      'is a negative value'.format(self._id, dim)
+                raise ValueError(msg)
+
+        self._pitch = pitch
+
+
     def set_universes(self, universes):
         super(HexLattice, self).set_universes(universes)
 
@@ -1058,7 +1073,78 @@ class HexLattice(Lattice):
         # The sub-lists are ordered from outermost ring to innermost ring.
         # The Universes within each sub-list are ordered from the "top" in a
         # clockwise fashion.
-        self.set_num_rings(universes.shape[0])
+
+        # Check to see if the given universes look like a 2D or a 3D array.
+        if isinstance(self._universes[0][0], Universe):
+            n_dims = 2
+
+        elif isinstance(self._universes[0][0][0], Universe):
+            n_dims = 3
+
+        else:
+            msg = 'HexLattice ID={0:d} does not appear to be either 2D or ' \
+                  '3D.  Make sure set_universes was given a two-deep or ' \
+                  'three-deep iterable of universes.'.format(self._id)
+            raise RuntimeError(msg)
+
+        # Set the number of axial positions.
+        if n_dims == 3:
+            self.set_num_axial(self._universes.shape[0])
+
+        else:
+            self._num_axial = None
+
+        # Set the number of rings and make sure this number is consistent for
+        # all axial positions.
+        if n_dims == 3:
+            self.set_num_rings(len(self._universes[0]))
+            for rings in self._universes:
+                if len(rings) != self._num_rings:
+                    msg = 'HexLattice ID={0:d} has an inconsistent number of ' \
+                          'rings per axial positon'.format(self._id)
+                    raise ValueError(msg)
+
+        else:
+            self.set_num_rings(self._universes.shape[0])
+
+        # Make sure there are the correct number of elements in each ring.
+        if n_dims == 3:
+            for axial_slice in self._universes:
+              # Check the center ring.
+              if len(axial_slice[-1]) != 1:
+                  msg = 'HexLattice ID={0:d} has the wrong number of ' \
+                        'elements in the innermost ring.  Only 1 element is ' \
+                        'allowed in the innermost ring.'.format(self._id)
+                  raise ValueError(msg)
+
+              # Check the outer rings.
+              for r in range(self._num_rings-1):
+                  if len(axial_slice[r]) != 6*(self._num_rings - 1 - r):
+                      msg = 'HexLattice ID={0:d} has the wrong number of ' \
+                            'elements in ring number {1:d} (counting from the '\
+                            'outermost ring).  This ring should have {2:d} ' \
+                            'elements.'.format(self._id, r,
+                            6*(self._num_rings - 1 - r))
+                      raise ValueError(msg)
+
+        else:
+          axial_slice = self._universes
+          # Check the center ring.
+          if len(axial_slice[-1]) != 1:
+              msg = 'HexLattice ID={0:d} has the wrong number of ' \
+                    'elements in the innermost ring.  Only 1 element is ' \
+                    'allowed in the innermost ring.'.format(self._id)
+              raise ValueError(msg)
+
+          # Check the outer rings.
+          for r in range(self._num_rings-1):
+              if len(axial_slice[r]) != 6*(self._num_rings - 1 - r):
+                  msg = 'HexLattice ID={0:d} has the wrong number of ' \
+                        'elements in ring number {1:d} (counting from the '\
+                        'outermost ring).  This ring should have {2:d} ' \
+                        'elements.'.format(self._id, r,
+                        6*(self._num_rings - 1 - r))
+                  raise ValueError(msg)
 
 
     def __repr__(self):
@@ -1081,16 +1167,12 @@ class HexLattice(Lattice):
 
         string += '{0: <16}\n'.format('\tUniverses')
 
-        # FIXME: This loop must be revised for hexagonal lattice ordering
-        # Lattice nested Universe IDs - column major for Fortran
-        for i, universe in enumerate(np.ravel(self._universes)):
-            string += '{0} '.format(universe._id)
+        if self._num_axial is not None:
+            slices = [self._repr_axial_slice(x) for x in self._universes]
+            string += '\n'.join(slices)
 
-            # Add a newline character every time we reach end of row of cells
-            if (i+1) % self._dimension[-1] == 0:
-                string += '\n'
-
-        string = string.rstrip('\n')
+        else:
+            string += self._repr_axial_slice(self._universes)
 
         return string
 
@@ -1098,18 +1180,34 @@ class HexLattice(Lattice):
     def create_xml_subelement(self, xml_element):
 
         # Determine if XML element already contains subelement for this Lattice
-        path = './lattice[@id=\'{0}\']'.format(self._id)
+        path = './hex_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 = \
-             super(HexLattice, self).create_xml_subelement(xml_element)
+        lattice_subelement = ET.Element("hex_lattice")
+        lattice_subelement.set("id", str(self._id))
+
+        # Export the Lattice cell pitch
+        if len(self._pitch) == 2:
+            pitch = ET.SubElement(lattice_subelement, "pitch")
+            pitch.text = '{0} {1}'.format(self._pitch[0], \
+                                          self._pitch[1])
+        else:
+            pitch = ET.SubElement(lattice_subelement, "pitch")
+            pitch.text = '{0}'.format(self._pitch[0])
+
+        # Export the Lattice outer Universe (if specified)
+        if self._outer is not None:
+            outer = ET.SubElement(lattice_subelement, "outer")
+            outer.text = '{0}'.format(self._outer._id)
 
         lattice_subelement.set("n_rings", str(self._num_rings))
-        lattice_subelement.set("n_axial", str(self._num_axial))
+
+        if self._num_axial is not None:
+            lattice_subelement.set("n_axial", str(self._num_axial))
 
         # Export Lattice cell center
         if len(self._center) == 3:
@@ -1122,54 +1220,156 @@ class HexLattice(Lattice):
             dimension.text = '{0} {1}'.format(self._center[0], \
                                               self._center[1])
 
-        # FIXME: The following must be revised for hexagonal lattice ordering
-        # Export the Lattice nested Universe IDs - column major for Fortran
-        universe_ids = '\n'
+        # Export the Lattice nested Universe IDs.
 
         # 3D Lattices
-        if len(self._dimension) == 3:
-            for z in range(self._dimension[2]):
-                for y in range(self._dimension[1]):
-                    for x in range(self._dimension[0]):
+        if self._num_axial is not None:
+            slices = []
+            for z in range(self._num_axial):
+                # Initialize the center universe.
+                universe = self._universes[z][-1][0]
+                universe.create_xml_subelement(xml_element)
 
-                        universe = self._universes[x][y][z]
-
-                        # Append Universe ID to the Lattice XML subelement
-                        universe_ids += '{0} '.format(universe._id)
-
-                        # Create XML subelement for this Universe
+                # Initialize the remaining universes.
+                for r in range(self._num_rings-1):
+                    for theta in range(6*(self._num_rings - 1 - r)):
+                        universe = self._universes[z][r][theta]
                         universe.create_xml_subelement(xml_element)
 
-                    # Add newline character when we reach end of row of cells
-                    universe_ids += '\n'
+                # Get a string representation of the universe IDs.
+                slices.append(self._repr_axial_slice(self._universes[z]))
 
-                # Add newline character when we reach end of row of cells
-                universe_ids += '\n'
+            # Collapse the list of axial slices into a single string.
+            universe_ids = '\n'.join(slices)
 
         # 2D Lattices
         else:
-            for y in range(self._dimension[1]):
-                for x in range(self._dimension[0]):
+            # Initialize the center universe.
+            universe = self._universes[-1][0]
+            universe.create_xml_subelement(xml_element)
 
-                    universe = self._universes[x][y]
-
-                    # Append Universe ID to Lattice XML subelement
-                    universe_ids += '{0} '.format(universe._id)
-
-                    # Create XML subelement for this Universe
+            # Initialize the remaining universes.
+            for r in range(self._num_rings-1):
+                for theta in range(2*(self._num_rings - r)):
+                    universe = self._universes[r][theta]
                     universe.create_xml_subelement(xml_element)
 
-                # Add newline character when we reach end of row of cells
-                universe_ids += '\n'
-
-        # Remove trailing newline character from Universe IDs string
-        universe_ids = universe_ids.rstrip('\n')
+            # Get a string representation of the universe IDs.
+            universe_ids = self._repr_axial_slice(self._universes)
 
         universes = ET.SubElement(lattice_subelement, "universes")
-        universes.text = universe_ids
+        universes.text = '\n' + 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)
+
+
+    def _repr_axial_slice(self, universes):
+        """Return string representation for the given 2D group of universes.
+
+        The 'universes' argument should be a list of lists of universes where
+        each sub-list represents a single ring.  The first list should be the
+        outer ring.
+        """
+
+        # Find the largest universe ID and count the number of digits so we can
+        # properly pad the output string later.
+        largest_id = max([max([univ._id for univ in ring])
+                          for ring in universes])
+        n_digits = len(str(largest_id))
+        pad = ' '*n_digits
+        id_form = '{: ^' + str(n_digits) + 'd}'
+
+        # Initialize the list for each row.
+        rows = [ [] for i in range(1 + 4 * (self._num_rings-1)) ]
+        middle = 2 * (self._num_rings - 1)
+
+        # Start with the degenerate first ring.
+        universe = universes[-1][0]
+        rows[middle] = [id_form.format(universe._id)]
+
+        # Add universes one ring at a time.
+        for r in range(1, self._num_rings):
+            # r_prime increments down while r increments up.
+            r_prime = self._num_rings - 1 - r
+            theta = 0
+            y = middle + 2*r
+
+            # Climb down the top-right.
+            for i in range(r):
+                # Add the universe.
+                universe = universes[r_prime][theta]
+                rows[y].append(id_form.format(universe._id))
+
+                # Translate the indices.
+                y -= 1
+                theta += 1
+
+            # Climb down the right.
+            for i in range(r):
+                # Add the universe.
+                universe = universes[r_prime][theta]
+                rows[y].append(id_form.format(universe._id))
+
+                # Translate the indices.
+                y -= 2
+                theta += 1
+
+            # Climb down the bottom-right.
+            for i in range(r):
+                # Add the universe.
+                universe = universes[r_prime][theta]
+                rows[y].append(id_form.format(universe._id))
+
+                # Translate the indices.
+                y -= 1
+                theta += 1
+
+            # Climb up the bottom-left.
+            for i in range(r):
+                # Add the universe.
+                universe = universes[r_prime][theta]
+                rows[y].insert(0, id_form.format(universe._id))
+
+                # Translate the indices.
+                y += 1
+                theta += 1
+
+            # Climb up the left.
+            for i in range(r):
+                # Add the universe.
+                universe = universes[r_prime][theta]
+                rows[y].insert(0, id_form.format(universe._id))
+
+                # Translate the indices.
+                y += 2
+                theta += 1
+
+            # Climb up the top-left.
+            for i in range(r):
+                # Add the universe.
+                universe = universes[r_prime][theta]
+                rows[y].insert(0, id_form.format(universe._id))
+
+                # Translate the indices.
+                y += 1
+                theta += 1
+
+        # Flip the rows and join each row into a single string.
+        rows = [pad.join(x) for x in rows[::-1]]
+
+        # Pad the beginning of the rows so they line up properly.
+        for y in range(self._num_rings - 1):
+            rows[y] = (self._num_rings - 1 - y)*pad + rows[y]
+            rows[-1 - y] = (self._num_rings - 1 - y)*pad + rows[-1 - y]
+
+        for y in range(self._num_rings % 2, self._num_rings, 2):
+            rows[middle + y] = pad + rows[middle + y]
+            if y != 0: rows[middle - y] = pad + rows[middle - y]
+
+        # Join the rows together and return the string.
+        universe_ids = '\n'.join(rows)
+        return universe_ids