diff --git a/openmc/mesh.py b/openmc/mesh.py
index ba07ec178b..a30b333fb3 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -246,6 +246,43 @@ class Mesh(IDManagerMixin):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate mesh from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.Mesh
+            Mesh generated from XML element
+
+        """
+        mesh_id = int(elem.get('id'))
+        mesh = cls(mesh_id)
+        mesh.type = elem.get('type')
+
+        dimension = elem.findtext('dimension')
+        if dimension is not None:
+            mesh.dimension = [int(x) for x in dimension.split()]
+
+        lower_left = elem.findtext('lower_left')
+        if lower_left is not None:
+            mesh.lower_left = [float(x) for x in lower_left.split()]
+
+        upper_right = elem.findtext('upper_right')
+        if upper_right is not None:
+            mesh.upper_right = [float(x) for x in upper_right.split()]
+
+        width = elem.findtext('width')
+        if width is not None:
+            mesh.width = [float(x) for x in width.split()]
+
+        return mesh
+
     def build_cells(self, bc=['reflective'] * 6):
         """Generates a lattice of universes with the same dimensionality
         as the mesh object.  The individual cells/universes produced
diff --git a/openmc/settings.py b/openmc/settings.py
index f348bc6cb0..a7ec086792 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -693,29 +693,29 @@ class Settings(object):
         elem = ET.SubElement(root, "run_mode")
         elem.text = self._run_mode
 
-    def _create_batches_subelement(self, run_mode_element):
+    def _create_batches_subelement(self, root):
         if self._batches is not None:
-            element = ET.SubElement(run_mode_element, "batches")
+            element = ET.SubElement(root, "batches")
             element.text = str(self._batches)
 
-    def _create_generations_per_batch_subelement(self, run_mode_element):
+    def _create_generations_per_batch_subelement(self, root):
         if self._generations_per_batch is not None:
-            element = ET.SubElement(run_mode_element, "generations_per_batch")
+            element = ET.SubElement(root, "generations_per_batch")
             element.text = str(self._generations_per_batch)
 
-    def _create_inactive_subelement(self, run_mode_element):
+    def _create_inactive_subelement(self, root):
         if self._inactive is not None:
-            element = ET.SubElement(run_mode_element, "inactive")
+            element = ET.SubElement(root, "inactive")
             element.text = str(self._inactive)
 
-    def _create_particles_subelement(self, run_mode_element):
+    def _create_particles_subelement(self, root):
         if self._particles is not None:
-            element = ET.SubElement(run_mode_element, "particles")
+            element = ET.SubElement(root, "particles")
             element.text = str(self._particles)
 
-    def _create_keff_trigger_subelement(self, run_mode_element):
+    def _create_keff_trigger_subelement(self, root):
         if self._keff_trigger is not None:
-            element = ET.SubElement(run_mode_element, "keff_trigger")
+            element = ET.SubElement(root, "keff_trigger")
 
             for key in self._keff_trigger:
                 subelement = ET.SubElement(element, key)
@@ -985,3 +985,238 @@ class Settings(object):
         # Write the XML Tree to the settings.xml file
         tree = ET.ElementTree(root_element)
         tree.write(str(p), xml_declaration=True, encoding='utf-8')
+
+    @classmethod
+    def from_xml(cls, path='settings.xml'):
+        """Generate settings from XML file
+
+        Parameters
+        ----------
+        path : str, optional
+            Path to settings XML file
+
+        Returns
+        -------
+        openmc.Settings
+            Settings object
+
+        """
+        tree = ET.parse(path)
+        root = tree.getroot()
+
+        settings = cls()
+
+        # Get the run mode
+        elem = root.find('run_mode')
+        if elem is not None:
+            settings.run_mode = elem.text
+
+        # Get number of particles
+        elem = root.find('particles')
+        if elem is not None:
+            settings.particles = int(elem.text)
+
+        # Get number of batches
+        elem = root.find('batches')
+        if elem is not None:
+            settings.batches = int(elem.text)
+
+        # Get number of inactive batches
+        elem = root.find('inactive')
+        if elem is not None:
+            settings.inactive = int(elem.text)
+
+        # Get number of generations per batch
+        elem = root.find('generations_per_batch')
+        if elem is not None:
+            settings.generations_per_batch = int(elem.text)
+
+        # Get keff trigger
+        elem = root.find('keff_trigger')
+        if elem is not None:
+            trigger = elem.findtext('type')
+            threshold = float(elem.findtext('threshold'))
+            settings.keff_trigger = {'type': trigger, 'threshold': threshold}
+
+        # Get the source
+        for elem in root.findall('source'):
+            settings.source.append(Source.from_xml_element(elem))
+
+        # Get the output
+        elem = root.find('output')
+        if elem is not None:
+            settings.output = {}
+            for entry in elem:
+                key = entry.tag
+                if key in ('summary', 'tallies'):
+                    value = entry.text == 'true'
+                else:
+                    value = entry.text
+                settings.output[key] = value
+
+        # Get the statepoint
+        elem = root.find('state_point')
+        if elem is not None:
+            batches = elem.findtext('batches')
+            if batches is not None:
+                settings.statepoint['batches'] = [int(x) for x in batches.split()]
+
+        # Get the sourcepoint
+        elem = root.find('source_point')
+        if elem is not None:
+            for entry in elem:
+                key = entry.tag
+                if key in ('separate', 'write', 'overwrite'):
+                    value = entry.text == 'true'
+                else:
+                    value = [int(x) for x in entry.text.split()]
+                settings.sourcepoint[key] = value
+
+        # Get confidence intervals
+        elem = root.find('confidence_intervals')
+        if elem is not None:
+            settings.confidence_intervals = elem.text == 'true'
+
+        # Get electron treatment
+        elem = root.find('electron_treatment')
+        if elem is not None:
+            settings.electron_treatment = elem.text
+
+        # Get energy mode
+        elem = root.find('energy_mode')
+        if elem is not None:
+            settings.energy_mode = elem.text
+
+        # Get max order
+        elem = root.find('max_order')
+        if elem is not None:
+            settings.max_order = int(elem.text)
+
+        # Get photon transport
+        elem = root.find('photon_transport')
+        if elem is not None:
+            settings.photon_transport = elem.text == 'true'
+
+        # Get probability tables
+        elem = root.find('ptables')
+        if elem is not None:
+            settings.ptables = elem.text == 'true'
+
+        # Get seed
+        elem = root.find('seed')
+        if elem is not None:
+            settings.seed = int(elem.text)
+
+        # Get survival biasing
+        elem = root.find('survival_biasing')
+        if elem is not None:
+            settings.survival_biasing = elem.text == 'true'
+
+        # Get cutoff
+        elem = root.find('cutoff')
+        if elem is not None:
+            settings.cutoff = {x.tag: float(x.text) for x in elem}
+
+        # Get entropy mesh
+        elem = root.find('entropy_mesh')
+        if elem is not None:
+            settings.entropy_mesh = Mesh.from_xml_element(elem)
+
+        # Get trigger
+        elem = root.find('trigger')
+        if elem is not None:
+            active = elem.find('active')
+            settings.trigger_active = active.text == 'true'
+            max_batches = elem.find('max_batches')
+            if max_batches is not None:
+                settings.trigger_max_batches = int(max_batches.text)
+            batch_interval = elem.find('batch_interval')
+            if batch_interval is not None:
+                settings.trigger_batch_interval = int(batch_interval.text)
+
+        # Get no reduce
+        elem = root.find('no_reduce')
+        if elem is not None:
+            settings.no_reduce = elem.text == 'true'
+
+        # Get verbosity
+        elem = root.find('verbosity')
+        if elem is not None:
+            settings.verbosity = int(elem.text)
+
+        # Get tabular legendre
+        elem = root.find('tabular_legendre')
+        if elem is not None:
+            enable = elem.findtext('eneable')
+            settings.tabular_legendre['enable'] = enable == 'true'
+            num_points = elem.findtext('num_points')
+            if num_points is not None:
+                settings.tabular_legendre['num_points'] = int(num_points)
+
+        # Get temperature
+        elem = root.findtext('temperature_default')
+        if elem is not None:
+            settings.temperature['default'] = float(elem)
+        elem = root.findtext('temperature_tolerance')
+        if elem is not None:
+            settings.temperature['tolerance'] = float(elem)
+        elem = root.findtext('temperature_method')
+        if elem is not None:
+            settings.temperature['method'] = elem
+        elem = root.findtext('temperature_range')
+        if elem is not None:
+            settings.temperature['range'] = [float(x) for x in elem.split()]
+        elem = root.findtext('temperature_multipole')
+        if elem is not None:
+            settings.temperature['multipole'] = elem == 'true'
+
+        # Get trace
+        elem = root.find('trace')
+        if elem is not None:
+            settings.trace = [int(x) for x in elem.text.split()]
+
+        # Get track
+        elem = root.find('track')
+        if elem is not None:
+            settings.track = [int(x) for x in elem.text.split()]
+
+        # Get UFS mesh
+        elem = root.find('ufs_mesh')
+        if elem is not None:
+            settings.ufs_mesh = Mesh.from_xml_element(elem)
+
+        # Get resonance scattering
+        elem = root.find('resonance_scattering')
+        if elem is not None:
+            for entry in elem:
+                key = entry.tag
+                if key == 'enable':
+                    value = entry.text == 'true'
+                elif key == 'method':
+                    value = entry.text
+                elif key == 'energy_min':
+                    value = float(entry.text)
+                elif key == 'energy_max':
+                    value = float(entry.text)
+                elif key == 'nuclides':
+                    value = entry.text.split()
+                settings.resonance_scattering[key] = value
+
+        # TODO: Get volume calculations
+
+        # Get fission neutrons
+        elem = root.find('create_fission_neutrons')
+        if elem is not None:
+            settings.create_fission_neutrons = elem.text == 'true'
+
+        # Get log grid bins
+        elem = root.find('log_grid_bins')
+        if elem is not None:
+            settings.log_grid_bins = int(elem.text)
+
+        # Get dagmc
+        elem = root.find('dagmc')
+        if elem is not None:
+            settings.dagmc = elem.text == 'true'
+
+        return settings
diff --git a/openmc/source.py b/openmc/source.py
index 6ec882ca6a..e278d80892 100644
--- a/openmc/source.py
+++ b/openmc/source.py
@@ -2,8 +2,11 @@ from numbers import Real
 import sys
 from xml.etree import ElementTree as ET
 
-from openmc.stats.univariate import Univariate
-from openmc.stats.multivariate import UnitSphere, Spatial
+from openmc.stats.univariate import (Univariate, Discrete, Uniform, Maxwell,
+                                     Watt, Normal, Muir, Tabular)
+from openmc.stats.multivariate import (UnitSphere, Spatial, PolarAzimuthal,
+                                       Isotropic, Monodirectional, Box, Point,
+                                       CartesianIndependent)
 import openmc.checkvalue as cv
 
 
@@ -137,3 +140,72 @@ class Source(object):
         if self.energy is not None:
             element.append(self.energy.to_xml_element('energy'))
         return element
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate source from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.Source
+            Source generated from XML element
+
+        """
+        source = cls()
+
+        strength = elem.find('strength')
+        if strength is not None:
+            source.strength = float(strength.text)
+
+        particle = elem.find('particle')
+        if particle is not None:
+            source.particle = particle.text
+
+        filename = elem.find('file')
+        if filename is not None:
+            source.file = filename.text
+
+        space = elem.find('space')
+        if space is not None:
+            space_type = space.get('type')
+            if space_type == 'cartesian':
+                source.space = CartesianIndependent.from_xml_element(space)
+            elif space_type == 'box' or space_type == 'fission':
+                source.space = Box.from_xml_element(space)
+            elif space_type == 'point':
+                source.space = Point.from_xml_element(space)
+
+        angle = elem.find('angle')
+        if angle is not None:
+            angle_type = angle.get('type')
+            if angle_type == 'mu-phi':
+                source.angle = PolarAzimuthal.from_xml_element(angle)
+            elif angle_type == 'isotropic':
+                source.angle = Isotropic.from_xml_element(angle)
+            elif angle_type == 'monodirectional':
+                source.angle = Monodirectional.from_xml_element(angle)
+
+        energy = elem.find('energy')
+        if energy is not None:
+            energy_type = energy.get('type')
+            if energy_type == 'discrete':
+                source.energy = Discrete.from_xml_element(energy)
+            elif energy_type == 'uniform':
+                source.energy = Uniform.from_xml_element(energy)
+            elif energy_type == 'maxwell':
+                source.energy = Maxwell.from_xml_element(energy)
+            elif energy_type == 'watt':
+                source.energy = Watt.from_xml_element(energy)
+            elif energy_type == 'normal':
+                source.energy = Normal.from_xml_element(energy)
+            elif energy_type == 'muir':
+                source.energy = Muir.from_xml_element(energy)
+            elif energy_type == 'tabular':
+                source.energy = Tabular.from_xml_element(energy)
+
+        return source
diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py
index ac788c3443..61122ada21 100644
--- a/openmc/stats/multivariate.py
+++ b/openmc/stats/multivariate.py
@@ -47,6 +47,11 @@ class UnitSphere(metaclass=ABCMeta):
     def to_xml_element(self):
         return ''
 
+    @classmethod
+    @abstractmethod
+    def from_xml_element(cls, elem):
+        pass
+
 
 class PolarAzimuthal(UnitSphere):
     """Angular distribution represented by polar and azimuthal angles
@@ -121,6 +126,29 @@ class PolarAzimuthal(UnitSphere):
         element.append(self.phi.to_xml_element('phi'))
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate angular distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.PolarAzimuthal
+            Angular distribution generated from XML element
+
+        """
+        mu_phi = cls()
+        params = elem.findtext('parameters')
+        if params is not None:
+            mu_phi.reference_uvw = [float(x) for x in params.split()]
+        mu_phi.mu = openmc.stats.Univariate.from_xml_element(elem.find('mu'))
+        mu_phi.phi = openmc.stats.Univariate.from_xml_element(elem.find('phi'))
+        return mu_phi
+
 
 class Isotropic(UnitSphere):
     """Isotropic angular distribution.
@@ -143,6 +171,23 @@ class Isotropic(UnitSphere):
         element.set("type", "isotropic")
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate isotropic distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Isotropic
+            Isotropic distribution generated from XML element
+
+        """
+        return cls()
+
 
 class Monodirectional(UnitSphere):
     """Monodirectional angular distribution.
@@ -178,6 +223,27 @@ class Monodirectional(UnitSphere):
             element.set("reference_uvw", ' '.join(map(str, self.reference_uvw)))
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate monodirectional distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Monodirectional
+            Monodirectional distribution generated from XML element
+
+        """
+        monodirectional = cls()
+        params = elem.findtext('parameters')
+        if params is not None:
+            monodirectional.reference_uvw = [float(x) for x in params.split()]
+        return monodirectional
+
 
 class Spatial(metaclass=ABCMeta):
     """Distribution of locations in three-dimensional Euclidean space.
@@ -193,6 +259,11 @@ class Spatial(metaclass=ABCMeta):
     def to_xml_element(self):
         return ''
 
+    @classmethod
+    @abstractmethod
+    def from_xml_element(cls, elem):
+        pass
+
 
 class CartesianIndependent(Spatial):
     """Spatial distribution with independent x, y, and z distributions.
@@ -270,6 +341,26 @@ class CartesianIndependent(Spatial):
         element.append(self.z.to_xml_element('z'))
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate spatial distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.CartesianIndependent
+            Spatial distribution generated from XML element
+
+        """
+        x = openmc.stats.Univariate.from_xml_element(elem.find('x'))
+        y = openmc.stats.Univariate.from_xml_element(elem.find('y'))
+        z = openmc.stats.Univariate.from_xml_element(elem.find('z'))
+        return cls(x, y, z)
+
 
 class Box(Spatial):
     """Uniform distribution of coordinates in a rectangular cuboid.
@@ -351,6 +442,27 @@ class Box(Spatial):
                       ' '.join(map(str, self.upper_right))
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate box distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Box
+            Box distribution generated from XML element
+
+        """
+        only_fissionable = elem.get('type') == 'fission'
+        params = [float(x) for x in elem.findtext('parameters').split()]
+        lower_left = params[:len(params)//2]
+        upper_right = paramx[len(params)//2:]
+        return cls(lower_left, upper_right, only_fissionable)
+
 
 class Point(Spatial):
     """Delta function in three dimensions.
@@ -398,3 +510,21 @@ class Point(Spatial):
         params = ET.SubElement(element, "parameters")
         params.text = ' '.join(map(str, self.xyz))
         return element
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate point distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Point
+            Point distribution generated from XML element
+
+        """
+        xyz = [float(x) for x in elem.findtext('parameters').split()]
+        return cls(xyz)
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index e61f216ef0..6a83cc780e 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -32,6 +32,11 @@ class Univariate(EqualityMixin, metaclass=ABCMeta):
     def __len__(self):
         return 0
 
+    @classmethod
+    @abstractmethod
+    def from_xml_element(cls, elem):
+        pass
+
 
 class Discrete(Univariate):
     """Distribution characterized by a probability mass function.
@@ -110,6 +115,26 @@ class Discrete(Univariate):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate discrete distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Discrete
+            Discrete distribution generated from XML element
+
+        """
+        params = [float(x) for x in elem.findtext('parameters').split()]
+        x = params[:len(params)//2]
+        p = params[len(params)//2:]
+        return cls(x, p)
+
 
 class Uniform(Univariate):
     """Distribution with constant probability over a finite interval [a,b]
@@ -181,6 +206,26 @@ class Uniform(Univariate):
         element.set("parameters", '{} {}'.format(self.a, self.b))
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate uniform distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Uniform
+            Uniform distribution generated from XML element
+
+        """
+        params = elem.findtext('parameters').split()
+        a = float(params[0])
+        b = float(params[1])
+        return cls(a, b)
+
 
 class Maxwell(Univariate):
     """Maxwellian distribution in energy.
@@ -237,6 +282,24 @@ class Maxwell(Univariate):
         element.set("parameters", str(self.theta))
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate Maxwellian distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Maxwell
+            Maxwellian distribution generated from XML element
+
+        """
+        theta = float(elem.findtext('parameters'))
+        return cls(theta)
+
 
 class Watt(Univariate):
     r"""Watt fission energy spectrum.
@@ -308,6 +371,27 @@ class Watt(Univariate):
         element.set("parameters", '{} {}'.format(self.a, self.b))
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate Watt distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Watt
+            Watt distribution generated from XML element
+
+        """
+        params = elem.findtext('parameters').split()
+        a = float(params[0])
+        b = float(params[1])
+        return watt(a, b)
+
+
 class Normal(Univariate):
     r"""Normally distributed sampling.
 
@@ -377,6 +461,27 @@ class Normal(Univariate):
         element.set("parameters", '{} {}'.format(self.mean_value, self.std_dev))
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate Normal distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Normal
+            Normal distribution generated from XML element
+
+        """
+        params = elem.findtext('parameters').split()
+        mean_value = float(params[0])
+        std_dev = float(params[1])
+        return cls(mean_value, std_dev)
+
+
 class Muir(Univariate):
     """Muir energy spectrum.
 
@@ -465,6 +570,27 @@ class Muir(Univariate):
         element.set("parameters", '{} {} {}'.format(self._e0, self._m_rat, self._kt))
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate Muir distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Muir
+            Muir distribution generated from XML element
+
+        """
+        params = elem.findtext('parameters').split()
+        e0 = float(params[0])
+        m_rat = float(params[1])
+        kt = float(params[2])
+        return muir(e0, m_rat, kt)
+
 
 class Tabular(Univariate):
     """Piecewise continuous probability distribution.
@@ -561,6 +687,27 @@ class Tabular(Univariate):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate tabular distribution from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.Tabular
+            Tabular distribution generated from XML element
+
+        """
+        interpolation = elem.get('interpolation')
+        params = [float(x) for x in elem.findtext('parameters').split()]
+        x = params[:len(params)//2]
+        p = paramx[len(params)//2:]
+        return cls(x, p, interpolation)
+
 
 class Legendre(Univariate):
     r"""Probability density given by a Legendre polynomial expansion
@@ -607,6 +754,10 @@ class Legendre(Univariate):
     def to_xml_element(self, element_name):
         raise NotImplementedError
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        raise NotImplementedError
+
 
 class Mixture(Univariate):
     """Probability distribution characterized by a mixture of random variables.
@@ -660,3 +811,7 @@ class Mixture(Univariate):
 
     def to_xml_element(self, element_name):
         raise NotImplementedError
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        raise NotImplementedError