diff --git a/openmc/_xml.py b/openmc/_xml.py
new file mode 100644
index 000000000..9c6c219e2
--- /dev/null
+++ b/openmc/_xml.py
@@ -0,0 +1,51 @@
+def clean_indentation(element, level=0, spaces_per_level=2):
+    """
+    copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint
+    it basically walks your tree and adds spaces and newlines so the tree is
+    printed in a nice way
+    """
+
+    i = "\n" + level*spaces_per_level*" "
+
+    if len(element):
+
+        if not element.text or not element.text.strip():
+            element.text = i + spaces_per_level*" "
+
+        if not element.tail or not element.tail.strip():
+            element.tail = i
+
+        for sub_element in element:
+            clean_indentation(sub_element, level+1, spaces_per_level)
+
+        if not sub_element.tail or not sub_element.tail.strip():
+            sub_element.tail = i
+
+    else:
+        if level and (not element.tail or not element.tail.strip()):
+            element.tail = i
+
+
+def get_text(elem, name, default=None):
+    """Retrieve text of an attribute or subelement.
+
+    Parameters
+    ----------
+    elem : xml.etree.ElementTree.Element
+        Element from which to search
+    name : str
+        Name of attribute/subelement
+    default : object
+        A defult value to return if matching attribute/subelement exists
+
+    Returns
+    -------
+    str
+        Text of attribute or subelement
+
+    """
+    if name in elem.attrib:
+        return elem.get(name, default)
+    else:
+        child = elem.find(name)
+        return child.text if child is not None else default
diff --git a/openmc/cell.py b/openmc/cell.py
index cbfd74b21..492bdeeb6 100644
--- a/openmc/cell.py
+++ b/openmc/cell.py
@@ -13,6 +13,7 @@ import openmc
 import openmc.checkvalue as cv
 from openmc.surface import Halfspace
 from openmc.region import Region, Intersection, Complement
+from openmc._xml import get_text
 from .mixin import IDManagerMixin
 
 
@@ -522,3 +523,61 @@ class Cell(IDManagerMixin):
             element.set("rotation", ' '.join(map(str, self.rotation)))
 
         return element
+
+    @classmethod
+    def from_xml_element(cls, elem, surfaces, materials, get_universe):
+        """Generate cell from XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            `<cell>` element
+        surfaces : dict
+            Dictionary mapping surface IDs to :class:`openmc.Surface` instances
+        materials : dict
+            Dictionary mapping material IDs to :class:`openmc.Material`
+            instances (defined in :math:`openmc.Geometry.from_xml`)
+        get_universe : function
+            Function returning universe (defined in
+            :meth:`openmc.Geometry.from_xml`)
+
+        Returns
+        -------
+        Cell
+            Cell instance
+
+        """
+        cell_id = int(get_text(elem, 'id'))
+        name = get_text(elem, 'name')
+        c = cls(cell_id, name)
+
+        # Assign material/distributed materials or fill
+        mat_text = get_text(elem, 'material')
+        if mat_text is not None:
+            mat_ids = mat_text.split()
+            if len(mat_ids) > 1:
+                c.fill = [materials[i] for i in mat_ids]
+            else:
+                c.fill = materials[mat_ids[0]]
+        else:
+            fill_id = int(get_text(elem, 'fill'))
+            c.fill = get_universe(fill_id)
+
+        # Assign region
+        region = get_text(elem, 'region')
+        if region is not None:
+            c.region = Region.from_expression(region, surfaces)
+
+        # Check for other attributes
+        t = get_text(elem, 'temperature')
+        if t is not None:
+            c.temperature = float(t)
+        for key in ('temperature', 'rotation', 'translation'):
+            value = get_text(elem, key)
+            if value is not None:
+                setattr(c, key, [float(x) for x in value.split()])
+
+        # Add this cell to appropriate universe
+        univ_id = int(get_text(elem, 'universe', 0))
+        get_universe(univ_id).add_cell(c)
+        return c
diff --git a/openmc/clean_xml.py b/openmc/clean_xml.py
deleted file mode 100644
index d1002070b..000000000
--- a/openmc/clean_xml.py
+++ /dev/null
@@ -1,26 +0,0 @@
-def clean_xml_indentation(element, level=0, spaces_per_level=2):
-    """
-    copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint
-    it basically walks your tree and adds spaces and newlines so the tree is
-    printed in a nice way
-    """
-
-    i = "\n" + level*spaces_per_level*" "
-
-    if len(element):
-
-        if not element.text or not element.text.strip():
-            element.text = i + spaces_per_level*" "
-
-        if not element.tail or not element.tail.strip():
-            element.tail = i
-
-        for sub_element in element:
-            clean_xml_indentation(sub_element, level+1, spaces_per_level)
-
-        if not sub_element.tail or not sub_element.tail.strip():
-            sub_element.tail = i
-
-    else:
-        if level and (not element.tail or not element.tail.strip()):
-            element.tail = i
diff --git a/openmc/cmfd.py b/openmc/cmfd.py
index 5444334b2..c156cd00a 100644
--- a/openmc/cmfd.py
+++ b/openmc/cmfd.py
@@ -15,7 +15,7 @@ from numbers import Real, Integral
 from xml.etree import ElementTree as ET
 import sys
 
-from openmc.clean_xml import clean_xml_indentation
+from openmc._xml import clean_indentation
 from openmc.checkvalue import (check_type, check_length, check_value,
                                check_greater_than, check_less_than)
 
@@ -512,7 +512,7 @@ class CMFD(object):
         self._create_write_matrices_subelement()
 
         # Clean the indentation in the file to be user-readable
-        clean_xml_indentation(self._cmfd_file)
+        clean_indentation(self._cmfd_file)
 
         # Write the XML Tree to the cmfd.xml file
         tree = ET.ElementTree(self._cmfd_file)
diff --git a/openmc/data/library.py b/openmc/data/library.py
index 58e5e4e16..8d821e49f 100644
--- a/openmc/data/library.py
+++ b/openmc/data/library.py
@@ -4,7 +4,7 @@ import xml.etree.ElementTree as ET
 import h5py
 
 from openmc.mixin import EqualityMixin
-from openmc.clean_xml import clean_xml_indentation
+from openmc._xml import clean_indentation
 from openmc.checkvalue import check_type
 
 
@@ -95,7 +95,7 @@ class DataLibrary(EqualityMixin):
             lib_element.set('type', library['type'])
 
         # Clean the indentation to be user-readable
-        clean_xml_indentation(root)
+        clean_indentation(root)
 
         # Write XML file
         tree = ET.ElementTree(root)
diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py
index 1826ca9ca..4635d99d1 100644
--- a/openmc/deplete/chain.py
+++ b/openmc/deplete/chain.py
@@ -23,7 +23,7 @@ except ImportError:
 import scipy.sparse as sp
 
 import openmc.data
-from openmc.clean_xml import clean_xml_indentation
+from openmc._xml import clean_indentation
 from .nuclide import Nuclide, DecayTuple, ReactionTuple
 
 
@@ -356,7 +356,7 @@ class Chain(object):
         if _have_lxml:
             tree.write(str(filename), encoding='utf-8', pretty_print=True)
         else:
-            clean_xml_indentation(root_elem)
+            clean_indentation(root_elem)
             tree.write(str(filename), encoding='utf-8')
 
     def form_matrix(self, rates):
diff --git a/openmc/geometry.py b/openmc/geometry.py
index 1ee93dfd6..4cbc60c24 100644
--- a/openmc/geometry.py
+++ b/openmc/geometry.py
@@ -1,10 +1,13 @@
-from collections import OrderedDict
+from collections import OrderedDict, defaultdict
 from collections.abc import Iterable
 from copy import deepcopy
+from pathlib import Path
 from xml.etree import ElementTree as ET
 
+import numpy as np
+
 import openmc
-from openmc.clean_xml import clean_xml_indentation
+import openmc._xml as xml
 from openmc.checkvalue import check_type
 
 
@@ -92,12 +95,104 @@ class Geometry(object):
             x.tag, int(x.get('id'))))
 
         # Clean the indentation in the file to be user-readable
-        clean_xml_indentation(root_element)
+        xml.clean_indentation(root_element)
 
         # Write the XML Tree to the geometry.xml file
         tree = ET.ElementTree(root_element)
         tree.write(path, xml_declaration=True, encoding='utf-8')
 
+    @classmethod
+    def from_xml(cls, path='geometry.xml', materials=None):
+        """Generate geometry from XML file
+
+        Parameters
+        ----------
+        path : str, optional
+            Path to geometry XML file
+        materials : openmc.Materials or None
+            Materials used to assign to cells. If None, an attempt is made to
+            generate it from the materials.xml file.
+
+        Returns
+        -------
+        openmc.Geometry
+            Geometry object
+
+        """
+        # Helper function for keeping a cache of Universe instances
+        universes = {}
+        def get_universe(univ_id):
+            if univ_id not in universes:
+                univ = openmc.Universe(univ_id)
+                universes[univ_id] = univ
+            return universes[univ_id]
+
+        tree = ET.parse(path)
+        root = tree.getroot()
+
+        # Get surfaces
+        surfaces = {}
+        periodic = {}
+        for surface in root.findall('surface'):
+            s = openmc.Surface.from_xml_element(surface)
+            surfaces[s.id] = s
+
+            # Check for periodic surface
+            other_id = xml.get_text(surface, 'periodic_surface_id')
+            if other_id is not None:
+                periodic[s.id] = int(other_id)
+
+        # Apply periodic surfaces
+        for s1, s2 in periodic.items():
+            surfaces[s1].periodic_surface = surfaces[s2]
+
+        # Dictionary that maps each universe to a list of cells/lattices that
+        # contain it (needed to determine which universe is the root)
+        child_of = defaultdict(list)
+
+        for elem in root.findall('lattice'):
+            lat = openmc.RectLattice.from_xml_element(elem, get_universe)
+            universes[lat.id] = lat
+            if lat.outer is not None:
+                child_of[lat.outer].append(lat)
+            for u in lat.universes.ravel():
+                child_of[u].append(lat)
+
+        for elem in root.findall('hex_lattice'):
+            lat = openmc.HexLattice.from_xml_element(elem, get_universe)
+            universes[lat.id] = lat
+            if lat.outer is not None:
+                child_of[lat.outer].append(lat)
+            if lat.ndim == 2:
+                for ring in lat.universes:
+                    for u in ring:
+                        child_of[u].append(lat)
+            else:
+                for axial_slice in lat.universes:
+                    for ring in axial_slice:
+                        for u in ring:
+                            child_of[u].append(lat)
+
+        # Create dictionary to easily look up materials
+        if materials is None:
+            filename = Path(path).parent / 'materials.xml'
+            materials = openmc.Materials.from_xml(str(filename))
+        mats = {str(m.id): m for m in materials}
+        mats['void'] = None
+
+        for elem in root.findall('cell'):
+            c = openmc.Cell.from_xml_element(elem, surfaces, mats, get_universe)
+            if c.fill_type in ('universe', 'lattice'):
+                child_of[c.fill].append(c)
+
+        # Determine which universe is the root by finding one which is not a
+        # child of any other object
+        for u in universes.values():
+            if not child_of[u]:
+                return cls(u)
+        else:
+            raise ValueError('Error determining root universe.')
+
     def find(self, point):
         """Find cells/universes/lattices which contain a given point
 
diff --git a/openmc/lattice.py b/openmc/lattice.py
index ef33247fc..8bf4e4757 100644
--- a/openmc/lattice.py
+++ b/openmc/lattice.py
@@ -10,6 +10,7 @@ import numpy as np
 
 import openmc.checkvalue as cv
 import openmc
+from openmc._xml import get_text
 from openmc.mixin import IDManagerMixin
 
 
@@ -768,6 +769,42 @@ class RectLattice(Lattice):
         # Append the XML subelement for this Lattice to the XML element
         xml_element.append(lattice_subelement)
 
+    @classmethod
+    def from_xml_element(cls, elem, get_universe):
+        """Generate rectangular lattice from XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            `<lattice>` element
+        get_universe : function
+            Function returning universe (defined in
+            :meth:`openmc.Geometry.from_xml`)
+
+        Returns
+        -------
+        RectLattice
+            Rectangular 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.pitch = [float(i) for i in get_text(elem, 'pitch').split()]
+        outer = get_text(elem, 'outer')
+        if outer is not None:
+            lat.outer = get_universe(int(outer))
+
+        # Get array of universes
+        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()])
+        uarray.shape = shape
+        lat.universes = uarray
+        return lat
+
 
 class HexLattice(Lattice):
     r"""A lattice consisting of hexagonal prisms.
@@ -1207,6 +1244,78 @@ class HexLattice(Lattice):
         # Append the XML subelement for this Lattice to the XML element
         xml_element.append(lattice_subelement)
 
+    @classmethod
+    def from_xml_element(cls, elem, get_universe):
+        """Generate hexagonal lattice from XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            `<hex_lattice>` element
+        get_universe : function
+            Function returning universe (defined in
+            :meth:`openmc.Geometry.from_xml`)
+
+        Returns
+        -------
+        HexLattice
+            Hexagonal lattice
+
+        """
+        lat_id = int(get_text(elem, 'id'))
+        name = get_text(elem, 'name')
+        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()]
+        outer = get_text(elem, 'outer')
+        if outer is not None:
+            lat.outer = get_universe(int(outer))
+
+        # Get nested lists of universes
+        lat._num_rings = n_rings = int(get_text(elem, 'n_rings'))
+        lat._num_axial = n_axial = int(get_text(elem, 'n_axial', 1))
+
+        # 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)]
+        if n_axial > 1:
+            univs = [deepcopy(univs) for i in range(n_axial)]
+
+        # Get flat array of universes numbers
+        uarray = np.array([get_universe(int(i)) for i in
+                            get_text(elem, 'universes').split()])
+
+        # Fill nested lists
+        j = 0
+        for z in range(n_axial):
+            # 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]
+
+                # 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
+                    if not lat.is_valid_index((x, alpha, z)):
+                        # Move to the right
+                        x += 2
+                        alpha -= 1
+                        if not lat.is_valid_index((x, alpha, z)):
+                            # Reached the bottom
+                            break
+                j += 1
+        lat.universes = univs
+        return lat
+
     def _repr_axial_slice(self, universes):
         """Return string representation for the given 2D group of universes.
 
diff --git a/openmc/material.py b/openmc/material.py
index bd31fbd23..d11f9de49 100644
--- a/openmc/material.py
+++ b/openmc/material.py
@@ -9,7 +9,7 @@ import numpy as np
 import openmc
 import openmc.data
 import openmc.checkvalue as cv
-from openmc.clean_xml import clean_xml_indentation
+from openmc._xml import clean_indentation
 from .mixin import IDManagerMixin
 
 
@@ -912,6 +912,56 @@ class Material(IDManagerMixin):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate material from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.Material
+            Material generated from XML element
+
+        """
+        mat_id = int(elem.get('id'))
+        mat = cls(mat_id)
+        mat.name = elem.get('name')
+        mat.temperature = elem.get('temperature')
+        mat.depletable = bool(elem.get('depletable'))
+
+        # Get each nuclide
+        for nuclide in elem.findall('nuclide'):
+            name = nuclide.attrib['name']
+            if 'ao' in nuclide.attrib:
+                mat.add_nuclide(name, float(nuclide.attrib['ao']))
+            elif 'wo' in nuclide.attrib:
+                mat.add_nuclide(name, float(nuclide.attrib['wo']), 'wo')
+
+        # Get each S(a,b) table
+        for sab in elem.findall('sab'):
+            fraction = float(sab.get('fraction', 1.0))
+            mat.add_s_alpha_beta(sab.get('name'), fraction)
+
+        # Get total material density
+        density = elem.find('density')
+        units = density.get('units')
+        if units == 'sum':
+            mat.set_density(units)
+        else:
+            value = float(density.get('value'))
+            mat.set_density(units, value)
+
+        # Check for isotropic scattering nuclides
+        isotropic = elem.find('isotropic')
+        if isotropic is not None:
+            mat.isotropic = isotropic.text.split()
+
+        return mat
+
 
 class Materials(cv.CheckedList):
     """Collection of Materials used for an OpenMC simulation.
@@ -1031,8 +1081,41 @@ class Materials(cv.CheckedList):
         self._create_material_subelements(root_element)
 
         # Clean the indentation in the file to be user-readable
-        clean_xml_indentation(root_element)
+        clean_indentation(root_element)
 
         # Write the XML Tree to the materials.xml file
         tree = ET.ElementTree(root_element)
         tree.write(path, xml_declaration=True, encoding='utf-8')
+
+    @classmethod
+    def from_xml(cls, path='materials.xml'):
+        """Generate materials collection from XML file
+
+        Parameters
+        ----------
+        path : str, optional
+            Path to materials XML file
+
+        Returns
+        -------
+        openmc.Materials
+            Materials collection
+
+        """
+        tree = ET.parse(path)
+        root = tree.getroot()
+
+        # Generate each material
+        materials = cls()
+        for material in root.findall('material'):
+            materials.append(Material.from_xml_element(material))
+
+        # Check for cross sections settings
+        xs = tree.find('cross_sections')
+        if xs is not None:
+            materials.cross_sections = xs.text
+        mpl = tree.find('multipole_library')
+        if mpl is not None:
+            materials.multipole_library = mpl.text
+
+        return materials
diff --git a/openmc/plots.py b/openmc/plots.py
index 8eff78b1d..f99392def 100644
--- a/openmc/plots.py
+++ b/openmc/plots.py
@@ -9,7 +9,7 @@ import numpy as np
 
 import openmc
 import openmc.checkvalue as cv
-from openmc.clean_xml import clean_xml_indentation
+from openmc._xml import clean_indentation
 from openmc.mixin import IDManagerMixin
 
 
@@ -816,7 +816,7 @@ class Plots(cv.CheckedList):
         self._create_plot_subelements()
 
         # Clean the indentation in the file to be user-readable
-        clean_xml_indentation(self._plots_file)
+        clean_indentation(self._plots_file)
 
         # Write the XML Tree to the plots.xml file
         tree = ET.ElementTree(self._plots_file)
diff --git a/openmc/settings.py b/openmc/settings.py
index 7629ea9af..efe41b2a1 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -6,7 +6,7 @@ import sys
 
 import numpy as np
 
-from openmc.clean_xml import clean_xml_indentation
+from openmc._xml import clean_indentation
 import openmc.checkvalue as cv
 from openmc import VolumeCalculation, Source, Mesh
 
@@ -994,7 +994,7 @@ class Settings(object):
         self._create_log_grid_bins_subelement(root_element)
 
         # Clean the indentation in the file to be user-readable
-        clean_xml_indentation(root_element)
+        clean_indentation(root_element)
 
         # Write the XML Tree to the settings.xml file
         tree = ET.ElementTree(root_element)
diff --git a/openmc/surface.py b/openmc/surface.py
index ad5053e37..826ac5cc8 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -59,7 +59,6 @@ class Surface(IDManagerMixin):
     def __init__(self, surface_id=None, boundary_type='transmission', name=''):
         self.id = surface_id
         self.name = name
-        self._type = ''
         self.boundary_type = boundary_type
 
         # A dictionary of the quadratic surface coefficients
@@ -67,10 +66,6 @@ class Surface(IDManagerMixin):
         # Value    - coefficient value
         self._coefficients = {}
 
-        # An ordered list of the coefficient names to export to XML in the
-        # proper order
-        self._coeff_keys = []
-
     def __neg__(self):
         return Halfspace(self, '-')
 
@@ -203,6 +198,49 @@ class Surface(IDManagerMixin):
 
         return element
 
+    @staticmethod
+    def from_xml_element(elem):
+        """Generate surface from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.Surface
+            Instance of a surface subclass
+
+        """
+
+        # Determine appropriate class
+        surf_type = elem.get('type')
+        surface_classes = {
+            'plane': Plane,
+            'x-plane': XPlane,
+            'y-plane': YPlane,
+            'z-plane': ZPlane,
+            'x-cylinder': XCylinder,
+            'y-cylinder': YCylinder,
+            'z-cylinder': ZCylinder,
+            'sphere': Sphere,
+            'x-cone': XCone,
+            'y-cone': YCone,
+            'z-cone': ZCone,
+            'quadric': Quadric,
+        }
+        cls = surface_classes[surf_type]
+
+        # Determine ID, boundary type, coefficients
+        kwargs = {}
+        kwargs['surface_id'] = int(elem.get('id'))
+        kwargs['boundary_type'] = elem.get('boundary', 'transmission')
+        coeffs = [float(x) for x in elem.get('coeffs').split()]
+        kwargs.update(dict(zip(cls._coeff_keys, coeffs)))
+
+        return cls(**kwargs)
+
     @staticmethod
     def from_hdf5(group):
         """Create surface from HDF5 group
@@ -324,12 +362,12 @@ class Plane(Surface):
 
     """
 
+    _type = 'plane'
+    _coeff_keys = ('A', 'B', 'C', 'D')
+
     def __init__(self, surface_id=None, boundary_type='transmission',
                  A=1., B=0., C=0., D=0., name=''):
         super().__init__(surface_id, boundary_type, name=name)
-
-        self._type = 'plane'
-        self._coeff_keys = ['A', 'B', 'C', 'D']
         self._periodic_surface = None
         self.a = A
         self.b = B
@@ -458,12 +496,12 @@ class XPlane(Plane):
 
     """
 
+    _type = 'x-plane'
+    _coeff_keys = ('x0',)
+
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=0., name=''):
         super().__init__(surface_id, boundary_type, name=name)
-
-        self._type = 'x-plane'
-        self._coeff_keys = ['x0']
         self.x0 = x0
 
     @property
@@ -563,13 +601,13 @@ class YPlane(Plane):
 
     """
 
+    _type = 'y-plane'
+    _coeff_keys = ('y0',)
+
     def __init__(self, surface_id=None, boundary_type='transmission',
                  y0=0., name=''):
         # Initialize YPlane class attributes
         super().__init__(surface_id, boundary_type, name=name)
-
-        self._type = 'y-plane'
-        self._coeff_keys = ['y0']
         self.y0 = y0
 
     @property
@@ -669,13 +707,13 @@ class ZPlane(Plane):
 
     """
 
+    _type = 'z-plane'
+    _coeff_keys = ('z0',)
+
     def __init__(self, surface_id=None, boundary_type='transmission',
                  z0=0., name=''):
         # Initialize ZPlane class attributes
         super().__init__(surface_id, boundary_type, name=name)
-
-        self._type = 'z-plane'
-        self._coeff_keys = ['z0']
         self.z0 = z0
 
     @property
@@ -774,8 +812,6 @@ class Cylinder(Surface, metaclass=ABCMeta):
     def __init__(self, surface_id=None, boundary_type='transmission',
                  R=1., name=''):
         super().__init__(surface_id, boundary_type, name=name)
-
-        self._coeff_keys = ['R']
         self.r = R
 
     @property
@@ -831,12 +867,12 @@ class XCylinder(Cylinder):
 
     """
 
+    _type = 'x-cylinder'
+    _coeff_keys = ('y0', 'z0', 'R')
+
     def __init__(self, surface_id=None, boundary_type='transmission',
                  y0=0., z0=0., R=1., name=''):
         super().__init__(surface_id, boundary_type, R, name=name)
-
-        self._type = 'x-cylinder'
-        self._coeff_keys = ['y0', 'z0', 'R']
         self.y0 = y0
         self.z0 = z0
 
@@ -953,12 +989,12 @@ class YCylinder(Cylinder):
 
     """
 
+    _type = 'y-cylinder'
+    _coeff_keys = ('x0', 'z0', 'R')
+
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=0., z0=0., R=1., name=''):
         super().__init__(surface_id, boundary_type, R, name=name)
-
-        self._type = 'y-cylinder'
-        self._coeff_keys = ['x0', 'z0', 'R']
         self.x0 = x0
         self.z0 = z0
 
@@ -1075,12 +1111,12 @@ class ZCylinder(Cylinder):
 
     """
 
+    _type = 'z-cylinder'
+    _coeff_keys = ('x0', 'y0', 'R')
+
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=0., y0=0., R=1., name=''):
         super().__init__(surface_id, boundary_type, R, name=name)
-
-        self._type = 'z-cylinder'
-        self._coeff_keys = ['x0', 'y0', 'R']
         self.x0 = x0
         self.y0 = y0
 
@@ -1201,12 +1237,12 @@ class Sphere(Surface):
 
     """
 
+    _type = 'sphere'
+    _coeff_keys = ('x0', 'y0', 'z0', 'R')
+
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=0., y0=0., z0=0., R=1., name=''):
         super().__init__(surface_id, boundary_type, name=name)
-
-        self._type = 'sphere'
-        self._coeff_keys = ['x0', 'y0', 'z0', 'R']
         self.x0 = x0
         self.y0 = y0
         self.z0 = z0
@@ -1348,11 +1384,12 @@ class Cone(Surface, metaclass=ABCMeta):
         Type of the surface
 
     """
+
+    _coeff_keys = ('x0', 'y0', 'z0', 'R2')
+
     def __init__(self, surface_id=None, boundary_type='transmission',
                  x0=0., y0=0., z0=0., R2=1., name=''):
         super().__init__(surface_id, boundary_type, name=name)
-
-        self._coeff_keys = ['x0', 'y0', 'z0', 'R2']
         self.x0 = x0
         self.y0 = y0
         self.z0 = z0
@@ -1443,12 +1480,7 @@ class XCone(Cone):
 
     """
 
-    def __init__(self, surface_id=None, boundary_type='transmission',
-                 x0=0., y0=0., z0=0., R2=1., name=''):
-        super().__init__(surface_id, boundary_type, x0, y0,
-                                    z0, R2, name=name)
-
-        self._type = 'x-cone'
+    _type = 'x-cone'
 
     def evaluate(self, point):
         """Evaluate the surface equation at a given point.
@@ -1519,12 +1551,7 @@ class YCone(Cone):
 
     """
 
-    def __init__(self, surface_id=None, boundary_type='transmission',
-                 x0=0., y0=0., z0=0., R2=1., name=''):
-        super().__init__(surface_id, boundary_type, x0, y0, z0,
-                                    R2, name=name)
-
-        self._type = 'y-cone'
+    _type = 'y-cone'
 
     def evaluate(self, point):
         """Evaluate the surface equation at a given point.
@@ -1595,12 +1622,7 @@ class ZCone(Cone):
 
     """
 
-    def __init__(self, surface_id=None, boundary_type='transmission',
-                 x0=0., y0=0., z0=0., R2=1., name=''):
-        super().__init__(surface_id, boundary_type, x0, y0, z0,
-                                    R2, name=name)
-
-        self._type = 'z-cone'
+    _type = 'z-cone'
 
     def evaluate(self, point):
         """Evaluate the surface equation at a given point.
@@ -1659,13 +1681,13 @@ class Quadric(Surface):
 
     """
 
+    _type = 'quadric'
+    _coeff_keys = ('a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k')
+
     def __init__(self, surface_id=None, boundary_type='transmission',
                  a=0., b=0., c=0., d=0., e=0., f=0., g=0.,
                  h=0., j=0., k=0., name=''):
         super().__init__(surface_id, boundary_type, name=name)
-
-        self._type = 'quadric'
-        self._coeff_keys = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k']
         self.a = a
         self.b = b
         self.c = c
diff --git a/openmc/tallies.py b/openmc/tallies.py
index 74f342903..1e4b142fe 100644
--- a/openmc/tallies.py
+++ b/openmc/tallies.py
@@ -15,7 +15,7 @@ import h5py
 
 import openmc
 import openmc.checkvalue as cv
-from openmc.clean_xml import clean_xml_indentation
+from openmc._xml import clean_indentation
 from .mixin import IDManagerMixin
 
 
@@ -3187,7 +3187,7 @@ class Tallies(cv.CheckedList):
         self._create_derivative_subelements(root_element)
 
         # Clean the indentation in the file to be user-readable
-        clean_xml_indentation(root_element)
+        clean_indentation(root_element)
 
         # Write the XML Tree to the tallies.xml file
         tree = ET.ElementTree(root_element)
diff --git a/tests/regression_tests/lattice_hex/geometry.xml b/tests/regression_tests/lattice_hex/geometry.xml
index cfe39849c..fa6a18ee1 100644
--- a/tests/regression_tests/lattice_hex/geometry.xml
+++ b/tests/regression_tests/lattice_hex/geometry.xml
@@ -7,42 +7,42 @@
 <geometry>
 
   <!-- Universe 011: Fuel pellet -->
-  <surface id="011" type="z-cylinder" coeffs="0.0 0.0 0.06"/>
-  <surface id="012" type="z-cylinder" coeffs="0.0 0.0 0.378"/>
-  <surface id="013" type="z-cylinder" coeffs="0.0 0.0 0.3865"/>
-  <surface id="014" type="z-cylinder" coeffs="0.0 0.0 0.455"/>
-  <surface id="015" type="sphere" coeffs="0.0 0.0 -2.5639 2.0"/>
-  <surface id="016" type="sphere" coeffs="0.0 0.0 2.5639 2.0"/>
+  <surface id="11" type="z-cylinder" coeffs="0.0 0.0 0.06"/>
+  <surface id="12" type="z-cylinder" coeffs="0.0 0.0 0.378"/>
+  <surface id="13" type="z-cylinder" coeffs="0.0 0.0 0.3865"/>
+  <surface id="14" type="z-cylinder" coeffs="0.0 0.0 0.455"/>
+  <surface id="15" type="sphere" coeffs="0.0 0.0 -2.5639 2.0"/>
+  <surface id="16" type="sphere" coeffs="0.0 0.0 2.5639 2.0"/>
 
-  <cell id="011" universe="011" material="void" region="-011"/>
-  <cell id="012" universe="011" material="13" region="011 -012 015 016"/>
-  <cell id="013" universe="011" material="void" region="011 -012 -015"/>
-  <cell id="014" universe="011" material="void" region="011 -012 -016"/>
-  <cell id="015" universe="011" material="void" region="012 -013"/>
-  <cell id="016" universe="011" material="21" region="013 -014"/>
-  <cell id="017" universe="011" material="01" region="014"/>
+  <cell id="11" universe="11" material="void" region="-11"/>
+  <cell id="12" universe="11" material="13" region="11 -12 15 16"/>
+  <cell id="13" universe="11" material="void" region="11 -12 -15"/>
+  <cell id="14" universe="11" material="void" region="11 -12 -16"/>
+  <cell id="15" universe="11" material="void" region="12 -13"/>
+  <cell id="16" universe="11" material="21" region="13 -14"/>
+  <cell id="17" universe="11" material="1" region="14"/>
 
 
   <!-- Universe 051: Central guide tube -->
-    <surface id="051" type="z-cylinder" coeffs="0.0 0.0 0.44"/>
-    <surface id="052" type="z-cylinder" coeffs="0.0 0.0 0.515"/>
-    <cell id="051" universe="051" material="01" region="-051"/>
-    <cell id="052" universe="051" material="21" region="051 -052"/>
-    <cell id="053" universe="051" material="01" region="052"/>
+    <surface id="51" type="z-cylinder" coeffs="0.0 0.0 0.44"/>
+    <surface id="52" type="z-cylinder" coeffs="0.0 0.0 0.515"/>
+    <cell id="51" universe="51" material="1" region="-51"/>
+    <cell id="52" universe="51" material="21" region="51 -52"/>
+    <cell id="53" universe="51" material="1" region="52"/>
 
 
   <!-- Universe 055: Cluster tube -->
     <surface id="155" type="z-cylinder" coeffs="0.0 0.0 0.55"/>
     <surface id="156" type="z-cylinder" coeffs="0.0 0.0 0.63"/>
-    <cell id="155" universe="055" material="01" region="-155"/>
-    <cell id="156" universe="055" material="22" region="155 -156"/>
-    <cell id="157" universe="055" material="01" region="156"/>
+    <cell id="155" universe="55" material="1" region="-155"/>
+    <cell id="156" universe="55" material="22" region="155 -156"/>
+    <cell id="157" universe="55" material="1" region="156"/>
 
 
 
 
   <!-- Lattice 101: Fuel assembly pins -->
-  <cell id="100" universe="100" material="01"/>
+  <cell id="100" universe="100" material="1"/>
   <hex_lattice id="101" outer="100" n_rings="11" n_axial="3"
       center="0.0 0.0 0.0" pitch="1.265 1.2">
     <universes>
diff --git a/tests/unit_tests/conftest.py b/tests/unit_tests/conftest.py
index 1434eaf3b..94fa6abc7 100644
--- a/tests/unit_tests/conftest.py
+++ b/tests/unit_tests/conftest.py
@@ -60,3 +60,54 @@ def cell_with_lattice():
     return ([inside_cyl, outside_cyl, main_cell],
             [m_inside[0], m_inside[1], m_inside[3], m_outside],
             univ, lattice)
+
+@pytest.fixture
+def mixed_lattice_model(uo2, water):
+    cyl = openmc.ZCylinder(R=0.4)
+    c1 = openmc.Cell(fill=uo2, region=-cyl)
+    c1.temperature = 600.0
+    c2 = openmc.Cell(fill=water, region=+cyl)
+    pin = openmc.Universe(cells=[c1, c2])
+
+    empty = openmc.Cell()
+    empty_univ = openmc.Universe(cells=[empty])
+
+    hex_lattice = openmc.HexLattice()
+    hex_lattice.center = (0.0, 0.0)
+    hex_lattice.pitch = (1.2, 10.0)
+    outer_ring = [pin]*6
+    inner_ring = [empty_univ]
+    axial_level = [outer_ring, inner_ring]
+    hex_lattice.universes = [axial_level]*3
+    hex_lattice.outer = empty_univ
+
+    cell_hex = openmc.Cell(fill=hex_lattice)
+    u = openmc.Universe(cells=[cell_hex])
+    rotated_cell_hex = openmc.Cell(fill=u)
+    rotated_cell_hex.rotation = (0., 0., 30.)
+    ur = openmc.Universe(cells=[rotated_cell_hex])
+
+    d = 6.0
+    rect_lattice = openmc.RectLattice()
+    rect_lattice.lower_left = (-d, -d)
+    rect_lattice.pitch = (d, d)
+    rect_lattice.outer = empty_univ
+    rect_lattice.universes = [
+        [ur, empty_univ],
+        [empty_univ, u]
+    ]
+
+    xmin = openmc.XPlane(x0=-d, boundary_type='periodic')
+    xmax = openmc.XPlane(x0=d, boundary_type='periodic')
+    xmin.periodic_surface = xmax
+    ymin = openmc.YPlane(y0=-d, boundary_type='periodic')
+    ymax = openmc.YPlane(y0=d, boundary_type='periodic')
+    main_cell = openmc.Cell(fill=rect_lattice,
+                            region=+xmin & -xmax & +ymin & -ymax)
+
+    # Create geometry and use unique material in each fuel cell
+    geometry = openmc.Geometry([main_cell])
+    geometry.determine_paths()
+    c1.fill = [water.clone() for i in range(c1.num_instances)]
+
+    return openmc.model.Model(geometry)
diff --git a/tests/unit_tests/test_geometry.py b/tests/unit_tests/test_geometry.py
index 93d2fa634..b5d6076bb 100644
--- a/tests/unit_tests/test_geometry.py
+++ b/tests/unit_tests/test_geometry.py
@@ -244,3 +244,14 @@ def test_determine_paths(cell_with_lattice):
     for i in range(4):
         assert geom.get_instances(cells[0].paths[i]) == i
         assert geom.get_instances(mats[-1].paths[i]) == i
+
+def test_from_xml(run_in_tmpdir, mixed_lattice_model):
+    # Export model
+    mixed_lattice_model.export_to_xml()
+
+    # Import geometry
+    geom = openmc.Geometry.from_xml()
+    assert isinstance(geom, openmc.Geometry)
+    ll, ur = geom.bounding_box
+    assert ll == pytest.approx((-6.0, -6.0, -np.inf))
+    assert ur == pytest.approx((6.0, 6.0, np.inf))
diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py
index b7b745408..7cedd8dee 100644
--- a/tests/unit_tests/test_material.py
+++ b/tests/unit_tests/test_material.py
@@ -182,3 +182,37 @@ def test_borated_water():
     # Test the density override
     m = openmc.model.borated_water(975, 566.5, 15.51, density=0.9)
     assert m.density == pytest.approx(0.9, 1e-3)
+
+
+def test_from_xml(run_in_tmpdir):
+    # Create a materials.xml file
+    m1 = openmc.Material(1, 'water')
+    m1.add_nuclide('H1', 1.0)
+    m1.add_nuclide('O16', 2.0)
+    m1.add_s_alpha_beta('c_H_in_H2O')
+    m1.set_density('g/cm3', 0.9)
+    m1.isotropic = ['H1']
+    m2 = openmc.Material(2, 'zirc')
+    m2.add_nuclide('Zr90', 1.0, 'wo')
+    m2.set_density('kg/m3', 10.0)
+    m3 = openmc.Material(3)
+    m3.add_nuclide('N14', 0.02)
+
+    mats = openmc.Materials([m1, m2, m3])
+    mats.cross_sections = 'fake_path.xml'
+    mats.multipole_library = 'fake_multipole/'
+    mats.export_to_xml()
+
+    # Regenerate materials from XML
+    mats = openmc.Materials.from_xml()
+    assert len(mats) == 3
+    m1 = mats[0]
+    assert m1.id == 1
+    assert m1.name == 'water'
+    assert m1.nuclides == [('H1', 1.0, 'ao'), ('O16', 2.0, 'ao')]
+    assert m1.isotropic == ['H1']
+    m2 = mats[1]
+    assert m2.nuclides == [('Zr90', 1.0, 'wo')]
+    assert m2.density == 10.0
+    assert m2.density_units == 'kg/m3'
+    assert mats[2].density_units == 'sum'