From f369a8ccac7280d62e9a3c90d01697ceacd8497b Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 20 Jan 2022 09:17:51 -0600
Subject: [PATCH 1/8] Add missing VolumeCalculation.from_xml_element method

---
 openmc/settings.py                |  7 +++++
 openmc/volume.py                  | 47 +++++++++++++++++++++++++++++++
 tests/unit_tests/test_settings.py |  7 +++++
 3 files changed, 61 insertions(+)

diff --git a/openmc/settings.py b/openmc/settings.py
index 1ae1597024..18d1fa4d9a 100644
--- a/openmc/settings.py
+++ b/openmc/settings.py
@@ -1247,6 +1247,12 @@ class Settings:
         for elem in root.findall('source'):
             self.source.append(Source.from_xml_element(elem))
 
+    def _volume_calcs_from_xml_element(self, root):
+        volume_elems = root.findall("volume_calc")
+        if volume_elems:
+            self.volume_calculations = [VolumeCalculation.from_xml_element(elem)
+                                        for elem in volume_elems]
+
     def _output_from_xml_element(self, root):
         elem = root.find('output')
         if elem is not None:
@@ -1585,6 +1591,7 @@ class Settings:
         settings._generations_per_batch_from_xml_element(root)
         settings._keff_trigger_from_xml_element(root)
         settings._source_from_xml_element(root)
+        settings._volume_calcs_from_xml_element(root)
         settings._output_from_xml_element(root)
         settings._statepoint_from_xml_element(root)
         settings._sourcepoint_from_xml_element(root)
diff --git a/openmc/volume.py b/openmc/volume.py
index 9acf58ca3b..67f4b7abe5 100644
--- a/openmc/volume.py
+++ b/openmc/volume.py
@@ -11,6 +11,7 @@ from uncertainties import ufloat
 
 import openmc
 import openmc.checkvalue as cv
+from openmc._xml import get_text
 
 _VERSION_VOLUME = 1
 
@@ -352,3 +353,49 @@ class VolumeCalculation:
             trigger_elem.set("type", self.trigger_type)
             trigger_elem.set("threshold", str(self.threshold))
         return element
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate volume calculation object from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.VolumeCalculation
+            Volume calculation object
+
+        """
+        domain_type = get_text(elem, "domain_type")
+        domain_ids = get_text(elem, "domain_ids").split()
+        ids = [int(x) for x in domain_ids]
+        samples = int(get_text(elem, "samples"))
+        lower_left = get_text(elem, "lower_left").split()
+        lower_left = tuple([float(x) for x in lower_left])
+        upper_right = get_text(elem, "upper_right").split()
+        upper_right = tuple([float(x) for x in upper_right])
+
+        # Instantiate some throw-away domains that are used by the constructor
+        # to assign IDs
+        with warnings.catch_warnings():
+            warnings.simplefilter('ignore', openmc.IDWarning)
+            if domain_type == 'cell':
+                domains = [openmc.Cell(uid) for uid in ids]
+            elif domain_type == 'material':
+                domains = [openmc.Material(uid) for uid in ids]
+            elif domain_type == 'universe':
+                domains = [openmc.Universe(uid) for uid in ids]
+
+        vol = cls(domains, samples, lower_left, upper_right)
+
+        # Check for trigger
+        trigger_elem = elem.find("threshold")
+        if trigger_elem is not None:
+            trigger_type = get_text(trigger_elem, "type")
+            threshold = float(get_text(trigger_elem, "threshold"))
+            vol.set_trigger(threshold, trigger_type)
+
+        return vol
diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py
index 832dc93989..e2c7259e75 100644
--- a/tests/unit_tests/test_settings.py
+++ b/tests/unit_tests/test_settings.py
@@ -112,3 +112,10 @@ def test_export_to_xml(run_in_tmpdir):
     assert not s.photon_transport
     assert s.electron_treatment == 'led'
     assert s.write_initial_source == True
+    assert len(s.volume_calculations) == 1
+    vol = s.volume_calculations[0]
+    assert vol.domain_type == 'cell'
+    assert len(vol.ids) == 1
+    assert vol.samples == 1000
+    assert vol.lower_left == (-10., -10., -10.)
+    assert vol.upper_right == (10., 10., 10.)

From 9df6820da1605ee8b1df3313896bd6c1a80d0466 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 20 Jan 2022 10:33:53 -0600
Subject: [PATCH 2/8] Implement from_xml_element methods for
 cylindrical/spherical mesh. Closes #1951

---
 openmc/mesh.py               | 64 ++++++++++++++++++++++++++++++++----
 openmc/stats/multivariate.py |  1 +
 openmc/weight_windows.py     |  2 +-
 3 files changed, 59 insertions(+), 8 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 2a16be1ad2..f12c5a3716 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -89,7 +89,7 @@ class MeshBase(IDManagerMixin, ABC):
             raise ValueError('Unrecognized mesh type: "' + mesh_type + '"')
 
     @classmethod
-    def from_xml(cls, elem):
+    def from_xml_element(cls, elem):
         """Generates a mesh from an XML element
 
         Parameters
@@ -109,6 +109,10 @@ class MeshBase(IDManagerMixin, ABC):
             return RegularMesh.from_xml_element(elem)
         elif mesh_type == 'rectilinear':
             return RectilinearMesh.from_xml_element(elem)
+        elif mesh_type == 'cylindrical':
+            return CylindricalMesh.from_xml_element(elem)
+        elif mesh_type == 'spherical':
+            return SphericalMesh.from_xml_element(elem)
         elif mesh_type == 'unstructured':
             return UnstructuredMesh.from_xml_element(elem)
         else:
@@ -697,10 +701,10 @@ class CylindricalMesh(MeshBase):
     n_dimension : int
         Number of mesh dimensions (always 3 for a CylindricalMesh).
     r_grid : Iterable of float
-        Mesh boundary points along the r-axis. 
+        Mesh boundary points along the r-axis.
         Requirement is r >= 0.
     phi_grid : Iterable of float
-        Mesh boundary points along the phi-axis. 
+        Mesh boundary points along the phi-axis.
         The default value is [0, 360], i.e. the full phi range.
     z_grid : Iterable of float
         Mesh boundary points along the z-axis.
@@ -822,6 +826,29 @@ class CylindricalMesh(MeshBase):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate a cylindrical mesh from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.CylindricalMesh
+            Cylindrical mesh object
+
+        """
+
+        mesh_id = int(get_text(elem, 'id'))
+        mesh = cls(mesh_id)
+        mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()]
+        mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()]
+        mesh.z_grid = [float(x) for x in get_text(elem, "z_grid").split()]
+        return mesh
+
     def calc_mesh_volumes(self):
         """Return Volumes for every mesh cell
 
@@ -986,6 +1013,29 @@ class SphericalMesh(MeshBase):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate a spherical mesh from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.SphericalMesh
+            Spherical mesh object
+
+        """
+
+        mesh_id = int(get_text(elem, 'id'))
+        mesh = cls(mesh_id)
+        mesh.r_grid = [float(x) for x in get_text(elem, "r_grid").split()]
+        mesh.theta_grid = [float(x) for x in get_text(elem, "theta_grid").split()]
+        mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()]
+        return mesh
+
     def calc_mesh_volumes(self):
         """Return Volumes for every mesh cell
 
@@ -1003,8 +1053,6 @@ class SphericalMesh(MeshBase):
         return np.multiply.outer(np.outer(V_r, V_t), V_p)
 
 
-
-
 class UnstructuredMesh(MeshBase):
     """A 3D unstructured mesh
 
@@ -1017,6 +1065,8 @@ class UnstructuredMesh(MeshBase):
     ----------
     filename : str
         Location of the unstructured mesh file
+    library : {'moab', 'libmesh'}
+        Mesh library used for the unstructured mesh tally
     mesh_id : int
         Unique identifier for the mesh
     name : str
@@ -1034,7 +1084,7 @@ class UnstructuredMesh(MeshBase):
         Name of the file containing the unstructured mesh
     length_multiplier: float
         Multiplicative factor to apply to mesh coordinates
-    library : str
+    library : {'moab', 'libmesh'}
         Mesh library used for the unstructured mesh tally
     output : bool
         Indicates whether or not automatic tally output should
@@ -1291,4 +1341,4 @@ class UnstructuredMesh(MeshBase):
         library = get_text(elem, 'library')
         length_multiplier = float(get_text(elem, 'length_multiplier', 1.0))
 
-        return cls(filename, library, mesh_id, '', length_multiplier)
\ No newline at end of file
+        return cls(filename, library, mesh_id, '', length_multiplier)
diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py
index f33427c34d..8635bcd71c 100644
--- a/openmc/stats/multivariate.py
+++ b/openmc/stats/multivariate.py
@@ -492,6 +492,7 @@ class SphericalIndependent(Spatial):
         origin = [float(x) for x in elem.get('origin').split()]
         return cls(r, theta, phi, origin=origin)
 
+
 class CylindricalIndependent(Spatial):
     r"""Spatial distribution represented in cylindrical coordinates.
 
diff --git a/openmc/weight_windows.py b/openmc/weight_windows.py
index ee8a4a2bf8..cb9fc51b81 100644
--- a/openmc/weight_windows.py
+++ b/openmc/weight_windows.py
@@ -270,7 +270,7 @@ class WeightWindows(IDManagerMixin):
         path = f"./mesh[@id='{mesh_id}']"
         mesh_elem = root.find(path)
         if mesh_elem is not None:
-            mesh = MeshBase.from_xml(mesh_elem)
+            mesh = MeshBase.from_xml_element(mesh_elem)
 
         # Read all other parameters
         lower_ww_bounds = [float(l) for l in get_text(elem, 'lower_ww_bounds').split()]

From 5dfc880d91e0596ab6c8e5625bf77ab0dc891734 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 20 Jan 2022 11:45:28 -0600
Subject: [PATCH 3/8] Use radians in CylindricalMesh and SphericalMesh

---
 openmc/mesh.py                                | 39 ++++++++++---------
 src/mesh.cpp                                  | 33 ++++++----------
 .../filter_mesh/inputs_true.dat               |  6 +--
 tests/regression_tests/filter_mesh/test.py    |  7 ++--
 4 files changed, 39 insertions(+), 46 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index f12c5a3716..4c973352b9 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -1,5 +1,6 @@
 from abc import ABC
 from collections.abc import Iterable
+from math import pi
 from numbers import Real, Integral
 import warnings
 from xml.etree import ElementTree as ET
@@ -705,7 +706,7 @@ class CylindricalMesh(MeshBase):
         Requirement is r >= 0.
     phi_grid : Iterable of float
         Mesh boundary points along the phi-axis.
-        The default value is [0, 360], i.e. the full phi range.
+        The default value is [0, 2π], i.e. the full phi range.
     z_grid : Iterable of float
         Mesh boundary points along the z-axis.
     indices : Iterable of tuple
@@ -718,7 +719,7 @@ class CylindricalMesh(MeshBase):
         super().__init__(mesh_id, name)
 
         self._r_grid = None
-        self._phi_grid = [0, 360]
+        self._phi_grid = [0.0, 2*pi]
         self._z_grid = None
 
     @property
@@ -761,12 +762,12 @@ class CylindricalMesh(MeshBase):
     @phi_grid.setter
     def phi_grid(self, grid):
         cv.check_type('mesh phi_grid', grid, Iterable, Real)
-        self._phi_grid = np.array(grid)
+        self._phi_grid = np.asarray(grid)
 
     @z_grid.setter
     def z_grid(self, grid):
         cv.check_type('mesh z_grid', grid, Iterable, Real)
-        self._z_grid = np.array(grid)
+        self._z_grid = np.asarray(grid)
 
     def __repr__(self):
         fmt = '{0: <16}{1}{2}\n'
@@ -796,7 +797,7 @@ class CylindricalMesh(MeshBase):
         # Read and assign mesh properties
         mesh = cls(mesh_id)
         mesh.r_grid = group['r_grid'][()]
-        mesh.phi_grid = 180 / np.pi * group['phi_grid'][()]
+        mesh.phi_grid = group['phi_grid'][()]
         mesh.z_grid = group['z_grid'][()]
 
         return mesh
@@ -859,9 +860,9 @@ class CylindricalMesh(MeshBase):
 
         """
 
-        V_r = np.diff(np.array(self.r_grid)**2 / 2)
-        V_p = np.diff(np.array(self.phi_grid) * np.pi / 180.0)
-        V_z = np.diff(np.array(self.z_grid))
+        V_r = np.diff(np.asarray(self.r_grid)**2 / 2)
+        V_p = np.diff(self.phi_grid)
+        V_z = np.diff(self.z_grid)
 
         return np.multiply.outer(np.outer(V_r, V_p), V_z)
 
@@ -891,10 +892,10 @@ class SphericalMesh(MeshBase):
         Requirement is r >= 0.
     theta_grid : Iterable of float
         Mesh boundary points along the theta-axis in degrees.
-        The default value is [0, 180], i.e. the full theta range.
+        The default value is [0, π], i.e. the full theta range.
     phi_grid : Iterable of float
         Mesh boundary points along the phi-axis in degrees.
-        The default value is [0, 360], i.e. the full phi range.
+        The default value is [0, 2π], i.e. the full phi range.
     indices : Iterable of tuple
         An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
         (2, 1, 1), ...]
@@ -905,8 +906,8 @@ class SphericalMesh(MeshBase):
         super().__init__(mesh_id, name)
 
         self._r_grid = None
-        self._theta_grid = [0, 180]
-        self._phi_grid = [0, 360]
+        self._theta_grid = [0, pi]
+        self._phi_grid = [0, 2*pi]
 
     @property
     def dimension(self):
@@ -948,12 +949,12 @@ class SphericalMesh(MeshBase):
     @theta_grid.setter
     def theta_grid(self, grid):
         cv.check_type('mesh theta_grid', grid, Iterable, Real)
-        self._theta_grid = np.array(grid)
+        self._theta_grid = np.asarray(grid)
 
     @phi_grid.setter
     def phi_grid(self, grid):
         cv.check_type('mesh phi_grid', grid, Iterable, Real)
-        self._phi_grid = np.array(grid)
+        self._phi_grid = np.asarray(grid)
 
     def __repr__(self):
         fmt = '{0: <16}{1}{2}\n'
@@ -983,8 +984,8 @@ class SphericalMesh(MeshBase):
         # Read and assign mesh properties
         mesh = cls(mesh_id)
         mesh.r_grid = group['r_grid'][()]
-        mesh.theta_grid = 180 / np.pi * group['theta_grid'][()]
-        mesh.phi_grid = 180 / np.pi * group['phi_grid'][()]
+        mesh.theta_grid = group['theta_grid'][()]
+        mesh.phi_grid = group['phi_grid'][()]
 
         return mesh
 
@@ -1046,9 +1047,9 @@ class SphericalMesh(MeshBase):
 
         """
 
-        V_r = np.diff(np.array(self.r_grid)**3 / 3)
-        V_t = np.diff(-np.cos(np.pi * np.array(self.theta_grid) / 180.0))
-        V_p = np.diff(np.array(self.phi_grid) * np.pi / 180)
+        V_r = np.diff(np.asarray(self.r_grid)**3 / 3)
+        V_t = np.diff(-np.cos(self.theta_grid))
+        V_p = np.diff(self.phi_grid)
 
         return np.multiply.outer(np.outer(V_r, V_t), V_p)
 
diff --git a/src/mesh.cpp b/src/mesh.cpp
index fa97318645..9d5da29e9b 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -1001,7 +1001,7 @@ double CylindricalMesh::find_r_crossing(
 double CylindricalMesh::find_phi_crossing(
   const Position& r, const Direction& u, double l, int shell) const
 {
-  // Phi grid is [0, 360], thus there is no real surface to cross
+  // Phi grid is [0, 2pi], thus there is no real surface to cross
   if (full_phi_ && (shape_[1] == 1))
     return INFTY;
 
@@ -1104,18 +1104,14 @@ int CylindricalMesh::set_grid()
                "cylindrical meshes must start at phi >= 0.");
     return OPENMC_E_INVALID_ARGUMENT;
   }
-  if (grid_[1].back() > 360) {
+  if (grid_[1].back() > 2.0 * PI) {
     set_errmsg("phi-grids for "
-               "cylindrical meshes must end with theta <= 360 degree.");
+               "cylindrical meshes must end with theta <= 2*pi.");
 
     return OPENMC_E_INVALID_ARGUMENT;
   }
 
-  full_phi_ = (grid_[1].front() == 0.0) && (grid_[1].back() == 360.0);
-
-  // Transform phi-grid from degrees to radians
-  std::transform(grid_[1].begin(), grid_[1].end(), grid_[1].begin(),
-    [](double d) { return M_PI * d / 180.0; });
+  full_phi_ = (grid_[1].front() == 0.0) && (grid_[1].back() == 2.0 * PI);
 
   lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()};
   upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()};
@@ -1225,7 +1221,7 @@ double SphericalMesh::find_r_crossing(
 double SphericalMesh::find_theta_crossing(
   const Position& r, const Direction& u, double l, int shell) const
 {
-  // Theta grid is [0, 180], thus there is no real surface to cross
+  // Theta grid is [0, pi], thus there is no real surface to cross
   if (full_theta_ && (shape_[1] == 1))
     return INFTY;
 
@@ -1287,7 +1283,7 @@ double SphericalMesh::find_theta_crossing(
 double SphericalMesh::find_phi_crossing(
   const Position& r, const Direction& u, double l, int shell) const
 {
-  // Phi grid is [0, 360], thus there is no real surface to cross
+  // Phi grid is [0, 2pi], thus there is no real surface to cross
   if (full_phi_ && (shape_[2] == 1))
     return INFTY;
 
@@ -1368,25 +1364,20 @@ int SphericalMesh::set_grid()
       return OPENMC_E_INVALID_ARGUMENT;
     }
   }
-  if (grid_[1].back() > 180) {
+  if (grid_[1].back() > PI) {
     set_errmsg("theta-grids for "
-               "spherical meshes must end with theta <= 180 degree.");
+               "spherical meshes must end with theta <= pi.");
 
     return OPENMC_E_INVALID_ARGUMENT;
   }
-  if (grid_[2].back() > 360) {
+  if (grid_[2].back() > 2 * PI) {
     set_errmsg("phi-grids for "
-               "spherical meshes must end with phi <= 180 degree.");
+               "spherical meshes must end with phi <= 2*pi.");
     return OPENMC_E_INVALID_ARGUMENT;
   }
 
-  full_theta_ = (grid_[1].front() == 0.0) && (grid_[1].back() == 180.0);
-  full_phi_ = (grid_[2].front() == 0.0) && (grid_[2].back() == 360.0);
-
-  // Transform theta- and phi-grid from degrees to radians
-  for (int i = 1; i < 3; i++)
-    std::transform(grid_[i].begin(), grid_[i].end(), grid_[i].begin(),
-      [](double d) { return M_PI * d / 180.0; });
+  full_theta_ = (grid_[1].front() == 0.0) && (grid_[1].back() == PI);
+  full_phi_ = (grid_[2].front() == 0.0) && (grid_[2].back() == 2 * PI);
 
   lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()};
   upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()};
diff --git a/tests/regression_tests/filter_mesh/inputs_true.dat b/tests/regression_tests/filter_mesh/inputs_true.dat
index e188768517..d0c481d604 100644
--- a/tests/regression_tests/filter_mesh/inputs_true.dat
+++ b/tests/regression_tests/filter_mesh/inputs_true.dat
@@ -55,13 +55,13 @@
   </mesh>
   <mesh id="5" type="cylindrical">
     <r_grid>0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5</r_grid>
-    <phi_grid>0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0</phi_grid>
+    <phi_grid>0.0 0.3490658503988659 0.6981317007977318 1.0471975511965976 1.3962634015954636 1.7453292519943295 2.0943951023931953 2.443460952792061 2.792526803190927 3.141592653589793 3.490658503988659 3.839724354387525 4.1887902047863905 4.537856055185257 4.886921905584122 5.235987755982989 5.585053606381854 5.93411945678072 6.283185307179586</phi_grid>
     <z_grid>-7.5 -6.5625 -5.625 -4.6875 -3.75 -2.8125 -1.875 -0.9375 0.0 0.9375 1.875 2.8125 3.75 4.6875 5.625 6.5625 7.5</z_grid>
   </mesh>
   <mesh id="6" type="spherical">
     <r_grid>0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5</r_grid>
-    <theta_grid>0.0 22.5 45.0 67.5 90.0 112.5 135.0 157.5 180.0</theta_grid>
-    <phi_grid>0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0</phi_grid>
+    <theta_grid>0.0 0.39269908169872414 0.7853981633974483 1.1780972450961724 1.5707963267948966 1.9634954084936207 2.356194490192345 2.748893571891069 3.141592653589793</theta_grid>
+    <phi_grid>0.0 0.3490658503988659 0.6981317007977318 1.0471975511965976 1.3962634015954636 1.7453292519943295 2.0943951023931953 2.443460952792061 2.792526803190927 3.141592653589793 3.490658503988659 3.839724354387525 4.1887902047863905 4.537856055185257 4.886921905584122 5.235987755982989 5.585053606381854 5.93411945678072 6.283185307179586</phi_grid>
   </mesh>
   <filter id="1" type="mesh">
     <bins>1</bins>
diff --git a/tests/regression_tests/filter_mesh/test.py b/tests/regression_tests/filter_mesh/test.py
index eddf0a9913..fc24b1b9bc 100644
--- a/tests/regression_tests/filter_mesh/test.py
+++ b/tests/regression_tests/filter_mesh/test.py
@@ -1,4 +1,5 @@
 import numpy as np
+from math import pi
 
 import openmc
 import pytest
@@ -53,13 +54,13 @@ def model():
 
     cyl_mesh = openmc.CylindricalMesh()
     cyl_mesh.r_grid = np.linspace(0, 7.5, 18)
-    cyl_mesh.phi_grid = np.linspace(0, 360, 19)
+    cyl_mesh.phi_grid = np.linspace(0, 2*pi, 19)
     cyl_mesh.z_grid = np.linspace(-7.5, 7.5, 17)
 
     sph_mesh = openmc.SphericalMesh()
     sph_mesh.r_grid = np.linspace(0, 7.5, 18)
-    sph_mesh.theta_grid = np.linspace(0, 180, 9)
-    sph_mesh.phi_grid = np.linspace(0, 360, 19)
+    sph_mesh.theta_grid = np.linspace(0, pi, 9)
+    sph_mesh.phi_grid = np.linspace(0, 2*pi, 19)
 
     # Create filters
     reg_filters = [

From 6bd33602cc1f81564227497759a71d939edcfde7 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 20 Jan 2022 17:51:16 -0600
Subject: [PATCH 4/8] Implement Tallies.from_xml classmethod

---
 openmc/filter.py           | 102 +++++++++++++++++++++++---------
 openmc/filter_expansion.py |  40 +++++++++++--
 openmc/model/model.py      |  13 +++--
 openmc/tallies.py          | 115 ++++++++++++++++++++++++++++++++++---
 openmc/tally_derivative.py |  21 +++++++
 openmc/trigger.py          |  38 ++++++++++--
 6 files changed, 281 insertions(+), 48 deletions(-)

diff --git a/openmc/filter.py b/openmc/filter.py
index 699fa1aef9..fbbc855c28 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -16,6 +16,7 @@ from .material import Material
 from .mixin import IDManagerMixin
 from .surface import Surface
 from .universe import UniverseBase
+from ._xml import get_text
 
 
 _FILTER_TYPES = (
@@ -231,9 +232,43 @@ class Filter(IDManagerMixin, metaclass=FilterMeta):
 
         subelement = ET.SubElement(element, 'bins')
         subelement.text = ' '.join(str(b) for b in self.bins)
-
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        """Generate a filter from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+        **kwargs
+            Keyword arguments (e.g., mesh information)
+
+        Returns
+        -------
+        openmc.Filter
+            Filter object
+
+        """
+        filter_type = elem.get('type')
+
+        # If the filter type matches this class's short_name, then
+        # there is no overriden from_xml_element method.
+        if filter_type == cls.short_name.lower():
+            # Get bins from element -- the default here works for any filters
+            # that just store a list of bins that can be represented as integers
+            filter_id = int(elem.get('id'))
+            bins = [int(x) for x in get_text(elem, 'bins').split()]
+            return cls(bins, filter_id=filter_id)
+
+        # Search through all subclasses and find the one matching the HDF5
+        # 'type'.  Call that class's from_hdf5 method.
+        for subclass in cls._recursive_subclasses():
+            if filter_type == subclass.short_name.lower():
+                return subclass.from_xml_element(elem, **kwargs)
+
+
     def can_merge(self, other):
         """Determine if filter can be merged with another.
 
@@ -622,6 +657,13 @@ class CellInstanceFilter(Filter):
         subelement.text = ' '.join(str(i) for i in self.bins.ravel())
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        filter_id = int(elem.get('id'))
+        bins = [int(x) for x in get_text(elem, 'bins').split()]
+        cell_instances = list(zip(bins[::2], bins[1::2]))
+        return cls(cell_instances, filter_id=filter_id)
+
 
 class SurfaceFilter(WithIDFilter):
     """Filters particles by surface crossing
@@ -661,8 +703,8 @@ class ParticleFilter(Filter):
 
     Attributes
     ----------
-    bins : Iterable of Integral
-        The Particles to tally
+    bins : iterable of str
+        The particles to tally
     id : int
         Unique identifier for the filter
     num_bins : Integral
@@ -698,6 +740,12 @@ class ParticleFilter(Filter):
         filter_id = int(group.name.split('/')[-1].lstrip('filter '))
         return cls(particles, filter_id=filter_id)
 
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        filter_id = int(elem.get('id'))
+        bins = get_text(elem, 'bins').split()
+        return cls(bins, filter_id=filter_id)
+
 
 class MeshFilter(Filter):
     """Bins tally event locations onto a regular, rectangular mesh.
@@ -877,6 +925,18 @@ class MeshFilter(Filter):
             element.set('translation', ' '.join(map(str, self.translation)))
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        mesh_id = int(get_text(elem, 'bins'))
+        mesh_obj = kwargs['meshes'][mesh_id]
+        filter_id = int(elem.get('id'))
+        out = cls(mesh_obj, filter_id=filter_id)
+
+        translation = elem.get('translation')
+        if translation:
+            out.translation = [float(x) for x in translation.split()]
+        return out
+
 
 class MeshSurfaceFilter(MeshFilter):
     """Filter events by surface crossings on a regular, rectangular mesh.
@@ -1019,35 +1079,12 @@ class CollisionFilter(Filter):
         self.bins = np.asarray(bins)
         self.id = filter_id
 
-    def __repr__(self):
-        string = type(self).__name__ + '\n'
-        string += '{: <16}=\t{}\n'.format('\tValues', self.bins)
-        string += '{: <16}=\t{}\n'.format('\tID', self.id)
-        return string
-
-    @Filter.bins.setter
-    def bins(self, bins):
-        Filter.bins.__set__(self, np.asarray(bins))
-
     def check_bins(self, bins):
         for x in bins:
             # Values should be integers
             cv.check_type('filter value', x, Integral)
             cv.check_greater_than('filter value', x, 0, equality=True)
 
-    def to_xml_element(self):
-        """Return XML Element representing the Filter.
-
-        Returns
-        -------
-        element : xml.etree.ElementTree.Element
-            XML element containing filter data
-
-        """
-        element = super().to_xml_element()
-        element[0].text = ' '.join(str(x) for x in self.bins)
-        return element
-
 
 class RealFilter(Filter):
     """Tally modifier that describes phase-space and other characteristics
@@ -1236,6 +1273,12 @@ class RealFilter(Filter):
         element[0].text = ' '.join(str(x) for x in self.values)
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        filter_id = int(elem.get('id'))
+        bins = [float(x) for x in get_text(elem, 'bins').split()]
+        return cls(bins, filter_id=filter_id)
+
 
 class EnergyFilter(RealFilter):
     """Bins tally events based on incident particle energy.
@@ -1969,6 +2012,13 @@ class EnergyFunctionFilter(Filter):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        filter_id = int(elem.get('id'))
+        energy = [float(x) for x in get_text(elem, 'energy').split()]
+        y = [float(x) for x in get_text(elem, 'y').split()]
+        return cls(energy, y, filter_id=filter_id)
+
     def can_merge(self, other):
         return False
 
diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py
index 9a915d18fa..1c07f58b71 100644
--- a/openmc/filter_expansion.py
+++ b/openmc/filter_expansion.py
@@ -46,6 +46,12 @@ class ExpansionFilter(Filter):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        filter_id = int(elem.get('id'))
+        order = int(elem.find('order').text)
+        return cls(order, filter_id=filter_id)
+
 
 class LegendreFilter(ExpansionFilter):
     r"""Score Legendre expansion moments up to specified order.
@@ -226,6 +232,15 @@ class SpatialLegendreFilter(ExpansionFilter):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        filter_id = int(elem.get('id'))
+        order = int(elem.find('order').text)
+        axis = elem.find('axis').text
+        minimum = float(elem.find('min').text)
+        maximum = float(elem.find('max').text)
+        return cls(order, axis, minimum, maximum, filter_id=filter_id)
+
 
 class SphericalHarmonicsFilter(ExpansionFilter):
     r"""Score spherical harmonic expansion moments up to specified order.
@@ -316,6 +331,14 @@ class SphericalHarmonicsFilter(ExpansionFilter):
         element.set('cosine', self.cosine)
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        filter_id = int(elem.get('id'))
+        order = int(elem.find('order').text)
+        filter = cls(order, filter_id=filter_id)
+        filter.cosine = elem.get('cosine')
+        return filter
+
 
 class ZernikeFilter(ExpansionFilter):
     r"""Score Zernike expansion moments in space up to specified order.
@@ -358,7 +381,7 @@ class ZernikeFilter(ExpansionFilter):
         x-coordinate of center of circle for normalization
     y : float
         y-coordinate of center of circle for normalization
-    r : int or None
+    r : float
         Radius of circle for normalization
 
     Attributes
@@ -369,7 +392,7 @@ class ZernikeFilter(ExpansionFilter):
         x-coordinate of center of circle for normalization
     y : float
         y-coordinate of center of circle for normalization
-    r : int or None
+    r : float
         Radius of circle for normalization
     id : int
         Unique identifier for the filter
@@ -464,6 +487,15 @@ class ZernikeFilter(ExpansionFilter):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        filter_id = int(elem.get('id'))
+        order = int(elem.find('order').text)
+        x = float(elem.find('x').text)
+        y = float(elem.find('y').text)
+        r = float(elem.find('r').text)
+        return cls(order, x, y, r, filter_id=filter_id)
+
 
 class ZernikeRadialFilter(ZernikeFilter):
     r"""Score the :math:`m = 0` (radial variation only) Zernike moments up to
@@ -499,7 +531,7 @@ class ZernikeRadialFilter(ZernikeFilter):
         x-coordinate of center of circle for normalization
     y : float
         y-coordinate of center of circle for normalization
-    r : int or None
+    r : float
         Radius of circle for normalization
 
     Attributes
@@ -510,7 +542,7 @@ class ZernikeRadialFilter(ZernikeFilter):
         x-coordinate of center of circle for normalization
     y : float
         y-coordinate of center of circle for normalization
-    r : int or None
+    r : float
         Radius of circle for normalization
     id : int
         Unique identifier for the filter
diff --git a/openmc/model/model.py b/openmc/model/model.py
index c03c151c6f..5c9fefda7d 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -204,10 +204,10 @@ class Model:
 
     @classmethod
     def from_xml(cls, geometry='geometry.xml', materials='materials.xml',
-                 settings='settings.xml'):
+                 settings='settings.xml', tallies='tallies.xml'):
         """Create model from existing XML files
-        When initializing this way, the user must manually load plots and
-        tallies.
+
+        When initializing this way, the user must manually load plots.
 
         Parameters
         ----------
@@ -217,6 +217,10 @@ class Model:
             Path to materials.xml file
         settings : str
             Path to settings.xml file
+        tallies : str
+            Path to tallies.xml file
+
+            .. versionadded:: 0.13.0
 
         Returns
         -------
@@ -227,7 +231,8 @@ class Model:
         materials = openmc.Materials.from_xml(materials)
         geometry = openmc.Geometry.from_xml(geometry, materials)
         settings = openmc.Settings.from_xml(settings)
-        return cls(geometry, materials, settings)
+        tallies = openmc.Tallies.from_xml(tallies)
+        return cls(geometry, materials, settings, tallies)
 
     def init_lib(self, threads=None, geometry_debug=False, restart_file=None,
                  tracks=False, output=True, event_based=None, intracomm=None):
diff --git a/openmc/tallies.py b/openmc/tallies.py
index d99eb41ff1..6755e805d0 100644
--- a/openmc/tallies.py
+++ b/openmc/tallies.py
@@ -10,12 +10,14 @@ from xml.etree import ElementTree as ET
 import h5py
 import numpy as np
 import pandas as pd
+from scipy.misc import derivative
 import scipy.sparse as sps
 
 import openmc
 import openmc.checkvalue as cv
-from ._xml import clean_indentation, reorder_attributes
+from ._xml import clean_indentation, reorder_attributes, get_text
 from .mixin import IDManagerMixin
+from .mesh import MeshBase
 
 
 # The tally arithmetic product types. The tensor product performs the full
@@ -844,13 +846,8 @@ class Tally(IDManagerMixin):
                   'not contain any scores'
             raise ValueError(msg)
 
-        else:
-            scores = ''
-            for score in self.scores:
-                scores += f'{score} '
-
-            subelement = ET.SubElement(element, "scores")
-            subelement.text = scores.rstrip(' ')
+        subelement = ET.SubElement(element, "scores")
+        subelement.text = ' '.join(str(x) for x in self.scores)
 
         # Tally estimator type
         if self.estimator is not None:
@@ -859,7 +856,7 @@ class Tally(IDManagerMixin):
 
         # Optional Triggers
         for trigger in self.triggers:
-            trigger.get_trigger_xml(element)
+            element.append(trigger.to_xml_element())
 
         # Optional derivatives
         if self.derivative is not None:
@@ -868,6 +865,60 @@ class Tally(IDManagerMixin):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem, **kwargs):
+        """Generate tally object from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.Tally
+            Tally object
+
+        """
+        tally_id = int(elem.get('id'))
+        name = elem.get('name', '')
+        tally = cls(tally_id=tally_id, name=name)
+
+        # Read filters
+        filters_elem = elem.find('filters')
+        if filters_elem is not None:
+            filter_ids = [int(x) for x in filters_elem.text.split()]
+            tally.filters = [kwargs['filters'][uid] for uid in filter_ids]
+
+        # Read nuclides
+        nuclides_elem = elem.find('nuclides')
+        if nuclides_elem is not None:
+            tally.nuclides = nuclides_elem.text.split()
+
+        # Read scores
+        scores_elem = elem.find('scores')
+        if scores_elem is not None:
+            tally.scores = scores_elem.text.split()
+
+        # Set estimator
+        estimator_elem = elem.find('estimator')
+        if estimator_elem is not None:
+            tally.estimator = estimator_elem.text
+
+        # Read triggers
+        tally.triggers = [
+            openmc.Trigger.from_xml_element(trigger_elem)
+            for trigger_elem in elem.findall('trigger')
+        ]
+
+        # Read tally derivative
+        deriv_elem = elem.find('derivative')
+        if deriv_elem is not None:
+            deriv_id = int(deriv_elem.text)
+            tally.derivative = kwargs['derivatives'][deriv_id]
+
+        return tally
+
     def contains_filter(self, filter_type):
         """Looks for a filter in the tally that matches a specified type
 
@@ -3143,3 +3194,49 @@ class Tallies(cv.CheckedList):
         reorder_attributes(root_element)  # TODO: Remove when support is Python 3.8+
         tree = ET.ElementTree(root_element)
         tree.write(str(p), xml_declaration=True, encoding='utf-8')
+
+    @classmethod
+    def from_xml(cls, path='tallies.xml'):
+        """Generate tallies from XML file
+
+        Parameters
+        ----------
+        path : str, optional
+            Path to tallies XML file
+
+        Returns
+        -------
+        openmc.Tallies
+            Tallies object
+
+        """
+        tree = ET.parse(path)
+        root = tree.getroot()
+
+        # Read mesh elements
+        meshes = {}
+        for elem in root.findall('mesh'):
+            mesh = MeshBase.from_xml_element(elem)
+            meshes[mesh.id] = mesh
+
+        # Read filter elements
+        filters = {}
+        for elem in root.findall('filter'):
+            filter = openmc.Filter.from_xml_element(elem, meshes=meshes)
+            filters[filter.id] = filter
+
+        # Read derivative elements
+        derivatives = {}
+        for elem in root.findall('derivative'):
+            deriv = openmc.TallyDerivative.from_xml_element(elem)
+            derivatives[deriv.id] = deriv
+
+        # Read tally elements
+        tallies = []
+        for elem in root.findall('tally'):
+            tally = openmc.Tally.from_xml_element(
+                elem, filters=filters, derivatives=derivatives
+            )
+            tallies.append(tally)
+
+        return cls(tallies)
diff --git a/openmc/tally_derivative.py b/openmc/tally_derivative.py
index 125946197b..05a27681d5 100644
--- a/openmc/tally_derivative.py
+++ b/openmc/tally_derivative.py
@@ -105,3 +105,24 @@ class TallyDerivative(EqualityMixin, IDManagerMixin):
         if self.variable == 'nuclide_density':
             element.set("nuclide", self.nuclide)
         return element
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate tally derivative from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.TallyDerivative
+            Tally derivative object
+
+        """
+        derivative_id = int(elem.get("id"))
+        variable = elem.get("variable")
+        material = int(elem.get("material"))
+        nuclide = elem.get("nuclide") if variable == "nuclide_density" else None
+        return cls(derivative_id, variable, material, nuclide)
diff --git a/openmc/trigger.py b/openmc/trigger.py
index 3b00ccd204..550f7ee312 100644
--- a/openmc/trigger.py
+++ b/openmc/trigger.py
@@ -74,7 +74,7 @@ class Trigger(EqualityMixin):
             if score not in self._scores:
                 self._scores.append(score)
 
-    def get_trigger_xml(self, element):
+    def to_xml_element(self):
         """Return XML representation of the trigger
 
         Returns
@@ -84,8 +84,36 @@ class Trigger(EqualityMixin):
 
         """
 
-        subelement = ET.SubElement(element, "trigger")
-        subelement.set("type", self._trigger_type)
-        subelement.set("threshold", str(self._threshold))
+        element = ET.Element("trigger")
+        element.set("type", self._trigger_type)
+        element.set("threshold", str(self._threshold))
         if len(self._scores) != 0:
-            subelement.set("scores", ' '.join(map(str, self._scores)))
+            element.set("scores", ' '.join(map(str, self._scores)))
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate trigger object from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.Trigger
+            Trigger object
+
+        """
+        # Generate trigger object
+        trigger_type = elem.get("type")
+        threshold = float(elem.get("threshold"))
+        trigger = cls(trigger_type, threshold)
+
+        # Add scores if present
+        scores = elem.get("scores")
+        if scores is not None:
+            trigger.scores = scores.split()
+
+        return trigger

From 5bb4b4bbdfb6c69917dad72a4f3c6509612d896b Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 20 Jan 2022 20:23:27 -0600
Subject: [PATCH 5/8] Implement Plots.from_xml classmethod

---
 openmc/model/model.py |  15 +++---
 openmc/plots.py       | 103 ++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 112 insertions(+), 6 deletions(-)

diff --git a/openmc/model/model.py b/openmc/model/model.py
index 5c9fefda7d..cf082407c0 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -5,7 +5,6 @@ import os
 from pathlib import Path
 from numbers import Integral
 from tempfile import NamedTemporaryFile
-import time
 
 import h5py
 
@@ -204,11 +203,10 @@ class Model:
 
     @classmethod
     def from_xml(cls, geometry='geometry.xml', materials='materials.xml',
-                 settings='settings.xml', tallies='tallies.xml'):
+                 settings='settings.xml', tallies='tallies.xml',
+                 plots='plots.xml'):
         """Create model from existing XML files
 
-        When initializing this way, the user must manually load plots.
-
         Parameters
         ----------
         geometry : str
@@ -220,6 +218,10 @@ class Model:
         tallies : str
             Path to tallies.xml file
 
+            .. versionadded:: 0.13.0
+        plots : str
+            Path to plots.xml file
+
             .. versionadded:: 0.13.0
 
         Returns
@@ -231,8 +233,9 @@ class Model:
         materials = openmc.Materials.from_xml(materials)
         geometry = openmc.Geometry.from_xml(geometry, materials)
         settings = openmc.Settings.from_xml(settings)
-        tallies = openmc.Tallies.from_xml(tallies)
-        return cls(geometry, materials, settings, tallies)
+        tallies = openmc.Tallies.from_xml(tallies) if Path(tallies).exists() else None
+        plots = openmc.Plots.from_xml(plots) if Path(plots).exists() else None
+        return cls(geometry, materials, settings, tallies, plots)
 
     def init_lib(self, threads=None, geometry_debug=False, restart_file=None,
                  tracks=False, output=True, event_based=None, intracomm=None):
diff --git a/openmc/plots.py b/openmc/plots.py
index 521be95ad2..6c5a5d4372 100644
--- a/openmc/plots.py
+++ b/openmc/plots.py
@@ -684,6 +684,85 @@ class Plot(IDManagerMixin):
 
         return element
 
+    @classmethod
+    def from_xml_element(cls, elem):
+        """Generate plot object from an XML element
+
+        Parameters
+        ----------
+        elem : xml.etree.ElementTree.Element
+            XML element
+
+        Returns
+        -------
+        openmc.Plot
+            Plot object
+
+        """
+        plot_id = int(elem.get("id"))
+        plot = cls(plot_id)
+        if "filename" in elem.keys():
+            plot.filename = elem.get("filename")
+        plot.color_by = elem.get("color_by")
+        plot.type = elem.get("type")
+        plot.basis = elem.get("basis")
+
+        # Helper function to get a tuple of values
+        def get_tuple(elem, name, dtype=int):
+            subelem = elem.find(name)
+            if subelem is not None:
+                return tuple([dtype(x) for x in subelem.text.split()])
+
+        plot.origin = get_tuple(elem, "origin", float)
+        plot.width = get_tuple(elem, "width", float)
+        plot.pixels = get_tuple(elem, "pixels")
+        plot._background = get_tuple(elem, "background")
+
+        # Set plot colors
+        colors = {}
+        for color_elem in elem.findall("color"):
+            uid = color_elem.get("id")
+            colors[uid] = tuple([int(x) for x in color_elem.get("rgb").split()])
+        # TODO: set colors (needs geometry information)
+
+        # Set masking information
+        mask_elem = elem.find("mask")
+        if mask_elem is not None:
+            mask_components = [int(x) for x in mask_elem.get("components").split()]
+            # TODO: set mask components (needs geometry information)
+            background = mask_elem.get("background")
+            if background is not None:
+                plot.mask_background = tuple([int(x) for x in background.split()])
+
+        # show overlaps
+        overlap_elem = elem.find("show_overlaps")
+        if overlap_elem is not None:
+            plot.show_overlaps = (overlap_elem.text in ('true', '1'))
+        overlap_color = get_tuple(elem, "overlap_color")
+        if overlap_color is not None:
+            plot.overlap_color = overlap_color
+
+        # Set universe level
+        level = elem.find("level")
+        if level is not None:
+            plot.level = int(level.text)
+
+        # Set meshlines
+        mesh_elem = elem.find("meshlines")
+        if mesh_elem is not None:
+            meshlines = {'type': mesh_elem.get('meshtype')}
+            if 'id' in mesh_elem.keys():
+                meshlines['id'] = int(mesh_elem.get('id'))
+            if 'linewidth' in mesh_elem.keys():
+                meshlines['linewidth'] = int(mesh_elem.get('linewidth'))
+            if 'color' in mesh_elem.keys():
+                meshlines['color'] = tuple(
+                    [int(x) for x in mesh_elem.get('color').split()]
+                )
+            plot.meshlines = meshlines
+
+        return plot
+
     def to_ipython_image(self, openmc_exec='openmc', cwd='.'):
         """Render plot as an image
 
@@ -848,3 +927,27 @@ class Plots(cv.CheckedList):
         reorder_attributes(self._plots_file)  # TODO: Remove when support is Python 3.8+
         tree = ET.ElementTree(self._plots_file)
         tree.write(str(p), xml_declaration=True, encoding='utf-8')
+
+    @classmethod
+    def from_xml(cls, path='plots.xml'):
+        """Generate plots collection from XML file
+
+        Parameters
+        ----------
+        path : str, optional
+            Path to plots XML file
+
+        Returns
+        -------
+        openmc.Plots
+            Plots collection
+
+        """
+        tree = ET.parse(path)
+        root = tree.getroot()
+
+        # Generate each plot
+        plots = cls()
+        for elem in root.findall('plot'):
+            plots.append(Plot.from_xml_element(elem))
+        return plots

From d57c14e06f3814d1d129244bff14bf4e85bead7b Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 20 Jan 2022 22:21:54 -0600
Subject: [PATCH 6/8] Add tests for XML roundtripping

---
 tests/unit_tests/test_filters.py | 41 ++++++++++++++++++++++++++++++--
 tests/unit_tests/test_plots.py   | 22 ++++++++++++++---
 tests/unit_tests/test_tallies.py | 40 +++++++++++++++++++++++++++++++
 3 files changed, 98 insertions(+), 5 deletions(-)
 create mode 100644 tests/unit_tests/test_tallies.py

diff --git a/tests/unit_tests/test_filters.py b/tests/unit_tests/test_filters.py
index ea94a600e3..49526ac673 100644
--- a/tests/unit_tests/test_filters.py
+++ b/tests/unit_tests/test_filters.py
@@ -1,5 +1,4 @@
-from math import sqrt, pi
-
+import numpy as np
 import openmc
 from pytest import fixture, approx
 
@@ -38,6 +37,11 @@ def test_cell_instance():
     assert all(x == c1.id for x in bins[:6:2])
     assert all(x == c2.id for x in bins[6::2])
 
+    # from_xml_element()
+    new_f = openmc.Filter.from_xml_element(elem)
+    assert new_f.id == f.id
+    assert np.all(new_f.bins == f.bins)
+
     # get_pandas_dataframe()
     df = f.get_pandas_dataframe(f.num_bins, 1)
     cells = df['cellinstance', 'cell']
@@ -61,6 +65,11 @@ def test_collision():
     assert elem.tag == 'filter'
     assert elem.attrib['type'] == 'collision'
 
+    # from_xml_element()
+    new_f = openmc.Filter.from_xml_element(elem)
+    assert new_f.id == f.id
+    assert np.all(new_f.bins == f.bins)
+
 
 def test_legendre():
     n = 5
@@ -79,6 +88,11 @@ def test_legendre():
     assert elem.attrib['type'] == 'legendre'
     assert elem.find('order').text == str(n)
 
+    # from_xml_element()
+    new_f = openmc.Filter.from_xml_element(elem)
+    assert new_f.id == f.id
+    assert new_f.bins, f.bins
+
 
 def test_spatial_legendre():
     n = 5
@@ -102,6 +116,12 @@ def test_spatial_legendre():
     assert elem.find('order').text == str(n)
     assert elem.find('axis').text == str(axis)
 
+    # from_xml_element()
+    new_f = openmc.Filter.from_xml_element(elem)
+    assert new_f.id == f.id
+    assert new_f.order == f.order
+    assert new_f.axis == f.axis
+
 
 def test_spherical_harmonics():
     n = 3
@@ -122,6 +142,12 @@ def test_spherical_harmonics():
     assert elem.attrib['cosine'] == f.cosine
     assert elem.find('order').text == str(n)
 
+    # from_xml_element()
+    new_f = openmc.Filter.from_xml_element(elem)
+    assert new_f.id == f.id
+    assert new_f.order == f.order
+    assert new_f.cosine == f.cosine
+
 
 def test_zernike():
     n = 4
@@ -140,6 +166,12 @@ def test_zernike():
     assert elem.attrib['type'] == 'zernike'
     assert elem.find('order').text == str(n)
 
+    # from_xml_element()
+    new_f = openmc.Filter.from_xml_element(elem)
+    for attr in ('id', 'order', 'x', 'y', 'r'):
+        assert getattr(new_f, attr) == getattr(f, attr)
+
+
 def test_zernike_radial():
     n = 4
     f = openmc.ZernikeRadialFilter(n, 0., 0., 1.)
@@ -157,6 +189,11 @@ def test_zernike_radial():
     assert elem.attrib['type'] == 'zernikeradial'
     assert elem.find('order').text == str(n)
 
+    # from_xml_element()
+    new_f = openmc.Filter.from_xml_element(elem)
+    for attr in ('id', 'order', 'x', 'y', 'r'):
+        assert getattr(new_f, attr) == getattr(f, attr)
+
 
 def test_first_moment(run_in_tmpdir, box_model):
     plain_tally = openmc.Tally()
diff --git a/tests/unit_tests/test_plots.py b/tests/unit_tests/test_plots.py
index f2d3e10148..b55dbe9be7 100644
--- a/tests/unit_tests/test_plots.py
+++ b/tests/unit_tests/test_plots.py
@@ -1,3 +1,4 @@
+import numpy as np
 import openmc
 import openmc.examples
 import pytest
@@ -12,8 +13,8 @@ def myplot():
     plot.filename = 'myplot'
     plot.type = 'slice'
     plot.basis = 'yz'
-    plot.background = (0, 0, 0)
     plot.background = 'black'
+    plot.background = (0, 0, 0)
 
     plot.color_by = 'material'
     m1, m2 = openmc.Material(), openmc.Material()
@@ -21,8 +22,8 @@ def myplot():
     plot.colors = {m1: 'green', m2: 'blue'}
 
     plot.mask_components = [openmc.Material()]
-    plot.mask_background = (255, 255, 255)
     plot.mask_background = 'white'
+    plot.mask_background = (255, 255, 255)
 
     plot.overlap_color = (255, 211, 0)
     plot.overlap_color = 'yellow'
@@ -70,7 +71,7 @@ def test_highlight_domains():
     plots.highlight_domains(model.geometry, mats)
 
 
-def test_to_xml_element(myplot):
+def test_xml_element(myplot):
     elem = myplot.to_xml_element()
     assert 'id' in elem.attrib
     assert 'color_by' in elem.attrib
@@ -80,10 +81,19 @@ def test_to_xml_element(myplot):
     assert elem.find('pixels') is not None
     assert elem.find('background').text == '0 0 0'
 
+    newplot = openmc.Plot.from_xml_element(elem)
+    attributes = ('id', 'color_by', 'filename', 'type', 'basis', 'level',
+                  'meshlines', 'show_overlaps', 'origin', 'width', 'pixels',
+                  'background', 'mask_background')
+    for attr in attributes:
+        assert getattr(newplot, attr) == getattr(myplot, attr), attr
+
 
 def test_plots(run_in_tmpdir):
     p1 = openmc.Plot(name='plot1')
+    p1.origin = (5., 5., 5.)
     p2 = openmc.Plot(name='plot2')
+    p2.origin = (-3., -3., -3.)
     plots = openmc.Plots([p1, p2])
     assert len(plots) == 2
 
@@ -92,3 +102,9 @@ def test_plots(run_in_tmpdir):
     assert len(plots) == 3
 
     plots.export_to_xml()
+
+    # from_xml
+    new_plots = openmc.Plots.from_xml()
+    assert len(plots)
+    assert plots[0].origin == p1.origin
+    assert plots[1].origin == p2.origin
diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py
new file mode 100644
index 0000000000..aeeb0612d5
--- /dev/null
+++ b/tests/unit_tests/test_tallies.py
@@ -0,0 +1,40 @@
+import numpy as np
+import openmc
+
+
+def test_xml_roundtrip(run_in_tmpdir):
+    # Create a tally with all possible gizmos
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = (-10., -10., -10.)
+    mesh.upper_right = (10., 10., 10.,)
+    mesh.dimension = (5, 5, 5)
+    mesh_filter = openmc.MeshFilter(mesh)
+    tally = openmc.Tally()
+    tally.filters = [mesh_filter]
+    tally.nuclides = ['U235', 'I135', 'Li6']
+    tally.scores = ['total', 'fission', 'heating']
+    tally.derivative = openmc.TallyDerivative(
+        variable='nuclide_density', material=1, nuclide='Li6'
+    )
+    tally.triggers = [openmc.Trigger('rel_err', 0.025)]
+    tally.triggers[0].scores = ['total', 'fission']
+    tallies = openmc.Tallies([tally])
+
+    # Roundtrip through XML and make sure we get what we started with
+    tallies.export_to_xml()
+    new_tallies = openmc.Tallies.from_xml()
+    assert len(new_tallies) == 1
+    new_tally = new_tallies[0]
+    assert new_tally.id == tally.id
+    assert len(new_tally.filters) == 1
+    assert isinstance(new_tally.filters[0], openmc.MeshFilter)
+    assert np.allclose(new_tally.filters[0].mesh.lower_left, mesh.lower_left)
+    assert new_tally.nuclides == tally.nuclides
+    assert new_tally.scores == tally.scores
+    assert new_tally.derivative.variable == tally.derivative.variable
+    assert new_tally.derivative.material == tally.derivative.material
+    assert new_tally.derivative.nuclide == tally.derivative.nuclide
+    assert len(new_tally.triggers) == 1
+    assert new_tally.triggers[0].trigger_type == tally.triggers[0].trigger_type
+    assert new_tally.triggers[0].threshold == tally.triggers[0].threshold
+    assert new_tally.triggers[0].scores == tally.triggers[0].scores

From 6bb6e2ec0cae0f562180018a4d6cd2ef73e5d6a2 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 21 Jan 2022 08:18:54 -0600
Subject: [PATCH 7/8] Fix failing unit tests

---
 openmc/lib/mesh.py             | 68 ++++++++++++++++++++++++++--------
 tests/unit_tests/test_lib.py   | 30 +++++++--------
 tests/unit_tests/test_model.py |  4 +-
 3 files changed, 69 insertions(+), 33 deletions(-)

diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 43f66d7a12..df64eee30e 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -285,6 +285,18 @@ class RectilinearMesh(Mesh):
         return (lower_left, upper_right, dimension, width)
 
     def set_grid(self, x_grid, y_grid, z_grid):
+        """Set grid values
+
+        Parameters
+        ----------
+        x_grid : iterable of float
+            Mesh boundary points along the x-axis.
+        y_grid : iterable of float
+            Mesh boundary points along the y-axis.
+        z_grid : iterable of float
+            Mesh boundary points along the z-axis.
+
+        """
         nx = len(x_grid)
         x_grid = (c_double*nx)(*x_grid)
         ny = len(y_grid)
@@ -371,15 +383,27 @@ class CylindricalMesh(Mesh):
 
         return (lower_left, upper_right, dimension, width)
 
-    def set_grid(self, x_grid, y_grid, z_grid):
-        nx = len(x_grid)
-        x_grid = (c_double*nx)(*x_grid)
-        ny = len(y_grid)
-        y_grid = (c_double*ny)(*y_grid)
+    def set_grid(self, r_grid, phi_grid, z_grid):
+        """Set grid values
+
+        Parameters
+        ----------
+        r_grid : iterable of float
+            Mesh boundary points along the r-axis
+        phi_grid : Iterable of float
+            Mesh boundary points along the phi-axis
+        z_grid : Iterable of float
+            Mesh boundary points along the z-axis
+
+        """
+        nr = len(r_grid)
+        r_grid = (c_double*nr)(*r_grid)
+        nphi = len(phi_grid)
+        phi_grid = (c_double*nphi)(*phi_grid)
         nz = len(z_grid)
         z_grid = (c_double*nz)(*z_grid)
-        _dll.openmc_cylindrical_mesh_set_grid(self._index, x_grid, nx, y_grid,
-                                              ny, z_grid, nz)                                              
+        _dll.openmc_cylindrical_mesh_set_grid(self._index, r_grid, nr, phi_grid,
+                                              nphi, z_grid, nz)
 
 class SphericalMesh(Mesh):
     """SphericalMesh stored internally.
@@ -457,15 +481,27 @@ class SphericalMesh(Mesh):
 
         return (lower_left, upper_right, dimension, width)
 
-    def set_grid(self, x_grid, y_grid, z_grid):
-        nx = len(x_grid)
-        x_grid = (c_double*nx)(*x_grid)
-        ny = len(y_grid)
-        y_grid = (c_double*ny)(*y_grid)
-        nz = len(z_grid)
-        z_grid = (c_double*nz)(*z_grid)
-        _dll.openmc_spherical_mesh_set_grid(self._index, x_grid, nx, y_grid,
-                                              ny, z_grid, nz)                                              
+    def set_grid(self, r_grid, theta_grid, phi_grid):
+        """Set grid values
+
+        Parameters
+        ----------
+        r_grid : iterable of float
+            Mesh boundary points along the r-axis
+        theta_grid : Iterable of float
+            Mesh boundary points along the theta-axis
+        phi_grid : Iterable of float
+            Mesh boundary points along the phi-axis
+
+        """
+        nr = len(r_grid)
+        r_grid = (c_double*nr)(*r_grid)
+        ntheta = len(theta_grid)
+        theta_grid = (c_double*ntheta)(*theta_grid)
+        nphi = len(phi_grid)
+        phi_grid = (c_double*nphi)(*phi_grid)
+        _dll.openmc_spherical_mesh_set_grid(self._index, r_grid, nr, theta_grid,
+                                              ntheta, phi_grid, nphi)
 
 
 _MESH_TYPE_MAP = {
diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py
index 3572e9863a..e14e04adb2 100644
--- a/tests/unit_tests/test_lib.py
+++ b/tests/unit_tests/test_lib.py
@@ -583,17 +583,17 @@ def test_rectilinear_mesh(lib_init):
 def test_cylindrical_mesh(lib_init):
     deg2rad = lambda deg: deg*np.pi/180
     mesh = openmc.lib.CylindricalMesh()
-    x_grid = [0., 5., 10.]
-    y_grid = [0., 10., 20.]
+    r_grid = [0., 5., 10.]
+    phi_grid = np.radians([0., 10., 20.])
     z_grid = [10., 20., 30.]
-    mesh.set_grid(x_grid, y_grid, z_grid)
+    mesh.set_grid(r_grid, phi_grid, z_grid)
     assert np.all(mesh.lower_left == (0., 0., 10.))
     assert np.all(mesh.upper_right == (10., deg2rad(20.), 30.))
     assert np.all(mesh.dimension == (2, 2, 2))
-    for i, diff_x in enumerate(np.diff(x_grid)):
-        for j, diff_y in enumerate(np.diff(y_grid)):
-            for k, diff_z in enumerate(np.diff(z_grid)):
-                assert np.all(mesh.width[i, j, k, :] == (5, deg2rad(10), 10))
+    for i, _ in enumerate(np.diff(r_grid)):
+        for j, _ in enumerate(np.diff(phi_grid)):
+            for k, _ in enumerate(np.diff(z_grid)):
+                assert np.allclose(mesh.width[i, j, k, :], (5, deg2rad(10), 10))
 
     with pytest.raises(exc.AllocationError):
         mesh2 = openmc.lib.CylindricalMesh(mesh.id)
@@ -614,17 +614,17 @@ def test_cylindrical_mesh(lib_init):
 def test_spherical_mesh(lib_init):
     deg2rad = lambda deg: deg*np.pi/180
     mesh = openmc.lib.SphericalMesh()
-    x_grid = [0., 5., 10.]
-    y_grid = [0., 10., 20.]
-    z_grid = [10., 20., 30.]
-    mesh.set_grid(x_grid, y_grid, z_grid)
+    r_grid = [0., 5., 10.]
+    theta_grid = np.radians([0., 10., 20.])
+    phi_grid = np.radians([10., 20., 30.])
+    mesh.set_grid(r_grid, theta_grid, phi_grid)
     assert np.all(mesh.lower_left == (0., 0., deg2rad(10.)))
     assert np.all(mesh.upper_right == (10., deg2rad(20.), deg2rad(30.)))
     assert np.all(mesh.dimension == (2, 2, 2))
-    for i, diff_x in enumerate(np.diff(x_grid)):
-        for j, diff_y in enumerate(np.diff(y_grid)):
-            for k, diff_z in enumerate(np.diff(z_grid)):
-                assert np.all(abs(mesh.width[i, j, k, :] - (5, deg2rad(10), deg2rad(10))) < 1e-16)
+    for i, _ in enumerate(np.diff(r_grid)):
+        for j, _ in enumerate(np.diff(theta_grid)):
+            for k, _ in enumerate(np.diff(phi_grid)):
+                assert np.allclose(mesh.width[i, j, k, :], (5, deg2rad(10), deg2rad(10)))
 
     with pytest.raises(exc.AllocationError):
         mesh2 = openmc.lib.SphericalMesh(mesh.id)
diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py
index d22485b3cd..80c13fe617 100644
--- a/tests/unit_tests/test_model.py
+++ b/tests/unit_tests/test_model.py
@@ -195,8 +195,8 @@ def test_from_xml(run_in_tmpdir, pin_model_attributes):
     keys = sorted(k for k in settings.__dict__.keys() if k not in no_test)
     for ref_k in keys:
         assert test_model.settings.__dict__[ref_k] == settings.__dict__[ref_k]
-    assert len(test_model.tallies) == 0
-    assert len(test_model.plots) == 0
+    assert len(test_model.tallies) == 1
+    assert len(test_model.plots) == 2
     assert test_model._materials_by_id == \
         {1: test_model.materials[0], 2: test_model.materials[1],
          3: test_model.materials[2]}

From ca58e1407388a66bc9890c524197cc44ca4d9e7c Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 23 Jan 2022 14:03:18 -0600
Subject: [PATCH 8/8] Apply suggestions from @pshriwise code review

Co-authored-by: Patrick Shriwise <pshriwise@gmail.com>
---
 openmc/filter.py   | 4 ++--
 openmc/lib/mesh.py | 6 +++---
 openmc/trigger.py  | 2 +-
 src/mesh.cpp       | 6 +++---
 4 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/openmc/filter.py b/openmc/filter.py
index fbbc855c28..65255bd3d4 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -254,7 +254,7 @@ class Filter(IDManagerMixin, metaclass=FilterMeta):
         filter_type = elem.get('type')
 
         # If the filter type matches this class's short_name, then
-        # there is no overriden from_xml_element method.
+        # there is no overriden from_xml_element method
         if filter_type == cls.short_name.lower():
             # Get bins from element -- the default here works for any filters
             # that just store a list of bins that can be represented as integers
@@ -263,7 +263,7 @@ class Filter(IDManagerMixin, metaclass=FilterMeta):
             return cls(bins, filter_id=filter_id)
 
         # Search through all subclasses and find the one matching the HDF5
-        # 'type'.  Call that class's from_hdf5 method.
+        # 'type'.  Call that class's from_hdf5 method
         for subclass in cls._recursive_subclasses():
             if filter_type == subclass.short_name.lower():
                 return subclass.from_xml_element(elem, **kwargs)
diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index df64eee30e..84571628c6 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -290,11 +290,11 @@ class RectilinearMesh(Mesh):
         Parameters
         ----------
         x_grid : iterable of float
-            Mesh boundary points along the x-axis.
+            Mesh boundary points along the x-axis
         y_grid : iterable of float
-            Mesh boundary points along the y-axis.
+            Mesh boundary points along the y-axis
         z_grid : iterable of float
-            Mesh boundary points along the z-axis.
+            Mesh boundary points along the z-axis
 
         """
         nx = len(x_grid)
diff --git a/openmc/trigger.py b/openmc/trigger.py
index 550f7ee312..c5151aa7e0 100644
--- a/openmc/trigger.py
+++ b/openmc/trigger.py
@@ -88,7 +88,7 @@ class Trigger(EqualityMixin):
         element.set("type", self._trigger_type)
         element.set("threshold", str(self._threshold))
         if len(self._scores) != 0:
-            element.set("scores", ' '.join(map(str, self._scores)))
+            element.set("scores", ' '.join(self._scores))
         return element
 
     @classmethod
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 9d5da29e9b..9310248c39 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -1001,7 +1001,7 @@ double CylindricalMesh::find_r_crossing(
 double CylindricalMesh::find_phi_crossing(
   const Position& r, const Direction& u, double l, int shell) const
 {
-  // Phi grid is [0, 2pi], thus there is no real surface to cross
+  // Phi grid is [0, 2π], thus there is no real surface to cross
   if (full_phi_ && (shape_[1] == 1))
     return INFTY;
 
@@ -1221,7 +1221,7 @@ double SphericalMesh::find_r_crossing(
 double SphericalMesh::find_theta_crossing(
   const Position& r, const Direction& u, double l, int shell) const
 {
-  // Theta grid is [0, pi], thus there is no real surface to cross
+  // Theta grid is [0, π], thus there is no real surface to cross
   if (full_theta_ && (shape_[1] == 1))
     return INFTY;
 
@@ -1283,7 +1283,7 @@ double SphericalMesh::find_theta_crossing(
 double SphericalMesh::find_phi_crossing(
   const Position& r, const Direction& u, double l, int shell) const
 {
-  // Phi grid is [0, 2pi], thus there is no real surface to cross
+  // Phi grid is [0, 2π], thus there is no real surface to cross
   if (full_phi_ && (shape_[2] == 1))
     return INFTY;