From aa24f27d4d514893c028dfaef319e08b445675fe Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Thu, 16 Jul 2026 15:47:49 -0600
Subject: [PATCH 01/13] Remove mesh labeling from MeshSurfaceFilter

---
 include/openmc/tallies/filter_meshsurface.h | 15 -------
 src/tallies/filter_meshsurface.cpp          | 48 +++------------------
 2 files changed, 5 insertions(+), 58 deletions(-)

diff --git a/include/openmc/tallies/filter_meshsurface.h b/include/openmc/tallies/filter_meshsurface.h
index 195995c699f..fc6915200cb 100644
--- a/include/openmc/tallies/filter_meshsurface.h
+++ b/include/openmc/tallies/filter_meshsurface.h
@@ -22,21 +22,6 @@ class MeshSurfaceFilter : public MeshFilter {
   // Accessors
 
   void set_mesh(int32_t mesh) override;
-
-  enum class MeshDir {
-    OUT_LEFT,   // x min
-    IN_LEFT,    // x min
-    OUT_RIGHT,  // x max
-    IN_RIGHT,   // x max
-    OUT_BACK,   // y min
-    IN_BACK,    // y min
-    OUT_FRONT,  // y max
-    IN_FRONT,   // y max
-    OUT_BOTTOM, // z min
-    IN_BOTTOM,  // z min
-    OUT_TOP,    // z max
-    IN_TOP      // z max
-  };
 };
 
 } // namespace openmc
diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp
index b26cd198b32..46626977bf0 100644
--- a/src/tallies/filter_meshsurface.cpp
+++ b/src/tallies/filter_meshsurface.cpp
@@ -28,52 +28,14 @@ std::string MeshSurfaceFilter::text_label(int bin) const
   auto& mesh = *model::meshes[mesh_];
   int n_dim = mesh.n_dimension_;
 
-  // Get flattend mesh index and surface index.
+  // Get flattened mesh index and surface index.
   int i_mesh = bin / (4 * n_dim);
-  MeshDir surf_dir = static_cast<MeshDir>(bin % (4 * n_dim));
+  int surf_index = bin % (4 * n_dim);
 
-  // Get mesh index part of label.
+  // Get mesh index part of label, then append the surface part. 
+  // The surface is labeled by the underlying mesh.
   std::string out = MeshFilter::text_label(i_mesh);
-
-  // Get surface part of label.
-  switch (surf_dir) {
-  case MeshDir::OUT_LEFT:
-    out += " Outgoing, x-min";
-    break;
-  case MeshDir::IN_LEFT:
-    out += " Incoming, x-min";
-    break;
-  case MeshDir::OUT_RIGHT:
-    out += " Outgoing, x-max";
-    break;
-  case MeshDir::IN_RIGHT:
-    out += " Incoming, x-max";
-    break;
-  case MeshDir::OUT_BACK:
-    out += " Outgoing, y-min";
-    break;
-  case MeshDir::IN_BACK:
-    out += " Incoming, y-min";
-    break;
-  case MeshDir::OUT_FRONT:
-    out += " Outgoing, y-max";
-    break;
-  case MeshDir::IN_FRONT:
-    out += " Incoming, y-max";
-    break;
-  case MeshDir::OUT_BOTTOM:
-    out += " Outgoing, z-min";
-    break;
-  case MeshDir::IN_BOTTOM:
-    out += " Incoming, z-min";
-    break;
-  case MeshDir::OUT_TOP:
-    out += " Outgoing, z-max";
-    break;
-  case MeshDir::IN_TOP:
-    out += " Incoming, z-max";
-    break;
-  }
+  out += mesh.surface_bin_label(surf_index);
 
   return out;
 }

From b416a938b64e69cb327d43186b5c211ac3c35357 Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Thu, 16 Jul 2026 16:04:30 -0600
Subject: [PATCH 02/13] Add the labeling to Mesh and its subclasses

---
 include/openmc/mesh.h | 10 ++++++++++
 src/mesh.cpp          | 29 +++++++++++++++++++++++++++++
 2 files changed, 39 insertions(+)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 0d8189caa1d..680c3b26d7c 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -236,6 +236,12 @@ class Mesh {
   //! \param[in] bin Mesh bin to generate a label for
   virtual std::string bin_label(int bin) const = 0;
 
+  //! Axis names (per dimension) used when labeling surface tally bins
+  virtual std::array<const char*, 3> axis_labels() const;
+
+  //! Build the surface component of a mesh surface tally bin label
+  std::string surface_bin_label(int surf_index) const;
+
   //! Get the volume of a mesh bin
   //
   //! \param[in] bin Bin to return the volume for
@@ -588,6 +594,8 @@ class CylindricalMesh : public PeriodicStructuredMesh {
 
   static const std::string mesh_type;
 
+  std::array<const char*, 3> axis_labels() const override;
+
   Position sample_element(const MeshIndex& ijk, uint64_t* seed) const override;
 
   MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
@@ -653,6 +661,8 @@ class SphericalMesh : public PeriodicStructuredMesh {
 
   static const std::string mesh_type;
 
+  std::array<const char*, 3> axis_labels() const override;
+
   Position sample_element(const MeshIndex& ijk, uint64_t* seed) const override;
 
   MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 89d1ff24b82..48c3acf40e5 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -430,6 +430,25 @@ vector<double> Mesh::volumes() const
   return volumes;
 }
 
+//! Default (Cartesian) axis labels used for surface bin labels.
+std::array<const char*, 3> Mesh::axis_labels() const
+{
+  return {"x", "y", "z"};
+}
+
+//! Build the surface component of a mesh surface tally bin label.
+//! surf_index: 0=out/min, 1=in/min, 2=out/max, 3=in/max
+std::string Mesh::surface_bin_label(int surf_index) const
+{
+  auto labels = this->axis_labels();
+  int dim = surf_index / 4;
+  int code = surf_index % 4;
+  bool incoming = (code == 1) || (code == 3);
+  bool max = (code == 2) || (code == 3);
+  return fmt::format(" {}, {}-{}", incoming ? "Incoming" : "Outgoing",
+    labels[dim], max ? "max" : "min");
+}
+
 void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
   int32_t* materials, double* volumes) const
 {
@@ -1815,6 +1834,11 @@ std::string CylindricalMesh::get_mesh_type() const
   return mesh_type;
 }
 
+std::array<const char*, 3> CylindricalMesh::axis_labels() const
+{
+  return {"r", "phi", "z"};
+}
+
 StructuredMesh::MeshIndex CylindricalMesh::get_indices(
   Position r, bool& in_mesh) const
 {
@@ -2108,6 +2132,11 @@ std::string SphericalMesh::get_mesh_type() const
   return mesh_type;
 }
 
+std::array<const char*, 3> SphericalMesh::axis_labels() const
+{
+  return {"r", "theta", "phi"};
+}
+
 StructuredMesh::MeshIndex SphericalMesh::get_indices(
   Position r, bool& in_mesh) const
 {

From 5f78fd9594a18230d407e67501ddd52bfd1b2d24 Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Thu, 16 Jul 2026 16:20:05 -0600
Subject: [PATCH 03/13] Propagate new names to openmc.filter and expose in
 openmc.mesh

---
 openmc/filter.py | 74 +++++++++++++++++++++++++-----------------------
 openmc/mesh.py   | 28 ++++++++++++++----
 2 files changed, 60 insertions(+), 42 deletions(-)

diff --git a/openmc/filter.py b/openmc/filter.py
index 87aeb70c3a7..b16f485df93 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -30,12 +30,19 @@
     'weight', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction',
 )
 
-_CURRENT_NAMES = (
-    'x-min out', 'x-min in', 'x-max out', 'x-max in',
-    'y-min out', 'y-min in', 'y-max out', 'y-max in',
-    'z-min out', 'z-min in', 'z-max out', 'z-max in'
-)
+def _mesh_current_names(mesh):
+    """Return surface-crossing bin names derived from a mesh's axis labels.
 
+    For each axis (in dimension order) the four crossing names are produced as
+    "<axis>-min out", "<axis>-min in", "<axis>-max out", "<axis>-max in". This
+    yields x/y/z names for Cartesian meshes, r/phi/z for cylindrical, and
+    r/theta/phi for spherical.
+    """
+    names = []
+    for axis in mesh.axis_labels:
+        names += [f'{axis}-min out', f'{axis}-min in',
+                  f'{axis}-max out', f'{axis}-max in']
+    return names
 
 
 class FilterMeta(ABCMeta):
@@ -996,7 +1003,7 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs):
         idx_start = 0 if isinstance(self.mesh, openmc.UnstructuredMesh) else 1
 
         # Generate a multi-index sub-column for each axis
-        for label, dim_size in zip(self.mesh._axis_labels, self.mesh.dimension):
+        for label, dim_size in zip(self.mesh.axis_labels, self.mesh.dimension):
             filter_dict[mesh_key, label] = _repeat_and_tile(
                 np.arange(idx_start, idx_start + dim_size), stride, data_size)
             stride *= dim_size
@@ -1273,7 +1280,8 @@ class MeshSurfaceFilter(MeshFilter):
         Unique identifier for the filter
     bins : list of tuple
         A list of mesh indices / surfaces for each filter bin, e.g. [(1, 1,
-        'x-min out'), (1, 1, 'x-min in'), ...]
+        'x-min out'), (1, 1, 'x-min in'), ...]. Surface names use the mesh's
+        axis labels (e.g. r/phi/z for a cylindrical mesh).
     num_bins : Integral
         The number of filter bins
 
@@ -1287,10 +1295,12 @@ def mesh(self, mesh):
         cv.check_type('filter mesh', mesh, openmc.MeshBase)
         self._mesh = mesh
 
-        # Take the product of mesh indices and current names
-        n_dim = mesh.n_dimension
+        # Take the product of mesh indices and surface-crossing names, using
+        # the mesh's own axis labels so cylindrical/spherical meshes get the
+        # correct names.
+        current_names = _mesh_current_names(mesh)
         self.bins = [mesh_tuple + (surf,) for mesh_tuple, surf in
-                     product(mesh.indices, _CURRENT_NAMES[:4*n_dim])]
+                     product(mesh.indices, current_names)]
 
     def get_pandas_dataframe(self, data_size, stride, **kwargs):
         """Builds a Pandas DataFrame for the Filter's bins.
@@ -1309,9 +1319,11 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs):
         Returns
         -------
         pandas.DataFrame
-            A Pandas DataFrame with three columns describing the x,y,z mesh
-            cell indices corresponding to each filter bin.  The number of rows
-            in the DataFrame is the same as the total number of bins in the
+            A Pandas DataFrame with columns describing the mesh cell indices
+            corresponding to each filter bin, plus a surface column. Column
+            names depend on the mesh type (e.g., x/y/z for RegularMesh, r/phi/z
+            for CylindricalMesh, r/theta/phi for SphericalMesh). The number of
+            rows in the DataFrame is the same as the total number of bins in the
             corresponding tally, with the filter bin appropriately tiled to map
             to the corresponding tally bins.
 
@@ -1329,34 +1341,24 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs):
         # Append mesh ID as outermost index of multi-index
         mesh_key = f'mesh {self.mesh.id}'
 
-        # Find mesh dimensions - use 3D indices for simplicity
-        n_surfs = 4 * len(self.mesh.dimension)
-        if len(self.mesh.dimension) == 3:
-            nx, ny, nz = self.mesh.dimension
-        elif len(self.mesh.dimension) == 2:
-            nx, ny = self.mesh.dimension
-            nz = 1
-        else:
-            nx = self.mesh.dimension
-            ny = nz = 1
-
-        # Generate multi-index sub-column for x-axis
-        filter_dict[mesh_key, 'x'] = _repeat_and_tile(
-            np.arange(1, nx + 1), n_surfs * stride, data_size)
+        # Number of surface-crossing bins per mesh element
+        dims = self.mesh.dimension
+        n_surfs = 4 * len(dims)
 
-        # Generate multi-index sub-column for y-axis
-        if len(self.mesh.dimension) > 1:
-            filter_dict[mesh_key, 'y'] = _repeat_and_tile(
-                np.arange(1, ny + 1), n_surfs * nx * stride, data_size)
+        # Surface-crossing names derived from the mesh's axis labels
+        current_names = _mesh_current_names(self.mesh)
 
-        # Generate multi-index sub-column for z-axis
-        if len(self.mesh.dimension) > 2:
-            filter_dict[mesh_key, 'z'] = _repeat_and_tile(
-                np.arange(1, nz + 1), n_surfs * nx * ny * stride, data_size)
+        # Generate a multi-index sub-column for each axis, using the mesh's own
+        # axis labels so curvilinear meshes are labeled correctly.
+        axis_stride = n_surfs * stride
+        for label, dim_size in zip(self.mesh.axis_labels, dims):
+            filter_dict[mesh_key, label] = _repeat_and_tile(
+                np.arange(1, dim_size + 1), axis_stride, data_size)
+            axis_stride *= dim_size
 
         # Generate multi-index sub-column for surface
         filter_dict[mesh_key, 'surf'] = _repeat_and_tile(
-            _CURRENT_NAMES[:n_surfs], stride, data_size)
+            current_names, stride, data_size)
 
         # Initialize a Pandas DataFrame from the mesh dictionary
         return pd.concat([df, pd.DataFrame(filter_dict)])
diff --git a/openmc/mesh.py b/openmc/mesh.py
index ff8e1457724..bcb3ad2538e 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -259,6 +259,12 @@ def indices(self):
     def n_elements(self):
         pass
 
+    @property
+    @abstractmethod
+    def axis_labels(self):
+        """tuple of str : Names of the mesh axes, one per dimension."""
+        pass
+
     def __repr__(self):
         string = type(self).__name__ + '\n'
         string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
@@ -538,7 +544,8 @@ def n_dimension(self):
 
     @property
     @abstractmethod
-    def _axis_labels(self):
+    def axis_labels(self):
+        """tuple of str : Names of the mesh axes, one per dimension."""
         pass
 
     @property
@@ -1042,6 +1049,9 @@ class RegularMesh(StructuredMesh):
         The number of mesh cells in each direction (x, y, z).
     n_dimension : int
         Number of mesh dimensions.
+    axis_labels : tuple of str
+        Names of the mesh axes ('x', 'y', 'z'), truncated to the mesh
+        dimensionality.
     lower_left : Iterable of float
         The lower-left corner of the structured mesh. If only two coordinate
         are given, it is assumed that the mesh is an x-y mesh.
@@ -1085,7 +1095,7 @@ def n_dimension(self):
             return None
 
     @property
-    def _axis_labels(self):
+    def axis_labels(self):
         return ('x', 'y', 'z')[:self.n_dimension]
 
     @property
@@ -1569,6 +1579,8 @@ class RectilinearMesh(StructuredMesh):
         The number of mesh cells in each direction (x, y, z).
     n_dimension : int
         Number of mesh dimensions (always 3 for a RectilinearMesh).
+    axis_labels : tuple of str
+        Names of the mesh axes ('x', 'y', 'z').
     x_grid : numpy.ndarray
         1-D array of mesh boundary points along the x-axis.
     y_grid : numpy.ndarray
@@ -1602,7 +1614,7 @@ def n_dimension(self):
         return 3
 
     @property
-    def _axis_labels(self):
+    def axis_labels(self):
         return ('x', 'y', 'z')
 
     @property
@@ -1870,6 +1882,8 @@ class CylindricalMesh(StructuredMesh):
         The number of mesh cells in each direction (r_grid, phi_grid, z_grid).
     n_dimension : int
         Number of mesh dimensions (always 3 for a CylindricalMesh).
+    axis_labels : tuple of str
+        Names of the mesh axes ('r', 'phi', 'z').
     r_grid : numpy.ndarray
         1-D array of mesh boundary points along the r-axis.
         Requirement is r >= 0.
@@ -1924,7 +1938,7 @@ def n_dimension(self):
         return 3
 
     @property
-    def _axis_labels(self):
+    def axis_labels(self):
         return ('r', 'phi', 'z')
 
     @property
@@ -2307,6 +2321,8 @@ class SphericalMesh(StructuredMesh):
         theta_grid, phi_grid).
     n_dimension : int
         Number of mesh dimensions (always 3 for a SphericalMesh).
+    axis_labels : tuple of str
+        Names of the mesh axes ('r', 'theta', 'phi').
     r_grid : numpy.ndarray
         1-D array of mesh boundary points along the r-axis.
         Requirement is r >= 0.
@@ -2361,7 +2377,7 @@ def n_dimension(self):
         return 3
 
     @property
-    def _axis_labels(self):
+    def axis_labels(self):
         return ('r', 'theta', 'phi')
 
     @property
@@ -2949,7 +2965,7 @@ def n_dimension(self):
         return 3
 
     @property
-    def _axis_labels(self):
+    def axis_labels(self):
         return ('element_index',)
 
     @property

From a9c046cdebffc209eab510c5766840334f17971c Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Thu, 16 Jul 2026 17:05:53 -0600
Subject: [PATCH 04/13] Add bin label unit test

---
 tests/unit_tests/test_mesh.py | 40 +++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 9b1469fc590..f63901f3883 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -267,6 +267,46 @@ def test_centroids():
     np.testing.assert_array_almost_equal(mesh.centroids[0, 0, 0], [x-5.0, y-5.0, z+5.0])
 
 
+@pytest.mark.parametrize('mesh_type', ('regular', 'rectilinear', 'cylindrical', 'spherical'))
+def test_axis_labels(mesh_type):
+    if mesh_type == 'regular':
+        mesh = openmc.RegularMesh()
+        mesh.lower_left = np.asarray([0.0]*3]
+        mesh.width = np.asarray([0.0]*3]
+        mesh.dimension = (1, 1, 1)
+    elif mesh_type == 'rectilinear':
+        mesh = openmc.RectilinearMesh()
+        mesh.x_grid = np.asarray([0.0, 1.0])
+        mesh.y_grid = np.asarray([0.0, 1.0])
+        mesh.z_grid = np.asarray([0.0, 1.0])
+    elif mesh_type == 'cylindrical':
+        r_grid = np.asarray([0.0, 1.0])
+        z_grid = np.asarray([0.0, 1.0])
+        p_grid = np.asarray([0.0, 1.0])
+        mesh = openmc.CylindricalMesh(r_grid=r_grid, z_grid=z_grid, phi_grid=p_grid)
+    elif mesh_type == 'spherical':
+        r_grid = np.asarray([0.0, 1.0])
+        t_grid = np.asarray([0.0, 1.0])
+        p_grid = np.asarray([0.0, 1.0])
+        mesh = openmc.SphericalMesh(r_grid=r_grid, theta_grid=t_grid, phi_grid=p_grid)
+    else:
+        raise ValueError(mesh_type)
+    
+    filt = openmc.MeshSurfaceFilter(mesh)
+    assert len(filt.bins) == 12
+    bin_names = [b[3] for b in filt.bins]
+    if mesh_type in {'regular', 'rectilinear'}:
+        assert bin_names[:8] == ['x-min out', 'x-min in', 'x-max out', 'x-max in', 
+                                 'y-min out', 'y-min in', 'y-max out', 'y-max in']
+    if mesh_type in {'regular', 'rectilinear', 'cylindrical'}:
+        assert bin_names[8:] == ['z-min out', 'z-min in', 'z-max out', 'z-max in']
+    if mesh_type in {'cylindrical', 'spherical'}:
+        assert bin_names[:4] == ['r-min out', 'r-min in', 'r-max out', 'r-max in']
+        assert bin_names[4:8] == ['phi-min out', 'phi-min in', 'phi-max out', 'phi-max in']
+    if mesh_type == 'spherical':
+        assert bin_names[8:] == ['theta-min out', 'theta-min in', 'theta-max out', 'theta-max in']
+   
+
 @pytest.mark.parametrize('mesh_type', ('regular', 'rectilinear', 'cylindrical', 'spherical'))
 def test_mesh_vertices(mesh_type):
 

From 6a255acfb173eb39cd2079ed372a6600e4d410e8 Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Thu, 16 Jul 2026 17:06:18 -0600
Subject: [PATCH 05/13] index surf

---
 openmc/filter.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/openmc/filter.py b/openmc/filter.py
index b16f485df93..674aa2803e9 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -1358,7 +1358,7 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs):
 
         # Generate multi-index sub-column for surface
         filter_dict[mesh_key, 'surf'] = _repeat_and_tile(
-            current_names, stride, data_size)
+            current_names[:n_surfs], stride, data_size)
 
         # Initialize a Pandas DataFrame from the mesh dictionary
         return pd.concat([df, pd.DataFrame(filter_dict)])

From 2fa461ffcbd87010ed1706cf867c8594c92636a9 Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Fri, 17 Jul 2026 09:51:58 -0600
Subject: [PATCH 06/13] Regression test for axis labels

---
 tests/regression_tests/filter_mesh/test.py | 43 ++++++++++++++++++++++
 1 file changed, 43 insertions(+)

diff --git a/tests/regression_tests/filter_mesh/test.py b/tests/regression_tests/filter_mesh/test.py
index 165ba2a0c09..4a9f1518919 100644
--- a/tests/regression_tests/filter_mesh/test.py
+++ b/tests/regression_tests/filter_mesh/test.py
@@ -7,6 +7,49 @@
 from tests.testing_harness import HashedPyAPITestHarness
 
 
+class MeshSurfaceFilterTest(HashedPyAPITestHarness):
+
+    def _compare_results(self):
+        tally_meshes = {}
+        with openmc.StatePoint(self.statepoint_name) as sp:
+            for tally in sp.tallies.values():
+                for filt in tally.filters:
+                    if isinstance(filt, openmc.MeshSurfaceFilter):
+                        tally_meshes[tally.id] = filt.mesh
+                        break
+
+        with open('tallies.out') as fh:
+            text = fh.read()
+
+        for tally_id, mesh in tally_meshes.items():
+            # Snip the relevant text.
+            header = f"TALLY {tally_id}"
+            assert header in text
+            text = text[text.find(header) + len(header):]
+            end = text.find("TALLY")
+            section = text if end == -1 else text[:end]
+            # Define the relevant labels and ensure they wrap.
+            nd = mesh.n_dimension
+            labels = [ln for ln in section.splitlines() if "Mesh Index" in ln]
+            labels = labels[:4*nd + 1]  # just the first 4, 8, or 12, plus one
+            expected = []
+            if isinstance(mesh, openmc.SphericalMesh):
+                axis_labels = ('r', 'phi', 'theta')
+            elif isinstance(mesh, openmc.CylindricalMesh):
+                axis_labels = ('r', 'phi', 'z')
+            else:
+                axis_labels = ('x', 'y', 'z')
+            for axis in axis_labels[:nd]:
+            expected += [f'Outgoing, {axis}-min', f'Incoming, {axis}-min',
+                         f'Outgoing, {axis}-max', f'Incoming, {axis}-max']
+        expected.append(expected[0])
+        assert len(labels) == len(expected)
+        for line, exp in zip(labels, expected):
+            assert exp in line, f"{exp} not in {line}"
+
+        return super()._compare_results()
+
+
 @pytest.fixture
 def model():
     model = openmc.model.Model()

From fd42b20de10fe1bd4dfc78789ca675b0272c1a26 Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Fri, 17 Jul 2026 09:59:20 -0600
Subject: [PATCH 07/13] Corrected spherical mesh coordinate order

---
 tests/regression_tests/filter_mesh/test.py | 2 +-
 tests/unit_tests/test_mesh.py              | 6 ++++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/tests/regression_tests/filter_mesh/test.py b/tests/regression_tests/filter_mesh/test.py
index 4a9f1518919..e53570ff487 100644
--- a/tests/regression_tests/filter_mesh/test.py
+++ b/tests/regression_tests/filter_mesh/test.py
@@ -34,7 +34,7 @@ def _compare_results(self):
             labels = labels[:4*nd + 1]  # just the first 4, 8, or 12, plus one
             expected = []
             if isinstance(mesh, openmc.SphericalMesh):
-                axis_labels = ('r', 'phi', 'theta')
+                axis_labels = ('r', 'theta', 'phi')
             elif isinstance(mesh, openmc.CylindricalMesh):
                 axis_labels = ('r', 'phi', 'z')
             else:
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index f63901f3883..12df74ae6af 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -302,9 +302,11 @@ def test_axis_labels(mesh_type):
         assert bin_names[8:] == ['z-min out', 'z-min in', 'z-max out', 'z-max in']
     if mesh_type in {'cylindrical', 'spherical'}:
         assert bin_names[:4] == ['r-min out', 'r-min in', 'r-max out', 'r-max in']
-        assert bin_names[4:8] == ['phi-min out', 'phi-min in', 'phi-max out', 'phi-max in']
     if mesh_type == 'spherical':
-        assert bin_names[8:] == ['theta-min out', 'theta-min in', 'theta-max out', 'theta-max in']
+        assert bin_names[4:] == ['theta-min out', 'theta-min in', 'theta-max out', 'theta-max in',
+                                 'phi-min out', 'phi-min in', 'phi-max out', 'phi-max in']
+    if mesh_type == 'cylindrical':
+        assert bin_names[4:8] == ['phi-min out', 'phi-min in', 'phi-max out', 'phi-max in']
    
 
 @pytest.mark.parametrize('mesh_type', ('regular', 'rectilinear', 'cylindrical', 'spherical'))

From 60a2ff1fc181c051e54d42f0d65b18da0871b31c Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Fri, 17 Jul 2026 10:01:16 -0600
Subject: [PATCH 08/13] This belongs in test_filter_mesh

---
 tests/unit_tests/test_filter_mesh.py | 42 ++++++++++++++++++++++++++++
 tests/unit_tests/test_mesh.py        | 42 ----------------------------
 2 files changed, 42 insertions(+), 42 deletions(-)

diff --git a/tests/unit_tests/test_filter_mesh.py b/tests/unit_tests/test_filter_mesh.py
index e26bea337bf..5f68de293d7 100644
--- a/tests/unit_tests/test_filter_mesh.py
+++ b/tests/unit_tests/test_filter_mesh.py
@@ -410,3 +410,45 @@ def test_mesh_filter_dataframe_rectilinear():
     assert (mesh_key, 'y') in df.columns
     assert (mesh_key, 'z') in df.columns
     assert len(df) == data_size
+
+
+@pytest.mark.parametrize('mesh_type', ('regular', 'rectilinear', 'cylindrical', 'spherical'))
+def test_axis_labels(mesh_type):
+    if mesh_type == 'regular':
+        mesh = openmc.RegularMesh()
+        mesh.lower_left = np.asarray([0.0]*3]
+        mesh.width = np.asarray([0.0]*3]
+        mesh.dimension = (1, 1, 1)
+    elif mesh_type == 'rectilinear':
+        mesh = openmc.RectilinearMesh()
+        mesh.x_grid = np.asarray([0.0, 1.0])
+        mesh.y_grid = np.asarray([0.0, 1.0])
+        mesh.z_grid = np.asarray([0.0, 1.0])
+    elif mesh_type == 'cylindrical':
+        r_grid = np.asarray([0.0, 1.0])
+        z_grid = np.asarray([0.0, 1.0])
+        p_grid = np.asarray([0.0, 1.0])
+        mesh = openmc.CylindricalMesh(r_grid=r_grid, z_grid=z_grid, phi_grid=p_grid)
+    elif mesh_type == 'spherical':
+        r_grid = np.asarray([0.0, 1.0])
+        t_grid = np.asarray([0.0, 1.0])
+        p_grid = np.asarray([0.0, 1.0])
+        mesh = openmc.SphericalMesh(r_grid=r_grid, theta_grid=t_grid, phi_grid=p_grid)
+    else:
+        raise ValueError(mesh_type)
+    
+    filt = openmc.MeshSurfaceFilter(mesh)
+    assert len(filt.bins) == 12
+    bin_names = [b[3] for b in filt.bins]
+    if mesh_type in {'regular', 'rectilinear'}:
+        assert bin_names[:8] == ['x-min out', 'x-min in', 'x-max out', 'x-max in', 
+                                 'y-min out', 'y-min in', 'y-max out', 'y-max in']
+    if mesh_type in {'regular', 'rectilinear', 'cylindrical'}:
+        assert bin_names[8:] == ['z-min out', 'z-min in', 'z-max out', 'z-max in']
+    if mesh_type in {'cylindrical', 'spherical'}:
+        assert bin_names[:4] == ['r-min out', 'r-min in', 'r-max out', 'r-max in']
+    if mesh_type == 'spherical':
+        assert bin_names[4:] == ['theta-min out', 'theta-min in', 'theta-max out', 'theta-max in',
+                                 'phi-min out', 'phi-min in', 'phi-max out', 'phi-max in']
+    if mesh_type == 'cylindrical':
+        assert bin_names[4:8] == ['phi-min out', 'phi-min in', 'phi-max out', 'phi-max in']
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 12df74ae6af..9b1469fc590 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -267,48 +267,6 @@ def test_centroids():
     np.testing.assert_array_almost_equal(mesh.centroids[0, 0, 0], [x-5.0, y-5.0, z+5.0])
 
 
-@pytest.mark.parametrize('mesh_type', ('regular', 'rectilinear', 'cylindrical', 'spherical'))
-def test_axis_labels(mesh_type):
-    if mesh_type == 'regular':
-        mesh = openmc.RegularMesh()
-        mesh.lower_left = np.asarray([0.0]*3]
-        mesh.width = np.asarray([0.0]*3]
-        mesh.dimension = (1, 1, 1)
-    elif mesh_type == 'rectilinear':
-        mesh = openmc.RectilinearMesh()
-        mesh.x_grid = np.asarray([0.0, 1.0])
-        mesh.y_grid = np.asarray([0.0, 1.0])
-        mesh.z_grid = np.asarray([0.0, 1.0])
-    elif mesh_type == 'cylindrical':
-        r_grid = np.asarray([0.0, 1.0])
-        z_grid = np.asarray([0.0, 1.0])
-        p_grid = np.asarray([0.0, 1.0])
-        mesh = openmc.CylindricalMesh(r_grid=r_grid, z_grid=z_grid, phi_grid=p_grid)
-    elif mesh_type == 'spherical':
-        r_grid = np.asarray([0.0, 1.0])
-        t_grid = np.asarray([0.0, 1.0])
-        p_grid = np.asarray([0.0, 1.0])
-        mesh = openmc.SphericalMesh(r_grid=r_grid, theta_grid=t_grid, phi_grid=p_grid)
-    else:
-        raise ValueError(mesh_type)
-    
-    filt = openmc.MeshSurfaceFilter(mesh)
-    assert len(filt.bins) == 12
-    bin_names = [b[3] for b in filt.bins]
-    if mesh_type in {'regular', 'rectilinear'}:
-        assert bin_names[:8] == ['x-min out', 'x-min in', 'x-max out', 'x-max in', 
-                                 'y-min out', 'y-min in', 'y-max out', 'y-max in']
-    if mesh_type in {'regular', 'rectilinear', 'cylindrical'}:
-        assert bin_names[8:] == ['z-min out', 'z-min in', 'z-max out', 'z-max in']
-    if mesh_type in {'cylindrical', 'spherical'}:
-        assert bin_names[:4] == ['r-min out', 'r-min in', 'r-max out', 'r-max in']
-    if mesh_type == 'spherical':
-        assert bin_names[4:] == ['theta-min out', 'theta-min in', 'theta-max out', 'theta-max in',
-                                 'phi-min out', 'phi-min in', 'phi-max out', 'phi-max in']
-    if mesh_type == 'cylindrical':
-        assert bin_names[4:8] == ['phi-min out', 'phi-min in', 'phi-max out', 'phi-max in']
-   
-
 @pytest.mark.parametrize('mesh_type', ('regular', 'rectilinear', 'cylindrical', 'spherical'))
 def test_mesh_vertices(mesh_type):
 

From 432298e7403570867e1c2d5c06f710a2ef928ccf Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Fri, 17 Jul 2026 10:04:56 -0600
Subject: [PATCH 09/13] clang-format

---
 src/mesh.cpp                       | 2 +-
 src/tallies/filter_meshsurface.cpp | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/mesh.cpp b/src/mesh.cpp
index 48c3acf40e5..7b5db5a2b52 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -3846,7 +3846,7 @@ void LibMesh::set_score_data(const std::string& var_name,
   unsigned int std_dev_num = variable_map_.at(std_dev_name);
 
   for (auto it = m_->local_elements_begin(); it != m_->local_elements_end();
-       it++) {
+    it++) {
     if (!(*it)->active()) {
       continue;
     }
diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp
index 46626977bf0..c7205aad344 100644
--- a/src/tallies/filter_meshsurface.cpp
+++ b/src/tallies/filter_meshsurface.cpp
@@ -32,7 +32,7 @@ std::string MeshSurfaceFilter::text_label(int bin) const
   int i_mesh = bin / (4 * n_dim);
   int surf_index = bin % (4 * n_dim);
 
-  // Get mesh index part of label, then append the surface part. 
+  // Get mesh index part of label, then append the surface part.
   // The surface is labeled by the underlying mesh.
   std::string out = MeshFilter::text_label(i_mesh);
   out += mesh.surface_bin_label(surf_index);

From 71cb56619a8b4dc99f6b4243426c079d1a38f3b7 Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Fri, 17 Jul 2026 10:15:40 -0600
Subject: [PATCH 10/13] clang-format==18.1.3

---
 src/mesh.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mesh.cpp b/src/mesh.cpp
index 7b5db5a2b52..48c3acf40e5 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -3846,7 +3846,7 @@ void LibMesh::set_score_data(const std::string& var_name,
   unsigned int std_dev_num = variable_map_.at(std_dev_name);
 
   for (auto it = m_->local_elements_begin(); it != m_->local_elements_end();
-    it++) {
+       it++) {
     if (!(*it)->active()) {
       continue;
     }

From b77bec0867c9e1f4255d2b68bb17fd1d62cbb79a Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Fri, 17 Jul 2026 10:41:25 -0600
Subject: [PATCH 11/13] commit local changes

---
 tests/regression_tests/filter_mesh/test.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/tests/regression_tests/filter_mesh/test.py b/tests/regression_tests/filter_mesh/test.py
index e53570ff487..6e0c08f6cec 100644
--- a/tests/regression_tests/filter_mesh/test.py
+++ b/tests/regression_tests/filter_mesh/test.py
@@ -40,9 +40,9 @@ def _compare_results(self):
             else:
                 axis_labels = ('x', 'y', 'z')
             for axis in axis_labels[:nd]:
-            expected += [f'Outgoing, {axis}-min', f'Incoming, {axis}-min',
-                         f'Outgoing, {axis}-max', f'Incoming, {axis}-max']
-        expected.append(expected[0])
+                expected += [f'Outgoing, {axis}-min', f'Incoming, {axis}-min',
+                             f'Outgoing, {axis}-max', f'Incoming, {axis}-max']
+            expected.append(expected[0])
         assert len(labels) == len(expected)
         for line, exp in zip(labels, expected):
             assert exp in line, f"{exp} not in {line}"
@@ -161,5 +161,5 @@ def model():
 
 
 def test_filter_mesh(model):
-    harness = HashedPyAPITestHarness('statepoint.5.h5', model)
+    harness = MeshSurfaceFilterTest('statepoint.5.h5', model)
     harness.main()

From eefe563f08e44c2cfa3b05dba59fcefed53da4c8 Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Fri, 17 Jul 2026 10:59:42 -0600
Subject: [PATCH 12/13] how did that even happen

---
 tests/unit_tests/test_filter_mesh.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/unit_tests/test_filter_mesh.py b/tests/unit_tests/test_filter_mesh.py
index 5f68de293d7..fa16f765eba 100644
--- a/tests/unit_tests/test_filter_mesh.py
+++ b/tests/unit_tests/test_filter_mesh.py
@@ -416,8 +416,8 @@ def test_mesh_filter_dataframe_rectilinear():
 def test_axis_labels(mesh_type):
     if mesh_type == 'regular':
         mesh = openmc.RegularMesh()
-        mesh.lower_left = np.asarray([0.0]*3]
-        mesh.width = np.asarray([0.0]*3]
+        mesh.lower_left = np.asarray([0.0]*3)
+        mesh.width = np.asarray([0.0]*3)
         mesh.dimension = (1, 1, 1)
     elif mesh_type == 'rectilinear':
         mesh = openmc.RectilinearMesh()

From 6ee22ed06141a3c6f8d97fde9877ae93914ab675 Mon Sep 17 00:00:00 2001
From: "Labossiere-Hickman, Travis James" <Travis.LabossiereHickman@inl.gov>
Date: Fri, 17 Jul 2026 13:23:48 -0600
Subject: [PATCH 13/13] Propagate attribute to openmc.mgxs

---
 openmc/mgxs/mdgxs.py | 2 +-
 openmc/mgxs/mgxs.py  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/openmc/mgxs/mdgxs.py b/openmc/mgxs/mdgxs.py
index 57cc955ba27..dcf9832eaa0 100644
--- a/openmc/mgxs/mdgxs.py
+++ b/openmc/mgxs/mdgxs.py
@@ -877,7 +877,7 @@ def get_pandas_dataframe(self, groups='all', nuclides='all',
         # energy groups such that data is from fast to thermal
         if self.domain_type == 'mesh':
             mesh_str = f'mesh {self.domain.id}'
-            mesh_cols = [(mesh_str, label) for label in self.domain._axis_labels]
+            mesh_cols = [(mesh_str, label) for label in self.domain.axis_labels]
             df.sort_values(by=mesh_cols + columns, inplace=True)
         else:
             df.sort_values(by=[self.domain_type] + columns, inplace=True)
diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index 4c13c75087c..7f0200c0976 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -2134,7 +2134,7 @@ def get_pandas_dataframe(self, groups='all', nuclides='all',
         # energy groups such that data is from fast to thermal
         if self.domain_type == 'mesh':
             mesh_str = f'mesh {self.domain.id}'
-            mesh_cols = [(mesh_str, label) for label in self.domain._axis_labels]
+            mesh_cols = [(mesh_str, label) for label in self.domain.axis_labels]
             df.sort_values(by=mesh_cols + columns, inplace=True)
         else:
             df.sort_values(by=[self.domain_type] + columns, inplace=True)