diff --git a/openmc/filter.py b/openmc/filter.py
index 8a54e05cee..f9f4f74766 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -376,7 +376,7 @@ class Filter(object):
             cell instance ID for 'distribcell' Filters. The bin is a 2-tuple of
             floats for 'energy' and 'energyout' filters corresponding to the
             energy boundaries of the bin of interest. The bin is an (x,y,z)
-            3-tuple for 'mesh' filters corresponding to the mesh cell
+            3-tuple for 'mesh' filters corresponding to the mesh cell of
             interest.
 
         Returns
diff --git a/openmc/mesh.py b/openmc/mesh.py
index 169e7705cc..9496a06511 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -6,6 +6,8 @@ import sys
 import numpy as np
 
 import openmc.checkvalue as cv
+import openmc.surface
+import openmc.cell
 
 
 if sys.version_info[0] >= 3:
@@ -181,6 +183,21 @@ class Mesh(object):
         string += '{0: <16}{1}{2}\n'.format('\tPixels', '=\t', self._width)
         return string
 
+    def cell_generator(self):
+        """Generator function to traverse through every [i,j,k] index
+        of the mesh.
+
+        """
+        if len(self.dimension) == 2:
+            for x in range(self.dimension[0]):
+                for y in range(self.dimension[1]):
+                    yield [x + 1, y + 1, 1]
+        else:
+            for x in range(self.dimension[0]):
+                for y in range(self.dimension[1]):
+                    for z in range(self.dimension[2]):
+                        yield [x + 1, y + 1, z + 1]
+
     def get_mesh_xml(self):
         """Return XML representation of the mesh
 
@@ -210,3 +227,90 @@ class Mesh(object):
             subelement.text = ' '.join(map(str, self._width))
 
         return element
+
+    def build_cells(self, bc=['reflective'] * 6):
+        """Generates a list of cells which mimic the mesh geometry.
+
+        Parameters
+        ----------
+        bc : iterable of {'reflective', 'periodic', or 'vacuum'}
+            Boundary conditions for each of the six faces of a parallelopiped:
+            North, East, South, West, Up, and Down.
+
+        Returns
+        -------
+        root_cell : openmc.Cell
+            The cell containing the lattice mimicking the mesh geometry.
+        cells : iterable of openmc.Cell
+            The list of cells within each lattice position mimicking the mesh
+            geometry.
+
+        """
+
+        cv.check_length('bc', bc, length_min=6, length_max=6)
+        for entry in bc:
+            cv.check_value('bc', entry, ['transmission', 'vacuum',
+                                         'reflective', 'periodic'])
+
+        if len(self.dimension) == 2:
+            twod = True
+        else:
+            twod = False
+
+        # Build enclosing cell
+        xplanes = [openmc.XPlane(x0=self.lower_left[0],
+                                 boundary_type=bc[3]),
+                   openmc.XPlane(x0=self.upper_right[0],
+                                 boundary_type=bc[1])]
+        yplanes = [openmc.YPlane(y0=self.lower_left[1],
+                                 boundary_type=bc[2]),
+                   openmc.YPlane(y0=self.upper_right[1],
+                                 boundary_type=bc[0])]
+        if twod:
+            zplanes = [openmc.ZPlane(z0=-1.E50, boundary_type=bc[5]),
+                       openmc.ZPlane(z0=+1.E50, boundary_type=bc[5])]
+        else:
+            zplanes = [openmc.ZPlane(z0=self.lower_left[2],
+                                     boundary_type=bc[5]),
+                       openmc.ZPlane(z0=self.upper_right[2],
+                                     boundary_type=bc[5])]
+        root_cell = openmc.Cell()
+        root_cell.region = ((+xplanes[0] & -xplanes[1]) &
+                            (+yplanes[0] & -yplanes[1]) &
+                            (+zplanes[0] & -zplanes[1]))
+
+        # Build our universes
+        universes = np.ndarray(self.dimension[::-1], dtype=np.object)
+        cells = []
+        for [i, j, k] in self.cell_generator():
+            if twod:
+                universes[j - 1, i - 1] = openmc.Universe()
+                cells.append(openmc.Cell())
+                universes[j - 1, i - 1].add_cells([cells[-1]])
+            else:
+                universes[k - 1, j - 1, i - 1] = openmc.Universe()
+                cells.append(openmc.Cell())
+                universes[k - 1, j - 1, i - 1].add_cells([cells[-1]])
+
+        lattice = openmc.RectLattice()
+        lattice.lower_left = self.lower_left
+
+        if self.width is not None:
+            lattice.pitch = self.width
+        else:
+            dx = ((self.upper_right[0] - self.lower_left[0]) /
+                  self.dimension[0])
+            dy = ((self.upper_right[1] - self.lower_left[1]) /
+                  self.dimension[1])
+            if twod:
+                lattice.pitch = [dx, dy]
+            else:
+                dz = ((self.upper_right[2] - self.lower_left[2]) /
+                      self.dimension[2])
+                lattice.pitch = [dx, dy, dz]
+        lattice.universes = universes
+
+        # Fill Cell with the Lattice
+        root_cell.fill = lattice
+
+        return root_cell, cells
diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py
index 84571ce236..c5e039cea5 100644
--- a/openmc/mgxs/library.py
+++ b/openmc/mgxs/library.py
@@ -756,14 +756,14 @@ class Library(object):
 
     def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro',
                    xs_id='1m', order=None, tabular_legendre=None,
-                   tabular_points=33):
+                   tabular_points=33, subdomain=None):
         """Generates an openmc.XSdata object describing a multi-group cross section
         data set for eventual combination in to an openmc.MGXSLibrary object
         (i.e., the library).
 
         Parameters
         ----------
-        domain : openmc.Material or openmc.Cell or openmc.Universe
+        domain : openmc.Material, openmc.Cell, openmc.Universe, or openmc.Mesh
             The domain for spatial homogenization
         xsdata_name : str
             Name to apply to the "xsdata" entry produced by this method
@@ -792,6 +792,10 @@ class Library(object):
             parameter is set to `True`.  In this case, this parameter sets the
             number of equally-spaced points in the domain of [-1,1] to be used
             in building the tabular distribution. Default is `33`.
+        subdomain : iterable of int
+            This parameter is not used unless using a mesh domain. In that
+            case, the subdomain is an [i,j,k] index (1-based indexing) of the
+            mesh cell of interest in the openmc.Mesh object.
 
         Returns
         -------
@@ -811,7 +815,7 @@ class Library(object):
         """
 
         cv.check_type('domain', domain, (openmc.Material, openmc.Cell,
-                                         openmc.Cell))
+                                         openmc.Universe, openmc.Mesh))
         cv.check_type('xsdata_name', xsdata_name, basestring)
         cv.check_type('nuclide', nuclide, basestring)
         cv.check_value('xs_type', xs_type, ['macro', 'micro'])
@@ -824,6 +828,9 @@ class Library(object):
                       (type(None), bool))
         if tabular_points is not None:
             cv.check_greater_than('tabular_points', tabular_points, 1)
+        if subdomain is not None:
+            cv.check_iterable_type('subdomain', subdomain, Integral,
+                                   max_depth=2)
 
         # Make sure statepoint has been loaded
         if self._sp_filename is None:
@@ -865,38 +872,43 @@ class Library(object):
             xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide])
         elif 'total' in self.mgxs_types:
             mymgxs = self.get_mgxs(domain, 'total')
-            xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide])
+            xsdata.set_total_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
+                                  subdomain=subdomain)
         if 'absorption' in self.mgxs_types:
             mymgxs = self.get_mgxs(domain, 'absorption')
             xsdata.set_absorption_mgxs(mymgxs, xs_type=xs_type,
-                                       nuclide=[nuclide])
+                                       nuclide=[nuclide], subdomain=subdomain)
         if 'fission' in self.mgxs_types:
             mymgxs = self.get_mgxs(domain, 'fission')
             xsdata.set_fission_mgxs(mymgxs, xs_type=xs_type,
-                                    nuclide=[nuclide])
+                                    nuclide=[nuclide], subdomain=subdomain)
         if 'kappa-fission' in self.mgxs_types:
             mymgxs = self.get_mgxs(domain, 'kappa-fission')
             xsdata.set_kappa_fission_mgxs(mymgxs, xs_type=xs_type,
-                                          nuclide=[nuclide])
+                                          nuclide=[nuclide],
+                                          subdomain=subdomain)
         # For chi and nu-fission we can either have only a nu-fission matrix
         # provided, or vectors of chi and nu-fission provided
         if 'nu-fission matrix' in self.mgxs_types:
             mymgxs = self.get_mgxs(domain, 'nu-fission matrix')
             xsdata.set_nu_fission_mgxs(mymgxs, xs_type=xs_type,
-                                       nuclide=[nuclide])
+                                       nuclide=[nuclide], subdomain=subdomain)
         else:
             if 'chi' in self.mgxs_types:
                 mymgxs = self.get_mgxs(domain, 'chi')
-                xsdata.set_chi_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide])
+                xsdata.set_chi_mgxs(mymgxs, xs_type=xs_type, nuclide=[nuclide],
+                                    subdomain=subdomain)
             if 'nu-fission' in self.mgxs_types:
                 mymgxs = self.get_mgxs(domain, 'nu-fission')
                 xsdata.set_nu_fission_mgxs(mymgxs, xs_type=xs_type,
-                                           nuclide=[nuclide])
+                                           nuclide=[nuclide],
+                                           subdomain=subdomain)
         # If multiplicity matrix is available, prefer that
         if 'multiplicity matrix' in self.mgxs_types:
             mymgxs = self.get_mgxs(domain, 'multiplicity matrix')
             xsdata.set_multiplicity_mgxs(mymgxs, xs_type=xs_type,
-                                         nuclide=[nuclide])
+                                         nuclide=[nuclide],
+                                         subdomain=subdomain)
             using_multiplicity = True
         # multiplicity wil fall back to using scatter and nu-scatter
         elif ((('scatter matrix' in self.mgxs_types) and
@@ -904,7 +916,8 @@ class Library(object):
             scatt_mgxs = self.get_mgxs(domain, 'scatter matrix')
             nuscatt_mgxs = self.get_mgxs(domain, 'nu-scatter matrix')
             xsdata.set_multiplicity_mgxs(nuscatt_mgxs, scatt_mgxs,
-                                         xs_type=xs_type, nuclide=[nuclide])
+                                         xs_type=xs_type, nuclide=[nuclide],
+                                         subdomain=subdomain)
             using_multiplicity = True
         else:
             using_multiplicity = False
@@ -912,12 +925,12 @@ class Library(object):
         if using_multiplicity:
             nuscatt_mgxs = self.get_mgxs(domain, 'nu-scatter matrix')
             xsdata.set_scatter_mgxs(nuscatt_mgxs, xs_type=xs_type,
-                                    nuclide=[nuclide])
+                                    nuclide=[nuclide], subdomain=subdomain)
         else:
             if 'nu-scatter matrix' in self.mgxs_types:
                 nuscatt_mgxs = self.get_mgxs(domain, 'nu-scatter matrix')
                 xsdata.set_scatter_mgxs(nuscatt_mgxs, xs_type=xs_type,
-                                        nuclide=[nuclide])
+                                        nuclide=[nuclide], subdomain=subdomain)
 
                 # Since we are not using multiplicity, then
                 # scattering multiplication (nu-scatter) must be
@@ -988,15 +1001,6 @@ class Library(object):
         cv.check_value('xs_type', xs_type, ['macro', 'micro'])
         if xsdata_names is not None:
             cv.check_iterable_type('xsdata_names', xsdata_names, basestring)
-        if xs_ids is not None:
-            if isinstance(xs_ids, basestring):
-                # If we only have a string lets convert it now to a list
-                # of strings.
-                xs_ids = [xs_ids for i in range(len(self.domains))]
-            else:
-                cv.check_iterable_type('xs_ids', xs_ids, basestring)
-        else:
-            xs_ids = ['1m' for i in range(len(self.domains))]
 
         # If gathering material-specific data, set the xs_type to macro
         if not self.by_nuclide:
@@ -1005,32 +1009,79 @@ class Library(object):
         # Initialize file
         mgxs_file = openmc.MGXSLibrary(self.energy_groups)
 
-        # Create the xsdata object and add it to the mgxs_file
-        for i, domain in enumerate(self.domains):
-            if self.by_nuclide:
-                nuclides = list(domain.get_all_nuclides().keys())
+        if xs_ids is not None:
+            if isinstance(xs_ids, basestring):
+                # If we only have a string lets convert it now to a list
+                # of strings.
+                xs_id_input = xs_ids
+                xs_ids = []
+                for domain in self.domains:
+                    if self.domain_type is 'mesh':
+                        for subdomain in domain.cell_generator():
+                            xs_ids.append(xs_id_input)
+                    else:
+                        xs_ids.append(xs_id_input)
             else:
-                nuclides = ['total']
-            for nuclide in nuclides:
-                # Build & add metadata to XSdata object
-                if xsdata_names is None:
-                    xsdata_name = 'set' + str(i + 1)
+                cv.check_iterable_type('xs_ids', xs_ids, basestring)
+        else:
+            xs_ids = []
+            for domain in self.domains:
+                if self.domain_type is 'mesh':
+                    for subdomain in domain.cell_generator():
+                        xs_ids.append('1m')
                 else:
-                    xsdata_name = xsdata_names[i]
-                if nuclide is not 'total':
-                    xsdata_name += '_' + nuclide
+                    xs_ids.append('1m')
 
-                xsdata = self.get_xsdata(domain, xsdata_name, nuclide=nuclide,
-                                         xs_type=xs_type, xs_id=xs_ids[i],
-                                         tabular_legendre=tabular_legendre,
-                                         tabular_points=tabular_points)
+        if self.domain_type == 'mesh':
+            # Create the xsdata objects and add to the mgxs_file
+            i = 0
+            for domain in self.domains:
+                for subdomain in domain.cell_generator():
+                    # Build & add metadata to XSdata object
+                    if xsdata_names is None:
+                        xsdata_name = 'set' + str(i + 1)
+                    else:
+                        xsdata_name = xsdata_names[i]
 
-                mgxs_file.add_xsdata(xsdata)
+                    # Create XSdata and Macroscopic for this domain
+                    xsdata = self.get_xsdata(domain, xsdata_name,
+                                             nuclide='total',
+                                             xs_type=xs_type, xs_id=xs_ids[i],
+                                             tabular_legendre=tabular_legendre,
+                                             tabular_points=tabular_points,
+                                             subdomain=subdomain)
+                    mgxs_file.add_xsdata(xsdata)
+                    i += 1
+
+        else:
+            # Create the xsdata object and add it to the mgxs_file
+            for i, domain in enumerate(self.domains):
+                if self.by_nuclide:
+                    nuclides = list(domain.get_all_nuclides().keys())
+                else:
+                    nuclides = ['total']
+                for nuclide in nuclides:
+                    # Build & add metadata to XSdata object
+                    if xsdata_names is None:
+                        xsdata_name = 'set' + str(i + 1)
+                    else:
+                        xsdata_name = xsdata_names[i]
+                    if nuclide is not 'total':
+                        xsdata_name += '_' + nuclide
+
+                    xsdata = self.get_xsdata(domain, xsdata_name,
+                                             nuclide=nuclide,
+                                             xs_type=xs_type, xs_id=xs_ids[i],
+                                             tabular_legendre=tabular_legendre,
+                                             tabular_points=tabular_points)
+
+                    mgxs_file.add_xsdata(xsdata)
 
         return mgxs_file
 
     def create_mg_mode(self, xsdata_names=None, xs_ids=None,
-                       tabular_legendre=None, tabular_points=33):
+                       tabular_legendre=None, tabular_points=33,
+                       bc=['reflective'] * 6):
         """Creates an openmc.MGXSLibrary object to contain the MGXS data for the
         Multi-Group mode of OpenMC as well as the associated openmc.Materials
         and openmc.Geometry objects. The created Geometry is the same as that
@@ -1060,6 +1111,9 @@ class Library(object):
             parameter is set to `True`.  In this case, this parameter sets the
             number of equally-spaced points in the domain of [-1,1] to be used
             in building the tabular distribution. Default is `33`.
+        bc : iterable of {'reflective', 'periodic', or 'vacuum'}
+            Boundary conditions for each of the six faces of a parallelopiped:
+            North, East, South, Wwest, Up, and Down.
 
         Returns
         -------
@@ -1092,60 +1146,121 @@ class Library(object):
 
         if xsdata_names is not None:
             cv.check_iterable_type('xsdata_names', xsdata_names, basestring)
-        if xs_ids is not None:
-            if isinstance(xs_ids, basestring):
-                # If we only have a string lets convert it now to a list
-                # of strings.
-                xs_ids = [xs_ids for i in range(len(self.domains))]
-            else:
-                cv.check_iterable_type('xs_ids', xs_ids, basestring)
-        else:
-            xs_ids = ['1m' for i in range(len(self.domains))]
         xs_type = 'macro'
 
         # Initialize MGXS File
         mgxs_file = openmc.MGXSLibrary(self.energy_groups)
 
-        # Create a copy of the Geometry to differentiate for these Macroscopics
-        geometry = copy.deepcopy(self.openmc_geometry)
-        materials = openmc.Materials()
-
-        # Get all Cells from the Geometry for differentiation
-        all_cells = geometry.get_all_material_cells()
-
-        # Create the xsdata object and add it to the mgxs_file
-        for i, domain in enumerate(self.domains):
-
-            # Build & add metadata to XSdata object
-            if xsdata_names is None:
-                xsdata_name = 'set' + str(i + 1)
+        # Set ID names
+        if xs_ids is not None:
+            if isinstance(xs_ids, basestring):
+                # If we only have a string lets convert it now to a list
+                # of strings.
+                xs_id_input = xs_ids
+                xs_ids = []
+                for domain in self.domains:
+                    if self.domain_type is 'mesh':
+                        for subdomain in domain.cell_generator():
+                            xs_ids.append(xs_id_input)
+                    else:
+                        xs_ids.append(xs_id_input)
             else:
-                xsdata_name = xsdata_names[i]
+                cv.check_iterable_type('xs_ids', xs_ids, basestring)
+        else:
+            xs_ids = []
+            for domain in self.domains:
+                if self.domain_type is 'mesh':
+                    for subdomain in domain.cell_generator():
+                        xs_ids.append('1m')
+                else:
+                    xs_ids.append('1m')
 
-            # Create XSdata and Macroscopic for this domain
-            xsdata = self.get_xsdata(domain, xsdata_name, nuclide='total',
-                                     xs_type=xs_type, xs_id=xs_ids[i],
-                                     tabular_legendre=tabular_legendre,
-                                     tabular_points=tabular_points)
-            mgxs_file.add_xsdata(xsdata)
-            macroscopic = openmc.Macroscopic(name=xsdata_name, xs=xs_ids[i])
+        if self.domain_type == 'mesh':
+            # Before continuing, we cant build a model automatically if
+            # the user provided multiple mesh domains for library generation.
+            # So check and complain as needed.
+            cv.check_length("domains", self.domains, 1, 1)
 
-            # Create Material and add to collection
-            material = openmc.Material(name=xsdata_name + '.' + xs_ids[i])
-            material.add_macroscopic(macroscopic)
-            materials.append(material)
+            # Create the xsdata objects and add to the mgxs_file
+            # and assign to materials and geometry as we go
+            materials = openmc.Materials()
+            root = openmc.Universe(name='root', universe_id=0)
+            # Add cells representative of the mesh with reflective BC
+            root_cell, cells = \
+                self.domains[0].build_cells(bc)
+            root.add_cell(root_cell)
+            geometry = openmc.Geometry()
+            geometry.root_universe = root
 
-            # Differentiate Geometry with new Material
-            if self.domain_type == 'material':
-                # Fill all appropriate Cells with new Material
-                for cell in all_cells:
-                    if cell.fill.id == domain.id:
-                        cell.fill = material
+            for i, subdomain in enumerate(self.domains[0].cell_generator()):
+                # Build & add metadata to XSdata object
+                if xsdata_names is None:
+                    xsdata_name = 'set' + str(i + 1)
+                else:
+                    xsdata_name = xsdata_names[i]
 
-            elif self.domain_type == 'cell':
-                for cell in all_cells:
-                    if cell.id == domain.id:
-                        cell.fill = material
+                # Create XSdata and Macroscopic for this domain
+                xsdata = self.get_xsdata(domain, xsdata_name,
+                                         nuclide='total',
+                                         xs_type=xs_type, xs_id=xs_ids[i],
+                                         tabular_legendre=tabular_legendre,
+                                         tabular_points=tabular_points,
+                                         subdomain=subdomain)
+                mgxs_file.add_xsdata(xsdata)
+                macroscopic = openmc.Macroscopic(name=xsdata_name,
+                                                 xs=xs_ids[i])
+
+                # Create Material and add to collection
+                material = openmc.Material(name=xsdata_name + '.' +
+                                           xs_ids[i])
+                material.add_macroscopic(macroscopic)
+                materials.append(material)
+
+                # Set the materials for each of the universes
+                cells[i].fill = materials[i]
+
+        else:
+            # Create a copy of the Geometry for these Macroscopics
+            geometry = copy.deepcopy(self.openmc_geometry)
+            materials = openmc.Materials()
+
+            # Get all Cells from the Geometry for differentiation
+            all_cells = geometry.get_all_material_cells()
+
+            # Create the xsdata object and add it to the mgxs_file
+            for i, domain in enumerate(self.domains):
+
+                # Build & add metadata to XSdata object
+                if xsdata_names is None:
+                    xsdata_name = 'set' + str(i + 1)
+                else:
+                    xsdata_name = xsdata_names[i]
+
+                # Create XSdata and Macroscopic for this domain
+                xsdata = self.get_xsdata(domain, xsdata_name, nuclide='total',
+                                         xs_type=xs_type, xs_id=xs_ids[i],
+                                         tabular_legendre=tabular_legendre,
+                                         tabular_points=tabular_points)
+                mgxs_file.add_xsdata(xsdata)
+                macroscopic = openmc.Macroscopic(name=xsdata_name,
+                                                 xs=xs_ids[i])
+
+                # Create Material and add to collection
+                material = openmc.Material(name=xsdata_name + '.' + xs_ids[i])
+                material.add_macroscopic(macroscopic)
+                materials.append(material)
+
+                # Differentiate Geometry with new Material
+                if self.domain_type == 'material':
+                    # Fill all appropriate Cells with new Material
+                    for cell in all_cells:
+                        if cell.fill.id == domain.id:
+                            cell.fill = material
+
+                elif self.domain_type == 'cell':
+                    for cell in all_cells:
+                        if cell.id == domain.id:
+                            cell.fill = material
 
         return mgxs_file, materials, geometry
 
diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index e1c2e62217..608f232b74 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -781,7 +781,8 @@ class MGXS(object):
 
         # Construct a collection of the domain filter bins
         if not isinstance(subdomains, basestring):
-            cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=3)
+            cv.check_iterable_type('subdomains', subdomains, Integral,
+                                   max_depth=3)
             for subdomain in subdomains:
                 filters.append(self.domain_type)
                 filter_bins.append((subdomain,))
diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py
index 6f8188b101..fca5f437cf 100644
--- a/openmc/mgxs_library.py
+++ b/openmc/mgxs_library.py
@@ -603,7 +603,8 @@ class XSdata(object):
         if np.sum(self._nu_fission) > 0.0:
             self._fissionable = True
 
-    def set_total_mgxs(self, total, nuclide='total', xs_type='macro'):
+    def set_total_mgxs(self, total, nuclide='total', xs_type='macro',
+                       subdomain=None):
         """This method allows for an openmc.mgxs.TotalXS or
         openmc.mgxs.TransportXS to be used to set the total cross section
         for this XSdata object.
@@ -619,6 +620,9 @@ class XSdata(object):
         xs_type: {'macro', 'micro'}
             Provide the macro or micro cross section in units of cm^-1 or
             barns. Defaults to 'macro'.
+        subdomain : iterable of int
+            If the MGXS contains a mesh domain type, the subdomain parameter
+            specifies which mesh cell ([i, j, k]) to use.
 
         See also
         --------
@@ -631,16 +635,21 @@ class XSdata(object):
                                     openmc.mgxs.TransportXS))
         check_value('energy_groups', total.energy_groups, [self.energy_groups])
         check_value('domain_type', total.domain_type,
-                    ['universe', 'cell', 'material'])
+                    ['universe', 'cell', 'material', 'mesh'])
+        if subdomain is None:
+            subdomain_val = None
+        else:
+            subdomain_val = [subdomain]
 
         if self.representation is 'isotropic':
-            self._total = total.get_xs(nuclides=nuclide, xs_type=xs_type)
+            self._total = total.get_xs(nuclides=nuclide, xs_type=xs_type,
+                                       subdomains=subdomain_val)
         elif self.representation is 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
     def set_absorption_mgxs(self, absorption, nuclide='total',
-                            xs_type='macro'):
+                            xs_type='macro', subdomain=None):
         """This method allows for an openmc.mgxs.AbsorptionXS
         to be used to set the absorption cross section for this XSdata object.
 
@@ -655,6 +664,9 @@ class XSdata(object):
         xs_type: {'macro', 'micro'}
             Provide the macro or micro cross section in units of cm^-1 or
             barns. Defaults to 'macro'.
+        subdomain : iterable of int
+            If the MGXS contains a mesh domain type, the subdomain parameter
+            specifies which mesh cell ([i, j, k]) to use.
 
         See also
         --------
@@ -667,16 +679,22 @@ class XSdata(object):
         check_value('energy_groups', absorption.energy_groups,
                     [self.energy_groups])
         check_value('domain_type', absorption.domain_type,
-                    ['universe', 'cell', 'material'])
+                    ['universe', 'cell', 'material', 'mesh'])
+        if subdomain is None:
+            subdomain_val = None
+        else:
+            subdomain_val = [subdomain]
 
         if self.representation is 'isotropic':
             self._absorption = absorption.get_xs(nuclides=nuclide,
-                                                 xs_type=xs_type)
+                                                 xs_type=xs_type,
+                                                 subdomains=subdomain_val)
         elif self.representation is 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
-    def set_fission_mgxs(self, fission, nuclide='total', xs_type='macro'):
+    def set_fission_mgxs(self, fission, nuclide='total', xs_type='macro',
+                         subdomain=None):
         """This method allows for an openmc.mgxs.FissionXS
         to be used to set the fission cross section for this XSdata object.
 
@@ -691,6 +709,9 @@ class XSdata(object):
         xs_type: {'macro', 'micro'}
             Provide the macro or micro cross section in units of cm^-1 or
             barns. Defaults to 'macro'.
+        subdomain : iterable of int
+            If the MGXS contains a mesh domain type, the subdomain parameter
+            specifies which mesh cell ([i, j, k]) to use.
 
         See also
         --------
@@ -703,17 +724,22 @@ class XSdata(object):
         check_value('energy_groups', fission.energy_groups,
                     [self.energy_groups])
         check_value('domain_type', fission.domain_type,
-                    ['universe', 'cell', 'material'])
+                    ['universe', 'cell', 'material', 'mesh'])
+        if subdomain is None:
+            subdomain_val = None
+        else:
+            subdomain_val = [subdomain]
 
         if self.representation is 'isotropic':
             self._fission = fission.get_xs(nuclides=nuclide,
-                                           xs_type=xs_type)
+                                           xs_type=xs_type,
+                                           subdomains=subdomain_val)
         elif self.representation is 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
     def set_nu_fission_mgxs(self, nu_fission, nuclide='total',
-                            xs_type='macro'):
+                            xs_type='macro', subdomain=None):
         """This method allows for an openmc.mgxs.NuFissionXS
         to be used to set the nu-fission cross section for this XSdata object.
 
@@ -728,6 +754,9 @@ class XSdata(object):
         xs_type: {'macro', 'micro'}
             Provide the macro or micro cross section in units of cm^-1 or
             barns. Defaults to 'macro'.
+        subdomain : iterable of int
+            If the MGXS contains a mesh domain type, the subdomain parameter
+            specifies which mesh cell ([i, j, k]) to use.
 
         See also
         --------
@@ -741,11 +770,16 @@ class XSdata(object):
         check_value('energy_groups', nu_fission.energy_groups,
                     [self.energy_groups])
         check_value('domain_type', nu_fission.domain_type,
-                    ['universe', 'cell', 'material'])
+                    ['universe', 'cell', 'material', 'mesh'])
+        if subdomain is None:
+            subdomain_val = None
+        else:
+            subdomain_val = [subdomain]
 
         if self.representation is 'isotropic':
             self._nu_fission = nu_fission.get_xs(nuclides=nuclide,
-                                                 xs_type=xs_type)
+                                                 xs_type=xs_type,
+                                                 subdomains=subdomain_val)
         elif self.representation is 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
@@ -759,7 +793,7 @@ class XSdata(object):
             self._fissionable = True
 
     def set_kappa_fission_mgxs(self, k_fission, nuclide='total',
-                               xs_type='macro'):
+                               xs_type='macro', subdomain=None):
         """This method allows for an openmc.mgxs.KappaFissionXS
         to be used to set the kappa-fission cross section for this XSdata
         object.
@@ -775,6 +809,9 @@ class XSdata(object):
         xs_type: {'macro', 'micro'}
             Provide the macro or micro cross section in units of cm^-1 or
             barns. Defaults to 'macro'.
+        subdomain : iterable of int
+            If the MGXS contains a mesh domain type, the subdomain parameter
+            specifies which mesh cell ([i, j, k]) to use.
 
         See also
         --------
@@ -787,16 +824,22 @@ class XSdata(object):
         check_value('energy_groups', k_fission.energy_groups,
                     [self.energy_groups])
         check_value('domain_type', k_fission.domain_type,
-                    ['universe', 'cell', 'material'])
+                    ['universe', 'cell', 'material', 'mesh'])
+        if subdomain is None:
+            subdomain_val = None
+        else:
+            subdomain_val = [subdomain]
 
         if self.representation is 'isotropic':
             self._kappa_fission = k_fission.get_xs(nuclides=nuclide,
-                                                   xs_type=xs_type)
+                                                   xs_type=xs_type,
+                                                   subdomains=subdomain_val)
         elif self.representation is 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
-    def set_chi_mgxs(self, chi, nuclide='total', xs_type='macro'):
+    def set_chi_mgxs(self, chi, nuclide='total', xs_type='macro',
+                     subdomain=None):
         """This method allows for an openmc.mgxs.Chi
         to be used to set chi for this XSdata object.
 
@@ -810,6 +853,9 @@ class XSdata(object):
         xs_type: {'macro', 'micro'}
             Provide the macro or micro cross section in units of cm^-1 or
             barns. Defaults to 'macro'.
+        subdomain : iterable of int
+            If the MGXS contains a mesh domain type, the subdomain parameter
+            specifies which mesh cell ([i, j, k]) to use.
 
         See also
         --------
@@ -827,11 +873,15 @@ class XSdata(object):
         check_type('chi', chi, openmc.mgxs.Chi)
         check_value('energy_groups', chi.energy_groups, [self.energy_groups])
         check_value('domain_type', chi.domain_type,
-                    ['universe', 'cell', 'material'])
+                    ['universe', 'cell', 'material', 'mesh'])
+        if subdomain is None:
+            subdomain_val = None
+        else:
+            subdomain_val = [subdomain]
 
         if self.representation is 'isotropic':
             self._chi = chi.get_xs(nuclides=nuclide,
-                                   xs_type=xs_type)
+                                   xs_type=xs_type, subdomains=subdomain_val)
         elif self.representation is 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
@@ -839,7 +889,8 @@ class XSdata(object):
         if self.use_chi is not None:
             self.use_chi = True
 
-    def set_scatter_mgxs(self, scatter, nuclide='total', xs_type='macro'):
+    def set_scatter_mgxs(self, scatter, nuclide='total', xs_type='macro',
+                         subdomain=None):
         """This method allows for an openmc.mgxs.ScatterMatrixXS
         to be used to set the scatter matrix cross section for this XSdata
         object.  If the XSdata.order attribute has not yet been set, then
@@ -856,6 +907,9 @@ class XSdata(object):
         xs_type: {'macro', 'micro'}
             Provide the macro or micro cross section in units of cm^-1 or
             barns. Defaults to 'macro'.
+        subdomain : iterable of int
+            If the MGXS contains a mesh domain type, the subdomain parameter
+            specifies which mesh cell ([i, j, k]) to use.
 
         See also
         --------
@@ -868,7 +922,7 @@ class XSdata(object):
         check_value('energy_groups', scatter.energy_groups,
                     [self.energy_groups])
         check_value('domain_type', scatter.domain_type,
-                    ['universe', 'cell', 'material'])
+                    ['universe', 'cell', 'material', 'mesh'])
 
         if self.scatt_type != 'legendre':
             msg = 'Anisotropic scattering representations other than ' \
@@ -886,22 +940,27 @@ class XSdata(object):
             check_value('legendre_order', scatter.legendre_order,
                         [self.order])
 
+        if subdomain is None:
+            subdomain_val = None
+        else:
+            subdomain_val = [subdomain]
+
         if self.representation is 'isotropic':
             # Get the scattering orders in the outermost dimension
             self._scatter = np.zeros((self.num_orders,
                                       self.energy_groups.num_groups,
                                       self.energy_groups.num_groups))
             for moment in range(self.num_orders):
-                self._scatter[moment, :, :] = scatter.get_xs(nuclides=nuclide,
-                                                             xs_type=xs_type,
-                                                             moment=moment)
+                self._scatter[moment, :, :] = \
+                    scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
+                                   moment=moment, subdomains=subdomain_val)
 
         elif self.representation is 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'
             raise ValueError(msg)
 
     def set_multiplicity_mgxs(self, nuscatter, scatter=None, nuclide='total',
-                              xs_type='macro'):
+                              xs_type='macro', subdomain=None):
         """This method allows for either the direct use of only an
         openmc.mgxs.MultiplicityMatrixXS OR an openmc.mgxs.NuScatterMatrixXS and
         openmc.mgxs.ScatterMatrixXS to be used to set the scattering
@@ -926,6 +985,9 @@ class XSdata(object):
         xs_type: {'macro', 'micro'}
             Provide the macro or micro cross section in units of cm^-1 or
             barns. Defaults to 'macro'.
+        subdomain : iterable of int
+            If the MGXS contains a mesh domain type, the subdomain parameter
+            specifies which mesh cell ([i, j, k]) to use.
 
         See also
         --------
@@ -939,7 +1001,7 @@ class XSdata(object):
         check_value('energy_groups', nuscatter.energy_groups,
                     [self.energy_groups])
         check_value('domain_type', nuscatter.domain_type,
-                    ['universe', 'cell', 'material'])
+                    ['universe', 'cell', 'material', 'mesh'])
         if scatter is not None:
             check_type('scatter', scatter, openmc.mgxs.ScatterMatrixXS)
             if isinstance(nuscatter, openmc.mgxs.MultiplicityMatrixXS):
@@ -950,16 +1012,22 @@ class XSdata(object):
             check_value('energy_groups', scatter.energy_groups,
                         [self.energy_groups])
             check_value('domain_type', scatter.domain_type,
-                        ['universe', 'cell', 'material'])
+                        ['universe', 'cell', 'material', 'mesh'])
+        if subdomain is None:
+            subdomain_val = None
+        else:
+            subdomain_val = [subdomain]
 
         if self.representation is 'isotropic':
             nuscatt = nuscatter.get_xs(nuclides=nuclide,
-                                       xs_type=xs_type, moment=0)
+                                       xs_type=xs_type, moment=0,
+                                       subdomains=subdomain_val)
             if isinstance(nuscatter, openmc.mgxs.MultiplicityMatrixXS):
                 self._multiplicity = nuscatt
             else:
                 scatt = scatter.get_xs(nuclides=nuclide,
-                                       xs_type=xs_type, moment=0)
+                                       xs_type=xs_type, moment=0,
+                                       subdomains=subdomain_val)
                 self._multiplicity = np.divide(nuscatt, scatt)
         elif self.representation is 'angle':
             msg = 'Angular-Dependent MGXS have not yet been implemented'