Merge branch 'develop' of https://github.com/openmc-dev/openmc into cmfd-added-funcs

.travis.yml

@@ -53,4 +53,5 @@ before_script:
 script:
   - ./tools/ci/travis-script.sh
 after_success:
-  - coveralls
+  - cpp-coveralls -i src -i include --exclude-pattern "/usr/*" --dump cpp_cov.json
+  - coveralls --merge=cpp_cov.json

CMakeLists.txt

@@ -95,6 +95,10 @@ if(optimize)
   list(REMOVE_ITEM cxxflags -O2)
   list(APPEND cxxflags -O3)
 endif()
+if(coverage)
+  list(APPEND cxxflags --coverage)
+  list(APPEND ldflags --coverage)
+endif()
 
 # Show flags being used
 message(STATUS "OpenMC C++ flags: ${cxxflags}")

MANIFEST.in

@@ -1,7 +1,10 @@
 include CMakeLists.txt
 include LICENSE
+include CODE_OF_CONDUCT.md
+include CONTRIBUTING.md
 include schemas.xml
 include pyproject.toml
+include pytest.ini
 include openmc/data/reconstruct.pyx
 include docs/source/_templates/layout.html
 include docs/sphinxext/LICENSE
@@ -33,4 +36,11 @@ recursive-include tests *.dat
 recursive-include tests *.h5
 recursive-include tests *.py
 recursive-include tests *.xml
+recursive-include vendor CMakeLists.txt
+recursive-include vendor *.cmake.in
+recursive-include vendor *.cc
+recursive-include vendor *.cpp
+recursive-include vendor *.hh
+recursive-include vendor *.hpp
 prune docs/build
+prune docs/source/pythonapi/generated/

docs/source/capi/index.rst

@@ -73,6 +73,17 @@ Functions
    :return: Return status (negative if an error occurred)
    :rtype: int
 
+.. c:function:: int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T)
+
+   Get the temperature of a cell
+
+   :param int32_t index: Index in the cells array
+   :param int32_t* instance: Which instance of the cell. If a null pointer is passed, the temperature
+                             of the first instance is returned.
+   :param double* T: temperature of the cell
+   :return: Return status (negative if an error occurred)
+   :rtype: int
+
 .. c:function:: int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices)
 
    Set the fill for a cell

docs/source/io_formats/nuclear_data.rst

@@ -165,6 +165,10 @@ Incident Photon Data
            - **J** (*double[][]*) -- Compton profile for each subshell in units
              of :math:`\hbar / (me^2)`
 
+**/<element>/heating/**
+
+:Datasets: - **xs** (*double[]*) -- Total heating cross section in [b-eV]
+
 **/<element>/incoherent/**
 
 :Datasets: - **xs** (*double[]*) -- Incoherent scattering cross section in [b]

docs/source/methods/neutron_physics.rst

@@ -86,7 +86,7 @@ Elastic Scattering
 ------------------
 
 Note that the multi-group mode makes no distinction between elastic or
-inelastic scattering reactions. The spceific multi-group scattering
+inelastic scattering reactions. The specific multi-group scattering
 implementation is discussed in the :ref:`multi-group-scatter` section.
 
 Elastic scattering refers to the process by which a neutron scatters off a

docs/source/usersguide/tallies.rst

@@ -256,9 +256,12 @@ The following tables show all valid scores:
     |inverse-velocity      |The flux-weighted inverse velocity where the       |
     |                      |velocity is in units of centimeters per second.    |
     +----------------------+---------------------------------------------------+
-    |heating               |Total neutron heating in units of eV per source    |
-    |                      |particle. This corresponds to MT=301 produced by   |
-    |                      |NJOY's HEATR module.                               |
+    |heating               |Total nuclear heating in units of eV per source    |
+    |                      |particle. For neutrons, this corresponds to MT=301 |
+    |                      |produced by NJOY's HEATR module while for photons, |
+    |                      |this is tallied from either direct photon energy   |
+    |                      |deposition (analog estimator) or pre-generated     |
+    |                      |photon heating number.                             |
     +----------------------+---------------------------------------------------+
     |kappa-fission         |The recoverable energy production rate due to      |
     |                      |fission. The recoverable energy is defined as the  |

include/openmc/capi.h

@@ -13,6 +13,7 @@ extern "C" {
   int openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n);
   int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n);
   int openmc_cell_get_id(int32_t index, int32_t* id);
+  int openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T);
   int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices);
   int openmc_cell_set_id(int32_t index, int32_t id);
   int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance);

include/openmc/constants.h

@@ -105,9 +105,10 @@ constexpr std::array<const char*, 39> SUBSHELLS =  {
   "Q1", "Q2", "Q3"
 };
 
-// Void material
+// Void material and nuclide
 // TODO: refactor and remove
 constexpr int MATERIAL_VOID {-1};
+constexpr int NUCLIDE_NONE  {-1};
 
 // ============================================================================
 // CROSS SECTION RELATED CONSTANTS
@@ -127,6 +128,7 @@ constexpr int TEMPERATURE_INTERPOLATION {2};
 
 // Reaction types
 // TODO: Convert to enum
+constexpr int REACTION_NONE {0};
 constexpr int TOTAL_XS {1};
 constexpr int ELASTIC {2};
 constexpr int N_NONELASTIC {3};
@@ -230,7 +232,7 @@ constexpr int N_XD    {204};
 constexpr int N_XT    {205};
 constexpr int N_X3HE  {206};
 constexpr int N_XA    {207};
-constexpr int HEATING {301};
+constexpr int NEUTRON_HEATING {301};
 constexpr int DAMAGE_ENERGY {444};
 constexpr int COHERENT {502};
 constexpr int INCOHERENT {504};
@@ -344,6 +346,7 @@ constexpr int ESTIMATOR_COLLISION   {3};
 // TODO: Convert to enum
 constexpr int EVENT_SURFACE {-2};
 constexpr int EVENT_LATTICE {-1};
+constexpr int EVENT_KILL    {0};
 constexpr int EVENT_SCATTER {1};
 constexpr int EVENT_ABSORB  {2};
 
@@ -366,6 +369,7 @@ constexpr int SCORE_INVERSE_VELOCITY   {-13}; // flux-weighted inverse velocity
 constexpr int SCORE_FISS_Q_PROMPT      {-14}; // prompt fission Q-value
 constexpr int SCORE_FISS_Q_RECOV       {-15}; // recoverable fission Q-value
 constexpr int SCORE_DECAY_RATE         {-16}; // delayed neutron precursor decay rate
+constexpr int SCORE_HEATING            {-17}; // nuclear heating (neutron or photon)
 
 // Tally map bin finding
 constexpr int NO_BIN_FOUND {-1};

include/openmc/particle.h

@@ -184,7 +184,7 @@ public:
   //! \param u Direction of the secondary particle
   //! \param E Energy of the secondary particle in [eV]
   //! \param type Particle type
-  void create_secondary(Direction u, double E, Type type) const;
+  void create_secondary(Direction u, double E, Type type);
 
   //! initialize from a source site
   //
@@ -261,6 +261,7 @@ public:
 
   // Post-collision physical data
   int n_bank_ {0};        //!< number of fission sites banked
+  int n_bank_second_ {0}; //!< number of secondary particles banked
   double wgt_bank_ {0.0}; //!< weight of fission sites banked
   int n_delayed_bank_[MAX_DELAYED_GROUPS];  //!< number of delayed fission
                                             //!< sites banked

include/openmc/photon.h

@@ -67,6 +67,7 @@ public:
   xt::xtensor<double, 1> pair_production_total_;
   xt::xtensor<double, 1> pair_production_electron_;
   xt::xtensor<double, 1> pair_production_nuclear_;
+  xt::xtensor<double, 1> heating_;
 
   // Form factors
   Tabulated1D incoherent_form_factor_;

include/openmc/string_utils.h

@@ -10,6 +10,8 @@ std::string& strtrim(std::string& s);
 
 char* strtrim(char* c_str);
 
+std::string to_element(const std::string& name);
+
 void to_lower(std::string& str);
 
 int word_count(std::string const& str);

openmc/capi/cell.py

@@ -24,6 +24,10 @@ _dll.openmc_cell_get_fill.argtypes = [
     c_int32, POINTER(c_int), POINTER(POINTER(c_int32)), POINTER(c_int32)]
 _dll.openmc_cell_get_fill.restype = c_int
 _dll.openmc_cell_get_fill.errcheck = _error_handler
+_dll.openmc_cell_get_temperature.argtypes = [
+    c_int32, POINTER(c_int32), POINTER(c_double)]
+_dll.openmc_cell_get_temperature.restype = c_int
+_dll.openmc_cell_get_temperature.errcheck = _error_handler
 _dll.openmc_cell_set_fill.argtypes = [
     c_int32, c_int, c_int32, POINTER(c_int32)]
 _dll.openmc_cell_set_fill.restype = c_int
@@ -128,6 +132,23 @@ class Cell(_FortranObjectWithID):
             indices = (c_int32*1)(-1)
             _dll.openmc_cell_set_fill(self._index, 1, 1, indices)
 
+    def get_temperature(self, instance=None):
+        """Get the temperature of a cell
+
+        Parameters
+        ----------
+        instance: int or None
+            Which instance of the cell
+
+        """
+
+        if instance is not None:
+            instance = c_int32(instance)
+
+        T = c_double()
+        _dll.openmc_cell_get_temperature(self._index, instance, T)
+        return T.value
+
     def set_temperature(self, T, instance=None):
         """Set the temperature of a cell
 
@@ -139,7 +160,11 @@ class Cell(_FortranObjectWithID):
             Which instance of the cell
 
         """
-        _dll.openmc_cell_set_temperature(self._index, T, c_int32(instance))
+
+        if instance is not None:
+            instance = c_int32(instance)
+
+        _dll.openmc_cell_set_temperature(self._index, T, instance)
 
 
 class _CellMapping(Mapping):

openmc/capi/core.py

@@ -243,18 +243,9 @@ def keff():
         Mean k-eigenvalue and standard deviation of the mean
 
     """
-    n = openmc.capi.num_realizations()
-    if n > 3:
-        # Use the combined estimator if there are enough realizations
-        k = (c_double*2)()
-        _dll.openmc_get_keff(k)
-        return tuple(k)
-    else:
-        # Otherwise, return the tracklength estimator
-        mean = c_double.in_dll(_dll, 'keff').value
-        std_dev = c_double.in_dll(_dll, 'keff_std').value \
-                  if n > 1 else np.inf
-        return (mean, std_dev)
+    k = (c_double*2)()
+    _dll.openmc_get_keff(k)
+    return tuple(k)
 
 
 def master():

openmc/capi/tally.py

@@ -89,7 +89,7 @@ _SCORES = {
     -5: 'absorption', -6: 'fission', -7: 'nu-fission', -8: 'kappa-fission',
     -9: 'current', -10: 'events', -11: 'delayed-nu-fission',
     -12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt',
-    -15: 'fission-q-recoverable', -16: 'decay-rate'
+    -15: 'fission-q-recoverable', -16: 'decay-rate', -17: 'heating'
 }
 _ESTIMATORS = {
     1: 'analog', 2: 'tracklength', 3: 'collision'

openmc/data/function.py

@@ -3,6 +3,7 @@ from collections.abc import Iterable, Callable
 from functools import reduce
 from itertools import zip_longest
 from numbers import Real, Integral
+from math import exp, log
 
 import numpy as np
 
@@ -153,13 +154,11 @@ class Tabulated1D(Function1D):
         self.y = np.asarray(y)
 
     def __call__(self, x):
-        # Check if input is array or scalar
-        if isinstance(x, Iterable):
-            iterable = True
-            x = np.array(x)
-        else:
-            iterable = False
-            x = np.array([x], dtype=float)
+        # Check if input is scalar
+        if not isinstance(x, Iterable):
+            return self._interpolate_scalar(x)
+
+        x = np.array(x)
 
         # Create output array
         y = np.zeros_like(x)
@@ -208,7 +207,46 @@ class Tabulated1D(Function1D):
         y[np.isclose(x, self.x[0], atol=1e-14)] = self.y[0]
         y[np.isclose(x, self.x[-1], atol=1e-14)] = self.y[-1]
 
-        return y if iterable else y[0]
+        return y
+
+    def _interpolate_scalar(self, x):
+        if x <= self._x[0]:
+            return self._y[0]
+        elif x >= self._x[-1]:
+            return self._y[-1]
+
+        # Get the index for interpolation
+        idx = np.searchsorted(self._x, x, side='right') - 1
+
+        # Loop over interpolation regions
+        for b, p in zip(self.breakpoints, self.interpolation):
+            if idx < b - 1:
+                break
+
+        xi = self._x[idx]       # low edge of the corresponding bin
+        xi1 = self._x[idx + 1]  # high edge of the corresponding bin
+        yi = self._y[idx]
+        yi1 = self._y[idx + 1]
+
+        if p == 1:
+            # Histogram
+            return yi
+
+        elif p == 2:
+            # Linear-linear
+            return yi + (x - xi)/(xi1 - xi)*(yi1 - yi)
+
+        elif p == 3:
+            # Linear-log
+            return yi + log(x/xi)/log(xi1/xi)*(yi1 - yi)
+
+        elif p == 4:
+            # Log-linear
+            return yi*exp((x - xi)/(xi1 - xi)*log(yi1/yi))
+
+        elif p == 5:
+            # Log-log
+            return yi*exp(log(x/xi)/log(xi1/xi)*log(yi1/yi))
 
     def __len__(self):
         return len(self.x)

openmc/data/photon.py

@@ -3,12 +3,14 @@ from collections.abc import Mapping, Callable
 from copy import deepcopy
 from io import StringIO
 from numbers import Integral, Real
+from math import pi, sqrt
 import os
 
 import h5py
 import numpy as np
 import pandas as pd
 from scipy.interpolate import CubicSpline
+from scipy.integrate import quad
 
 from openmc.mixin import EqualityMixin
 import openmc.checkvalue as cv
@@ -19,6 +21,14 @@ from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record
 from .function import Tabulated1D
 
 
+# Constants
+MASS_ELECTRON_EV = 0.5109989461e6 # Electron mass energy
+PLANCK_C = 1.2398419739062977e4 # Planck's constant times c in eV-Angstroms
+FINE_STRUCTURE = 137.035999139 # Inverse fine structure constant
+CM_PER_ANGSTROM = 1.0e-8
+# classical electron radius in cm
+R0 = CM_PER_ANGSTROM * PLANCK_C / (2.0 * pi * FINE_STRUCTURE * MASS_ELECTRON_EV)
+
 # Electron subshell labels
 _SUBSHELLS = [None, 'K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5',
               'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'O1', 'O2', 'O3',
@@ -33,6 +43,7 @@ _REACTION_NAME = {
     516: ('Total pair production', 'pair_production_total'),
     517: ('Pair production, nuclear field', 'pair_production_nuclear'),
     522: ('Photoelectric absorption', 'photoelectric'),
+    525: ('Heating', 'heating'),
     526: ('Electro-atomic scattering', 'electro_atomic_scat'),
     527: ('Electro-atomic bremsstrahlung', 'electro_atomic_brem'),
     528: ('Electro-atomic excitation', 'electro_atomic_excit'),
@@ -150,6 +161,7 @@ class AtomicRelaxation(EqualityMixin):
         self.binding_energy = binding_energy
         self.num_electrons = num_electrons
         self.transitions = transitions
+        self._e_fluorescence = {}
 
     @property
     def binding_energy(self):
@@ -379,6 +391,49 @@ class AtomicRelaxation(EqualityMixin):
                  _SUBSHELLS, range(len(_SUBSHELLS)))
             group.create_dataset('transitions', data=df.values.astype(float))
 
+    def energy_fluorescence(self, shell):
+        """Compute expected energy of fluorescent photons for the shell
+
+        Parameters
+        ----------
+        shell : str
+            The subshell to compute
+
+        Returns
+        -------
+        float
+            Energy of fluorescent photons
+
+        """
+
+        if shell not in self.binding_energy:
+            raise KeyError('Invalid shell {}.'.format(shell))
+
+        if shell in self._e_fluorescence:
+            # Already computed
+            return self._e_fluorescence[shell]
+        e = 0.0
+        if shell not in self.transitions or self.transitions[shell].empty:
+            e = self.binding_energy[shell]
+        else:
+            df = self.transitions[shell]
+            for primary, secondary, energy, prob in df.itertuples(index=False):
+                e_row = 0.0
+                if secondary is None:
+                    # Fluorescent photon release in radiative transition
+                    e_row += energy
+                else:
+                    # Fill the hole left by auger electron
+                    e_row += self.energy_fluorescence(secondary)
+
+                # Fill the photoelectron hole
+                e_row += self.energy_fluorescence(primary)
+
+                # Expected fluorescent photon energy
+                e += e_row * prob
+
+        self._e_fluorescence[shell] = e
+        return e
 
 class IncidentPhoton(EqualityMixin):
     r"""Photon interaction data.
@@ -499,9 +554,14 @@ class IncidentPhoton(EqualityMixin):
 
         # Read each reaction
         data = cls(Z)
-        for mt in (502, 504, 515, 522):
+        for mt in (502, 504, 515, 522, 525):
             data.reactions[mt] = PhotonReaction.from_ace(ace, mt)
 
+        # Get heating cross sections [eV-barn] from factors [eV per collision]
+        # by multiplying with total xs
+        data.reactions[525].xs.y *= sum([data.reactions[mt].xs.y for mt in
+                                         (502, 504, 515, 522)])
+
         # Compton profiles
         n_shell = ace.nxs[5]
         if n_shell != 0:
@@ -631,6 +691,9 @@ class IncidentPhoton(EqualityMixin):
         # Add bremsstrahlung DCS data
         data._add_bremsstrahlung()
 
+        # Add heating cross sections
+        data._compute_heating()
+
         return data
 
     @classmethod
@@ -751,7 +814,7 @@ class IncidentPhoton(EqualityMixin):
         designators = []
         for mt, rx in self.reactions.items():
             name, key = _REACTION_NAME[mt]
-            if mt in [502, 504, 515, 517, 522]:
+            if mt in (502, 504, 515, 517, 522, 525):
                 sub_group = group.create_group(key)
             elif mt >= 534 and mt <= 572:
                 # Subshell
@@ -858,6 +921,80 @@ class IncidentPhoton(EqualityMixin):
         self.bremsstrahlung['photon_energy'] = _BREMSSTRAHLUNG['photon_energy']
         self.bremsstrahlung.update(_BREMSSTRAHLUNG[self.atomic_number])
 
+    def _compute_heating(self):
+        r"""Compute heating cross sections (KERMA)
+
+        Photon energy is deposited as energy loss in three reactions:
+        incoherent scattering, pair production and photoelectric effect.
+        The point-wise heating cross section is calculated as:
+
+        .. math::
+            \sigma_{Hx}(E) &= (E - \overline{E}_x(E)) \cdot \sigma_x(E), x \in \left\{I, PP, PE \right\}
+
+            \overline{E}_I(E) &= \frac {\int E' \sigma_I (E,E',\mu) d\mu} {\int \sigma_I (E,E',\mu) d\mu}
+
+            \overline{E}_{PP} &= 2 m_e c^2 = 1.022 \times 10^6 eV
+
+            \overline{E}_{PE} &= E(\text{fluorescent photons})
+
+        The differential cross section representation for incoherent
+        scattering can be found in the theory manual.
+
+        """
+
+        # Determine a union energy grid
+        energy = np.array([])
+        for mt in (504, 515, 517, 522):
+            if mt in self:
+                energy = np.union1d(energy, self[mt].xs.x)
+
+        heating_xs = np.zeros_like(energy)
+
+        # Incoherent scattering
+        if 504 in self:
+            rx = self[504]
+
+            def dsigma_dmu(mu, E):
+                k = E / MASS_ELECTRON_EV
+                krat = 1.0 / (1.0 + k * (1.0 - mu))
+                x = E * sqrt(0.5 * (1.0 - mu)) / PLANCK_C
+                return pi * R0*R0 * krat*krat * (krat + 1/krat +
+                       mu*mu - 1.0) * rx.scattering_factor(x)
+
+            def eout_dsigma_dmu(mu, E):
+                Eout = E / (1.0 + E / MASS_ELECTRON_EV * (1.0 - mu))
+                return Eout * dsigma_dmu(mu, E)
+
+            def eout_average(E):
+                integral_sigma = quad(dsigma_dmu, -1.0, 1.0,
+                                      args=(E,), epsabs=0.0, epsrel=1e-3)[0]
+                integral_sigma_e = quad(eout_dsigma_dmu, -1.0, 1.0,
+                                        args=(E,), epsabs=0.0, epsrel=1e-3)[0]
+                return integral_sigma_e / integral_sigma
+
+            e_out = np.vectorize(eout_average)(energy)
+            heating_xs += (energy - e_out) * rx.xs(energy)
+
+        # Pair production, electron field
+        if 515 in self:
+            heating_xs += (energy - 2*MASS_ELECTRON_EV)*self[515].xs(energy)
+
+        # Pair production, nuclear field
+        if 517 in self:
+            heating_xs += (energy - 2*MASS_ELECTRON_EV)*self[517].xs(energy)
+
+        # Photoelectric effect
+        if 522 in self:
+            # Account for fluorescent photons
+            for mt, rx in self.reactions.items():
+                if mt >= 534 and mt <= 572:
+                    shell = _REACTION_NAME[mt][1]
+                    e_f = self.atomic_relaxation.energy_fluorescence(shell)
+                    heating_xs += (energy - e_f) * rx.xs(energy)
+
+        heat_rx = PhotonReaction(525)
+        heat_rx.xs = Tabulated1D(energy, heating_xs, [energy.size], [5])
+        self.reactions[525] = heat_rx
 
 class PhotonReaction(EqualityMixin):
     """Photon-induced reaction
@@ -972,14 +1109,21 @@ class PhotonReaction(EqualityMixin):
         elif mt == 522:
             # Photoelectric
             idx = ace.jxs[1] + 3*n
+        elif mt == 525:
+            # Heating
+            idx = ace.jxs[5]
         else:
             raise ValueError('ACE photoatomic cross sections do not have '
                              'data for MT={}.'.format(mt))
 
         # Store cross section
         xs = ace.xss[idx : idx+n].copy()
-        nonzero = (xs != 0.0)
-        xs[nonzero] = np.exp(xs[nonzero])
+        if mt == 525:
+            # Get heating factors in [eV per collision]
+            xs *= EV_PER_MEV
+        else:
+            nonzero = (xs != 0.0)
+            xs[nonzero] = np.exp(xs[nonzero])
         rx.xs = Tabulated1D(energy, xs, [n], [5])
 
         # Get form factors for incoherent/coherent scattering

openmc/mesh.py

@@ -7,6 +7,7 @@ import numpy as np
 
 import openmc.checkvalue as cv
 import openmc
+from openmc._xml import get_text
 from openmc.mixin import EqualityMixin, IDManagerMixin
 
 
@@ -230,8 +231,9 @@ class Mesh(IDManagerMixin):
         element.set("id", str(self._id))
         element.set("type", self._type)
 
-        subelement = ET.SubElement(element, "dimension")
-        subelement.text = ' '.join(map(str, self._dimension))
+        if self._dimension is not None:
+            subelement = ET.SubElement(element, "dimension")
+            subelement.text = ' '.join(map(str, self._dimension))
 
         subelement = ET.SubElement(element, "lower_left")
         subelement.text = ' '.join(map(str, self._lower_left))
@@ -246,6 +248,46 @@ 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(get_text(elem, 'id'))
+        mesh = cls(mesh_id)
+
+        mesh_type = get_text(elem, 'type')
+        if mesh_type is not None:
+            mesh.type = mesh_type
+
+        dimension = get_text(elem, 'dimension')
+        if dimension is not None:
+            mesh.dimension = [int(x) for x in dimension.split()]
+
+        lower_left = get_text(elem, 'lower_left')
+        if lower_left is not None:
+            mesh.lower_left = [float(x) for x in lower_left.split()]
+
+        upper_right = get_text(elem, 'upper_right')
+        if upper_right is not None:
+            mesh.upper_right = [float(x) for x in upper_right.split()]
+
+        width = get_text(elem, '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

openmc/plots.py

@@ -439,7 +439,7 @@ class Plot(IDManagerMixin):
         string += '{: <16}=\t{}\n'.format('\tBasis', self._basis)
         string += '{: <16}=\t{}\n'.format('\tWidth', self._width)
         string += '{: <16}=\t{}\n'.format('\tOrigin', self._origin)
-        string += '{: <16}=\t{}\n'.format('\tPixels', self._origin)
+        string += '{: <16}=\t{}\n'.format('\tPixels', self._pixels)
         string += '{: <16}=\t{}\n'.format('\tColor by', self._color_by)
         string += '{: <16}=\t{}\n'.format('\tBackground', self._background)
         string += '{: <16}=\t{}\n'.format('\tMask components',

openmc/settings.py

@@ -7,7 +7,7 @@ import sys
 
 import numpy as np
 
-from openmc._xml import clean_indentation
+from openmc._xml import clean_indentation, get_text
 import openmc.checkvalue as cv
 from openmc import VolumeCalculation, Source, Mesh
 
@@ -174,7 +174,6 @@ class Settings(object):
         self._source = cv.CheckedList(Source, 'source distributions')
 
         self._confidence_intervals = None
-        self._cross_sections = None
         self._electron_treatment = None
         self._photon_transport = None
         self._ptables = None
@@ -552,7 +551,8 @@ class Settings(object):
     @entropy_mesh.setter
     def entropy_mesh(self, entropy):
         cv.check_type('entropy mesh', entropy, Mesh)
-        cv.check_length('entropy mesh dimension', entropy.dimension, 3)
+        if entropy.dimension:
+            cv.check_length('entropy mesh dimension', entropy.dimension, 3)
         cv.check_length('entropy mesh lower-left corner', entropy.lower_left, 3)
         cv.check_length('entropy mesh upper-right corner', entropy.upper_right, 3)
         self._entropy_mesh = entropy
@@ -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)
@@ -927,6 +927,237 @@ class Settings(object):
             elem = ET.SubElement(root, "dagmc")
             elem.text = str(self._dagmc).lower()
 
+    def _eigenvalue_from_xml_element(self, root):
+        elem = root.find('eigenvalue')
+        if elem is not None:
+            self._run_mode_from_xml_element(elem)
+            self._particles_from_xml_element(elem)
+            self._batches_from_xml_element(elem)
+            self._inactive_from_xml_element(elem)
+            self._generations_per_batch_from_xml_element(elem)
+
+    def _run_mode_from_xml_element(self, root):
+        text = get_text(root, 'run_mode')
+        if text is not None:
+            self.run_mode = text
+
+    def _particles_from_xml_element(self, root):
+        text = get_text(root, 'particles')
+        if text is not None:
+            self.particles = int(text)
+
+    def _batches_from_xml_element(self, root):
+        text = get_text(root, 'batches')
+        if text is not None:
+            self.batches = int(text)
+
+    def _inactive_from_xml_element(self, root):
+        text = get_text(root, 'inactive')
+        if text is not None:
+            self.inactive = int(text)
+
+    def _generations_per_batch_from_xml_element(self, root):
+        text = get_text(root, 'generations_per_batch')
+        if text is not None:
+            self.generations_per_batch = int(text)
+
+    def _keff_trigger_from_xml_element(self, root):
+        elem = root.find('keff_trigger')
+        if elem is not None:
+            trigger = get_text(elem, 'type')
+            threshold = float(get_text(elem, 'threshold'))
+            self.keff_trigger = {'type': trigger, 'threshold': threshold}
+
+    def _source_from_xml_element(self, root):
+        for elem in root.findall('source'):
+            self.source.append(Source.from_xml_element(elem))
+
+    def _output_from_xml_element(self, root):
+        elem = root.find('output')
+        if elem is not None:
+            self.output = {}
+            for key in ('summary', 'tallies', 'path'):
+                value = get_text(elem, key)
+                if value is not None:
+                    if key in ('summary', 'tallies'):
+                        value = value in ('true', '1')
+                self.output[key] = value
+
+    def _statepoint_from_xml_element(self, root):
+        elem = root.find('state_point')
+        if elem is not None:
+            text = get_text(elem, 'batches')
+            if text is not None:
+                self.statepoint['batches'] = [int(x) for x in text.split()]
+
+    def _sourcepoint_from_xml_element(self, root):
+        elem = root.find('source_point')
+        if elem is not None:
+            for key in ('separate', 'write', 'overwrite_latest', 'batches'):
+                value = get_text(elem, key)
+                if value is not None:
+                    if key in ('separate', 'write'):
+                        value = value in ('true', '1')
+                    elif key == 'overwrite_latest':
+                        value = value in ('true', '1')
+                        key = 'overwrite'
+                    else:
+                        value = [int(x) for x in value.split()]
+                    self.sourcepoint[key] = value
+
+    def _confidence_intervals_from_xml_element(self, root):
+        text = get_text(root, 'confidence_intervals')
+        if text is not None:
+            self.confidence_intervals = text in ('true', '1')
+
+    def _electron_treatment_from_xml_element(self, root):
+        text = get_text(root, 'electron_treatment')
+        if text is not None:
+            self.electron_treatment = text
+
+    def _energy_mode_from_xml_element(self, root):
+        text = get_text(root, 'energy_mode')
+        if text is not None:
+            self.energy_mode = text
+
+    def _max_order_from_xml_element(self, root):
+        text = get_text(root, 'max_order')
+        if text is not None:
+            self.max_order = int(text)
+
+    def _photon_transport_from_xml_element(self, root):
+        text = get_text(root, 'photon_transport')
+        if text is not None:
+            self.photon_transport = text in ('true', '1')
+
+    def _ptables_from_xml_element(self, root):
+        text = get_text(root, 'ptables')
+        if text is not None:
+            self.ptables = text in ('true', '1')
+
+    def _seed_from_xml_element(self, root):
+        text = get_text(root, 'seed')
+        if text is not None:
+            self.seed = int(text)
+
+    def _survival_biasing_from_xml_element(self, root):
+        text = get_text(root, 'survival_biasing')
+        if text is not None:
+            self.survival_biasing = text in ('true', '1')
+
+    def _cutoff_from_xml_element(self, root):
+        elem = root.find('cutoff')
+        if elem is not None:
+            self.cutoff = {}
+            for key in ('energy_neutron', 'energy_photon', 'energy_electron',
+                        'energy_positron', 'weight', 'weight_avg'):
+                value = get_text(elem, key)
+                if value is not None:
+                    self.cutoff[key] = float(value)
+
+    def _entropy_mesh_from_xml_element(self, root):
+        text = get_text(root, 'entropy_mesh')
+        if text is not None:
+            path = "./mesh[@id='{}']".format(int(text))
+            elem = root.find(path)
+            if elem is not None:
+                self.entropy_mesh = Mesh.from_xml_element(elem)
+
+    def _trigger_from_xml_element(self, root):
+        elem = root.find('trigger')
+        if elem is not None:
+            self.trigger_active = get_text(elem, 'active') in ('true', '1')
+            text = get_text(elem, 'max_batches')
+            if text is not None:
+                self.trigger_max_batches = int(text)
+            text = get_text(elem, 'batch_interval')
+            if text is not None:
+                self.trigger_batch_interval = int(text)
+
+    def _no_reduce_from_xml_element(self, root):
+        text = get_text(root, 'no_reduce')
+        if text is not None:
+            self.no_reduce = text in ('true', '1')
+
+    def _verbosity_from_xml_element(self, root):
+        text = get_text(root, 'verbosity')
+        if text is not None:
+            self.verbosity = int(text)
+
+    def _tabular_legendre_from_xml_element(self, root):
+        elem = root.find('tabular_legendre')
+        if elem is not None:
+            text = get_text(elem, 'enable')
+            self.tabular_legendre['enable'] = text in ('true', '1')
+            text = get_text(elem, 'num_points')
+            if text is not None:
+                self.tabular_legendre['num_points'] = int(text)
+
+    def _temperature_from_xml_element(self, root):
+        text = get_text(root, 'temperature_default')
+        if text is not None:
+            self.temperature['default'] = float(text)
+        text = get_text(root, 'temperature_tolerance')
+        if text is not None:
+            self.temperature['tolerance'] = float(text)
+        text = get_text(root, 'temperature_method')
+        if text is not None:
+            self.temperature['method'] = text
+        text = get_text(root, 'temperature_range')
+        if text is not None:
+            self.temperature['range'] = [float(x) for x in text.split()]
+        text = get_text(root, 'temperature_multipole')
+        if text is not None:
+            self.temperature['multipole'] = text in ('true', '1')
+
+    def _trace_from_xml_element(self, root):
+        text = get_text(root, 'trace')
+        if text is not None:
+            self.trace = [int(x) for x in text.split()]
+
+    def _track_from_xml_element(self, root):
+        text = get_text(root, 'track')
+        if text is not None:
+            self.track = [int(x) for x in text.split()]
+
+    def _ufs_mesh_from_xml_element(self, root):
+        text = get_text(root, 'ufs_mesh')
+        if text is not None:
+            path = "./mesh[@id='{}']".format(int(text))
+            elem = root.find(path)
+            if elem is not None:
+                self.ufs_mesh = Mesh.from_xml_element(elem)
+
+    def _resonance_scattering_from_xml_element(self, root):
+        elem = root.find('resonance_scattering')
+        if elem is not None:
+            keys = ('enable', 'method', 'energy_min', 'energy_max', 'nuclides')
+            for key in keys:
+                value = get_text(elem, key)
+                if value is not None:
+                    if key == 'enable':
+                        value = value in ('true', '1')
+                    elif key in ('energy_min', 'energy_max'):
+                        value = float(value)
+                    elif key == 'nuclides':
+                        value = value.split()
+                    self.resonance_scattering[key] = value
+
+    def _create_fission_neutrons_from_xml_element(self, root):
+        text = get_text(root, 'create_fission_neutrons')
+        if text is not None:
+            self.create_fission_neutrons = text in ('true', '1')
+
+    def _log_grid_bins_from_xml_element(self, root):
+        text = get_text(root, 'log_grid_bins')
+        if text is not None:
+            self.log_grid_bins = int(text)
+
+    def _dagmc_from_xml_element(self, root):
+        text = get_text(root, 'dagmc')
+        if text is not None:
+            self.dagmc = text in ('true', '1')
+
     def export_to_xml(self, path='settings.xml'):
         """Export simulation settings to an XML file.
 
@@ -985,3 +1216,60 @@ 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()
+        settings._eigenvalue_from_xml_element(root)
+        settings._run_mode_from_xml_element(root)
+        settings._particles_from_xml_element(root)
+        settings._batches_from_xml_element(root)
+        settings._inactive_from_xml_element(root)
+        settings._generations_per_batch_from_xml_element(root)
+        settings._keff_trigger_from_xml_element(root)
+        settings._source_from_xml_element(root)
+        settings._output_from_xml_element(root)
+        settings._statepoint_from_xml_element(root)
+        settings._sourcepoint_from_xml_element(root)
+        settings._confidence_intervals_from_xml_element(root)
+        settings._electron_treatment_from_xml_element(root)
+        settings._energy_mode_from_xml_element(root)
+        settings._max_order_from_xml_element(root)
+        settings._photon_transport_from_xml_element(root)
+        settings._ptables_from_xml_element(root)
+        settings._seed_from_xml_element(root)
+        settings._survival_biasing_from_xml_element(root)
+        settings._cutoff_from_xml_element(root)
+        settings._entropy_mesh_from_xml_element(root)
+        settings._trigger_from_xml_element(root)
+        settings._no_reduce_from_xml_element(root)
+        settings._verbosity_from_xml_element(root)
+        settings._tabular_legendre_from_xml_element(root)
+        settings._temperature_from_xml_element(root)
+        settings._trace_from_xml_element(root)
+        settings._track_from_xml_element(root)
+        settings._ufs_mesh_from_xml_element(root)
+        settings._resonance_scattering_from_xml_element(root)
+        settings._create_fission_neutrons_from_xml_element(root)
+        settings._log_grid_bins_from_xml_element(root)
+        settings._dagmc_from_xml_element(root)
+
+        # TODO: Get volume calculations
+
+        return settings

openmc/source.py

@@ -2,6 +2,7 @@ from numbers import Real
 import sys
 from xml.etree import ElementTree as ET
 
+from openmc._xml import get_text
 from openmc.stats.univariate import Univariate
 from openmc.stats.multivariate import UnitSphere, Spatial
 import openmc.checkvalue as cv
@@ -137,3 +138,46 @@ 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 = get_text(elem, 'strength')
+        if strength is not None:
+            source.strength = float(strength)
+
+        particle = get_text(elem, 'particle')
+        if particle is not None:
+            source.particle = particle
+
+        filename = get_text(elem, 'file')
+        if filename is not None:
+            source.file = filename
+
+        space = elem.find('space')
+        if space is not None:
+            source.space = Spatial.from_xml_element(space)
+
+        angle = elem.find('angle')
+        if angle is not None:
+            source.angle = UnitSphere.from_xml_element(angle)
+
+        energy = elem.find('energy')
+        if energy is not None:
+            source.energy = Univariate.from_xml_element(energy)
+
+        return source

openmc/stats/multivariate.py

@@ -8,6 +8,7 @@ from xml.etree import ElementTree as ET
 import numpy as np
 
 import openmc.checkvalue as cv
+from openmc._xml import get_text
 from openmc.stats.univariate import Univariate, Uniform
 
 
@@ -47,6 +48,17 @@ class UnitSphere(metaclass=ABCMeta):
     def to_xml_element(self):
         return ''
 
+    @classmethod
+    @abstractmethod
+    def from_xml_element(cls, elem):
+        distribution = get_text(elem, 'type')
+        if distribution == 'mu-phi':
+            return PolarAzimuthal.from_xml_element(elem)
+        elif distribution == 'isotropic':
+            return Isotropic.from_xml_element(elem)
+        elif distribution == 'monodirectional':
+            return Monodirectional.from_xml_element(elem)
+
 
 class PolarAzimuthal(UnitSphere):
     """Angular distribution represented by polar and azimuthal angles
@@ -121,6 +133,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 = get_text(elem, 'parameters')
+        if params is not None:
+            mu_phi.reference_uvw = [float(x) for x in params.split()]
+        mu_phi.mu = Univariate.from_xml_element(elem.find('mu'))
+        mu_phi.phi = Univariate.from_xml_element(elem.find('phi'))
+        return mu_phi
+
 
 class Isotropic(UnitSphere):
     """Isotropic angular distribution.
@@ -143,6 +178,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 +230,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 = get_text(elem, '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 +266,17 @@ class Spatial(metaclass=ABCMeta):
     def to_xml_element(self):
         return ''
 
+    @classmethod
+    @abstractmethod
+    def from_xml_element(cls, elem):
+        distribution = get_text(elem, 'type')
+        if distribution == 'cartesian':
+            return CartesianIndependent.from_xml_element(elem)
+        elif distribution == 'box' or distribution == 'fission':
+            return Box.from_xml_element(elem)
+        elif distribution == 'point':
+            return Point.from_xml_element(elem)
+
 
 class CartesianIndependent(Spatial):
     """Spatial distribution with independent x, y, and z distributions.
@@ -270,6 +354,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 = Univariate.from_xml_element(elem.find('x'))
+        y = Univariate.from_xml_element(elem.find('y'))
+        z = 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 +455,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 = get_text(elem, 'type') == 'fission'
+        params = [float(x) for x in get_text(elem, 'parameters').split()]
+        lower_left = params[:len(params)//2]
+        upper_right = params[len(params)//2:]
+        return cls(lower_left, upper_right, only_fissionable)
+
 
 class Point(Spatial):
     """Delta function in three dimensions.
@@ -398,3 +523,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 get_text(elem, 'parameters').split()]
+        return cls(xyz)

openmc/stats/univariate.py

@@ -7,6 +7,7 @@ from xml.etree import ElementTree as ET
 import numpy as np
 
 import openmc.checkvalue as cv
+from openmc._xml import get_text
 from openmc.mixin import EqualityMixin
 
 
@@ -32,6 +33,29 @@ class Univariate(EqualityMixin, metaclass=ABCMeta):
     def __len__(self):
         return 0
 
+    @classmethod
+    @abstractmethod
+    def from_xml_element(cls, elem):
+        distribution = get_text(elem, 'type')
+        if distribution == 'discrete':
+            return Discrete.from_xml_element(elem)
+        elif distribution == 'uniform':
+            return Uniform.from_xml_element(elem)
+        elif distribution == 'maxwell':
+            return Maxwell.from_xml_element(elem)
+        elif distribution == 'watt':
+            return Watt.from_xml_element(elem)
+        elif distribution == 'normal':
+            return Normal.from_xml_element(elem)
+        elif distribution == 'muir':
+            return Muir.from_xml_element(elem)
+        elif distribution == 'tabular':
+            return Tabular.from_xml_element(elem)
+        elif distribution == 'legendre':
+            return Legendre.from_xml_element(elem)
+        elif distribution == 'mixture':
+            return Mixture.from_xml_element(elem)
+
 
 class Discrete(Univariate):
     """Distribution characterized by a probability mass function.
@@ -110,6 +134,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 get_text(elem, '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 +225,24 @@ 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 = get_text(elem, 'parameters').split()
+        return cls(*map(float, params))
+
 
 class Maxwell(Univariate):
     """Maxwellian distribution in energy.
@@ -237,6 +299,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(get_text(elem, 'parameters'))
+        return cls(theta)
+
 
 class Watt(Univariate):
     r"""Watt fission energy spectrum.
@@ -308,6 +388,25 @@ 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 = get_text(elem, 'parameters').split()
+        return cls(*map(float, params))
+
+
 class Normal(Univariate):
     r"""Normally distributed sampling.
 
@@ -377,6 +476,25 @@ 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 = get_text(elem, 'parameters').split()
+        return cls(*map(float, params))
+
+
 class Muir(Univariate):
     """Muir energy spectrum.
 
@@ -465,6 +583,24 @@ 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 = get_text(elem, 'parameters').split()
+        return cls(*map(float, params))
+
 
 class Tabular(Univariate):
     """Piecewise continuous probability distribution.
@@ -561,6 +697,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 = get_text(elem, 'interpolation')
+        params = [float(x) for x in get_text(elem, 'parameters').split()]
+        x = params[:len(params)//2]
+        p = params[len(params)//2:]
+        return cls(x, p, interpolation)
+
 
 class Legendre(Univariate):
     r"""Probability density given by a Legendre polynomial expansion
@@ -607,6 +764,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 +821,7 @@ class Mixture(Univariate):
 
     def to_xml_element(self, element_name):
         raise NotImplementedError
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        raise NotImplementedError

setup.py

@@ -39,7 +39,13 @@ kwargs = {
     'author': 'The OpenMC Development Team',
     'author_email': 'openmc-dev@googlegroups.com',
     'description': 'OpenMC',
-    'url': 'https://github.com/openmc-dev/openmc',
+    'url': 'https://openmc.org',
+    'download_url': 'https://github.com/openmc-dev/openmc/releases',
+    'project_urls': {
+        'Issue Tracker': 'https://github.com/openmc-dev/openmc/issues',
+        'Documentation': 'https://openmc.readthedocs.io',
+        'Source Code': 'https://github.com/openmc-dev/openmc',
+    },
     'classifiers': [
         'Development Status :: 4 - Beta',
         'Intended Audience :: Developers',
@@ -48,6 +54,7 @@ kwargs = {
         'License :: OSI Approved :: MIT License',
         'Natural Language :: English',
         'Topic :: Scientific/Engineering'
+        'Programming Language :: C++',
         'Programming Language :: Python :: 3',
         'Programming Language :: Python :: 3.4',
         'Programming Language :: Python :: 3.5',
@@ -55,13 +62,12 @@ kwargs = {
         'Programming Language :: Python :: 3.7',
     ],
 
-    # Required dependencies
+    # Dependencies
+    'python_requires': '>=3.4',
     'install_requires': [
         'numpy>=1.9', 'h5py', 'scipy', 'ipython', 'matplotlib',
         'pandas', 'lxml', 'uncertainties'
     ],
-
-    # Optional dependencies
     'extras_require': {
         'test': ['pytest', 'pytest-cov'],
         'vtk': ['vtk'],

src/cell.cpp

@@ -922,6 +922,36 @@ openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance)
   return 0;
 }
 
+extern "C" int
+openmc_cell_get_temperature(int32_t index, const int32_t* instance, double* T)
+{
+  if (index < 0 || index >= model::cells.size()) {
+    strcpy(openmc_err_msg, "Index in cells array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  Cell& c {*model::cells[index]};
+
+  if (c.sqrtkT_.size() < 1) {
+    strcpy(openmc_err_msg, "Cell temperature has not yet been set.");
+    return OPENMC_E_UNASSIGNED;
+  }
+
+  if (instance) {
+    if (*instance >= 0 && *instance < c.n_instances_) {
+      double sqrtkT = c.sqrtkT_.size() == 1 ? c.sqrtkT_[0] : c.sqrtkT_[*instance];
+      *T = sqrtkT * sqrtkT / K_BOLTZMANN;
+    } else {
+      strcpy(openmc_err_msg, "Distribcell instance is out of bounds.");
+      return OPENMC_E_OUT_OF_BOUNDS;
+    }
+  } else {
+    *T = c.sqrtkT_[0] * c.sqrtkT_[0] / K_BOLTZMANN;
+  }
+
+  return 0;
+}
+
 //! Return the index in the cells array of a cell with a given ID
 extern "C" int
 openmc_get_cell_index(int32_t id, int32_t* index)

src/cross_sections.cpp

@@ -239,8 +239,7 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
       already_read.insert(name);
 
       // Check if elemental data has been read, if needed
-      int pos = name.find_first_of("0123456789");
-      std::string element = name.substr(0, pos);
+      std::string element = to_element(name);
       if (settings::photon_transport) {
         if (already_read.find(element) == already_read.end()) {
           // Read photon interaction data from HDF5 photon library

src/eigenvalue.cpp

@@ -29,6 +29,7 @@
 #include <cmath> // for sqrt, abs, pow
 #include <iterator> // for back_inserter
 #include <string>
+#include <limits> //for infinity
 
 namespace openmc {
 
@@ -388,10 +389,15 @@ int openmc_get_keff(double* k_combined)
   k_combined[0] = 0.0;
   k_combined[1] = 0.0;
 
-  // Make sure we have at least four realizations. Notice that at the end,
+  // Special case for n <=3. Notice that at the end,
   // there is a N-3 term in a denominator.
   if (simulation::n_realizations <= 3) {
-    return -1;
+    k_combined[0] = simulation::keff;
+    k_combined[1] = simulation::keff_std;
+    if (simulation::n_realizations <=1) {
+      k_combined[1] = std::numeric_limits<double>::infinity();
+    }
+    return 0;
   }
 
   // Initialize variables

src/material.cpp

@@ -227,8 +227,7 @@ Material::Material(pugi::xml_node node)
     // If the corresponding element hasn't been encountered yet and photon
     // transport will be used, we need to add its symbol to the element_dict
     if (settings::photon_transport) {
-      int pos = name.find_first_of("0123456789");
-      std::string element = name.substr(0, pos);
+      std::string element = to_element(name);
 
       // Make sure photon cross section data is available
       LibraryKey key {Library::Type::photon, element};

src/output.cpp

@@ -617,6 +617,7 @@ const std::unordered_map<int, const char*> score_names = {
   {SCORE_FISS_Q_PROMPT,      "Prompt fission power"},
   {SCORE_FISS_Q_RECOV,       "Recoverable fission power"},
   {SCORE_CURRENT,            "Current"},
+  {SCORE_HEATING,            "Heating"},
 };
 
 //! Create an ASCII output file showing all tally results.

src/particle.cpp

@@ -73,7 +73,7 @@ Particle::clear()
 }
 
 void
-Particle::create_secondary(Direction u, double E, Type type) const
+Particle::create_secondary(Direction u, double E, Type type)
 {
   simulation::secondary_bank.emplace_back();
 
@@ -83,6 +83,8 @@ Particle::create_secondary(Direction u, double E, Type type) const
   bank.r = this->r();
   bank.u = u;
   bank.E = settings::run_CE ? E : g_;
+
+  n_bank_second_ += 1;
 }
 
 void
@@ -157,6 +159,11 @@ Particle::transport()
     u_last_ = this->u();
     r_last_ = this->r();
 
+    // Reset event variables
+    event_ = EVENT_KILL;
+    event_nuclide_ = NUCLIDE_NONE;
+    event_mt_ = REACTION_NONE;
+
     // If the cell hasn't been determined based on the particle's location,
     // initiate a search for the current cell. This generally happens at the
     // beginning of the history and again for any secondary particles
@@ -309,6 +316,7 @@ Particle::transport()
 
       // Reset banked weight during collision
       n_bank_ = 0;
+      n_bank_second_ = 0;
       wgt_bank_ = 0.0;
       for (int& v : n_delayed_bank_) v = 0;
 

src/photon.cpp

@@ -96,6 +96,15 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
   read_dataset(rgroup, "xs", photoelectric_total_);
   close_group(rgroup);
 
+  // Read heating
+  if (object_exists(group, "heating")) {
+    rgroup = open_group(group, "heating");
+    read_dataset(rgroup, "xs", heating_);
+    close_group(rgroup);
+  } else {
+    heating_ = xt::zeros_like(energy_);
+  }
+
   // Read subshell photoionization cross section and atomic relaxation data
   rgroup = open_group(group, "subshells");
   std::vector<std::string> designators;
@@ -280,6 +289,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
     xt::log(photoelectric_total_), -500.0);
   pair_production_total_ = xt::where(pair_production_total_ > 0.0,
     xt::log(pair_production_total_), -500.0);
+  heating_ = xt::where(heating_ > 0.0, xt::log(heating_), -500.0);
 }
 
 void PhotonInteraction::compton_scatter(double alpha, bool doppler,

src/physics.cpp

@@ -17,6 +17,7 @@
 #include "openmc/search.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
+#include "openmc/string_utils.h"
 #include "openmc/thermal.h"
 #include "openmc/tallies/tally.h"
 
@@ -61,12 +62,16 @@ void collision(Particle* p)
   // Display information about collision
   if (settings::verbosity >= 10 || simulation::trace) {
     std::stringstream msg;
-    if (p->type_ == Particle::Type::neutron) {
+    if (p->event_ == EVENT_KILL) {
+      msg << "    Killed. Energy = " << p->E_ << " eV.";
+    } else if (p->type_ == Particle::Type::neutron) {
       msg << "    " << reaction_name(p->event_mt_) << " with " <<
         data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV.";
-    } else {
+    } else if (p->type_ == Particle::Type::photon) {
       msg << "    " << reaction_name(p->event_mt_) << " with " <<
-        data::elements[p->event_nuclide_].name_ << ". Energy = " << p->E_ << " eV.";
+        to_element(data::nuclides[p->event_nuclide_]->name_) << ". Energy = " << p->E_ << " eV.";
+    } else {
+      msg << "    Disappeared. Energy = " << p->E_ << " eV.";
     }
     write_message(msg, 1);
   }
@@ -209,7 +214,6 @@ void sample_photon_reaction(Particle* p)
 
   // Sample element within material
   int i_element = sample_element(p);
-  p->event_nuclide_ = i_element;
   const auto& micro {p->photon_xs_[i_element]};
   const auto& element {data::elements[i_element]};
 
@@ -226,6 +230,7 @@ void sample_photon_reaction(Particle* p)
   if (prob > cutoff) {
     double mu = element.rayleigh_scatter(alpha);
     p->u() = rotate_angle(p->u(), mu, nullptr);
+    p->event_ = EVENT_SCATTER;
     p->event_mt_ = COHERENT;
     return;
   }
@@ -268,6 +273,7 @@ void sample_photon_reaction(Particle* p)
     phi += PI;
     p->E_ = alpha_out*MASS_ELECTRON_EV;
     p->u() = rotate_angle(p->u(), mu, &phi);
+    p->event_ = EVENT_SCATTER;
     p->event_mt_ = INCOHERENT;
     return;
   }
@@ -319,6 +325,7 @@ void sample_photon_reaction(Particle* p)
         // Allow electrons to fill orbital and produce auger electrons
         // and fluorescent photons
         element.atomic_relaxation(shell, *p);
+        p->event_ = EVENT_ABSORB;
         p->event_mt_ = 533 + shell.index_subshell;
         p->alive_ = false;
         p->E_ = 0.0;
@@ -344,6 +351,7 @@ void sample_photon_reaction(Particle* p)
     u = rotate_angle(p->u(), mu_positron, nullptr);
     p->create_secondary(u, E_positron, Particle::Type::positron);
 
+    p->event_ = EVENT_ABSORB;
     p->event_mt_ = PAIR_PROD;
     p->alive_ = false;
     p->E_ = 0.0;
@@ -361,6 +369,7 @@ void sample_electron_reaction(Particle* p)
 
   p->E_ = 0.0;
   p->alive_ = false;
+  p->event_ = EVENT_ABSORB;
 }
 
 void sample_positron_reaction(Particle* p)
@@ -386,6 +395,7 @@ void sample_positron_reaction(Particle* p)
 
   p->E_ = 0.0;
   p->alive_ = false;
+  p->event_ = EVENT_ABSORB;
 }
 
 int sample_nuclide(const Particle* p)
@@ -432,7 +442,12 @@ int sample_element(Particle* p)
 
     // Increment probability to compare to cutoff
     prob += sigma;
-    if (prob > cutoff) return i_element;
+    if (prob > cutoff) {
+      // Save which nuclide particle had collision with for tally purpose
+      p->event_nuclide_ = mat->nuclide_[i];
+
+      return i_element;
+    }
   }
 
   // If we made it here, no element was sampled

src/reaction.cpp

@@ -121,6 +121,8 @@ std::string reaction_name(int mt)
     return  "fission-q-prompt";
   } else if (mt == SCORE_FISS_Q_RECOV) {
     return  "fission-q-recoverable";
+  } else if (mt == SCORE_HEATING) {
+    return  "heating";
 
   // Normal ENDF-based reactions
   } else if (mt == TOTAL_XS) {

src/settings.cpp

@@ -401,7 +401,7 @@ void read_settings_xml()
     SourceDistribution source {
       UPtrSpace{new SpatialPoint({0.0, 0.0, 0.0})},
       UPtrAngle{new Isotropic()},
-      UPtrDist{new Watt(0.988, 2.249e-6)}
+      UPtrDist{new Watt(0.988e6, 2.249e-6)}
     };
     model::external_sources.push_back(std::move(source));
   }

src/string_utils.cpp

@@ -26,6 +26,12 @@ char* strtrim(char* c_str)
 }
 
 
+std::string to_element(const std::string& name) {
+  int pos = name.find_first_of("0123456789");
+  return name.substr(0, pos);
+}
+
+
 void to_lower(std::string& str)
 {
   for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);

src/tallies/tally.cpp

@@ -9,6 +9,7 @@
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/particle.h"
+#include "openmc/reaction.h"
 #include "openmc/reaction_product.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
@@ -109,6 +110,9 @@ score_str_to_int(std::string score_str)
   if (score_str == "fission-q-recoverable")
     return SCORE_FISS_Q_RECOV;
 
+  if (score_str == "heating")
+    return SCORE_HEATING;
+
   if (score_str == "current")
     return SCORE_CURRENT;
 
@@ -207,8 +211,6 @@ score_str_to_int(std::string score_str)
     return N_X3HE;
   if (score_str == "(n,Xa)" || score_str == "He4-production")
     return N_XA;
-  if (score_str == "heating")
-    return HEATING;
   if (score_str == "damage-energy")
     return DAMAGE_ENERGY;
 
@@ -771,7 +773,7 @@ void read_tallies_xml()
           case SCORE_PROMPT_NU_FISSION:
           case SCORE_DECAY_RATE:
             warning("Particle filter is not used with photon transport"
-              " on and " + std::to_string(score) + " score.");
+              " on and " + reaction_name(score) + " score.");
             break;
           }
         }

src/tallies/tally_scoring.cpp

@@ -7,10 +7,12 @@
 #include "openmc/material.h"
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
+#include "openmc/photon.h"
 #include "openmc/reaction_product.h"
 #include "openmc/search.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
+#include "openmc/string_utils.h"
 #include "openmc/tallies/derivative.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/filter_delayedgroup.h"
@@ -1149,6 +1151,156 @@ score_general_ce(Particle* p, int i_tally, int start_index,
       break;
 
 
+    case SCORE_HEATING:
+      score = 0.;
+      if (p->type_ == Particle::Type::neutron) {
+        if (tally.estimator_ == ESTIMATOR_ANALOG) {
+          // All events score to a heating tally bin. We actually use a
+          // collision estimator in place of an analog one since there is no
+          // reaction-wise heating cross section
+          if (settings::survival_biasing) {
+            // We need to account for the fact that some weight was already
+            // absorbed
+            score = p->wgt_last_ + p->wgt_absorb_;
+          } else {
+            score = p->wgt_last_;
+          }
+          if (i_nuclide >= 0) {
+            // Calculate nuclide heating cross section
+            double macro_heating = 0.;
+            const auto& nuc {*data::nuclides[i_nuclide]};
+            auto m = nuc.reaction_index_[NEUTRON_HEATING];
+            if (m == C_NONE) continue;
+            const auto& rxn {*nuc.reactions_[m]};
+            auto i_temp = p->neutron_xs_[i_nuclide].index_temp;
+            if (i_temp >= 0) { // Can be false due to multipole
+              auto i_grid = p->neutron_xs_[i_nuclide].index_grid;
+              auto f = p->neutron_xs_[i_nuclide].interp_factor;
+              const auto& xs {rxn.xs_[i_temp]};
+              if (i_grid >= xs.threshold) {
+                macro_heating = ((1.0 - f) * xs.value[i_grid-xs.threshold]
+                  + f * xs.value[i_grid-xs.threshold+1]);
+              }
+            }
+            score *= macro_heating * flux / p->neutron_xs_[i_nuclide].total;
+          } else {
+            if (p->material_ != MATERIAL_VOID) {
+              // Calculate material heating cross section
+              double macro_heating = 0.;
+              const Material& material {*model::materials[p->material_]};
+              for (auto i = 0; i < material.nuclide_.size(); ++i) {
+                auto j_nuclide = material.nuclide_[i];
+                auto atom_density = material.atom_density_(i);
+                const auto& nuc {*data::nuclides[j_nuclide]};
+                auto m = nuc.reaction_index_[NEUTRON_HEATING];
+                if (m == C_NONE) continue;
+                const auto& rxn {*nuc.reactions_[m]};
+                auto i_temp = p->neutron_xs_[j_nuclide].index_temp;
+                if (i_temp >= 0) { // Can be false due to multipole
+                  auto i_grid = p->neutron_xs_[j_nuclide].index_grid;
+                  auto f = p->neutron_xs_[j_nuclide].interp_factor;
+                  const auto& xs {rxn.xs_[i_temp]};
+                  if (i_grid >= xs.threshold) {
+                    macro_heating += ((1.0 - f) * xs.value[i_grid-xs.threshold]
+                      + f * xs.value[i_grid-xs.threshold+1]) * atom_density;
+                  }
+                }
+              }
+              score *= macro_heating * flux / p->macro_xs_.total;
+            } else {
+              score = 0.;
+            }
+          }
+        } else {
+          // Calculate neutron heating cross section on-the-fly
+          if (i_nuclide >= 0) {
+            const auto& nuc {*data::nuclides[i_nuclide]};
+            auto m = nuc.reaction_index_[NEUTRON_HEATING];
+            if (m == C_NONE) continue;
+            const auto& rxn {*nuc.reactions_[m]};
+            auto i_temp = p->neutron_xs_[i_nuclide].index_temp;
+            if (i_temp >= 0) { // Can be false due to multipole
+              auto i_grid = p->neutron_xs_[i_nuclide].index_grid;
+              auto f = p->neutron_xs_[i_nuclide].interp_factor;
+              const auto& xs {rxn.xs_[i_temp]};
+              if (i_grid >= xs.threshold) {
+                score = ((1.0 - f) * xs.value[i_grid-xs.threshold]
+                  + f * xs.value[i_grid-xs.threshold+1]) * atom_density * flux;
+              }
+            }
+          } else {
+            if (p->material_ != MATERIAL_VOID) {
+              const Material& material {*model::materials[p->material_]};
+              for (auto i = 0; i < material.nuclide_.size(); ++i) {
+                auto j_nuclide = material.nuclide_[i];
+                auto atom_density = material.atom_density_(i);
+                const auto& nuc {*data::nuclides[j_nuclide]};
+                auto m = nuc.reaction_index_[NEUTRON_HEATING];
+                if (m == C_NONE) continue;
+                const auto& rxn {*nuc.reactions_[m]};
+                auto i_temp = p->neutron_xs_[j_nuclide].index_temp;
+                if (i_temp >= 0) { // Can be false due to multipole
+                  auto i_grid = p->neutron_xs_[j_nuclide].index_grid;
+                  auto f = p->neutron_xs_[j_nuclide].interp_factor;
+                  const auto& xs {rxn.xs_[i_temp]};
+                  if (i_grid >= xs.threshold) {
+                    score += ((1.0 - f) * xs.value[i_grid-xs.threshold]
+                      + f * xs.value[i_grid-xs.threshold+1]) * atom_density
+                      * flux;
+                  }
+                }
+              }
+            }
+          }
+        }
+      } else if (p->type_ == Particle::Type::photon) {
+        if (tally.estimator_ == ESTIMATOR_ANALOG) {
+          // Score direct energy deposition in the collision
+          score = E - p->E_;
+          // We need to substract the energy of the secondary particles since
+          // they will be transported individually later
+          for (auto i = 0; i < p->n_bank_second_; ++i) {
+            auto i_bank = simulation::secondary_bank.size() - p->n_bank_second_ + i;
+            const auto& bank = simulation::secondary_bank[i_bank];
+            if (bank.particle == Particle::Type::photon ||
+                bank.particle == Particle::Type::neutron) {
+              score -= bank.E;
+            } else if (bank.particle == Particle::Type::positron) {
+              // Annihilation of the positron will produce two new photons
+              score -= 2*MASS_ELECTRON_EV;
+            }
+          }
+          score *= p->wgt_last_ * flux;
+        } else {
+          // Calculate photon heating cross section on-the-fly
+          if (i_nuclide >= 0) {
+            // Find the element corresponding to the nuclide
+            auto name = data::nuclides[i_nuclide]->name_;
+            std::string element = to_element(name);
+            int i_element = data::element_map[element];
+            auto& heating {data::elements[i_element].heating_};
+            auto i_grid = p->photon_xs_[i_element].index_grid;
+            auto f = p->photon_xs_[i_element].interp_factor;
+            score = std::exp(heating(i_grid) + f * (heating(i_grid+1) -
+              heating(i_grid))) * atom_density * flux;
+          } else {
+            if (p->material_ != MATERIAL_VOID) {
+              const Material& material {*model::materials[p->material_]};
+              for (auto i = 0; i < material.nuclide_.size(); ++i) {
+                auto i_element = material.element_[i];
+                auto atom_density = material.atom_density_(i);
+                auto& heating {data::elements[i_element].heating_};
+                auto i_grid = p->photon_xs_[i_element].index_grid;
+                auto f = p->photon_xs_[i_element].interp_factor;
+                score += std::exp(heating(i_grid) + f * (heating(i_grid+1) -
+                  heating(i_grid))) * atom_density * flux;
+              }
+            }
+          }
+        }
+      }
+      break;
+
     default:
       if (tally.estimator_ == ESTIMATOR_ANALOG) {
         // Any other score is assumed to be a MT number. Thus, we just need
@@ -1979,12 +2131,11 @@ void score_analog_tally_ce(Particle* p)
           auto i_nuclide = tally.nuclides_[i];
 
           // Tally this event in the present nuclide bin if that bin represents
-          // the event nuclide or the total material.  Note that the i_nuclide
-          // and flux arguments for score_general are not used for analog
-          // tallies.
+          // the event nuclide or the total material.  Note that the atomic
+          // density argument for score_general is not used for analog tallies.
           if (i_nuclide == p->event_nuclide_ || i_nuclide == -1)
             score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index,
-              -1, -1., filter_weight);
+              i_nuclide, -1., filter_weight);
         }
 
       } else {

src/thermal.cpp

@@ -199,10 +199,10 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
 
   // Calculate S(a,b) inelastic scattering cross section
   auto& xs = sab.inelastic_sigma_;
-  *inelastic = (1.0 - f) * xs[i_grid] + f * xs[i_grid + 1];
+  *inelastic = xs[i_grid] + f * (xs[i_grid + 1] - xs[i_grid]);
 
   // Check for elastic data
-  if (E < sab.threshold_elastic_) {
+  if (!sab.elastic_e_in_.empty()) {
     // Determine whether elastic scattering is given in the coherent or
     // incoherent approximation. For coherent, the cross section is
     // represented as P/E whereas for incoherent, it is simply P
@@ -231,7 +231,7 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
 
       // Calculate S(a,b) elastic scattering cross section
       auto& xs = sab.elastic_P_;
-      *elastic = (1.0 - f) * xs[i_grid] + f * xs[i_grid + 1];
+      *elastic = xs[i_grid] + f*(xs[i_grid + 1] - xs[i_grid]);
     }
   } else {
     // No elastic data

tests/regression_tests/asymmetric_lattice/results_true.dat

@@ -1 +1 @@
-4b75e203d06d0fc1b4c4dfcb8c180d6f3df8fa2bc44e9775b59bbfd8f7a3785f956f9a9f301526f69c0f8a963ce2e553e0c62c7f6c696656ce1d47b415af6076
+4401f503237c94e9d9cfc9f60e0269d5ae5bb67be3225e18c5510ed08616482964e2962a06268751f66a455fac3ddd5faf91555638dfb56fcd09eee60219edff

tests/regression_tests/cmfd_feed_2g/results_true.dat

@@ -1,5 +1,5 @@
 k-combined:
-1.038883E+00 1.017026E-02
+1.038883E+00 1.017030E-02
 tally 1:
 1.167304E+02
 1.362680E+03
@@ -74,7 +74,7 @@ tally 3:
 9.717439E+01
 4.724254E+02
 8.435691E-01
-3.666151E-02
+3.666152E-02
 6.192881E+01
 1.924016E+02
 0.000000E+00
@@ -82,7 +82,7 @@ tally 3:
 1.796382E-02
 3.190855E-05
 4.343238E+00
-9.514039E-01
+9.514040E-01
 3.504683E+00
 6.157615E-01
 0.000000E+00
@@ -383,7 +383,7 @@ cmfd balance
 1.08402E-03
 1.09177E-03
 5.45977E-04
-4.45554E-04
+4.45555E-04
 4.01147E-04
 3.71025E-04
 3.57715E-04

tests/regression_tests/cmfd_feed_ng/results_true.dat

@@ -1,5 +1,5 @@
 k-combined:
-1.029540E+00 1.765356E-02
+1.029540E+00 1.765357E-02
 tally 1:
 1.144958E+02
 1.311468E+03
@@ -93,7 +93,7 @@ tally 3:
 0.000000E+00
 6.939496E+01
 3.014077E+02
-6.294738E-01
+6.294737E-01
 2.532146E-02
 4.991526E+01
 1.566280E+02
@@ -129,7 +129,7 @@ tally 3:
 0.000000E+00
 6.885265E+01
 2.965191E+02
-6.537576E-01
+6.537575E-01
 2.783733E-02
 4.958341E+01
 1.546656E+02

tests/regression_tests/filter_distribcell/case-1/results_true.dat

@@ -1,5 +1,5 @@
 k-combined:
-0.000000E+00 0.000000E+00
+5.497140E-02          INF
 tally 1:
 1.548980E-02
 2.399339E-04

tests/regression_tests/filter_distribcell/case-2/results_true.dat

@@ -1,5 +1,5 @@
 k-combined:
-0.000000E+00 0.000000E+00
+2.149726E-02          INF
 tally 1:
 7.588170E-03
 5.758032E-05

tests/regression_tests/filter_distribcell/case-3/results_true.dat

@@ -1 +1 @@
-11755ecac8355b5e79384f5d72974e618b6f95500c0a5718c01bb340c5d1c8ceafc33de65b9e94b7e9d78f16ec209f8a0fdf6a531eab5430657688d0125db7ef
+2ee0162762999f71ad2178936509fd9e054928023a9e0e90c078cede3ecf6583267069486adf15f1b02188333d1bfe1fc41b341bc00aab53feddd1395441f1f8

tests/regression_tests/filter_distribcell/case-4/results_true.dat

@@ -1,5 +1,5 @@
 k-combined:
-0.000000E+00 0.000000E+00
+1.121246E-01          INF
 tally 1:
 2.265319E-02
 5.131668E-04

tests/regression_tests/lattice_hex_coincident/inputs_true.dat

@@ -39,7 +39,7 @@
     <nuclide ao="0.0049817" name="U235" />
   </material>
   <material id="14">
-    <density units="atom/b-cm" value="0.087742" />
+    <density units="atom/b-cm" value="0.017742" />
     <nuclide ao="1.0" name="C0" />
     <sab name="c_Graphite" />
   </material>
@@ -67,8 +67,8 @@
 <settings>
   <run_mode>eigenvalue</run_mode>
   <particles>1000</particles>
-  <batches>10</batches>
-  <inactive>5</inactive>
+  <batches>5</batches>
+  <inactive>2</inactive>
   <source strength="1.0">
     <space type="box">
       <parameters>-0.9899494936611666 -0.9899494936611666 0.0 0.9899494936611666 0.9899494936611666 10.0</parameters>

tests/regression_tests/lattice_hex_coincident/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-1.741370E+00 1.384609E-03
+1.910374E+00 2.685991E-02

tests/regression_tests/lattice_hex_coincident/test.py

@@ -1,4 +1,5 @@
-import numpy as np
+from math import sqrt
+
 import openmc
 
 from tests.testing_harness import PyAPITestHarness
@@ -13,7 +14,7 @@ class HexLatticeCoincidentTestHarness(PyAPITestHarness):
         materials.append(fuel_mat)
 
         matrix = openmc.Material()
-        matrix.set_density('atom/b-cm', 8.7742E-02)
+        matrix.set_density('atom/b-cm', 1.7742E-02)
         matrix.add_element('C', 1.0, 'ao')
         matrix.add_s_alpha_beta('c_Graphite')
         materials.append(matrix)
@@ -94,7 +95,7 @@ class HexLatticeCoincidentTestHarness(PyAPITestHarness):
         coolant_univ.add_cells(coolant_channel)
 
         half_width = assembly_pitch # cm
-        edge_length = (2./np.sqrt(3.0)) * half_width
+        edge_length = (2./sqrt(3.0)) * half_width
 
         inf_mat = openmc.Cell()
         inf_mat.fill = matrix
@@ -132,19 +133,19 @@ class HexLatticeCoincidentTestHarness(PyAPITestHarness):
         settings.run_mode = 'eigenvalue'
 
         source = openmc.Source()
-        corner_dist = np.sqrt(2) * pin_rad
+        corner_dist = sqrt(2) * pin_rad
         ll = [-corner_dist, -corner_dist, 0.0]
         ur = [corner_dist, corner_dist, 10.0]
         source.space = openmc.stats.Box(ll, ur)
         source.strength = 1.0
         settings.source = source
         settings.output = {'summary' : False}
-        settings.batches = 10
-        settings.inactive = 5
+        settings.batches = 5
+        settings.inactive = 2
         settings.particles = 1000
         settings.seed = 22
         settings.export_to_xml()
 
 def test_lattice_hex_coincident_surf():
-    harness = HexLatticeCoincidentTestHarness('statepoint.10.h5')
+    harness = HexLatticeCoincidentTestHarness('statepoint.5.h5')
     harness.main()

tests/regression_tests/photon_production/inputs_true.dat

@@ -49,17 +49,20 @@
   </tally>
   <tally id="2">
     <filters>2</filters>
-    <scores>total</scores>
+    <nuclides>Al27 total</nuclides>
+    <scores>total heating</scores>
     <estimator>tracklength</estimator>
   </tally>
   <tally id="3">
     <filters>2</filters>
-    <scores>total</scores>
+    <nuclides>Al27 total</nuclides>
+    <scores>total heating</scores>
     <estimator>collision</estimator>
   </tally>
   <tally id="4">
     <filters>2</filters>
-    <scores>total</scores>
+    <nuclides>Al27 total</nuclides>
+    <scores>total heating</scores>
     <estimator>analog</estimator>
   </tally>
 </tallies>

tests/regression_tests/photon_production/results_true.dat

@@ -2,11 +2,29 @@ tally 1:
 9.403000E-01
 8.841641E-01
 tally 2:
+1.819886E-01
+3.311985E-02
+1.960159E+05
+3.842224E+10
 8.281718E-01
 6.858685E-01
+1.960159E+05
+3.842224E+10
 tally 3:
+1.799069E-01
+3.236650E-02
+1.945225E+05
+3.783901E+10
 8.242000E-01
 6.793056E-01
+1.945225E+05
+3.783901E+10
 tally 4:
+2.308000E-01
+5.326864E-02
+1.988563E+05
+3.954384E+10
 8.242000E-01
 6.793056E-01
+1.988612E+05
+3.954578E+10

tests/regression_tests/photon_production/test.py

@@ -51,20 +51,23 @@ class SourceTestHarness(PyAPITestHarness):
         current_tally = openmc.Tally()
         current_tally.filters = [surface_filter, particle_filter]
         current_tally.scores = ['current']
-        total_tally_tracklength = openmc.Tally()
-        total_tally_tracklength.filters = [particle_filter]
-        total_tally_tracklength.scores = ['total']
-        total_tally_tracklength.estimator = 'tracklength'
-        total_tally_collision = openmc.Tally()
-        total_tally_collision.filters = [particle_filter]
-        total_tally_collision.scores = ['total']
-        total_tally_collision.estimator = 'collision'
-        total_tally_analog = openmc.Tally()
-        total_tally_analog.filters = [particle_filter]
-        total_tally_analog.scores = ['total']
-        total_tally_analog.estimator = 'analog'
-        tallies = openmc.Tallies([current_tally, total_tally_tracklength,
-                                  total_tally_collision, total_tally_analog])
+        tally_tracklength = openmc.Tally()
+        tally_tracklength.filters = [particle_filter]
+        tally_tracklength.scores = ['total', 'heating']
+        tally_tracklength.nuclides = ['Al27', 'total']
+        tally_tracklength.estimator = 'tracklength'
+        tally_collision = openmc.Tally()
+        tally_collision.filters = [particle_filter]
+        tally_collision.scores = ['total', 'heating']
+        tally_collision.nuclides = ['Al27', 'total']
+        tally_collision.estimator = 'collision'
+        tally_analog = openmc.Tally()
+        tally_analog.filters = [particle_filter]
+        tally_analog.scores = ['total', 'heating']
+        tally_analog.nuclides = ['Al27', 'total']
+        tally_analog.estimator = 'analog'
+        tallies = openmc.Tallies([current_tally, tally_tracklength,
+                                  tally_collision, tally_analog])
         tallies.export_to_xml()
 
     def _get_results(self):

tests/regression_tests/salphabeta/inputs_true.dat

@@ -49,7 +49,7 @@
 <?xml version='1.0' encoding='utf-8'?>
 <settings>
   <run_mode>eigenvalue</run_mode>
-  <particles>1000</particles>
+  <particles>400</particles>
   <batches>5</batches>
   <inactive>0</inactive>
   <source strength="1.0">

tests/regression_tests/salphabeta/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-8.403447E-01 2.461538E-02
+8.474822E-01 1.767966E-02

tests/regression_tests/salphabeta/test.py

@@ -66,7 +66,7 @@ def make_model():
     # Settings
     model.settings.batches = 5
     model.settings.inactive = 0
-    model.settings.particles = 1000
+    model.settings.particles = 400
     model.settings.source = openmc.Source(space=openmc.stats.Box(
         [-4, -4, -4], [4, 4, 4]))
 

tests/regression_tests/sourcepoint_batch/results_true.dat

@@ -1,3 +1,3 @@
 k-combined:
-0.000000E+00 0.000000E+00
+2.976389E-01 3.770725E-03
 1.892327E+00 -3.385257E+00 6.702632E-01

tests/regression_tests/triso/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-1.707485E+00 9.795497E-02
+1.701412E+00 3.180881E-02

tests/regression_tests/void/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-1.062505E+00 2.674375E-02
+9.612556E-01 1.990135E-02

tests/unit_tests/test_capi.py

@@ -75,6 +75,14 @@ def test_cell(capi_init):
     assert str(cell) == 'Cell[0]'
 
 
+def test_cell_temperature(capi_init):
+    cell = openmc.capi.cells[1]
+    cell.set_temperature(100.0, 0)
+    assert cell.get_temperature(0) == 100.0
+    cell.set_temperature(200)
+    assert cell.get_temperature() == 200.0
+
+
 def test_new_cell(capi_init):
     with pytest.raises(exc.AllocationError):
         openmc.capi.Cell(1)

tests/unit_tests/test_settings.py

@@ -19,11 +19,12 @@ def test_export_to_xml(run_in_tmpdir):
                      'write': True, 'overwrite': True}
     s.statepoint = {'batches': [50, 150, 500, 1000]}
     s.confidence_intervals = True
-    s.cross_sections = '/path/to/cross_sections.xml'
     s.ptables = True
     s.seed = 17
     s.survival_biasing = True
-    s.cutoff = {'weight': 0.25, 'weight_avg': 0.5, 'energy': 1.0e-5}
+    s.cutoff = {'weight': 0.25, 'weight_avg': 0.5, 'energy_neutron': 1.0e-5,
+                'energy_photon': 1000.0, 'energy_electron': 1.0e-5,
+                'energy_positron': 1.0e-5}
     mesh = openmc.Mesh()
     mesh.lower_left = (-10., -10., -10.)
     mesh.upper_right = (10., 10., 10.)
@@ -47,6 +48,59 @@ def test_export_to_xml(run_in_tmpdir):
         upper_right = (10., 10., 10.))
     s.create_fission_neutrons = True
     s.log_grid_bins = 2000
+    s.photon_transport = False
+    s.electron_treatment = 'led'
+    s.dagmc = False
 
     # Make sure exporting XML works
     s.export_to_xml()
+
+    # Generate settings from XML
+    s = openmc.Settings.from_xml()
+    assert s.run_mode == 'fixed source'
+    assert s.batches == 1000
+    assert s.generations_per_batch == 10
+    assert s.inactive == 100
+    assert s.particles == 1000000
+    assert s.keff_trigger == {'type': 'std_dev', 'threshold': 0.001}
+    assert s.energy_mode == 'continuous-energy'
+    assert s.max_order == 5
+    assert isinstance(s.source[0], openmc.Source)
+    assert isinstance(s.source[0].space, openmc.stats.Point)
+    assert s.output == {'summary': True, 'tallies': False, 'path': 'here'}
+    assert s.verbosity == 7
+    assert s.sourcepoint == {'batches': [50, 150, 500, 1000], 'separate': True,
+                             'write': True, 'overwrite': True}
+    assert s.statepoint == {'batches': [50, 150, 500, 1000]}
+    assert s.confidence_intervals
+    assert s.ptables
+    assert s.seed == 17
+    assert s.survival_biasing
+    assert s.cutoff == {'weight': 0.25, 'weight_avg': 0.5,
+                        'energy_neutron': 1.0e-5, 'energy_photon': 1000.0,
+                        'energy_electron': 1.0e-5, 'energy_positron': 1.0e-5}
+    assert isinstance(s.entropy_mesh, openmc.Mesh)
+    assert s.entropy_mesh.lower_left == [-10., -10., -10.]
+    assert s.entropy_mesh.upper_right == [10., 10., 10.]
+    assert s.entropy_mesh.dimension == [5, 5, 5]
+    assert s.trigger_active
+    assert s.trigger_max_batches == 10000
+    assert s.trigger_batch_interval == 50
+    assert not s.no_reduce
+    assert s.tabular_legendre == {'enable': True, 'num_points': 50}
+    assert s.temperature == {'default': 293.6, 'method': 'interpolation',
+                             'multipole': True, 'range': [200., 1000.]}
+    assert s.trace == [10, 1, 20]
+    assert s.track == [1, 1, 1, 2, 1, 1]
+    assert isinstance(s.ufs_mesh, openmc.Mesh)
+    assert s.ufs_mesh.lower_left == [-10., -10., -10.]
+    assert s.ufs_mesh.upper_right == [10., 10., 10.]
+    assert s.ufs_mesh.dimension == [5, 5, 5]
+    assert s.resonance_scattering == {'enable': True, 'method': 'rvs',
+                                      'energy_min': 1.0, 'energy_max': 1000.0,
+                                      'nuclides': ['U235', 'U238', 'Pu239']}
+    assert s.create_fission_neutrons
+    assert s.log_grid_bins == 2000
+    assert not s.photon_transport
+    assert s.electron_treatment == 'led'
+    assert not s.dagmc

tests/unit_tests/test_source.py

@@ -11,7 +11,6 @@ def test_source():
     assert src.space == space
     assert src.angle == angle
     assert src.energy == energy
-    assert src.strength == 1.0
 
     elem = src.to_xml_element()
     assert 'strength' in elem.attrib
@@ -19,6 +18,13 @@ def test_source():
     assert elem.find('angle') is not None
     assert elem.find('energy') is not None
 
+    src = openmc.Source.from_xml_element(elem)
+    assert isinstance(src.angle, openmc.stats.Isotropic)
+    assert src.space.xyz == [0.0, 0.0, 0.0]
+    assert src.energy.x == [1.0e6]
+    assert src.energy.p == [1.0]
+    assert src.strength == 1.0
+
 
 def test_source_file():
     filename = 'source.h5'

tests/unit_tests/test_stats.py

@@ -10,10 +10,15 @@ def test_discrete():
     x = [0.0, 1.0, 10.0]
     p = [0.3, 0.2, 0.5]
     d = openmc.stats.Discrete(x, p)
+    elem = d.to_xml_element('distribution')
+
+    d = openmc.stats.Discrete.from_xml_element(elem)
     assert d.x == x
     assert d.p == p
     assert len(d) == len(x)
-    d.to_xml_element('distribution')
+
+    d = openmc.stats.Univariate.from_xml_element(elem)
+    assert isinstance(d, openmc.stats.Discrete)
 
     # Single point
     d2 = openmc.stats.Discrete(1e6, 1.0)
@@ -25,6 +30,9 @@ def test_discrete():
 def test_uniform():
     a, b = 10.0, 20.0
     d = openmc.stats.Uniform(a, b)
+    elem = d.to_xml_element('distribution')
+
+    d = openmc.stats.Uniform.from_xml_element(elem)
     assert d.a == a
     assert d.b == b
     assert len(d) == 2
@@ -34,35 +42,39 @@ def test_uniform():
     assert t.p == [1/(b-a), 1/(b-a)]
     assert t.interpolation == 'histogram'
 
-    d.to_xml_element('distribution')
-
 
 def test_maxwell():
     theta = 1.2895e6
     d = openmc.stats.Maxwell(theta)
+    elem = d.to_xml_element('distribution')
+
+    d = openmc.stats.Maxwell.from_xml_element(elem)
     assert d.theta == theta
     assert len(d) == 1
-    d.to_xml_element('distribution')
 
 
 def test_watt():
     a, b = 0.965e6, 2.29e-6
     d = openmc.stats.Watt(a, b)
+    elem = d.to_xml_element('distribution')
+
+    d = openmc.stats.Watt.from_xml_element(elem)
     assert d.a == a
     assert d.b == b
     assert len(d) == 2
-    d.to_xml_element('distribution')
 
 
 def test_tabular():
     x = [0.0, 5.0, 7.0]
     p = [0.1, 0.2, 0.05]
     d = openmc.stats.Tabular(x, p, 'linear-linear')
+    elem = d.to_xml_element('distribution')
+
+    d = openmc.stats.Tabular.from_xml_element(elem)
     assert d.x == x
     assert d.p == p
     assert d.interpolation == 'linear-linear'
     assert len(d) == len(x)
-    d.to_xml_element('distribution')
 
 
 def test_legendre():
@@ -115,6 +127,15 @@ def test_polar_azimuthal():
     assert elem.find('mu') is not None
     assert elem.find('phi') is not None
 
+    d = openmc.stats.PolarAzimuthal.from_xml_element(elem)
+    assert d.mu.x == [1.]
+    assert d.mu.p == [1.]
+    assert d.phi.x == [0.]
+    assert d.phi.p == [1.]
+
+    d = openmc.stats.UnitSphere.from_xml_element(elem)
+    assert isinstance(d, openmc.stats.PolarAzimuthal)
+
 
 def test_isotropic():
     d = openmc.stats.Isotropic()
@@ -122,24 +143,25 @@ def test_isotropic():
     assert elem.tag == 'angle'
     assert elem.attrib['type'] == 'isotropic'
 
+    d = openmc.stats.Isotropic.from_xml_element(elem)
+    assert isinstance(d, openmc.stats.Isotropic)
+
 
 def test_monodirectional():
     d = openmc.stats.Monodirectional((1., 0., 0.))
-    assert d.reference_uvw == pytest.approx((1., 0., 0.))
-
     elem = d.to_xml_element()
     assert elem.tag == 'angle'
     assert elem.attrib['type'] == 'monodirectional'
 
+    d = openmc.stats.Monodirectional.from_xml_element(elem)
+    assert d.reference_uvw == pytest.approx((1., 0., 0.))
+
 
 def test_cartesian():
     x = openmc.stats.Uniform(-10., 10.)
     y = openmc.stats.Uniform(-10., 10.)
     z = openmc.stats.Uniform(0., 20.)
     d = openmc.stats.CartesianIndependent(x, y, z)
-    assert d.x == x
-    assert d.y == y
-    assert d.z == z
 
     elem = d.to_xml_element()
     assert elem.tag == 'space'
@@ -147,55 +169,75 @@ def test_cartesian():
     assert elem.find('x') is not None
     assert elem.find('y') is not None
 
+    d = openmc.stats.CartesianIndependent.from_xml_element(elem)
+    assert d.x == x
+    assert d.y == y
+    assert d.z == z
+
+    d = openmc.stats.Spatial.from_xml_element(elem)
+    assert isinstance(d, openmc.stats.CartesianIndependent)
+
 
 def test_box():
     lower_left = (-10., -10., -10.)
     upper_right = (10., 10., 10.)
     d = openmc.stats.Box(lower_left, upper_right)
-    assert d.lower_left == pytest.approx(lower_left)
-    assert d.upper_right == pytest.approx(upper_right)
-    assert not d.only_fissionable
 
     elem = d.to_xml_element()
     assert elem.tag == 'space'
     assert elem.attrib['type'] == 'box'
     assert elem.find('parameters') is not None
 
+    d = openmc.stats.Box.from_xml_element(elem)
+    assert d.lower_left == pytest.approx(lower_left)
+    assert d.upper_right == pytest.approx(upper_right)
+    assert not d.only_fissionable
+
     # only fissionable parameter
     d2 = openmc.stats.Box(lower_left, upper_right, True)
     assert d2.only_fissionable
     elem = d2.to_xml_element()
     assert elem.attrib['type'] == 'fission'
+    d = openmc.stats.Spatial.from_xml_element(elem)
+    assert isinstance(d, openmc.stats.Box)
 
 
 def test_point():
     p = (-4., 2., 10.)
     d = openmc.stats.Point(p)
-    assert d.xyz == pytest.approx(p)
 
     elem = d.to_xml_element()
     assert elem.tag == 'space'
     assert elem.attrib['type'] == 'point'
     assert elem.find('parameters') is not None
 
+    d = openmc.stats.Point.from_xml_element(elem)
+    assert d.xyz == pytest.approx(p)
+
 def test_normal():
     mean = 10.0
     std_dev = 2.0
     d = openmc.stats.Normal(mean,std_dev)
+
+    elem = d.to_xml_element('distribution')
+    assert elem.attrib['type'] == 'normal'
+
+    d = openmc.stats.Normal.from_xml_element(elem)
     assert d.mean_value == pytest.approx(mean)
     assert d.std_dev == pytest.approx(std_dev)
     assert len(d) == 2
-    elem = d.to_xml_element('distribution')
-    assert elem.attrib['type'] == 'normal'
 
 def test_muir():
     mean = 10.0
     mass = 5.0
     temp = 20000.
     d = openmc.stats.Muir(mean,mass,temp)
+
+    elem = d.to_xml_element('energy')
+    assert elem.attrib['type'] == 'muir'
+
+    d = openmc.stats.Muir.from_xml_element(elem)
     assert d.e0 == pytest.approx(mean)
     assert d.m_rat == pytest.approx(mass)
     assert d.kt == pytest.approx(temp)
     assert len(d) == 3
-    elem = d.to_xml_element('energy')
-    assert elem.attrib['type'] == 'muir'

tools/ci/download-xs.sh

@@ -3,7 +3,7 @@ set -ex
 
 # Download HDF5 data
 if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
-    wget -q -O - https://anl.box.com/shared/static/pzutl4i2717yypv12l78l7fn5nmg6grs.xz | tar -C $HOME -xJ
+    wget -q -O - https://anl.box.com/shared/static/u1g3n8iai0u1n5f6ev3pg2j3ff941bqa.xz | tar -C $HOME -xJ
 fi
 
 # Download ENDF/B-VII.1 distribution

tools/ci/travis-install.py

@@ -48,6 +48,9 @@ def install(omp=False, mpi=False, phdf5=False, dagmc=False):
     if dagmc:
         cmake_cmd.append('-Ddagmc=ON')
 
+    # Build in coverage mode for coverage testing
+    cmake_cmd.append('-Dcoverage=on')
+
     # Build and install
     cmake_cmd.append('..')
     print(' '.join(cmake_cmd))

tools/ci/travis-install.sh

@@ -25,5 +25,8 @@ python tools/ci/travis-install.py
 # Install Python API in editable mode
 pip install -e .[test,vtk]
 
-# For uploading to coveralls
+# For coverage testing of the C++ source files
+pip install cpp-coveralls
+
+# For coverage testing of the Python source files
 pip install coveralls