Repository: openmc-dev/openmc

Issue:
CylinderSector  doesn't properly generate planes for off origin centers
## Bug Description
If you create an `openmc.model.CylinderSector` that is not located at the origin, the planes that define the sector do not appear to intersect the center. 

Reading the code it seems like this issue is [this line and the next few](https://github.com/openmc-dev/openmc/blob/1702b4554bfb1a7f19cc2b9f1318edeb0f7ad0d2/openmc/model/surface_composite.py#L136). It seems to me that the assumption with these planes is that they are going through the origin and they are parallel to the Z-axis. I tried to read through #2027 to better understand the logic, and there doesn't seem to be any discussion of this. @yardasol, am I missing something?


## Steps to Reproduce
<!--Steps to reproduce the behavior (input file, or modifications to an existing input file, etc.)-->

1. Create an off-origin cylinder sector

``` python
import openmc
center = (10, 0, 0)
sector = openmc.model.CylinderSector(1.0, 1.1, 90, 95, center[0:2])
```

2. If we try to plot this we get errors, likely due to the null set region (related to #3015):

``` python
cell = openmc.Cell(region=-sector)
cell.plot()
```
This leads to errors  shown in #3015.

3. What is going on with these surfaces?

``` python
>  print(sector.inner_cyl, sector.outer_cyl, sector.plane1, sector.plane2)
Surface
        ID                 =    16
        Name               =
        Type               =    z-cylinder
        Boundary           =    transmission
        Coefficients
x0                  =   10
y0                  =   0
r                   =   1.0
 Surface
        ID                 =    17
        Name               =
        Type               =    z-cylinder
        Boundary           =    transmission
        Coefficients
x0                  =   10
y0                  =   0
r                   =   1.1
 Surface
        ID                 =    18
        Name               =
        Type               =    plane
        Boundary           =    transmission
        Coefficients
a                   =   0.10000000000000009
b                   =   -6.123233995736777e-18
c                   =   0.0
d                   =   0.0
 Surface
        ID                 =    19
        Name               =
        Type               =    plane
        Boundary           =    transmission
        Coefficients
a                   =   0.09961946980917458
b                   =   0.008715574274765825
c                   =   0.0
d                   =   0.0
```

The cylinders are properly located. However it looks like the planes go through the origin. This can be easily confirmed:

``` python
In [22]: half = -sector.plane1

In [23]: epsilon = 1e-6

In [24]: (epsilon, 0, 0) in half
Out[24]: False

In [25]: (-epsilon, 0, 0) in half
Out[25]: True
```

Therefore, the plane does go through the origin. 


## Environment
* openmc: 0.14.0


Discussion:
MicahGale: Don't the planes also need to be rotated if the axis isn't `z`?

yardasol: > Don't the planes also need to be rotated if the axis isn't `z`?

@MicahGale I believe the planes are rotated. It's a little subtle, but the way we handle this in this funtion is by constructing the planes directly from a set of points. Those points will have the name numerical values, and which value goes to (x,y,z) depends on the axis we want the planes to be parallel to. Lines 159-161 shuffle the points accordingly. I added a comment to make this more clear:

```python
        # Reorder the points to correspond to the correct central axis
        for p in points:
            p[:] = p[coord_map]
```
MicahGale: Oh I missed; very clever. I see also that you do test for different orientations. Thanks for clarifying.

Patch:
diff --git a/openmc/cell.py b/openmc/cell.py
index 94fac8413..fe3939bbe 100644
--- a/openmc/cell.py
+++ b/openmc/cell.py
@@ -426,15 +426,13 @@ class Cell(IDManagerMixin):
             instances
 
         """
+        if memo is None:
+            memo = set()
+        elif self in memo:
+            return {}
+        memo.add(self)
 
         cells = {}
-
-        if memo and self in memo:
-            return cells
-
-        if memo is not None:
-            memo.add(self)
-
         if self.fill_type in ('universe', 'lattice'):
             cells.update(self.fill.get_all_cells(memo))
 
@@ -465,7 +463,7 @@ class Cell(IDManagerMixin):
 
         return materials
 
-    def get_all_universes(self):
+    def get_all_universes(self, memo=None):
         """Return all universes that are contained within this one if any of
         its cells are filled with a universe or lattice.
 
@@ -476,18 +474,22 @@ class Cell(IDManagerMixin):
             :class:`Universe` instances
 
         """
+        if memo is None:
+            memo = set()
+        if self in memo:
+            return {}
+        memo.add(self)
 
         universes = {}
-
         if self.fill_type == 'universe':
             universes[self.fill.id] = self.fill
-            universes.update(self.fill.get_all_universes())
+            universes.update(self.fill.get_all_universes(memo))
         elif self.fill_type == 'lattice':
-            universes.update(self.fill.get_all_universes())
+            universes.update(self.fill.get_all_universes(memo))
 
         return universes
 
-    def clone(self, clone_materials=True, clone_regions=True,  memo=None):
+    def clone(self, clone_materials=True, clone_regions=True, memo=None):
         """Create a copy of this cell with a new unique ID, and clones
         the cell's region and fill.
 
@@ -678,10 +680,11 @@ class Cell(IDManagerMixin):
             # thus far.
             def create_surface_elements(node, element, memo=None):
                 if isinstance(node, Halfspace):
-                    if memo and node.surface in memo:
+                    if memo is None:
+                        memo = set()
+                    elif node.surface in memo:
                         return
-                    if memo is not None:
-                        memo.add(node.surface)
+                    memo.add(node.surface)
                     xml_element.append(node.surface.to_xml_element())
 
                 elif isinstance(node, Complement):
diff --git a/openmc/checkvalue.py b/openmc/checkvalue.py
index 840306486..42df3e5ef 100644
--- a/openmc/checkvalue.py
+++ b/openmc/checkvalue.py
@@ -171,6 +171,29 @@ def check_length(name, value, length_min, length_max=None):
         raise ValueError(msg)
 
 
+def check_increasing(name: str, value, equality: bool = False):
+    """Ensure that a list's elements are strictly or loosely increasing.
+
+    Parameters
+    ----------
+    name : str
+        Description of value being checked
+    value : iterable
+        Object to check if increasing
+    equality : bool, optional
+        Whether equality is allowed. Defaults to False.
+
+    """
+    if equality:
+        if not np.all(np.diff(value) >= 0.0):
+            raise ValueError(f'Unable to set "{name}" to "{value}" since its '
+                             'elements must be increasing.')
+    elif not equality:
+        if not np.all(np.diff(value) > 0.0):
+            raise ValueError(f'Unable to set "{name}" to "{value}" since its '
+                             'elements must be strictly increasing.')
+
+
 def check_value(name, value, accepted_values):
     """Ensure that an object's value is contained in a set of acceptable values.
 
diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py
index f5841d9c9..6df171d3a 100644
--- a/openmc/data/thermal.py
+++ b/openmc/data/thermal.py
@@ -698,8 +698,8 @@ class ThermalScattering(EqualityMixin):
                 # add an outgoing energy 0 eV that has a PDF of 0 (and of
                 # course, a CDF of 0 as well).
                 if eout_i.c[0] > 0.:
-                    eout_i.x = np.insert(eout_i.x, 0, 0.)
-                    eout_i.p = np.insert(eout_i.p, 0, 0.)
+                    eout_i._x = np.insert(eout_i.x, 0, 0.)
+                    eout_i._p = np.insert(eout_i.p, 0, 0.)
                     eout_i.c = np.insert(eout_i.c, 0, 0.)
 
                     # For this added outgoing energy (of 0 eV) we add a set of
diff --git a/openmc/geometry.py b/openmc/geometry.py
index 175cef2bf..db927c46e 100644
--- a/openmc/geometry.py
+++ b/openmc/geometry.py
@@ -134,7 +134,7 @@ class Geometry:
 
         # Create XML representation
         element = ET.Element("geometry")
-        self.root_universe.create_xml_subelement(element, memo=set())
+        self.root_universe.create_xml_subelement(element)
 
         # Sort the elements in the file
         element[:] = sorted(element, key=lambda x: (
@@ -373,7 +373,7 @@ class Geometry:
 
         """
         if self.root_universe is not None:
-            return self.root_universe.get_all_cells(memo=set())
+            return self.root_universe.get_all_cells()
         else:
             return {}
 
@@ -417,7 +417,7 @@ class Geometry:
 
         """
         if self.root_universe is not None:
-            return self.root_universe.get_all_materials(memo=set())
+            return self.root_universe.get_all_materials()
         else:
             return {}
 
diff --git a/openmc/lattice.py b/openmc/lattice.py
index 40b97f6cf..f7f9da8ea 100644
--- a/openmc/lattice.py
+++ b/openmc/lattice.py
@@ -170,11 +170,11 @@ class Lattice(IDManagerMixin, ABC):
         """
         cells = {}
 
-        if memo and self in memo:
+        if memo is None:
+            memo = set()
+        elif self in memo:
             return cells
-
-        if memo is not None:
-            memo.add(self)
+        memo.add(self)
 
         unique_universes = self.get_unique_universes()
 
@@ -194,6 +194,9 @@ class Lattice(IDManagerMixin, ABC):
 
         """
 
+        if memo is None:
+            memo = set()
+
         materials = {}
 
         # Append all Cells in each Cell in the Universe to the dictionary
@@ -203,7 +206,7 @@ class Lattice(IDManagerMixin, ABC):
 
         return materials
 
-    def get_all_universes(self):
+    def get_all_universes(self, memo=None):
         """Return all universes that are contained within the lattice
 
         Returns
@@ -213,11 +216,16 @@ class Lattice(IDManagerMixin, ABC):
             :class:`Universe` instances
 
         """
-
         # Initialize a dictionary of all Universes contained by the Lattice
         # in each nested Universe level
         all_universes = {}
 
+        if memo is None:
+            memo = set()
+        elif self in memo:
+            return all_universes
+        memo.add(self)
+
         # Get all unique Universes contained in each of the lattice cells
         unique_universes = self.get_unique_universes()
 
@@ -226,7 +234,7 @@ class Lattice(IDManagerMixin, ABC):
 
         # Append all Universes containing each cell to the dictionary
         for universe in unique_universes.values():
-            all_universes.update(universe.get_all_universes())
+            all_universes.update(universe.get_all_universes(memo))
 
         return all_universes
 
@@ -846,10 +854,11 @@ class RectLattice(Lattice):
 
         """
         # If the element already contains the Lattice subelement, then return
-        if memo and self in memo:
+        if memo is None:
+            memo = set()
+        elif self in memo:
             return
-        if memo is not None:
-            memo.add(self)
+        memo.add(self)
 
         # Make sure universes have been assigned
         if self.universes is None:
@@ -1417,10 +1426,11 @@ class HexLattice(Lattice):
 
     def create_xml_subelement(self, xml_element, memo=None):
         # If this subelement has already been written, return
-        if memo and self in memo:
+        if memo is None:
+            memo = set()
+        elif self in memo:
             return
-        if memo is not None:
-            memo.add(self)
+        memo.add(self)
 
         lattice_subelement = ET.Element("hex_lattice")
         lattice_subelement.set("id", str(self._id))
@@ -1564,6 +1574,7 @@ class HexLattice(Lattice):
                             alpha -= 1
                             if not lat.is_valid_index((x, alpha, z)):
                                 # Reached the bottom
+                                j += 1
                                 break
                     j += 1
             else:
@@ -1583,6 +1594,7 @@ class HexLattice(Lattice):
 
                         # Check if we've reached the bottom
                         if y == -n_rings:
+                            j += 1
                             break
 
                         while not lat.is_valid_index((alpha, y, z)):
diff --git a/openmc/mesh.py b/openmc/mesh.py
index 8777da325..48263c205 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -1820,6 +1820,8 @@ class SphericalMesh(StructuredMesh):
     @r_grid.setter
     def r_grid(self, grid):
         cv.check_type('mesh r_grid', grid, Iterable, Real)
+        cv.check_length('mesh r_grid', grid, 2)
+        cv.check_increasing('mesh r_grid', grid)
         self._r_grid = np.asarray(grid, dtype=float)
 
     @property
@@ -1829,6 +1831,8 @@ class SphericalMesh(StructuredMesh):
     @theta_grid.setter
     def theta_grid(self, grid):
         cv.check_type('mesh theta_grid', grid, Iterable, Real)
+        cv.check_length('mesh theta_grid', grid, 2)
+        cv.check_increasing('mesh theta_grid', grid)
         self._theta_grid = np.asarray(grid, dtype=float)
 
     @property
@@ -1838,6 +1842,8 @@ class SphericalMesh(StructuredMesh):
     @phi_grid.setter
     def phi_grid(self, grid):
         cv.check_type('mesh phi_grid', grid, Iterable, Real)
+        cv.check_length('mesh phi_grid', grid, 2)
+        cv.check_increasing('mesh phi_grid', grid)
         self._phi_grid = np.asarray(grid, dtype=float)
 
     @property
diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 6f9123311..1fa4234fb 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -90,7 +90,7 @@ class CylinderSector(CompositeSurface):
         counterclockwise direction with respect to the first basis axis
         (+y, +z, or +x). Must be greater than :attr:`theta1`.
     center : iterable of float
-       Coordinate for central axes of cylinders in the (y, z), (z, x), or (x, y)
+       Coordinate for central axes of cylinders in the (y, z), (x, z), or (x, y)
        basis. Defaults to (0,0).
     axis : {'x', 'y', 'z'}
         Central axis of the cylinders defining the inner and outer surfaces of
@@ -133,13 +133,13 @@ class CylinderSector(CompositeSurface):
         phi2 = pi / 180 * theta2
 
         # Coords for axis-perpendicular planes
-        p1 = np.array([0., 0., 1.])
+        p1 = np.array([center[0], center[1], 1.])
 
-        p2_plane1 = np.array([r1 * cos(phi1), r1 * sin(phi1), 0.])
-        p3_plane1 = np.array([r2 * cos(phi1), r2 * sin(phi1), 0.])
+        p2_plane1 = np.array([r1 * cos(phi1) + center[0], r1 * sin(phi1) + center[1], 0.])
+        p3_plane1 = np.array([r2 * cos(phi1) + center[0], r2 * sin(phi1) + center[1], 0.])
 
-        p2_plane2 = np.array([r1 * cos(phi2), r1 * sin(phi2), 0.])
-        p3_plane2 = np.array([r2 * cos(phi2), r2 * sin(phi2), 0.])
+        p2_plane2 = np.array([r1 * cos(phi2) + center[0], r1 * sin(phi2)+ center[1], 0.])
+        p3_plane2 = np.array([r2 * cos(phi2) + center[0], r2 * sin(phi2)+ center[1], 0.])
 
         points = [p1, p2_plane1, p3_plane1, p2_plane2, p3_plane2]
         if axis == 'z':
@@ -147,7 +147,7 @@ class CylinderSector(CompositeSurface):
             self.inner_cyl = openmc.ZCylinder(*center, r1, **kwargs)
             self.outer_cyl = openmc.ZCylinder(*center, r2, **kwargs)
         elif axis == 'y':
-            coord_map = [1, 2, 0]
+            coord_map = [0, 2, 1]
             self.inner_cyl = openmc.YCylinder(*center, r1, **kwargs)
             self.outer_cyl = openmc.YCylinder(*center, r2, **kwargs)
         elif axis == 'x':
@@ -155,6 +155,7 @@ class CylinderSector(CompositeSurface):
             self.inner_cyl = openmc.XCylinder(*center, r1, **kwargs)
             self.outer_cyl = openmc.XCylinder(*center, r2, **kwargs)
 
+        # Reorder the points to correspond to the correct central axis
         for p in points:
             p[:] = p[coord_map]
 
@@ -192,8 +193,8 @@ class CylinderSector(CompositeSurface):
             with respect to the first basis axis (+y, +z, or +x). Note that
             negative values translate to an offset in the clockwise direction.
         center : iterable of float
-            Coordinate for central axes of cylinders in the (y, z), (z, x), or (x, y)
-            basis. Defaults to (0,0).
+            Coordinate for central axes of cylinders in the (y, z), (x, z), or
+            (x, y) basis. Defaults to (0,0).
         axis : {'x', 'y', 'z'}
             Central axis of the cylinders defining the inner and outer surfaces
             of the sector. Defaults to 'z'.
@@ -216,10 +217,16 @@ class CylinderSector(CompositeSurface):
         return cls(r1, r2, theta1, theta2, center=center, axis=axis, **kwargs)
 
     def __neg__(self):
-        return -self.outer_cyl & +self.inner_cyl & -self.plane1 & +self.plane2
+        if isinstance(self.inner_cyl, openmc.YCylinder):
+            return -self.outer_cyl & +self.inner_cyl & +self.plane1 & -self.plane2
+        else:
+            return -self.outer_cyl & +self.inner_cyl & -self.plane1 & +self.plane2
 
     def __pos__(self):
-        return +self.outer_cyl | -self.inner_cyl | +self.plane1 | -self.plane2
+        if isinstance(self.inner_cyl, openmc.YCylinder):
+            return +self.outer_cyl | -self.inner_cyl | -self.plane1 | +self.plane2
+        else:
+            return +self.outer_cyl | -self.inner_cyl | +self.plane1 | -self.plane2
 
 
 class IsogonalOctagon(CompositeSurface):
@@ -284,11 +291,11 @@ class IsogonalOctagon(CompositeSurface):
         c1, c2 = center
 
         # Coords for axis-perpendicular planes
-        ctop = c1 + r1
-        cbottom = c1 - r1
+        cright = c1 + r1
+        cleft = c1 - r1
 
-        cright = c2 + r1
-        cleft = c2 - r1
+        ctop = c2 + r1
+        cbottom = c2 - r1
 
         # Side lengths
         if r2 > r1 * sqrt(2):
@@ -309,7 +316,16 @@ class IsogonalOctagon(CompositeSurface):
         p2_lr = np.array([L_basis_ax, -r1, 0.])
         p3_lr = np.array([L_basis_ax, -r1, 1.])
 
-        points = [p1_ur, p2_ur, p3_ur, p1_lr, p2_lr, p3_lr]
+        p1_ll = -p1_ur
+        p2_ll = -p2_ur
+        p3_ll = -p3_ur
+
+        p1_ul = -p1_lr
+        p2_ul = -p2_lr
+        p3_ul = -p3_lr
+
+        points = [p1_ur, p2_ur, p3_ur, p1_lr, p2_lr, p3_lr,
+                  p1_ll, p2_ll, p3_ll, p1_ul, p2_ul, p3_ul]
 
         # Orientation specific variables
         if axis == 'z':
@@ -331,17 +347,19 @@ class IsogonalOctagon(CompositeSurface):
             self.right = openmc.YPlane(cright, **kwargs)
             self.left = openmc.YPlane(cleft, **kwargs)
 
-        # Put our coordinates in (x,y,z) order
+        # Put our coordinates in (x,y,z) order and add the offset
         for p in points:
+            p[0] += c1
+            p[1] += c2
             p[:] = p[coord_map]
 
         self.upper_right = openmc.Plane.from_points(p1_ur, p2_ur, p3_ur,
                                                     **kwargs)
         self.lower_right = openmc.Plane.from_points(p1_lr, p2_lr, p3_lr,
                                                     **kwargs)
-        self.lower_left = openmc.Plane.from_points(-p1_ur, -p2_ur, -p3_ur,
+        self.lower_left = openmc.Plane.from_points(p1_ll, p2_ll, p3_ll,
                                                    **kwargs)
-        self.upper_left = openmc.Plane.from_points(-p1_lr, -p2_lr, -p3_lr,
+        self.upper_left = openmc.Plane.from_points(p1_ul, p2_ul, p3_ul,
                                                    **kwargs)
 
     def __neg__(self):
diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py
index a0e86c3f8..f44bce67c 100644
--- a/openmc/stats/univariate.py
+++ b/openmc/stats/univariate.py
@@ -837,10 +837,11 @@ class Tabular(Univariate):
     x : Iterable of float
         Tabulated values of the random variable
     p : Iterable of float
-        Tabulated probabilities
+        Tabulated probabilities. For histogram interpolation, if the length of
+        `p` is the same as `x`, the last value is ignored.
     interpolation : {'histogram', 'linear-linear', 'linear-log', 'log-linear', 'log-log'}, optional
-        Indicate whether the density function is constant between tabulated
-        points or linearly-interpolated. Defaults to 'linear-linear'.
+        Indicates how the density function is interpolated between tabulated
+        points. Defaults to 'linear-linear'.
     ignore_negative : bool
         Ignore negative probabilities
 
@@ -850,9 +851,9 @@ class Tabular(Univariate):
         Tabulated values of the random variable
     p : numpy.ndarray
         Tabulated probabilities
-    interpolation : {'histogram', 'linear-linear', 'linear-log', 'log-linear', 'log-log'}, optional
-        Indicate whether the density function is constant between tabulated
-        points or linearly-interpolated.
+    interpolation : {'histogram', 'linear-linear', 'linear-log', 'log-linear', 'log-log'}
+        Indicates how the density function is interpolated between tabulated
+        points. Defaults to 'linear-linear'.
 
     """
 
@@ -863,35 +864,38 @@ class Tabular(Univariate):
             interpolation: str = 'linear-linear',
             ignore_negative: bool = False
         ):
-        self._ignore_negative = ignore_negative
-        self.x = x
-        self.p = p
         self.interpolation = interpolation
 
+        cv.check_type('tabulated values', x, Iterable, Real)
+        cv.check_type('tabulated probabilities', p, Iterable, Real)
+
+        x = np.array(x, dtype=float)
+        p = np.array(p, dtype=float)
+
+        if p.size > x.size:
+            raise ValueError('Number of probabilities exceeds number of table values.')
+        if self.interpolation != 'histogram' and x.size != p.size:
+            raise ValueError(f'Tabulated values ({x.size}) and probabilities '
+                             f'({p.size}) should have the same length')
+
+        if not ignore_negative:
+            for pk in p:
+                cv.check_greater_than('tabulated probability', pk, 0.0, True)
+
+        self._x = x
+        self._p = p
+
     def __len__(self):
-        return len(self.x)
+        return self.p.size
 
     @property
     def x(self):
         return self._x
 
-    @x.setter
-    def x(self, x):
-        cv.check_type('tabulated values', x, Iterable, Real)
-        self._x = np.array(x, dtype=float)
-
     @property
     def p(self):
         return self._p
 
-    @p.setter
-    def p(self, p):
-        cv.check_type('tabulated probabilities', p, Iterable, Real)
-        if not self._ignore_negative:
-            for pk in p:
-                cv.check_greater_than('tabulated probability', pk, 0.0, True)
-        self._p = np.array(p, dtype=float)
-
     @property
     def interpolation(self):
         return self._interpolation
@@ -907,7 +911,7 @@ class Tabular(Univariate):
         p = self.p
 
         if self.interpolation == 'histogram':
-            c[1:] = p[:-1] * np.diff(x)
+            c[1:] = p[:x.size-1] * np.diff(x)
         elif self.interpolation == 'linear-linear':
             c[1:] = 0.5 * (p[:-1] + p[1:]) * np.diff(x)
         else:
@@ -938,7 +942,7 @@ class Tabular(Univariate):
         elif self.interpolation == 'histogram':
             x_l = self.x[:-1]
             x_r = self.x[1:]
-            p_l = self.p[:-1]
+            p_l = self.p[:self.x.size-1]
             mean = (0.5 * (x_l + x_r) * (x_r - x_l) * p_l).sum()
         else:
             raise NotImplementedError('Can only compute mean for tabular '
@@ -952,7 +956,7 @@ class Tabular(Univariate):
 
     def normalize(self):
         """Normalize the probabilities stored on the distribution"""
-        self.p /= self.cdf().max()
+        self._p /= self.cdf().max()
 
     def sample(self, n_samples: int = 1, seed: typing.Optional[int] = None):
         rng = np.random.RandomState(seed)
@@ -1070,7 +1074,7 @@ class Tabular(Univariate):
             Integral of tabular distrbution
         """
         if self.interpolation == 'histogram':
-            return np.sum(np.diff(self.x) * self.p[:-1])
+            return np.sum(np.diff(self.x) * self.p[:self.x.size-1])
         elif self.interpolation == 'linear-linear':
             return trapezoid(self.p, self.x)
         else:
@@ -1346,7 +1350,7 @@ def combine_distributions(
     # Apply probabilites to continuous distributions
     for i in cont_index:
         dist = dist_list[i]
-        dist.p *= probs[i]
+        dist._p *= probs[i]
 
     if discrete_index:
         # Create combined discrete distribution
diff --git a/openmc/tallies.py b/openmc/tallies.py
index 04ae00b6a..f39ffa078 100644
--- a/openmc/tallies.py
+++ b/openmc/tallies.py
@@ -3209,6 +3209,7 @@ class Tallies(cv.CheckedList):
     def to_xml_element(self, memo=None):
         """Creates a 'tallies' element to be written to an XML file.
         """
+        memo = memo if memo is not None else set()
         element = ET.Element("tallies")
         self._create_mesh_subelements(element, memo)
         self._create_filter_subelements(element)
diff --git a/openmc/universe.py b/openmc/universe.py
index 2e86ab05a..9fab9ae51 100644
--- a/openmc/universe.py
+++ b/openmc/universe.py
@@ -90,7 +90,7 @@ class UniverseBase(ABC, IDManagerMixin):
         else:
             raise ValueError('No volume information found for this universe.')
 
-    def get_all_universes(self):
+    def get_all_universes(self, memo=None):
         """Return all universes that are contained within this one.
 
         Returns
@@ -100,10 +100,16 @@ class UniverseBase(ABC, IDManagerMixin):
             :class:`Universe` instances
 
         """
+        if memo is None:
+            memo = set()
+        elif self in memo:
+            return {}
+        memo.add(self)
+
         # Append all Universes within each Cell to the dictionary
         universes = {}
         for cell in self.get_all_cells().values():
-            universes.update(cell.get_all_universes())
+            universes.update(cell.get_all_universes(memo))
 
         return universes
 
@@ -639,15 +645,14 @@ class Universe(UniverseBase):
 
         """
 
-        cells = {}
-
-        if memo and self in memo:
-            return cells
-
-        if memo is not None:
-            memo.add(self)
+        if memo is None:
+            memo = set()
+        elif self in memo:
+            return {}
+        memo.add(self)
 
         # Add this Universe's cells to the dictionary
+        cells = {}
         cells.update(self._cells)
 
         # Append all Cells in each Cell in the Universe to the dictionary
@@ -667,6 +672,9 @@ class Universe(UniverseBase):
 
         """
 
+        if memo is None:
+            memo = set()
+
         materials = {}
 
         # Append all Cells in each Cell in the Universe to the dictionary
@@ -677,15 +685,17 @@ class Universe(UniverseBase):
         return materials
 
     def create_xml_subelement(self, xml_element, memo=None):
+        if memo is None:
+            memo = set()
+
         # Iterate over all Cells
         for cell in self._cells.values():
 
             # If the cell was already written, move on
-            if memo and cell in memo:
+            if cell in memo:
                 continue
 
-            if memo is not None:
-                memo.add(cell)
+            memo.add(cell)
 
             # Create XML subelement for this Cell
             cell_element = cell.create_xml_subelement(xml_element, memo)
@@ -928,11 +938,13 @@ class DAGMCUniverse(UniverseBase):
         return self._n_geom_elements('surface')
 
     def create_xml_subelement(self, xml_element, memo=None):
-        if memo and self in memo:
+        if memo is None:
+            memo = set()
+
+        if self in memo:
             return
 
-        if memo is not None:
-            memo.add(self)
+        memo.add(self)
 
         # Set xml element values
         dagmc_element = ET.Element('dagmc_universe')


Test patch:
diff --git a/tests/unit_tests/test_lattice.py b/tests/unit_tests/test_lattice.py
index 6fa32760e..d72d9c5c3 100644
--- a/tests/unit_tests/test_lattice.py
+++ b/tests/unit_tests/test_lattice.py
@@ -377,3 +377,46 @@ def test_unset_universes():
     hex_lattice.pitch = (1.,)
     with pytest.raises(ValueError):
         hex_lattice.create_xml_subelement(elem)
+
+
+@pytest.mark.parametrize("orientation", ['x', 'y'])
+def test_hex_lattice_roundtrip(orientation):
+    openmc.reset_auto_ids()
+
+    # ensure that the lattice universes are the same on all axial levels
+    def check_lattice_universes(og_lattice, xml_lattice):
+        for axial_og, axial_rt in zip(og_lattice.universes, xml_lattice.universes):
+            for ring_og, ring_rt in zip(axial_og, axial_rt):
+                assert [u.id for u in ring_og] == [u.id for u in ring_rt]
+
+    latt = openmc.HexLattice()
+    latt.pitch = (1.0, 1.0)
+    latt.center = (0.0, 0.0, 0.0)
+    latt.orientation = orientation
+
+    # fill the lattice with universes in increasing order and repeat for
+    # the second actial level
+    lvl_one_univs = [openmc.Universe(cells=[openmc.Cell()]) for _ in range(19)]
+    lvl_one_univs = [lvl_one_univs[-12:], lvl_one_univs[1:7], lvl_one_univs[:1]]
+    latt.universes = [lvl_one_univs, lvl_one_univs]
+
+    geom = openmc.Geometry([openmc.Cell(fill=latt)])
+    geom.export_to_xml()
+
+    xml_geom = openmc.Geometry.from_xml(materials=openmc.Materials())
+
+    xml_latt = xml_geom.get_all_lattices()[latt.id]
+
+    check_lattice_universes(latt, xml_latt)
+
+    # same test but with unique universes for each axial level
+    lvl_two_univs = [openmc.Universe(cells=[openmc.Cell()]) for _ in range(19)]
+    lvl_two_univs = [lvl_two_univs[-12:], lvl_two_univs[1:7], lvl_two_univs[:1]]
+    latt.universes = [lvl_one_univs, lvl_two_univs]
+
+    geom.export_to_xml()
+
+    xml_geom = openmc.Geometry.from_xml(materials=openmc.Materials())
+    xml_latt = xml_geom.get_all_lattices()[latt.id]
+
+    check_lattice_universes(latt, xml_latt)
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index c49930650..5bacd5fc6 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -132,6 +132,27 @@ def test_SphericalMesh_initiation():
     mesh.r_grid = (0, 11)
     assert (mesh.r_grid == np.array([0., 11.])).all()
 
+    # test invalid r_grid values
+    with pytest.raises(ValueError):
+        openmc.SphericalMesh(r_grid=[1, 1])
+
+    with pytest.raises(ValueError):
+        openmc.SphericalMesh(r_grid=[0])
+
+    # test invalid theta_grid values
+    with pytest.raises(ValueError):
+        openmc.SphericalMesh(r_grid=[1, 2], theta_grid=[1, 1])
+
+    with pytest.raises(ValueError):
+        openmc.SphericalMesh(r_grid=[1, 2], theta_grid=[0])
+
+    # test invalid phi_grid values
+    with pytest.raises(ValueError):
+        openmc.SphericalMesh(r_grid=[1, 2], phi_grid=[1, 1])
+
+    with pytest.raises(ValueError):
+        openmc.SphericalMesh(r_grid=[1, 2], phi_grid=[0])
+
     # waffles and pancakes are unfortunately not valid radii
     with pytest.raises(TypeError):
         openmc.SphericalMesh(('🧇', '🥞'))
diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py
index 761f26ab3..f5b467d03 100644
--- a/tests/unit_tests/test_stats.py
+++ b/tests/unit_tests/test_stats.py
@@ -92,7 +92,7 @@ def test_clip_discrete():
 
     with pytest.raises(ValueError):
         d.clip(-1.)
-    
+
     with pytest.raises(ValueError):
         d.clip(5)
 
@@ -188,8 +188,8 @@ def test_watt():
 
 
 def test_tabular():
-    x = np.array([0.0, 5.0, 7.0])
-    p = np.array([10.0, 20.0, 5.0])
+    x = np.array([0.0, 5.0, 7.0, 10.0])
+    p = np.array([10.0, 20.0, 5.0, 6.0])
     d = openmc.stats.Tabular(x, p, 'linear-linear')
     elem = d.to_xml_element('distribution')
 
@@ -217,6 +217,23 @@ def test_tabular():
     d.normalize()
     assert d.integral() == pytest.approx(1.0)
 
+    # ensure that passing a set of probabilities shorter than x works
+    # for histogram interpolation
+    d = openmc.stats.Tabular(x, p[:-1], interpolation='histogram')
+    d.cdf()
+    d.mean()
+    assert_sample_mean(d.sample(n_samples), d.mean())
+
+    # passing a shorter probability set should raise an error for linear-linear
+    with pytest.raises(ValueError):
+        d = openmc.stats.Tabular(x, p[:-1], interpolation='linear-linear')
+        d.cdf()
+
+    # Use probabilities of correct length for linear-linear interpolation and
+    # call the CDF method
+    d = openmc.stats.Tabular(x, p, interpolation='linear-linear')
+    d.cdf()
+
 
 def test_legendre():
     # Pu239 elastic scattering at 100 keV
diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py
index 19221212e..62ec18b32 100644
--- a/tests/unit_tests/test_surface_composite.py
+++ b/tests/unit_tests/test_surface_composite.py
@@ -158,18 +158,23 @@ def test_cone_one_sided(axis, point_pos, point_neg, ll_true):
 
 
 @pytest.mark.parametrize(
-    "axis, indices", [
-        ("X", [0, 1, 2]),
-        ("Y", [1, 2, 0]),
-        ("Z", [2, 0, 1]),
+    "axis, indices, center", [
+        ("X", [2, 0, 1], (0., 0.)),
+        ("Y", [0, 2, 1], (0., 0.)),
+        ("Z", [0, 1, 2], (0., 0.)),
+        ("X", [2, 0, 1], (10., 5.)),
+        ("Y", [0, 2, 1], (10., 5.)),
+        ("Z", [0, 1, 2], (10., 5.)),
+
     ]
 )
-def test_cylinder_sector(axis, indices):
+def test_cylinder_sector(axis, indices, center):
+    c1, c2 = center
     r1, r2 = 0.5, 1.5
     d = (r2 - r1) / 2
     phi1 = -60.
     phi2 = 60
-    s = openmc.model.CylinderSector(r1, r2, phi1, phi2,
+    s = openmc.model.CylinderSector(r1, r2, phi1, phi2, center=center,
                                     axis=axis.lower())
     assert isinstance(s.outer_cyl, getattr(openmc, axis + "Cylinder"))
     assert isinstance(s.inner_cyl, getattr(openmc, axis + "Cylinder"))
@@ -189,16 +194,18 @@ def test_cylinder_sector(axis, indices):
     assert np.all(np.isinf(ll))
     assert np.all(np.isinf(ur))
     ll, ur = (-s).bounding_box
-    assert ll == pytest.approx(np.roll([-np.inf, -r2, -r2], indices[0]))
-    assert ur == pytest.approx(np.roll([np.inf, r2, r2], indices[0]))
+    test_point_ll = np.array([-r2 + c1, -r2 + c2, -np.inf])
+    assert ll == pytest.approx(test_point_ll[indices])
+    test_point_ur = np.array([r2 + c1, r2 + c2, np.inf])
+    assert ur == pytest.approx(test_point_ur[indices])
 
     # __contains__ on associated half-spaces
-    point_pos = np.roll([0, r2 + 1, 0], indices[0])
-    assert point_pos in +s
-    assert point_pos not in -s
-    point_neg = np.roll([0, r1 + d, r1 + d], indices[0])
-    assert point_neg in -s
-    assert point_neg not in +s
+    point_pos = np.array([0 + c1, r2 + 1 + c2, 0])
+    assert point_pos[indices] in +s
+    assert point_pos[indices] not in -s
+    point_neg = np.array([r1 + d + c1, r1 + d + c2, 0])
+    assert point_neg[indices] in -s
+    assert point_neg[indices] not in +s
 
     # translate method
     t = uniform(-5.0, 5.0)
@@ -399,7 +406,7 @@ def test_polygon():
     with pytest.raises(ValueError):
         openmc.model.Polygon(rz_points)
 
-    # Test "M" shaped polygon 
+    # Test "M" shaped polygon
     points = np.array([[8.5151581, -17.988337],
                        [10.381711000000001, -17.988337],
                        [12.744357, -24.288728000000003],