diff --git a/.gitignore b/.gitignore
index 383232296..f73309646 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,6 +9,7 @@
 *.pyc
 
 # Python distribution
+.settings/
 dist/
 openmc.egg-info/
 
diff --git a/docs/source/io_formats/geometry.rst b/docs/source/io_formats/geometry.rst
index 94511842d..617ac6e87 100644
--- a/docs/source/io_formats/geometry.rst
+++ b/docs/source/io_formats/geometry.rst
@@ -318,6 +318,13 @@ the following attributes or sub-elements:
 
     *Default*: None
 
+  :orientation:
+    The orientation of the hexagonal lattice. The string "x" indicates that two
+    sides of the lattice are parallel to the x-axis, whereas the string "y"
+    indicates that two sides are parallel to the y-axis.
+
+    *Default*: "y"
+
   :center:
     The coordinates of the center of the lattice. If the lattice does not have
     axial sections then only the x- and y-coordinates are specified.
diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h
index 25ab4c034..cd73b9640 100644
--- a/include/openmc/lattice.h
+++ b/include/openmc/lattice.h
@@ -27,6 +27,7 @@ enum class LatticeType {
   rect, hex
 };
 
+
 //==============================================================================
 // Global variables
 //==============================================================================
@@ -265,8 +266,20 @@ public:
   void to_hdf5_inner(hid_t group_id) const;
 
 private:
+  enum class Orientation {
+      y, //!< Flat side of lattice parallel to y-axis
+      x  //!< Flat side of lattice parallel to x-axis
+  };
+
+  //! Fill universes_ vector for 'y' orientation
+  void fill_lattice_y(const std::vector<std::string>& univ_words);
+
+  //! Fill universes_ vector for 'x' orientation
+  void fill_lattice_x(const std::vector<std::string>& univ_words);
+
   int n_rings_;                   //!< Number of radial tile positions
   int n_axial_;                   //!< Number of axial tile positions
+  Orientation orientation_;       //!< Orientation of lattice
   Position center_;               //!< Global center of lattice
   std::array<double, 2> pitch_;   //!< Lattice tile width and height
 };
diff --git a/openmc/lattice.py b/openmc/lattice.py
index fb0a7a1bd..4776222da 100644
--- a/openmc/lattice.py
+++ b/openmc/lattice.py
@@ -36,7 +36,7 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta):
     outer : openmc.Universe
         A universe to fill all space outside the lattice
     universes : Iterable of Iterable of openmc.Universe
-        A two- or three-dimensional list/array of universes filling each element
+        A two-or three-dimensional list/array of universes filling each element
         of the lattice
 
     """
@@ -141,6 +141,10 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta):
             center = group['center'][()]
             pitch = group['pitch'][()]
             outer = group['outer'][()]
+            if 'orientation' in group:
+                orientation = group['orientation'][()].decode()
+            else:
+                orientation = "y"
 
             universe_ids = group['universes'][()]
 
@@ -148,67 +152,124 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta):
             lattice = openmc.HexLattice(lattice_id, name)
             lattice.center = center
             lattice.pitch = pitch
-
+            lattice.orientation = orientation
             # If the Universe specified outer the Lattice is not void
             if outer >= 0:
                 lattice.outer = universes[outer]
+            if orientation == "y":
+                # Build array of Universe pointers for the Lattice.  Note that
+                # we need to convert between the HDF5's square array of
+                # (x, alpha, z) to the Python API's format of a ragged nested
+                # list of (z, ring, theta).
+                uarray = []
+                for z in range(n_axial):
+                    # Add a list for this axial level.
+                    uarray.append([])
+                    x = n_rings - 1
+                    a = 2*n_rings - 2
+                    for r in range(n_rings - 1, 0, -1):
+                        # Add a list for this ring.
+                        uarray[-1].append([])
 
-            # Build array of Universe pointers for the Lattice.  Note that
-            # we need to convert between the HDF5's square array of
-            # (x, alpha, z) to the Python API's format of a ragged nested
-            # list of (z, ring, theta).
-            uarray = []
-            for z in range(n_axial):
-                # Add a list for this axial level.
-                uarray.append([])
-                x = n_rings - 1
-                a = 2*n_rings - 2
-                for r in range(n_rings - 1, 0, -1):
-                    # Add a list for this ring.
-                    uarray[-1].append([])
+                        # Climb down the top-right.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, a, x])
+                            x += 1
+                            a -= 1
 
-                    # Climb down the top-right.
-                    for i in range(r):
-                        uarray[-1][-1].append(universe_ids[z, a, x])
-                        x += 1
+                        # Climb down the right.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, a, x])
+                            a -= 1
+
+                        # Climb down the bottom-right.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, a, x])
+                            x -= 1
+
+                        # Climb up the bottom-left.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, a, x])
+                            x -= 1
+                            a += 1
+
+                        # Climb up the left.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, a, x])
+                            a += 1
+
+                        # Climb up the top-left.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, a, x])
+                            x += 1
+
+                        # Move down to the next ring.
                         a -= 1
 
-                    # Climb down the right.
-                    for i in range(r):
-                        uarray[-1][-1].append(universe_ids[z, a, x])
+                        # Convert the ids into Universe objects.
+                        uarray[-1][-1] = [universes[u_id]
+                                          for u_id in uarray[-1][-1]]
+
+                    # Handle the degenerate center ring separately.
+                    u_id = universe_ids[z, a, x]
+                    uarray[-1].append([universes[u_id]])
+            else:
+                # Build array of Universe pointers for the Lattice.  Note that
+                # we need to convert between the HDF5's square array of
+                # (alpha, y, z) to the Python API's format of a ragged nested
+                # list of (z, ring, theta).
+                uarray = []
+                for z in range(n_axial):
+                    # Add a list for this axial level.
+                    uarray.append([])
+                    a = 2*n_rings - 2
+                    y = n_rings - 1
+                    for r in range(n_rings - 1, 0, -1):
+                        # Add a list for this ring.
+                        uarray[-1].append([])
+
+                        # Climb down the bottom-right.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, y, a])
+                            y -= 1
+
+                        # Climb across the bottom.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, y, a])
+                            a -= 1
+
+                        # Climb up the bottom-left.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, y, a])
+                            a -= 1
+                            y +=1
+
+                        # Climb up the top-left.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, y, a])
+                            y += 1
+
+                        # Climb across the top.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, y, a])
+                            a += 1
+
+                        # Climb down the top-right.
+                        for i in range(r):
+                            uarray[-1][-1].append(universe_ids[z, y, a])
+                            a += 1
+                            y -= 1
+
+                        # Move down to the next ring.
                         a -= 1
 
-                    # Climb down the bottom-right.
-                    for i in range(r):
-                        uarray[-1][-1].append(universe_ids[z, a, x])
-                        x -= 1
+                        # Convert the ids into Universe objects.
+                        uarray[-1][-1] = [universes[u_id]
+                                          for u_id in uarray[-1][-1]]
 
-                    # Climb up the bottom-left.
-                    for i in range(r):
-                        uarray[-1][-1].append(universe_ids[z, a, x])
-                        x -= 1
-                        a += 1
-
-                    # Climb up the left.
-                    for i in range(r):
-                        uarray[-1][-1].append(universe_ids[z, a, x])
-                        a += 1
-
-                    # Climb up the top-left.
-                    for i in range(r):
-                        uarray[-1][-1].append(universe_ids[z, a, x])
-                        x += 1
-
-                    # Move down to the next ring.
-                    a -= 1
-
-                    # Convert the ids into Universe objects.
-                    uarray[-1][-1] = [universes[u_id]
-                                      for u_id in uarray[-1][-1]]
-
-                # Handle the degenerate center ring separately.
-                u_id = universe_ids[z, a, x]
-                uarray[-1].append([universes[u_id]])
+                    # Handle the degenerate center ring separately.
+                    u_id = universe_ids[z, y, a]
+                    uarray[-1].append([universes[u_id]])
 
             # Add the universes to the lattice.
             if len(pitch) == 2:
@@ -346,7 +407,9 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta):
             Lattice element indices. For a rectangular lattice, the indices are
             given in the :math:`(x,y)` or :math:`(x,y,z)` coordinate system. For
             hexagonal lattices, they are given in the :math:`x,\alpha` or
-            :math:`x,\alpha,z` coordinate systems.
+            :math:`x,\alpha,z` coordinate systems for "y" orientations and
+            :math:`\alpha,y` or :math:`\alpha,y,z` coordinate systems for "x"
+            orientations.
 
         Returns
         -------
@@ -646,7 +709,8 @@ class RectLattice(Lattice):
         return (x, y, z)
 
     def get_universe_index(self, idx):
-        """Return index in the universes array corresponding to a lattice element index
+        """Return index in the universes array corresponding
+        to a lattice element index
 
         Parameters
         ----------
@@ -790,7 +854,8 @@ class RectLattice(Lattice):
         lat_id = int(get_text(elem, 'id'))
         name = get_text(elem, 'name')
         lat = cls(lat_id, name)
-        lat.lower_left = [float(i) for i in get_text(elem, 'lower_left').split()]
+        lat.lower_left = [float(i)
+                          for i in get_text(elem, 'lower_left').split()]
         lat.pitch = [float(i) for i in get_text(elem, 'pitch').split()]
         outer = get_text(elem, 'outer')
         if outer is not None:
@@ -800,7 +865,7 @@ class RectLattice(Lattice):
         dimension = get_text(elem, 'dimension').split()
         shape = np.array(dimension, dtype=int)[::-1]
         uarray = np.array([get_universe(int(i)) for i in
-                            get_text(elem, 'universes').split()])
+                           get_text(elem, 'universes').split()])
         uarray.shape = shape
         lat.universes = uarray
         return lat
@@ -815,10 +880,15 @@ class HexLattice(Lattice):
 
     Most methods for this class use a natural indexing scheme wherein elements
     are assigned an index corresponding to their position relative to skewed
-    :math:`(x,\alpha,z)` axes as described fully in
-    :ref:`hexagonal_indexing`. However, note that when universes are assigned to
-    lattice elements using the :attr:`HexLattice.universes` property, the array
-    indices do not correspond to natural indices.
+    :math:`(x,\alpha,z)` or :math:`(\alpha,y,z)` bases, depending on the lattice
+    orientation, as described fully in :ref:`hexagonal_indexing`. However, note
+    that when universes are assigned to lattice elements using the
+    :attr:`HexLattice.universes` property, the array indices do not correspond
+    to natural indices.
+
+    .. versionchanged:: 0.11
+        The orientation of the lattice can now be changed with the
+        :attr:`orientation` attribute.
 
     Parameters
     ----------
@@ -855,6 +925,9 @@ class HexLattice(Lattice):
         possible, where z is the axial index, r is in the ring index (starting
         from the outermost ring), and i is the index with a ring starting from
         the top and proceeding clockwise.
+    orientation : {'x', 'y'}
+        str by default 'y' orientation of main lattice diagonal another option
+        - 'x'
     num_rings : int
         Number of radial ring positions in the xy-plane
     num_axial : int
@@ -869,11 +942,14 @@ class HexLattice(Lattice):
         self._num_rings = None
         self._num_axial = None
         self._center = None
+        self._orientation = 'y'
 
     def __repr__(self):
         string = 'HexLattice\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('\tOrientation', '=\t',
+                                            self._orientation)
         string += '{0: <16}{1}{2}\n'.format('\t# Rings', '=\t', self._num_rings)
         string += '{0: <16}{1}{2}\n'.format('\t# Axial', '=\t', self._num_axial)
         string += '{0: <16}{1}{2}\n'.format('\tCenter', '=\t',
@@ -902,6 +978,10 @@ class HexLattice(Lattice):
     def num_rings(self):
         return self._num_rings
 
+    @property
+    def orientation(self):
+        return self._orientation
+
     @property
     def num_axial(self):
         return self._num_axial
@@ -955,6 +1035,11 @@ class HexLattice(Lattice):
         cv.check_length('lattice center', center, 2, 3)
         self._center = center
 
+    @orientation.setter
+    def orientation(self, orientation):
+        cv.check_value('orientation', orientation.lower(), ('x', 'y'))
+        self._orientation = orientation.lower()
+
     @Lattice.pitch.setter
     def pitch(self, pitch):
         cv.check_type('lattice pitch', pitch, Iterable, Real)
@@ -1000,7 +1085,7 @@ class HexLattice(Lattice):
                 # 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 ' \
+                          'elements in the innermost ring. Only 1 element is ' \
                           'allowed in the innermost ring.'.format(self._id)
                     raise ValueError(msg)
 
@@ -1009,7 +1094,7 @@ class HexLattice(Lattice):
                     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} ' \
+                              'outermost ring). This ring should have {2:d} ' \
                               'elements.'.format(self._id, r,
                                                  6*(self._num_rings - 1 - r))
                         raise ValueError(msg)
@@ -1045,7 +1130,7 @@ class HexLattice(Lattice):
         -------
         3-tuple of int
             Indices of corresponding lattice element in :math:`(x,\alpha,z)`
-            bases
+            or :math:`(\alpha,y,z)` bases
         numpy.ndarray
             Carestian coordinates of the point in the corresponding lattice
             element coordinate system
@@ -1059,15 +1144,21 @@ class HexLattice(Lattice):
         else:
             z = point[2] - self.center[2]
             iz = floor(z/self.pitch[1] + 0.5*self.num_axial)
-        alpha = y - x/sqrt(3.)
-        ix = floor(x/(sqrt(0.75) * self.pitch[0]))
-        ia = floor(alpha/self.pitch[0])
-
+        if self._orientation == 'x':
+            alpha = y - x*sqrt(3.)
+            i1 = floor(-alpha/(sqrt(3.0) * self.pitch[0]))
+            i2 = floor(y/(sqrt(0.75) * self.pitch[0]))
+        else:
+            alpha = y - x/sqrt(3.)
+            i1 = floor(x/(sqrt(0.75) * self.pitch[0]))
+            i2 = floor(alpha/self.pitch[0])
         # Check four lattice elements to see which one is closest based on local
         # coordinates
+        indices = [(i1, i2, iz), (i1 + 1, i2, iz), (i1, i2 + 1, iz),
+                   (i1 + 1, i2 + 1, iz)]
         d_min = np.inf
-        for idx in [(ix, ia, iz), (ix + 1, ia, iz), (ix, ia + 1, iz),
-                    (ix + 1, ia + 1, iz)]:
+
+        for idx in indices:
             p = self.get_local_coordinates(point, idx)
             d = p[0]**2 + p[1]**2
             if d < d_min:
@@ -1085,7 +1176,8 @@ class HexLattice(Lattice):
         point : Iterable of float
             Cartesian coordinates of point
         idx : Iterable of int
-            Indices of lattice element in :math:`(x,\alpha,z)` bases
+            Indices of lattice element in :math:`(x,\alpha,z)`
+            or :math:`(\alpha,y,z)` bases
 
         Returns
         -------
@@ -1094,8 +1186,17 @@ class HexLattice(Lattice):
             system
 
         """
-        x = point[0] - (self.center[0] + sqrt(0.75)*self.pitch[0]*idx[0])
-        y = point[1] - (self.center[1] + (0.5*idx[0] + idx[1])*self.pitch[0])
+        if self._orientation == 'x':
+            x = point[0] - (self.center[0] + self.pitch[0]*idx[0] +
+                            0.5*self.pitch[0]*idx[1])
+            y = point[1] - (self.center[1] +
+                            sqrt(0.75)*self.pitch[0]*idx[1])
+        else:
+            x = point[0] - (self.center[0]
+                            + sqrt(0.75)*self.pitch[0]*idx[0])
+            y = point[1] - (self.center[1]
+                            + (0.5*idx[0] + idx[1])*self.pitch[0])
+
         if self._num_axial is None:
             z = point[2]
         else:
@@ -1104,18 +1205,21 @@ class HexLattice(Lattice):
         return (x, y, z)
 
     def get_universe_index(self, idx):
-        r"""Return index in the universes array corresponding to a lattice element index
+        r"""Return index in the universes array corresponding
+        to a lattice element index
 
         Parameters
         ----------
         idx : Iterable of int
             Lattice element indices in the :math:`(x,\alpha,z)` coordinate
-            system
+            system in 'y' orientation case, or indices in the
+            :math:`(\alpha,y,z)` coordinate system in 'x' one
 
         Returns
         -------
-        2- or 3-tuple of int
-            Indices used when setting the :attr:`HexLattice.universes` property
+            2- or 3-tuple of int
+            Indices used when setting the :attr:`HexLattice.universes`
+            property
 
         """
 
@@ -1138,6 +1242,9 @@ class HexLattice(Lattice):
             else:
                 i_within = 5*g - z
 
+        if self._orientation == 'x' and g > 0:
+            i_within = (i_within + 5*g) % (6*g)
+
         if self.num_axial is None:
             return (i_ring, i_within)
         else:
@@ -1149,8 +1256,8 @@ class HexLattice(Lattice):
         Parameters
         ----------
         idx : Iterable of int
-            Lattice element indices in the :math:`(x,\alpha,z)` coordinate
-            system
+            Lattice element indices in the both :math:`(x,\alpha,z)`
+            and :math:`(\alpha,y,z)` coordinate system
 
         Returns
         -------
@@ -1193,6 +1300,9 @@ class HexLattice(Lattice):
             self._outer.create_xml_subelement(xml_element)
 
         lattice_subelement.set("n_rings", str(self._num_rings))
+        # If orientation is "x" export it to XML
+        if self._orientation == 'x':
+            lattice_subelement.set("orientation", "x")
 
         if self._num_axial is not None:
             lattice_subelement.set("n_axial", str(self._num_axial))
@@ -1267,6 +1377,7 @@ class HexLattice(Lattice):
         lat = cls(lat_id, name)
         lat.center = [float(i) for i in get_text(elem, 'center').split()]
         lat.pitch = [float(i) for i in get_text(elem, 'pitch').split()]
+        lat.orientation = get_text(elem, 'orientation', 'y')
         outer = get_text(elem, 'outer')
         if outer is not None:
             lat.outer = get_universe(int(outer))
@@ -1277,13 +1388,13 @@ class HexLattice(Lattice):
 
         # Create empty nested lists for one axial level
         univs = [[None for _ in range(max(6*(n_rings - 1 - r), 1))]
-                    for r in range(n_rings)]
+                 for r in range(n_rings)]
         if n_axial > 1:
             univs = [deepcopy(univs) for i in range(n_axial)]
 
-        # Get flat array of universes numbers
+        # Get flat array of universes
         uarray = np.array([get_universe(int(i)) for i in
-                            get_text(elem, 'universes').split()])
+                           get_text(elem, 'universes').split()])
 
         # Fill nested lists
         j = 0
@@ -1291,34 +1402,174 @@ class HexLattice(Lattice):
             # Get list for a single axial level
             axial_level = univs[z] if n_axial > 1 else univs
 
-            # Start iterating from top
-            x, alpha = 0, n_rings - 1
-            while True:
-                # Set entry in list based on (x,alpha,z) coordinates
-                _, i_ring, i_within = lat.get_universe_index((x, alpha, z))
-                axial_level[i_ring][i_within] = uarray[j]
+            if lat.orientation == 'y':
+                # Start iterating from top
+                x, alpha = 0, n_rings - 1
+                while True:
+                    # Set entry in list based on (x,alpha,z) coordinates
+                    _, i_ring, i_within = lat.get_universe_index((x, alpha, z))
+                    axial_level[i_ring][i_within] = uarray[j]
 
-                # Move to the right
-                x += 2
-                alpha -= 1
-                if not lat.is_valid_index((x, alpha, z)):
-                    # Move down in y direction
-                    alpha += x - 1
-                    x = 1 - x
+                    # Move to the right
+                    x += 2
+                    alpha -= 1
                     if not lat.is_valid_index((x, alpha, z)):
-                        # Move to the right
-                        x += 2
-                        alpha -= 1
+                        # Move down in y direction
+                        alpha += x - 1
+                        x = 1 - x
                         if not lat.is_valid_index((x, alpha, z)):
-                            # Reached the bottom
+                            # Move to the right
+                            x += 2
+                            alpha -= 1
+                            if not lat.is_valid_index((x, alpha, z)):
+                                # Reached the bottom
+                                break
+                    j += 1
+            else:
+                # Start iterating from top
+                alpha, y = 1 - n_rings, n_rings - 1
+                while True:
+                    # Set entry in list based on (alpha,y,z) coordinates
+                    _, i_ring, i_within = lat.get_universe_index((alpha, y, z))
+                    axial_level[i_ring][i_within] = uarray[j]
+
+                    # Move to the right
+                    alpha += 1
+                    if not lat.is_valid_index((alpha, y, z)):
+                        # Move down to next row
+                        alpha = 1 - n_rings
+                        y -= 1
+
+                        # Check if we've reached the bottom
+                        if y == -n_rings:
                             break
-                j += 1
+
+                        while not lat.is_valid_index((alpha, y, z)):
+                            # Move to the right
+                            alpha += 1
+                    j += 1
+
         lat.universes = univs
         return lat
 
     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.
+        """
+        if self._orientation == 'x':
+            return self._repr_axial_slice_x(universes)
+        else:
+            return self._repr_axial_slice_y(universes)
+
+    def _repr_axial_slice_x(self, universes):
+        """Return string representation for the given 2D group of universes
+        in 'x' orientation case.
+
+        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(2*self._num_rings - 1)]
+        middle = 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
+
+            # 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 left across the bottom
+            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.
+                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 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
+
+            # Climb right across the top
+            for i in range(r):
+                # Add the universe.
+                universe = universes[r_prime][theta]
+                rows[y].append(id_form.format(universe._id))
+
+                # Translate the indices.
+                theta += 1
+
+            # 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
+
+        # Flip the rows and join each row into a single string.
+        rows = [pad.join(x) for x in rows]
+
+        # 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]
+
+        # Join the rows together and return the string.
+        universe_ids = '\n'.join(rows)
+        return universe_ids
+
+    def _repr_axial_slice_y(self, universes):
+        """Return string representation for the given 2D group of universes in
+        'y' orientation case..
+
         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.
@@ -1425,7 +1676,7 @@ class HexLattice(Lattice):
         return universe_ids
 
     @staticmethod
-    def show_indices(num_rings):
+    def _show_indices_y(num_rings):
         """Return a diagram of the hexagonal lattice layout with indices.
 
         This method can be used to show the proper indices to be used when
@@ -1527,3 +1778,124 @@ class HexLattice(Lattice):
 
         # Join the rows together and return the string.
         return '\n'.join(rows)
+
+    @staticmethod
+    def _show_indices_x(num_rings):
+        """Return a diagram of the hexagonal lattice with x orientation
+        layout with indices.
+
+        This method can be used to show the proper indices to be used when
+        setting the :attr:`HexLattice.universes` property. For example,running
+        this method with num_rings=3 will return the similar diagram::
+
+                      (0, 8)      (0, 9)      (0,10)
+
+                (0, 7)      (1, 4)      (1, 5)      (0,11)
+
+          (0, 6)      (1, 3)      (2, 0)      (1, 0)      (0, 0)
+
+                (0, 5)      (1, 2)      (1, 1)      (0, 1)
+
+                      (0, 4)      (0, 3)      (0, 2)
+
+        Parameters
+        ----------
+        num_rings : int
+            Number of rings in the hexagonal lattice
+
+        Returns
+        -------
+        str
+            Diagram of the hexagonal lattice showing indices in OX orientation
+
+        """
+
+        # Find the largest string and count the number of digits so we can
+        # properly pad the output string later
+        largest_index = 6*(num_rings - 1)
+        n_digits_index = len(str(largest_index))
+        n_digits_ring = len(str(num_rings - 1))
+        str_form = '({{:{}}},{{:{}}})'.format(n_digits_ring, n_digits_index)
+        pad = ' '*(n_digits_index + n_digits_ring + 3)
+
+        # Initialize the list for each row.
+        rows = [[] for i in range(2*num_rings - 1)]
+        middle = num_rings - 1
+
+        # Start with the degenerate first ring.
+        rows[middle] = [str_form.format(num_rings - 1, 0)]
+
+        # Add universes one ring at a time.
+        for r in range(1, num_rings):
+            # r_prime increments down while r increments up.
+            r_prime = num_rings - 1 - r
+            theta = 0
+            y = middle
+
+            for i in range(r):
+                # Climb down the bottom-right
+                rows[y].append(str_form.format(r_prime, theta))
+                y += 1
+                theta += 1
+
+            for i in range(r):
+                # Climb left across the bottom
+                rows[y].insert(0, str_form.format(r_prime, theta))
+                theta += 1
+
+            for i in range(r):
+                # Climb up the bottom-left
+                rows[y].insert(0, str_form.format(r_prime, theta))
+                y -= 1
+                theta += 1
+
+            for i in range(r):
+                # Climb up the top-left
+                rows[y].insert(0, str_form.format(r_prime, theta))
+                y -= 1
+                theta += 1
+
+            for i in range(r):
+                # Climb right across the top
+                rows[y].append(str_form.format(r_prime, theta))
+                theta += 1
+
+            for i in range(r):
+                # Climb down the top-right
+                rows[y].append(str_form.format(r_prime, theta))
+                y += 1
+                theta += 1
+
+        # Flip the rows and join each row into a single string.
+        rows = [pad.join(x) for x in rows]
+
+        # Pad the beginning of the rows so they line up properly.
+        for y in range(num_rings - 1):
+            rows[y] = (num_rings - 1 - y)*pad + rows[y]
+            rows[-1 - y] = (num_rings - 1 - y)*pad + rows[-1 - y]
+
+        # Join the rows together and return the string.
+        return '\n\n'.join(rows)
+
+    @staticmethod
+    def show_indices(num_rings, orientation="y"):
+        """Return a diagram of the hexagonal lattice layout with indices.
+
+        Parameters
+        ----------
+        num_rings : int
+            Number of rings in the hexagonal lattice
+        orientation : {"x", "y"}
+            Orientation of the hexagonal lattice
+
+        Returns
+        -------
+        str
+            Diagram of the hexagonal lattice showing indices
+
+        """
+
+        if orientation == 'x':
+            return HexLattice._show_indices_x(num_rings)
+        else:
+            return HexLattice._show_indices_y(num_rings)
diff --git a/src/lattice.cpp b/src/lattice.cpp
index 42e15013d..6f10ed358 100644
--- a/src/lattice.cpp
+++ b/src/lattice.cpp
@@ -440,6 +440,21 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
     is_3d_ = false;
   }
 
+  // Read the lattice orientation.  Default to 'y'.
+  if (check_for_node(lat_node, "orientation")) {
+    std::string orientation = get_node_value(lat_node, "orientation");
+    if (orientation == "y") {
+      orientation_ = Orientation::y;
+    } else if (orientation == "x") {
+      orientation_ = Orientation::x;
+    } else {
+      fatal_error("Unrecognized orientation '" + orientation
+                  + "' for lattice " + std::to_string(id_));
+    }
+  } else {
+    orientation_ = Orientation::y;
+  }
+
   // Read the lattice center.
   std::string center_str {get_node_value(lat_node, "center")};
   std::vector<std::string> center_words {split(center_str)};
@@ -482,14 +497,79 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
 
   // Parse the universes.
   // Universes in hexagonal lattices are stored in a manner that represents
-  // a skewed coordinate system: (x, alpha) rather than (x, y).  There is
+  // a skewed coordinate system: (x, alpha) in case of 'y' orientation
+  // and (alpha,y) in 'x' one rather than (x, y).  There is
   // no obvious, direct relationship between the order of universes in the
   // input and the order that they will be stored in the skewed array so
   // the following code walks a set of index values across the skewed array
   // in a manner that matches the input order.  Note that i_x = 0, i_a = 0
-  // corresponds to the center of the hexagonal lattice.
-
+  // or i_a = 0, i_y = 0 corresponds to the center of the hexagonal lattice.
   universes_.resize((2*n_rings_-1) * (2*n_rings_-1) * n_axial_, C_NONE);
+  if (orientation_ == Orientation::y) {
+    fill_lattice_y(univ_words);
+  } else {
+    fill_lattice_x(univ_words);
+  }
+}
+
+//==============================================================================
+
+void
+HexLattice::fill_lattice_x(const std::vector<std::string>& univ_words)
+{
+  int input_index = 0;
+  for (int m = 0; m < n_axial_; m++) {
+    // Initialize lattice indecies.
+    int i_a = -(n_rings_ - 1);
+    int i_y = n_rings_ - 1;
+
+    // Map upper region of hexagonal lattice which is found in the
+    // first n_rings-1 rows of the input.
+    for (int k = 0; k < n_rings_-1; k++) {
+
+      // Iterate over the input columns.
+      for (int j = 0; j < k+n_rings_; j++) {
+        int indx = (2*n_rings_-1)*(2*n_rings_-1) * m
+                    + (2*n_rings_-1) * (i_y+n_rings_-1)
+                    + (i_a+n_rings_-1);
+        universes_[indx] = std::stoi(univ_words[input_index]);
+        input_index++;
+        // Move to the next right neighbour cell
+        i_a += 1;
+      }
+
+      // Return the lattice index to the start of the current row.
+      i_a = -(n_rings_ - 1);
+      i_y -= 1;
+    }
+
+    // Map the lower region from the centerline of cart to down side
+    for (int k = 0; k < n_rings_; k++) {
+      // Walk the index to the lower-right neighbor of the last row start.
+      i_a = -(n_rings_ - 1) + k;
+
+      // Iterate over the input columns.
+      for (int j = 0; j < 2*n_rings_-k-1; j++) {
+        int indx = (2*n_rings_-1)*(2*n_rings_-1) * m
+                    + (2*n_rings_-1) * (i_y+n_rings_-1)
+                    + (i_a+n_rings_-1);
+        universes_[indx] = std::stoi(univ_words[input_index]);
+        input_index++;
+        // Move to the next right neighbour cell
+        i_a += 1;
+      }
+
+      // Return lattice index to start of current row.
+      i_y -= 1;
+    }
+  }
+}
+
+//==============================================================================
+
+void
+HexLattice::fill_lattice_y(const std::vector<std::string>& univ_words)
+{
   int input_index = 0;
   for (int m = 0; m < n_axial_; m++) {
     // Initialize lattice indecies.
@@ -610,9 +690,33 @@ std::pair<double, std::array<int, 3>>
 HexLattice::distance(Position r, Direction u, const std::array<int, 3>& i_xyz)
 const
 {
-  // Compute the direction on the hexagonal basis.
-  double beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0;
-  double gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0;
+  // Short description of the direction vectors used here.  The beta, gamma, and
+  // delta vectors point towards the flat sides of each hexagonal tile.
+  // Y - orientation:
+  //   basis0 = (1, 0)
+  //   basis1 = (-1/sqrt(3), 1)   = +120 degrees from basis0
+  //   beta   = (sqrt(3)/2, 1/2)  = +30 degrees from basis0
+  //   gamma  = (sqrt(3)/2, -1/2) = -60 degrees from beta
+  //   delta  = (0, 1)            = +60 degrees from beta
+  // X - orientation:
+  //   basis0 = (1/sqrt(3), -1)
+  //   basis1 = (0, 1)            = +120 degrees from basis0
+  //   beta   = (1, 0)            = +30 degrees from basis0
+  //   gamma  = (1/2, -sqrt(3)/2) = -60 degrees from beta
+  //   delta  = (1/2, sqrt(3)/2)  = +60 degrees from beta
+  // The z-axis is considered separately.
+  double beta_dir;
+  double gamma_dir;
+  double delta_dir;
+  if (orientation_ == Orientation::y) {
+    beta_dir = u.x * std::sqrt(3.0) / 2.0  + u.y / 2.0;
+    gamma_dir = u.x * std::sqrt(3.0) / 2.0  - u.y / 2.0;
+    delta_dir = u.y;
+  } else {
+    beta_dir = u.x;
+    gamma_dir = u.x / 2.0  - u.y * std::sqrt(3.0) / 2.0;
+    delta_dir = u.x / 2.0  + u.y * std::sqrt(3.0) / 2.0;
+  }
 
   // Note that hexagonal lattice distance calculations are performed
   // using the particle's coordinates relative to the neighbor lattice
@@ -620,7 +724,7 @@ const
   // because there is significant disagreement between neighboring cells
   // on where the lattice boundary is due to finite precision issues.
 
-  // Upper-right and lower-left sides.
+  // beta direction
   double d {INFTY};
   std::array<int, 3> lattice_trans;
   double edge = -copysign(0.5*pitch_[0], beta_dir);  // Oncoming edge
@@ -632,7 +736,12 @@ const
     const std::array<int, 3> i_xyz_t {i_xyz[0]-1, i_xyz[1], i_xyz[2]};
     r_t = get_local_position(r, i_xyz_t);
   }
-  double beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0;
+  double beta;
+  if (orientation_ == Orientation::y) {
+    beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0;
+  } else {
+    beta = r_t.x;
+  }
   if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) {
     d = (edge - beta) / beta_dir;
     if (beta_dir > 0) {
@@ -642,7 +751,7 @@ const
     }
   }
 
-  // Lower-right and upper-left sides.
+  // gamma direction
   edge = -copysign(0.5*pitch_[0], gamma_dir);
   if (gamma_dir > 0) {
     const std::array<int, 3> i_xyz_t {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]};
@@ -651,7 +760,12 @@ const
     const std::array<int, 3> i_xyz_t {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]};
     r_t = get_local_position(r, i_xyz_t);
   }
-  double gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0;
+  double gamma;
+  if (orientation_ == Orientation::y) {
+    gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0;
+  } else {
+    gamma = r_t.x  / 2.0 - r_t.y * std::sqrt(3.0) / 2.0;
+  }
   if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) {
     double this_d = (edge - gamma) / gamma_dir;
     if (this_d < d) {
@@ -664,19 +778,25 @@ const
     }
   }
 
-  // Upper and lower sides.
-  edge = -copysign(0.5*pitch_[0], u.y);
-  if (u.y > 0) {
+  // delta direction
+  edge = -copysign(0.5*pitch_[0], delta_dir);
+  if (delta_dir > 0) {
     const std::array<int, 3> i_xyz_t {i_xyz[0], i_xyz[1]+1, i_xyz[2]};
     r_t = get_local_position(r, i_xyz_t);
   } else {
     const std::array<int, 3> i_xyz_t {i_xyz[0], i_xyz[1]-1, i_xyz[2]};
     r_t = get_local_position(r, i_xyz_t);
   }
-  if ((std::abs(r_t.y - edge) > FP_PRECISION) && u.y != 0) {
-    double this_d = (edge - r_t.y) / u.y;
+  double delta;
+  if (orientation_ == Orientation::y) {
+    delta =  r_t.y;
+  } else {
+    delta = r_t.x  / 2.0 + r_t.y * std::sqrt(3.0) / 2.0;
+  }
+  if ((std::abs(delta - edge) > FP_PRECISION) && delta_dir != 0) {
+    double this_d = (edge - delta) / delta_dir;
     if (this_d < d) {
-      if (u.y > 0) {
+      if (delta_dir > 0) {
         lattice_trans = {0, 1, 0};
       } else {
         lattice_trans = {0, -1, 0};
@@ -727,16 +847,25 @@ HexLattice::get_indices(Position r, Direction u) const
     }
   }
 
-  // Convert coordinates into skewed bases.  The (x, alpha) basis is used to
-  // find the index of the global coordinates to within 4 cells.
-  double alpha = r_o.y - r_o.x / std::sqrt(3.0);
-  int ix = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0]));
-  int ia = std::floor(alpha / pitch_[0]);
+  int i1, i2;
+  if (orientation_ == Orientation::y) {
+    // Convert coordinates into skewed bases.  The (x, alpha) basis is used to
+    // find the index of the global coordinates to within 4 cells.
+    double alpha = r_o.y - r_o.x / std::sqrt(3.0);
+    i1 = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0]));
+    i2 = std::floor(alpha / pitch_[0]);
+  } else {
+    // Convert coordinates into skewed bases.  The (alpha, y) basis is used to
+    // find the index of the global coordinates to within 4 cells.
+    double alpha = r_o.y - r_o.x * std::sqrt(3.0);
+    i1 = std::floor(-alpha / (std::sqrt(3.0) * pitch_[0]));
+    i2 = std::floor(r_o.y / (0.5*std::sqrt(3.0) * pitch_[0]));
+  }
 
-  // Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but
+  // Add offset to indices (the center cell is (i1, i2) = (0, 0) but
   // the array is offset so that the indices never go below 0).
-  ix += n_rings_-1;
-  ia += n_rings_-1;
+  i1 += n_rings_-1;
+  i2 += n_rings_-1;
 
   // Calculate the (squared) distance between the particle and the centers of
   // the four possible cells.  Regular hexagonal tiles form a Voronoi
@@ -755,14 +884,14 @@ HexLattice::get_indices(Position r, Direction u) const
   // is kept (i.e. the cell with the lowest dot product as the vectors will be
   // completely opposed if the particle is moving directly toward the center of
   // the cell).
-  int ix_chg {};
-  int ia_chg {};
+  int i1_chg {};
+  int i2_chg {};
   double d_min {INFTY};
   double dp_min {INFTY};
   for (int i = 0; i < 2; i++) {
     for (int j = 0; j < 2; j++) {
       // get local coordinates
-      const std::array<int, 3> i_xyz {ix + j, ia + i, 0};
+      const std::array<int, 3> i_xyz {i1 + j, i2 + i, 0};
       Position r_t = get_local_position(r, i_xyz);
       // calculate distance
       double d = r_t.x*r_t.x + r_t.y*r_t.y;
@@ -777,18 +906,18 @@ HexLattice::get_indices(Position r, Direction u) const
         if (on_boundary && dp > dp_min) continue;
         // update values
         d_min = d;
-        ix_chg = j;
-        ia_chg = i;
+        i1_chg = j;
+        i2_chg = i;
         dp_min = dp;
       }
     }
   }
 
   // update outgoing indices
-  ix += ix_chg;
-  ia += ia_chg;
+  i1 += i1_chg;
+  i2 += i2_chg;
 
-  return {ix, ia, iz};
+  return {i1, i2, iz};
 }
 
 //==============================================================================
@@ -797,15 +926,26 @@ Position
 HexLattice::get_local_position(Position r, const std::array<int, 3> i_xyz)
 const
 {
-  // x_l = x_g - (center + pitch_x*cos(30)*index_x)
-  r.x -= center_.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0];
-  // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y)
-  r.y -= (center_.y + (i_xyz[1] - n_rings_ + 1) * pitch_[0]
-          + (i_xyz[0] - n_rings_ + 1) * pitch_[0] / 2.0);
-  if (is_3d_) {
-    r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1];
+  if (orientation_ == Orientation::y) {
+    // x_l = x_g - (center + pitch_x*cos(30)*index_x)
+    r.x -= center_.x
+           + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings_ + 1) * pitch_[0];
+    // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y)
+    r.y -= (center_.y + (i_xyz[1] - n_rings_ + 1) * pitch_[0]
+            + (i_xyz[0] - n_rings_ + 1) * pitch_[0] / 2.0);
+  } else {
+    // x_l = x_g - (center + pitch_x*index_a + pitch_y*sin(30)*index_y)
+    r.x -= (center_.x + (i_xyz[0] - n_rings_ + 1) * pitch_[0]
+            + (i_xyz[1] - n_rings_ + 1) * pitch_[0] / 2.0);
+    // y_l = y_g - (center + pitch_y*cos(30)*index_y)
+    r.y -= center_.y
+           + std::sqrt(3.0)/2.0 * (i_xyz[1] - n_rings_ + 1) * pitch_[0];
   }
 
+  if (is_3d_) {
+      r.z -= center_.z - (0.5 * n_axial_ - i_xyz[2] - 0.5) * pitch_[1];
+    }
+
   return r;
 }
 
@@ -863,6 +1003,11 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const
   write_string(lat_group, "type", "hexagonal", false);
   write_dataset(lat_group, "n_rings", n_rings_);
   write_dataset(lat_group, "n_axial", n_axial_);
+  if (orientation_ == Orientation::y) {
+    write_string(lat_group, "orientation", "y", false);
+  } else {
+    write_string(lat_group, "orientation", "x", false);
+  }
   if (is_3d_) {
     write_dataset(lat_group, "pitch", pitch_);
     write_dataset(lat_group, "center", center_);
diff --git a/src/relaxng/geometry.rnc b/src/relaxng/geometry.rnc
index c9c53f6d6..e3c88c445 100644
--- a/src/relaxng/geometry.rnc
+++ b/src/relaxng/geometry.rnc
@@ -48,6 +48,7 @@ element geometry {
     (element n_axial { xsd:int } | attribute n_axial { xsd:int })? &
     (element center { list { xsd:double+ } } | attribute center { list { xsd:double+ } }) &
     (element pitch { list { xsd:double+ } } | attribute pitch { list { xsd:double+ } }) &
+    (element orientation { ( "x" | "y" ) } | attribute orientation { ( "x" | "y" ) })? &
     (element universes { list { xsd:int+ } } | attribute universes { list { xsd:int+ } }) &
     (element outer { xsd:int } | attribute outer { xsd:int })?
   }*
diff --git a/src/relaxng/geometry.rng b/src/relaxng/geometry.rng
index d2c1f635e..56bf38580 100644
--- a/src/relaxng/geometry.rng
+++ b/src/relaxng/geometry.rng
@@ -398,6 +398,22 @@
               </list>
             </attribute>
           </choice>
+          <optional>
+            <choice>
+              <element name="orientation">
+                <choice>
+                  <value>x</value>
+                  <value>y</value>
+                </choice>
+              </element>
+              <attribute name="orientation">
+                <choice>
+                  <value>x</value>
+                  <value>y</value>
+                </choice>
+              </attribute>
+            </choice>
+          </optional>
           <choice>
             <element name="universes">
               <list>
diff --git a/tests/regression_tests/lattice_hex_x/__init__.py b/tests/regression_tests/lattice_hex_x/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/regression_tests/lattice_hex_x/inputs_true.dat b/tests/regression_tests/lattice_hex_x/inputs_true.dat
new file mode 100644
index 000000000..b528a97b7
--- /dev/null
+++ b/tests/regression_tests/lattice_hex_x/inputs_true.dat
@@ -0,0 +1,123 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="1" material="1" region="-1" universe="1" />
+  <cell id="2" material="4" region="-2 1" universe="1" />
+  <cell id="3" material="2" region="2" universe="1" />
+  <cell id="4" material="2" region="-3" universe="2" />
+  <cell id="5" material="4" region="-4 3" universe="2" />
+  <cell id="6" material="2" region="4" universe="2" />
+  <cell id="7" material="3" region="-5" universe="3" />
+  <cell id="8" material="4" region="-6 5" universe="3" />
+  <cell id="9" material="2" region="-7 6" universe="3" />
+  <cell id="10" material="4" region="-8 7" universe="3" />
+  <cell id="11" material="2" region="8" universe="3" />
+  <cell id="12" material="2" universe="4" />
+  <cell fill="9" id="13" name="container assembly cell" region="-11 12 -13 14 15 -16 9 -10" universe="5" />
+  <hex_lattice id="9" n_axial="2" n_rings="11" name="regular fuel assembly" orientation="x">
+    <pitch>1.235 5.0</pitch>
+    <outer>4</outer>
+    <center>0.0 0.0 5.0</center>
+    <universes>
+          1 1 1 1 1 1 1 1 1 1 1
+         1 1 1 1 1 1 1 1 1 1 1 1
+        1 1 1 1 1 1 1 1 1 1 1 1 1
+       1 1 1 1 1 1 1 1 1 1 1 1 1 1
+      1 1 1 1 1 1 1 3 1 1 1 1 1 1 1
+     1 1 1 1 1 3 1 1 1 1 3 1 1 1 1 1
+    1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+   1 1 1 1 3 1 1 1 1 3 1 1 1 3 1 1 1 1
+  1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1
+1 1 1 1 1 3 1 1 1 1 2 1 1 1 1 3 1 1 1 1 1
+ 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1
+  1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1
+   1 1 1 1 3 1 1 1 3 1 1 1 1 3 1 1 1 1
+    1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+     1 1 1 1 1 3 1 1 1 1 3 1 1 1 1 1
+      1 1 1 1 1 1 1 3 1 1 1 1 1 1 1
+       1 1 1 1 1 1 1 1 1 1 1 1 1 1
+        1 1 1 1 1 1 1 1 1 1 1 1 1
+         1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1
+         1 1 1 1 1 1 1 1 1 1 1 1
+        1 1 1 1 1 1 1 1 1 1 1 1 1
+       1 1 1 1 1 1 1 1 1 1 1 1 1 1
+      1 1 1 1 1 1 1 3 1 1 1 1 1 1 1
+     1 1 1 1 1 3 1 1 1 1 3 1 1 1 1 1
+    1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+   1 1 1 1 3 1 1 1 1 3 1 1 1 3 1 1 1 1
+  1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1
+ 1 1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1
+1 1 1 1 1 3 1 1 1 1 2 1 1 1 1 3 1 1 1 1 1
+ 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 1 1 1 1
+  1 1 1 1 1 1 1 1 1 1 1 3 1 1 1 1 1 1 1
+   1 1 1 1 3 1 1 1 3 1 1 1 1 3 1 1 1 1
+    1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
+     1 1 1 1 1 3 1 1 1 1 3 1 1 1 1 1
+      1 1 1 1 1 1 1 3 1 1 1 1 1 1 1
+       1 1 1 1 1 1 1 1 1 1 1 1 1 1
+        1 1 1 1 1 1 1 1 1 1 1 1 1
+         1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1</universes>
+  </hex_lattice>
+  <surface coeffs="0.0 0.0 0.386" id="1" type="z-cylinder" />
+  <surface coeffs="0.0 0.0 0.4582" id="2" type="z-cylinder" />
+  <surface coeffs="0.0 0.0 0.45" id="3" type="z-cylinder" />
+  <surface coeffs="0.0 0.0 0.5177" id="4" type="z-cylinder" />
+  <surface coeffs="0.0 0.0 0.35" id="5" type="z-cylinder" />
+  <surface coeffs="0.0 0.0 0.41" id="6" type="z-cylinder" />
+  <surface coeffs="0.0 0.0 0.545" id="7" type="z-cylinder" />
+  <surface coeffs="0.0 0.0 0.6323" id="8" type="z-cylinder" />
+  <surface boundary="reflective" coeffs="0.0" id="9" type="z-plane" />
+  <surface boundary="reflective" coeffs="10.0" id="10" type="z-plane" />
+  <surface boundary="reflective" coeffs="11.8" id="11" type="y-plane" />
+  <surface boundary="reflective" coeffs="-11.8" id="12" type="y-plane" />
+  <surface boundary="reflective" coeffs="1.7320508075688772 1.0 0.0 23.6" id="13" type="plane" />
+  <surface boundary="reflective" coeffs="-1.7320508075688772 1.0 0.0 -23.6" id="14" type="plane" />
+  <surface boundary="reflective" coeffs="1.7320508075688772 1.0 0.0 -23.6" id="15" type="plane" />
+  <surface boundary="reflective" coeffs="-1.7320508075688772 1.0 0.0 23.6" id="16" type="plane" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material depletable="true" id="1" name="UO2">
+    <density units="sum" />
+    <nuclide ao="0.0008737" name="U235" />
+    <nuclide ao="0.018744" name="U238" />
+    <nuclide ao="0.039235" name="O16" />
+  </material>
+  <material id="2" name="borated H2O">
+    <density units="sum" />
+    <nuclide ao="0.06694" name="H1" />
+    <nuclide ao="0.03347" name="O16" />
+    <nuclide ao="6.6262e-06" name="B10" />
+    <nuclide ao="2.6839e-05" name="B11" />
+  </material>
+  <material id="3" name="pellet B4C">
+    <density units="sum" />
+    <nuclide ao="0.01966" name="C0" />
+    <nuclide ao="4.7344e-06" name="B11" />
+    <nuclide ao="1.9177e-05" name="B10" />
+  </material>
+  <material id="4" name="Zirc4">
+    <density units="sum" />
+    <nuclide ao="0.021763349999999997" name="Zr90" />
+    <nuclide ao="0.00474606" name="Zr91" />
+    <nuclide ao="0.00725445" name="Zr92" />
+    <nuclide ao="0.00735174" name="Zr94" />
+    <nuclide ao="0.0011844" name="Zr96" />
+  </material>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>eigenvalue</run_mode>
+  <particles>1000</particles>
+  <batches>10</batches>
+  <inactive>5</inactive>
+  <source strength="1.0">
+    <space type="box">
+      <parameters>-13.62546635287517 -13.62546635287517 0.0 13.62546635287517 13.62546635287517 10.0</parameters>
+    </space>
+  </source>
+  <seed>22</seed>
+</settings>
diff --git a/tests/regression_tests/lattice_hex_x/results_true.dat b/tests/regression_tests/lattice_hex_x/results_true.dat
new file mode 100644
index 000000000..a742aa0d5
--- /dev/null
+++ b/tests/regression_tests/lattice_hex_x/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+1.355663E+00 2.896562E-02
diff --git a/tests/regression_tests/lattice_hex_x/test.py b/tests/regression_tests/lattice_hex_x/test.py
new file mode 100644
index 000000000..538f790da
--- /dev/null
+++ b/tests/regression_tests/lattice_hex_x/test.py
@@ -0,0 +1,206 @@
+from tests.testing_harness import PyAPITestHarness
+import openmc
+import numpy as np
+
+
+class HexLatticeOXTestHarness(PyAPITestHarness):
+
+    def _build_inputs(self):
+        materials = openmc.Materials()
+
+        fuel_mat = openmc.Material(material_id=1, name="UO2")
+        fuel_mat.set_density('sum')
+        fuel_mat.add_nuclide('U235', 0.87370e-03)
+        fuel_mat.add_nuclide('U238', 1.87440e-02)
+        fuel_mat.add_nuclide('O16', 3.92350e-02)
+        materials.append(fuel_mat)
+
+        coolant = openmc.Material(material_id=2, name="borated H2O")
+        coolant.set_density('sum')
+        coolant.add_nuclide('H1', 0.06694)
+        coolant.add_nuclide('O16', 0.03347)
+        coolant.add_nuclide('B10', 6.6262e-6)
+        coolant.add_nuclide('B11', 2.6839e-5)
+        materials.append(coolant)
+
+        absorber = openmc.Material(material_id=3, name="pellet B4C")
+        absorber.set_density('sum')
+        absorber.add_nuclide('C0', 0.01966)
+        absorber.add_nuclide('B11', 4.7344e-6)
+        absorber.add_nuclide('B10', 1.9177e-5)
+        materials.append(absorber)
+
+        zirc = openmc.Material(material_id=4, name="Zirc4")
+        zirc.set_density('sum')
+        zirc.add_element('Zr', 4.23e-2)
+        materials.append(zirc)
+
+        materials.export_to_xml()
+
+        # Geometry #
+
+        pin_rad = 0.7 # cm
+        assembly_pitch = 1.235 # cm
+        hexagonal_pitch = 23.6 # cm
+        length = 10.0 # cm
+
+        # Fuel pin surfaces
+
+        cylfuelin = openmc.ZCylinder(surface_id=1, r=0.386)
+        cylfuelout = openmc.ZCylinder(surface_id=2, r=0.4582)
+
+        # Fuel cells
+
+        infcell = openmc.Cell(cell_id=1)
+        infcell.region = -cylfuelin
+        infcell.fill = fuel_mat
+
+        clfcell = openmc.Cell(cell_id=2)
+        clfcell.region = -cylfuelout & +cylfuelin
+        clfcell.fill = zirc
+
+        outfcell = openmc.Cell(cell_id=3)
+        outfcell.region = +cylfuelout
+        outfcell.fill = coolant
+
+        # Fuel universe
+
+        fuel_ch_univ = openmc.Universe(universe_id=1, name="Fuel channel",
+                                       cells=[infcell, clfcell, outfcell])
+
+        # Central tube surfaces
+
+        cyltubein = openmc.ZCylinder(surface_id=3, r=0.45)
+        cyltubeout = openmc.ZCylinder(surface_id=4, r=0.5177)
+
+        # Central tube cells
+
+        inctcell = openmc.Cell(cell_id=4)
+        inctcell.region = -cyltubein
+        inctcell.fill = coolant
+
+        clctcell = openmc.Cell(cell_id=5)
+        clctcell.region = -cyltubeout & +cyltubein
+        clctcell.fill = zirc
+
+        outctcell = openmc.Cell(cell_id=6)
+        outctcell.region = +cyltubeout
+        outctcell.fill = coolant
+
+        # Central tubel universe
+
+        tube_ch_univ = openmc.Universe(universe_id=2,
+                                       name="Central tube channel",
+                                       cells=[inctcell, clctcell, outctcell])
+
+        # Absorber tube surfaces
+
+        cylabsin = openmc.ZCylinder(surface_id=5, r=0.35)
+        cylabsout = openmc.ZCylinder(surface_id=6, r=0.41)
+        cylabsclin = openmc.ZCylinder(surface_id=7, r=0.545)
+        cylabsclout = openmc.ZCylinder(surface_id=8, r=0.6323)
+
+        # Absorber tube cells
+
+        inabscell = openmc.Cell(cell_id=7)
+        inabscell.region = -cylabsin
+        inabscell.fill = absorber
+
+        clabscell = openmc.Cell(cell_id=8)
+        clabscell.region = -cylabsout & +cylabsin
+        clabscell.fill = zirc
+
+        interabscell = openmc.Cell(cell_id=9)
+        interabscell.region = -cylabsclin & +cylabsout
+        interabscell.fill = coolant
+
+        clatcell = openmc.Cell(cell_id=10)
+        clatcell.region = -cylabsclout & +cylabsclin
+        clatcell.fill = zirc
+
+        outabscell = openmc.Cell(cell_id=11)
+        outabscell.region = +cylabsclout
+        outabscell.fill = coolant
+
+        # Absorber tube universe
+
+        abs_ch_univ = openmc.Universe(universe_id=3,
+                                      name="Central tube channel",
+                                      cells=[inabscell, clabscell,
+                                             interabscell,
+                                             clatcell, outabscell])
+        # Assembly surfaces
+
+        edge_length = (1./np.sqrt(3.0)) * hexagonal_pitch
+        fuel_bottom = openmc.ZPlane(surface_id=9, z0=0.0,
+                                    boundary_type='reflective')
+        fuel_top = openmc.ZPlane(surface_id=10, z0=length,
+                                 boundary_type='reflective')
+
+        # a hex surface for the core to go inside of
+
+        hexprism = openmc.model.get_hexagonal_prism(edge_length=edge_length,
+                                                    origin=(0.0,  0.0),
+                                                    boundary_type='reflective',
+                                                    orientation='x')
+        region = hexprism & +fuel_bottom & -fuel_top
+
+        inf_mat = openmc.Cell(cell_id=12)
+        inf_mat.fill = coolant
+        inf_mat_univ = openmc.Universe(universe_id=4, cells=[inf_mat])
+
+        # Fill lattice by channels
+
+        nring = 11
+        universes = []
+        for ring in range(nring - 1, -1, -1):
+            arr = []
+            arr.append(fuel_ch_univ)
+            for cell in range(ring * 6 - 1):
+                arr.append(fuel_ch_univ)
+            universes.append(arr)
+        universes[-1] = [tube_ch_univ]
+        channels = [(7, 2), (7, 5), (7, 8), (7, 11), (7, 14), (7, 17), (5, 0),
+                    (4, 3), (5, 5), (4, 9), (5, 10), (4, 15), (5, 15),
+                    (4, 21), (5, 20), (4, 27), (5, 25), (4, 33)]
+        for i, j in channels:
+            universes[i][j] = abs_ch_univ
+        lattice = openmc.HexLattice(name="regular fuel assembly")
+        lattice.orientation = "x"
+        lattice.center = (0., 0., length/2.0)
+        lattice.pitch = (assembly_pitch, length/2.0)
+        lattice.universes = 2*[universes]
+        lattice.outer = inf_mat_univ
+
+        assembly_cell = openmc.Cell(cell_id=13,
+                                    name="container assembly cell")
+        assembly_cell.region = region
+        assembly_cell.fill = lattice
+
+        root_univ = openmc.Universe(universe_id=5, name="root universe",
+                                    cells=[assembly_cell])
+
+        geom = openmc.Geometry(root_univ)
+        geom.export_to_xml()
+
+        # Settings #
+
+        settings = openmc.Settings()
+        settings.run_mode = 'eigenvalue'
+
+        source = openmc.Source()
+        ll = [-edge_length, -edge_length, 0.0]
+        ur = [edge_length, edge_length, 10.0]
+        source.space = openmc.stats.Box(ll, ur)
+        source.strength = 1.0
+        settings.source = source
+        settings.batches = 10
+        settings.inactive = 5
+        settings.particles = 1000
+        settings.seed = 22
+        settings.export_to_xml()
+
+
+def test_lattice_hex_ox_surf():
+    harness = HexLatticeOXTestHarness('statepoint.10.h5')
+    harness.main()
diff --git a/tests/unit_tests/test_lattice.py b/tests/unit_tests/test_lattice.py
index fe86e2a78..c4a351776 100644
--- a/tests/unit_tests/test_lattice.py
+++ b/tests/unit_tests/test_lattice.py
@@ -202,6 +202,13 @@ def test_get_universe(rlat2, rlat3, hlat2, hlat3):
     assert hlat2.get_universe((1, 0)) == u1
     assert hlat2.get_universe((-2, 2)) == u1
 
+    hlat2.orientation = 'x'
+    assert hlat2.get_universe((2, 0)) == u2
+    assert hlat2.get_universe((1, 0)) == u2
+    assert hlat2.get_universe((1, 1)) == u1
+    assert hlat2.get_universe((-1, 1)) == u1
+    hlat2.orientation = 'y'
+
     u1, u2, u3, outer = hlat3.univs
     assert hlat3.get_universe((0, 0, 0)) == u2
     assert hlat3.get_universe((0, 0, 1)) == u3
@@ -342,3 +349,5 @@ def test_show_indices():
     for i in range(1, 11):
         lines = openmc.HexLattice.show_indices(i).split('\n')
         assert len(lines) == 4*i - 3
+        lines_x = openmc.HexLattice.show_indices(i, 'x').split('\n')
+        assert len(lines) == 4*i - 3