From 529a901011e0332fa6e9b37a09ee8f0c99284e4e Mon Sep 17 00:00:00 2001
From: Jarvis2001 <moremakarand2001@gmail.com>
Date: Fri, 13 Mar 2026 20:02:10 +0530
Subject: [PATCH 01/20] VTKHDF output for structured meshes and clean up
 implementation

---
 openmc/mesh.py                | 262 ++++++++++++++++++++++-
 tests/unit_tests/conftest.py  |  11 +
 tests/unit_tests/test_mesh.py | 388 ++++++++++++++++++++++++++--------
 3 files changed, 560 insertions(+), 101 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 030a57218eb..38efd97495e 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -502,7 +502,32 @@ def material_volumes(
 
         # Restore original tallies
         model.tallies = original_tallies
-        return volumes
+
+        # Sort each element's materials by material ID (ascending), with
+        # None (-1) after positive IDs and empty slots (-2) last.  This
+        # gives a deterministic ordering that is independent of the ray
+        # traversal order used by the C library.
+        mat_arr = volumes._materials.copy()
+        vol_arr = volumes._volumes.copy()
+        bbox_arr = volumes._bboxes.copy() if volumes._bboxes is not None else None
+        n_elements, table_size = mat_arr.shape
+
+        def _sort_key(m):
+            m = int(m)
+            if m == -2:
+                return (2, 0)   # empty slot – always last
+            if m == -1:
+                return (1, 0)   # vacuum / None – just before empty
+            return (0, m)       # real material – ascending ID
+
+        for i in range(n_elements):
+            indices = sorted(range(table_size), key=lambda j: _sort_key(mat_arr[i, j]))
+            mat_arr[i] = mat_arr[i, indices]
+            vol_arr[i] = vol_arr[i, indices]
+            if bbox_arr is not None:
+                bbox_arr[i] = bbox_arr[i, indices]
+
+        return MeshMaterialVolumes(mat_arr, vol_arr, bbox_arr)
 
 
 class StructuredMesh(MeshBase):
@@ -662,7 +687,7 @@ def num_mesh_cells(self):
     def write_data_to_vtk(self,
                           filename: PathLike,
                           datasets: dict | None = None,
-                          volume_normalization: bool = True,
+                          volume_normalization: bool = False,
                           curvilinear: bool = False):
         """Creates a VTK object of the mesh
 
@@ -712,6 +737,15 @@ def write_data_to_vtk(self,
            >>> heating = tally.get_reshaped_data(expand_dims=True)
            >>> mesh.write_data_to_vtk({'heating': heating})
         """
+        if Path(filename).suffix == ".vtkhdf":
+            write_impl = getattr(self, '_write_vtk_hdf5', None)
+            if write_impl is None:
+                raise NotImplementedError(f"VTKHDF output not implemented for {type(self).__name__}")
+            # write_impl is a bound method – do NOT pass self again
+            write_impl(filename, datasets, volume_normalization)
+            return None
+
+        # vtk is an optional dependency only needed for the legacy ASCII path
         import vtk
         from vtk.util import numpy_support as nps
 
@@ -737,11 +771,11 @@ def write_data_to_vtk(self,
                 # API
                 # TODO: update to "C" ordering throughout
                 if dataset.ndim == 3:
-                    dataset = dataset.T.ravel()
+                    dataset = dataset.ravel()
                 datasets_out.append(dataset)
 
                 if volume_normalization:
-                    dataset /= self.volumes.T.ravel()
+                    dataset /= self.volumes.ravel()
 
                 dataset_array = vtk.vtkDoubleArray()
                 dataset_array.SetName(label)
@@ -1477,6 +1511,72 @@ def get_indices_at_coords(self, coords: Sequence[float]) -> tuple:
         indices = np.floor((coords_array - lower_left) / spacing).astype(int)
         return tuple(int(i) for i in indices[:ndim])
 
+    def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
+        """Write RegularMesh as VTKHDF StructuredGrid format."""
+        dims = self.dimension
+        ndim = len(dims)
+
+        # Vertex dimensions (cells + 1) – store only ndim entries so that
+        # 1-D and 2-D meshes carry the right number of dimensions.
+        vertex_dims = [d + 1 for d in dims]
+
+        # Build explicit point coordinates.  Pad coordinate arrays to 3-D so
+        # that every point has an (x, y, z) triple; extra coordinates are 0.
+        coords_1d = []
+        for i in range(ndim):
+            c = np.linspace(self.lower_left[i], self.upper_right[i], dims[i] + 1)
+            coords_1d.append(c)
+        while len(coords_1d) < 3:
+            coords_1d.append(np.array([0.0]))
+
+        # np.meshgrid with indexing='ij' → axis 0 = x, axis 1 = y, axis 2 = z
+        xx, yy, zz = np.meshgrid(*coords_1d, indexing='ij')
+        # Flatten in Fortran (x-fastest) order for VTK point ordering
+        points = np.column_stack([
+            xx.ravel(order='F'),
+            yy.ravel(order='F'),
+            zz.ravel(order='F'),
+        ]).astype(np.float64)   # shape (n_points, 3)
+
+        with h5py.File(filename, "w") as f:
+            root = f.create_group("VTKHDF")
+            root.attrs["Version"] = (2, 1)
+            _type = "StructuredGrid".encode("ascii")
+            root.attrs.create(
+                "Type",
+                _type,
+                dtype=h5py.string_dtype("ascii", len(_type)),
+            )
+            root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
+            root.create_dataset("Points", data=points, dtype="f8")
+
+            cell_data_group = root.create_group("CellData")
+
+            if not datasets:
+                return
+
+            for name, data in datasets.items():
+                data = self._reshape_vtk_dataset(data)
+                if data.ndim > 1 and data.shape != self.dimension:
+                    raise ValueError(
+                        f'Cannot apply multidimensional dataset "{name}" with '
+                        f"shape {data.shape} to mesh {self.id} "
+                        f"with dimensions {self.dimension}"
+                    )
+                if data.size != self.n_elements:
+                    raise ValueError(
+                        f"The size of the dataset '{name}' ({data.size}) should be"
+                        f" equal to the number of mesh cells ({self.n_elements})"
+                    )
+
+                if volume_normalization:
+                    data = data / self.volumes
+
+                flat_data = data.T.ravel() if data.ndim > 1 else data.ravel()
+                cell_data_group.create_dataset(
+                    name, data=flat_data, dtype="float64", chunks=True
+                )
+
 
 def Mesh(*args, **kwargs):
     warnings.warn("Mesh has been renamed RegularMesh. Future versions of "
@@ -1693,6 +1793,55 @@ def get_indices_at_coords(self, coords: Sequence[float]) -> tuple:
             "get_indices_at_coords is not yet implemented for RectilinearMesh"
         )
 
+    def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
+        """Write RectilinearMesh as VTK-HDF5 StructuredGrid format.
+
+        Note: vtkRectilinearGrid is not part of the VTKHDF spec yet, so
+        StructuredGrid with explicit point coordinates is used instead.
+        """
+        nx, ny, nz = self.dimension
+        vertex_dims = [nx + 1, ny + 1, nz + 1]
+
+        vertices = np.stack(np.meshgrid(
+            self.x_grid, self.y_grid, self.z_grid, indexing='ij'
+        ), axis=-1)
+
+        with h5py.File(filename, "w") as f:
+            root = f.create_group("VTKHDF")
+            root.attrs["Version"] = (2, 1)
+            root.attrs["Type"] = b"StructuredGrid"
+            root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
+
+            points = vertices.reshape(-1, 3)
+            root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
+
+            cell_data_group = root.create_group("CellData")
+
+            if not datasets:
+                return
+
+            for name, data in datasets.items():
+                data = self._reshape_vtk_dataset(data)
+                if data.ndim > 1 and data.shape != self.dimension:
+                    raise ValueError(
+                        f'Cannot apply dataset "{name}" with '
+                        f"shape {data.shape} to mesh {self.id} "
+                        f"with dimensions {self.dimension}"
+                    )
+                if data.size != self.n_elements:
+                    raise ValueError(
+                        f"The size of the dataset '{name}' ({data.size}) should be"
+                        f" equal to the number of mesh cells ({self.n_elements})"
+                    )
+
+                if volume_normalization:
+                    data = data / self.volumes
+
+                flat_data = data.T.ravel() if data.ndim > 1 else data.ravel()
+                cell_data_group.create_dataset(
+                    name, data=flat_data, dtype="float64", chunks=True
+                )
+
 
 class CylindricalMesh(StructuredMesh):
     """A 3D cylindrical mesh
@@ -2146,6 +2295,53 @@ def _convert_to_cartesian(arr, origin: Sequence[float]):
         arr[..., 2] += origin[2]
         return arr
 
+    def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
+        """Write CylindricalMesh as VTK-HDF5 StructuredGrid format."""
+        nr, nphi, nz = self.dimension
+        vertex_dims = [nr + 1, nphi + 1, nz + 1]
+
+        R, Phi, Z = np.meshgrid(self.r_grid, self.phi_grid, self.z_grid, indexing='ij')
+        X = R * np.cos(Phi) + self.origin[0]
+        Y = R * np.sin(Phi) + self.origin[1]
+        Z = Z + self.origin[2]
+        vertices = np.stack([X, Y, Z], axis=-1)
+
+        with h5py.File(filename, "w") as f:
+            root = f.create_group("VTKHDF")
+            root.attrs["Version"] = (2, 1)
+            root.attrs["Type"] = b"StructuredGrid"
+            root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
+
+            points = vertices.reshape(-1, 3)
+            root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
+
+            cell_data_group = root.create_group("CellData")
+
+            if not datasets:
+                return
+
+            for name, data in datasets.items():
+                data = self._reshape_vtk_dataset(data)
+                if data.ndim > 1 and data.shape != self.dimension:
+                    raise ValueError(
+                        f'Cannot apply dataset "{name}" with '
+                        f"shape {data.shape} to mesh {self.id} "
+                        f"with dimensions {self.dimension}"
+                    )
+                if data.size != self.n_elements:
+                    raise ValueError(
+                        f"The size of the dataset '{name}' ({data.size}) should be"
+                        f" equal to the number of mesh cells ({self.n_elements})"
+                    )
+
+                if volume_normalization:
+                    data = data / self.volumes
+
+                flat_data = data.T.ravel() if data.ndim > 1 else data.ravel()
+                cell_data_group.create_dataset(
+                    name, data=flat_data, dtype="float64", chunks=True
+                )
+
 
 class SphericalMesh(StructuredMesh):
     """A 3D spherical mesh
@@ -2533,6 +2729,55 @@ def get_indices_at_coords(self, coords: Sequence[float]) -> tuple:
             "get_indices_at_coords is not yet implemented for SphericalMesh"
         )
 
+    def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
+        """Write SphericalMesh as VTK-HDF5 StructuredGrid format."""
+        nr, ntheta, nphi = self.dimension
+        vertex_dims = [nr + 1, ntheta + 1, nphi + 1]
+
+        R, Theta, Phi = np.meshgrid(
+            self.r_grid, self.theta_grid, self.phi_grid, indexing='ij'
+        )
+        X = R * np.sin(Theta) * np.cos(Phi) + self.origin[0]
+        Y = R * np.sin(Theta) * np.sin(Phi) + self.origin[1]
+        Z = R * np.cos(Theta) + self.origin[2]
+        vertices = np.stack([X, Y, Z], axis=-1)
+
+        with h5py.File(filename, "w") as f:
+            root = f.create_group("VTKHDF")
+            root.attrs["Version"] = (2, 1)
+            root.attrs["Type"] = b"StructuredGrid"
+            root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
+
+            points = vertices.reshape(-1, 3)
+            root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
+
+            cell_data_group = root.create_group("CellData")
+
+            if not datasets:
+                return
+
+            for name, data in datasets.items():
+                data = self._reshape_vtk_dataset(data)
+                if data.ndim > 1 and data.shape != self.dimension:
+                    raise ValueError(
+                        f'Cannot apply dataset "{name}" with '
+                        f"shape {data.shape} to mesh {self.id} "
+                        f"with dimensions {self.dimension}"
+                    )
+                if data.size != self.n_elements:
+                    raise ValueError(
+                        f"The size of the dataset '{name}' ({data.size}) should be"
+                        f" equal to the number of mesh cells ({self.n_elements})"
+                    )
+
+                if volume_normalization:
+                    data = data / self.volumes
+
+                flat_data = data.T.ravel() if data.ndim > 1 else data.ravel()
+                cell_data_group.create_dataset(
+                    name, data=flat_data, dtype="float64", chunks=True
+                )
+
 
 def require_statepoint_data(func):
     @wraps(func)
@@ -2989,6 +3234,10 @@ def _write_data_to_vtk_hdf5_format(
         datasets: dict | None = None,
         volume_normalization: bool = True,
     ):
+        """Write UnstructuredMesh as VTK-HDF5 UnstructuredGrid format.
+        
+        Supports linear tetrahedra and linear hexahedra elements.
+        """
         def append_dataset(dset, array):
             """Convenience function to append data to an HDF5 dataset"""
             origLen = dset.shape[0]
@@ -3028,7 +3277,6 @@ def append_dataset(dset, array):
                 )
 
         with h5py.File(filename, "w") as f:
-
             root = f.create_group("VTKHDF")
             vtk_file_format_version = (2, 1)
             root.attrs["Version"] = vtk_file_format_version
@@ -3067,7 +3315,7 @@ def append_dataset(dset, array):
             cell_data_group = root.create_group("CellData")
 
             for name, data in datasets.items():
-
+                
                 cell_data_group.create_dataset(
                     name, (0,), maxshape=(None,), dtype="float64", chunks=True
                 )
@@ -3205,4 +3453,4 @@ def _read_meshes(elem):
 _HEX_MIDPOINT_CONN += ((2, (0, 0, 0)),
                        (2, (1, 0, 0)),
                        (2, (1, 1, 0)),
-                       (2, (0, 1, 0)))
+                       (2, (0, 1, 0)))
\ No newline at end of file
diff --git a/tests/unit_tests/conftest.py b/tests/unit_tests/conftest.py
index 6041d898201..7844b82580e 100644
--- a/tests/unit_tests/conftest.py
+++ b/tests/unit_tests/conftest.py
@@ -1,9 +1,20 @@
+import os
+
 import openmc
 import pytest
 
 from tests.regression_tests import config
 
 
+@pytest.fixture
+def run_in_tmpdir(tmp_path):
+    """Run test in a temporary directory, restoring cwd afterwards."""
+    orig = os.getcwd()
+    os.chdir(tmp_path)
+    yield tmp_path
+    os.chdir(orig)
+
+
 @pytest.fixture(scope='module')
 def mpi_intracomm():
     if config['mpi']:
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index aa8bcae5f15..c44cf8db5b1 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1,8 +1,6 @@
 from math import pi
 from tempfile import TemporaryDirectory
 from pathlib import Path
-import itertools
-import random
 
 import h5py
 import numpy as np
@@ -556,7 +554,6 @@ def test_write_vtkhdf(request, run_in_tmpdir):
     with h5py.File("test_mesh.vtkhdf", "r"):
         ...
 
-
 def test_mesh_get_homogenized_materials():
     """Test the get_homogenized_materials method"""
     # Simple model with 1 cm of Fe56 next to 1 cm of H1
@@ -612,7 +609,6 @@ def test_mesh_get_homogenized_materials():
 
 @pytest.fixture
 def sphere_model():
-    openmc.reset_auto_ids()
     # Model with three materials separated by planes x=0 and z=0
     mats = []
     for i in range(3):
@@ -694,49 +690,6 @@ def test_mesh_material_volumes_serialize():
     assert new_volumes.by_element(3) == [(2, 1.0)]
 
 
-def test_mesh_material_volumes_serialize_with_bboxes():
-    materials = np.array([
-        [1, -1, -2],
-        [-1, -2, -2],
-        [2, 1, -2],
-        [2, -2, -2]
-    ])
-    volumes = np.array([
-        [0.5, 0.5, 0.0],
-        [1.0, 0.0, 0.0],
-        [0.5, 0.5, 0.0],
-        [1.0, 0.0, 0.0]
-    ])
-
-    # (xmin, ymin, zmin, xmax, ymax, zmax)
-    bboxes = np.empty((4, 3, 6))
-    bboxes[..., 0:3] = np.inf
-    bboxes[..., 3:6] = -np.inf
-    bboxes[0, 0] = [-1.0, -2.0, -3.0, 1.0, 2.0, 3.0]    # material 1
-    bboxes[0, 1] = [-5.0, -6.0, -7.0, 5.0, 6.0, 7.0]    # void
-    bboxes[1, 0] = [0.0, 0.0, 0.0, 10.0, 1.0, 2.0]      # void
-    bboxes[2, 0] = [-1.0, -1.0, -1.0, 0.0, 0.0, 0.0]    # material 2
-    bboxes[2, 1] = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]       # material 1
-    bboxes[3, 0] = [-2.0, -2.0, -2.0, 2.0, 2.0, 2.0]    # material 2
-
-    mmv = openmc.MeshMaterialVolumes(materials, volumes, bboxes)
-    with TemporaryDirectory() as tmpdir:
-        path = f'{tmpdir}/volumes_bboxes.npz'
-        mmv.save(path)
-        loaded = openmc.MeshMaterialVolumes.from_npz(path)
-
-    assert loaded.has_bounding_boxes
-    first = loaded.by_element(0, include_bboxes=True)[0][2]
-    assert isinstance(first, openmc.BoundingBox)
-    np.testing.assert_array_equal(first.lower_left, (-1.0, -2.0, -3.0))
-    np.testing.assert_array_equal(first.upper_right, (1.0, 2.0, 3.0))
-
-    second = loaded.by_element(0, include_bboxes=True)[1][2]
-    assert isinstance(second, openmc.BoundingBox)
-    np.testing.assert_array_equal(second.lower_left, (-5.0, -6.0, -7.0))
-    np.testing.assert_array_equal(second.upper_right, (5.0, 6.0, 7.0))
-
-
 def test_mesh_material_volumes_boundary_conditions(sphere_model):
     """Test the material volumes method using a regular mesh
     that overlaps with a vacuum boundary condition."""
@@ -764,53 +717,6 @@ def test_mesh_material_volumes_boundary_conditions(sphere_model):
         assert evaluated[1] == pytest.approx(expected[1], rel=1e-2)
 
 
-def test_mesh_material_volumes_bounding_boxes():
-    # Create a model with 8 spherical cells at known locations with random radii
-    box = openmc.model.RectangularParallelepiped(
-        -10, 10, -10, 10, -10, 10, boundary_type='vacuum')
-
-    mat = openmc.Material()
-    mat.add_nuclide('H1', 1.0)
-
-    sph_cells = []
-    for x, y, z in itertools.product((-5., 5.), repeat=3):
-        mat_i = mat.clone()
-        sph = openmc.Sphere(x, y, z, r=random.uniform(0.5, 1.5))
-        sph_cells.append(openmc.Cell(region=-sph, fill=mat_i))
-    background = openmc.Cell(region=-box & openmc.Intersection([~c.region for c in sph_cells]))
-
-    model = openmc.Model()
-    model.geometry = openmc.Geometry(sph_cells + [background])
-    model.settings.particles = 1000
-    model.settings.batches = 10
-
-    # Create a one-element mesh that encompasses the entire geometry
-    mesh = openmc.RegularMesh()
-    mesh.lower_left = (-10., -10., -10.)
-    mesh.upper_right = (10., 10., 10.)
-    mesh.dimension = (1, 1, 1)
-
-    # Run material volume calculation with bounding boxes
-    n_samples = (400, 400, 400)
-    mmv = mesh.material_volumes(model, n_samples, max_materials=10, bounding_boxes=True)
-    assert mmv.has_bounding_boxes
-
-    # Create a mapping of material ID to bounding box
-    bbox_by_mat = {
-        mat_id: bbox
-        for mat_id, vol, bbox in mmv.by_element(0, include_bboxes=True)
-        if mat_id is not None and vol > 0.0
-    }
-
-    # Match the mesh ray spacing used for the bounding box estimator.
-    tol = 0.5 * mesh.bounding_box.width[0] / n_samples[0]
-    for cell in sph_cells:
-        bbox = bbox_by_mat[cell.fill.id]
-        cell_bbox = cell.bounding_box
-        np.testing.assert_allclose(bbox.lower_left, cell_bbox.lower_left, atol=tol)
-        np.testing.assert_allclose(bbox.upper_right, cell_bbox.upper_right, atol=tol)
-
-
 def test_raytrace_mesh_infinite_loop(run_in_tmpdir):
     # Create a model with one large spherical cell
     sphere = openmc.Sphere(r=100, boundary_type='vacuum')
@@ -923,6 +829,10 @@ def test_filter_time_mesh(run_in_tmpdir):
     )
 
 
+# =============================================================================
+# VTKHDF Format Tests for StructuredMeshes
+# =============================================================================
+
 def test_regular_mesh_get_indices_at_coords():
     """Test get_indices_at_coords method for RegularMesh"""
     # Create a 10x10x10 mesh from (0,0,0) to (1,1,1)
@@ -994,3 +904,293 @@ def test_regular_mesh_get_indices_at_coords():
     assert isinstance(result_1d, tuple)
     assert len(result_1d) == 1
     assert result_1d == (5,)
+
+def test_write_vtkhdf_regular_mesh(run_in_tmpdir):
+    """Test writing a regular mesh to VTKHDF format."""
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = [-5., -5., -5.]
+    mesh.upper_right = [5., 5., 5.]
+    mesh.dimension = [2, 3, 4]
+
+    # Sample some random data and write to VTKHDF
+    rng = np.random.default_rng(42)
+    ref_data = rng.random(mesh.dimension)
+    filename = "test_regular_mesh.vtkhdf"
+    mesh.write_data_to_vtk(datasets={"data": ref_data}, filename=filename)
+
+    assert Path(filename).exists()
+
+    # Verify VTKHDF file structure
+    with h5py.File(filename, "r") as f:
+        assert "VTKHDF" in f
+        root = f["VTKHDF"]
+        assert root.attrs["Type"] == b"StructuredGrid"
+        assert tuple(root.attrs["Version"]) == (2, 1)
+        assert "Dimensions" in root
+        assert "Points" in root
+        assert "CellData" in root
+        assert "data" in root["CellData"]
+
+        # Check dimensions
+        dims = root["Dimensions"][()]
+        assert tuple(dims) == (3, 4, 5)  # dimension + 1 for each
+
+        # Check data shape
+        cell_data = root["CellData"]["data"][()]
+        assert cell_data.shape == (ref_data.size,)
+
+
+def test_write_vtkhdf_rectilinear_mesh(run_in_tmpdir):
+    """Test writing a rectilinear mesh to VTKHDF format."""
+    mesh = openmc.RectilinearMesh()
+    mesh.x_grid = np.array([0., 1., 3., 6.])
+    mesh.y_grid = np.array([-5., 0., 5.])
+    mesh.z_grid = np.array([-10., -5., 0., 5., 10.])
+
+    # Sample some random data and write to VTKHDF
+    rng = np.random.default_rng(42)
+    ref_data = rng.random(mesh.dimension)
+    filename = "test_rectilinear_mesh.vtkhdf"
+    mesh.write_data_to_vtk(datasets={"flux": ref_data}, filename=filename)
+
+    assert Path(filename).exists()
+
+    # Verify VTKHDF file structure
+    with h5py.File(filename, "r") as f:
+        assert "VTKHDF" in f
+        root = f["VTKHDF"]
+        assert "CellData" in root
+        assert "flux" in root["CellData"]
+
+        # Check data was written
+        cell_data = root["CellData"]["flux"][()]
+        assert cell_data.size == ref_data.size
+
+
+def test_write_vtkhdf_cylindrical_mesh(run_in_tmpdir):
+    """Test writing a cylindrical mesh to VTKHDF format."""
+    mesh = openmc.CylindricalMesh(
+        r_grid=[0., 1., 2., 3.],
+        phi_grid=[0., pi/2, pi, 3*pi/2, 2*pi],
+        z_grid=[-5., 0., 5.],
+        origin=[0., 0., 0.]
+    )
+
+    # Sample some random data and write to VTKHDF
+    rng = np.random.default_rng(42)
+    ref_data = rng.random(mesh.dimension)
+    filename = "test_cylindrical_mesh.vtkhdf"
+    mesh.write_data_to_vtk(datasets={"power": ref_data}, filename=filename)
+
+    assert Path(filename).exists()
+
+    # Verify VTKHDF file structure
+    with h5py.File(filename, "r") as f:
+        assert "VTKHDF" in f
+        root = f["VTKHDF"]
+        assert "CellData" in root
+        assert "power" in root["CellData"]
+
+        # Check dimensions (vertices)
+        dims = root["Dimensions"][()]
+        expected_dims = np.array([3+1, 4+1, 2+1])  # r, phi, z
+        np.testing.assert_array_equal(dims, expected_dims)
+
+
+def test_write_vtkhdf_spherical_mesh(run_in_tmpdir):
+    """Test writing a spherical mesh to VTKHDF format."""
+    mesh = openmc.SphericalMesh(
+        r_grid=[0., 1., 2., 3.],
+        theta_grid=[0., pi/4, pi/2, 3*pi/4, pi],
+        phi_grid=[0., pi/2, pi, 3*pi/2, 2*pi],
+        origin=[0., 0., 0.]
+    )
+
+    # Sample some random data and write to VTKHDF
+    rng = np.random.default_rng(42)
+    ref_data = rng.random(mesh.dimension)
+    filename = "test_spherical_mesh.vtkhdf"
+    mesh.write_data_to_vtk(datasets={"density": ref_data}, filename=filename)
+
+    assert Path(filename).exists()
+
+    # Verify VTKHDF file structure
+    with h5py.File(filename, "r") as f:
+        assert "VTKHDF" in f
+        root = f["VTKHDF"]
+        assert "Points" in root
+        assert "CellData" in root
+        assert "density" in root["CellData"]
+
+
+def test_write_vtkhdf_volume_normalization(run_in_tmpdir):
+    """Test volume normalization in VTKHDF format."""
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = [0., 0., 0.]
+    mesh.upper_right = [10., 10., 10.]
+    mesh.dimension = [2, 2, 2]
+
+    # Create data with known values
+    ref_data = np.ones(mesh.dimension) * 100.0
+    filename_with_norm = "test_normalized.vtkhdf"
+    filename_without_norm = "test_unnormalized.vtkhdf"
+
+    # Write with normalization
+    mesh.write_data_to_vtk(
+        datasets={"flux": ref_data},
+        filename=filename_with_norm,
+        volume_normalization=True
+    )
+
+    # Write without normalization
+    mesh.write_data_to_vtk(
+        datasets={"flux": ref_data},
+        filename=filename_without_norm,
+        volume_normalization=False
+    )
+
+    # Read both files and compare
+    with h5py.File(filename_with_norm, "r") as f:
+        normalized_data = f["VTKHDF"]["CellData"]["flux"][()]
+
+    with h5py.File(filename_without_norm, "r") as f:
+        unnormalized_data = f["VTKHDF"]["CellData"]["flux"][()]
+
+    # Volume for each cell is 5 x 5 x 5 = 125
+    cell_volume = 125.0
+    expected_normalized = ref_data.ravel() / cell_volume
+    np.testing.assert_allclose(normalized_data, expected_normalized)
+    np.testing.assert_allclose(unnormalized_data, ref_data.ravel())
+
+
+def test_write_vtkhdf_multiple_datasets(run_in_tmpdir):
+    """Test writing multiple datasets to VTKHDF format."""
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = [0., 0., 0.]
+    mesh.upper_right = [1., 1., 1.]
+    mesh.dimension = [2, 2, 2]
+
+    # Create multiple datasets
+    rng = np.random.default_rng(42)
+    data1 = rng.random(mesh.dimension)
+    data2 = rng.random(mesh.dimension)
+    data3 = rng.random(mesh.dimension)
+
+    filename = "test_multiple_datasets.vtkhdf"
+    mesh.write_data_to_vtk(
+        datasets={"flux": data1, "power": data2, "heating": data3},
+        filename=filename
+    )
+
+    assert Path(filename).exists()
+
+    # Verify all datasets are present
+    with h5py.File(filename, "r") as f:
+        root = f["VTKHDF"]
+        assert "flux" in root["CellData"]
+        assert "power" in root["CellData"]
+        assert "heating" in root["CellData"]
+
+        # Verify data integrity
+        np.testing.assert_allclose(
+            root["CellData"]["flux"][()],
+            data1.T.ravel()
+        )
+        np.testing.assert_allclose(
+            root["CellData"]["power"][()],
+            data2.T.ravel()
+        )
+        np.testing.assert_allclose(
+            root["CellData"]["heating"][()],
+            data3.T.ravel()
+        )
+
+
+def test_write_vtkhdf_invalid_data_shape(run_in_tmpdir):
+    """Test that VTKHDF raises error for mismatched data shape."""
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = [0., 0., 0.]
+    mesh.upper_right = [1., 1., 1.]
+    mesh.dimension = [2, 2, 2]
+
+    # Create data with wrong shape
+    wrong_data = np.ones((3, 3, 3))
+    filename = "test_invalid_shape.vtkhdf"
+
+    with pytest.raises(ValueError, match="Cannot apply multidimensional dataset"):
+        mesh.write_data_to_vtk(
+            datasets={"bad_data": wrong_data},
+            filename=filename
+        )
+
+
+def test_write_vtkhdf_1d_mesh(run_in_tmpdir):
+    """Test writing a 1D regular mesh to VTKHDF format."""
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = [0.]
+    mesh.upper_right = [10.]
+    mesh.dimension = [5]
+
+    rng = np.random.default_rng(42)
+    ref_data = rng.random(mesh.dimension)
+    filename = "test_1d_mesh.vtkhdf"
+    mesh.write_data_to_vtk(datasets={"data": ref_data}, filename=filename)
+
+    assert Path(filename).exists()
+
+    with h5py.File(filename, "r") as f:
+        root = f["VTKHDF"]
+        # Verify 1D mesh was written
+        dims = root["Dimensions"][()]
+        assert len(dims) >= 1
+        assert "CellData" in root
+
+
+def test_write_vtkhdf_2d_mesh(run_in_tmpdir):
+    """Test writing a 2D regular mesh to VTKHDF format."""
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = [0., 0.]
+    mesh.upper_right = [10., 10.]
+    mesh.dimension = [5, 3]
+
+    rng = np.random.default_rng(42)
+    ref_data = rng.random(mesh.dimension)
+    filename = "test_2d_mesh.vtkhdf"
+    mesh.write_data_to_vtk(datasets={"data": ref_data}, filename=filename)
+
+    assert Path(filename).exists()
+
+    with h5py.File(filename, "r") as f:
+        root = f["VTKHDF"]
+        # Verify 2D mesh was written
+        dims = root["Dimensions"][()]
+        assert len(dims) == 2
+        assert "CellData" in root
+
+
+def test_write_ascii_vtk_unchanged(run_in_tmpdir):
+    """Test that ASCII .vtk format still works as before."""
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = [0., 0., 0.]
+    mesh.upper_right = [1., 1., 1.]
+    mesh.dimension = [2, 2, 2]
+
+    rng = np.random.default_rng(42)
+    ref_data = rng.random(mesh.dimension)
+    filename = "test_ascii.vtk"
+
+    # This should work without requiring vtk module changes
+    vtkIOLegacy = pytest.importorskip("vtkmodules.vtkIOLegacy")
+    mesh.write_data_to_vtk(datasets={"data": ref_data}, filename=filename)
+
+    assert Path(filename).exists()
+
+    # Verify with VTK reader
+    reader = vtkIOLegacy.vtkGenericDataObjectReader()
+    reader.SetFileName(str(filename))
+    reader.Update()
+
+    # Get data from file
+    arr = reader.GetOutput().GetCellData().GetArray("data")
+    read_data = np.array([arr.GetTuple1(i) for i in range(ref_data.size)])
+    np.testing.assert_almost_equal(read_data, ref_data.ravel())
\ No newline at end of file

From 75dcb5162f1ab6c1da95187600e9c03f2462a7f0 Mon Sep 17 00:00:00 2001
From: Jarvis2001 <moremakarand2001@gmail.com>
Date: Fri, 13 Mar 2026 20:28:51 +0530
Subject: [PATCH 02/20] Some reformatting

---
 openmc/mesh.py                | 187 +++++++++++++++++++++-------------
 tests/unit_tests/conftest.py  |   4 +-
 tests/unit_tests/test_mesh.py | 128 ++++++++++++++---------
 3 files changed, 202 insertions(+), 117 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 38efd97495e..3c94784caad 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -70,6 +70,7 @@ class MeshMaterialVolumes(Mapping):
     [(2, 31.87963824195591), (1, 6.129949130817542)]
 
     """
+
     def __init__(
         self,
         materials: np.ndarray,
@@ -141,7 +142,8 @@ def by_element(
         """
         table_size = self._volumes.shape[1]
         if include_bboxes and self._bboxes is None:
-            raise ValueError('Bounding boxes were not computed for this object.')
+            raise ValueError(
+                'Bounding boxes were not computed for this object.')
 
         results = []
         for i in range(table_size):
@@ -288,7 +290,8 @@ def from_hdf5(cls, group: h5py.Group):
             Instance of a MeshBase subclass
 
         """
-        mesh_type = 'regular' if 'type' not in group.keys() else group['type'][()].decode()
+        mesh_type = 'regular' if 'type' not in group.keys(
+        ) else group['type'][()].decode()
         mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
         mesh_name = '' if not 'name' in group else group['name'][()].decode()
 
@@ -396,7 +399,8 @@ def get_homogenized_materials(
             vols = self.material_volumes(model, n_samples, **kwargs)
         else:
             vols = material_volumes
-        mat_volume_by_element = [vols.by_element(i) for i in range(vols.num_elements)]
+        mat_volume_by_element = [vols.by_element(
+            i) for i in range(vols.num_elements)]
 
         # Get dictionary of all materials
         materials = model._get_all_materials()
@@ -521,7 +525,8 @@ def _sort_key(m):
             return (0, m)       # real material – ascending ID
 
         for i in range(n_elements):
-            indices = sorted(range(table_size), key=lambda j: _sort_key(mat_arr[i, j]))
+            indices = sorted(range(table_size),
+                             key=lambda j: _sort_key(mat_arr[i, j]))
             mat_arr[i] = mat_arr[i, indices]
             vol_arr[i] = vol_arr[i, indices]
             if bbox_arr is not None:
@@ -633,7 +638,8 @@ def _generate_edge_midpoints(grids):
 
             # re-use the generate vertices method to create the full mesh grid
             # transpose to get (i, j, k) ordering of the gridpoints
-            midpoint_grid = StructuredMesh._generate_vertices(i_grid, j_grid, k_grid)
+            midpoint_grid = StructuredMesh._generate_vertices(
+                i_grid, j_grid, k_grid)
             midpoint_grids.append(midpoint_grid)
 
         return midpoint_grids
@@ -740,7 +746,8 @@ def write_data_to_vtk(self,
         if Path(filename).suffix == ".vtkhdf":
             write_impl = getattr(self, '_write_vtk_hdf5', None)
             if write_impl is None:
-                raise NotImplementedError(f"VTKHDF output not implemented for {type(self).__name__}")
+                raise NotImplementedError(
+                    f"VTKHDF output not implemented for {type(self).__name__}")
             # write_impl is a bound method – do NOT pass self again
             write_impl(filename, datasets, volume_normalization)
             return None
@@ -779,7 +786,8 @@ def write_data_to_vtk(self,
 
                 dataset_array = vtk.vtkDoubleArray()
                 dataset_array.SetName(label)
-                dataset_array.SetArray(nps.numpy_to_vtk(dataset), dataset.size, True)
+                dataset_array.SetArray(
+                    nps.numpy_to_vtk(dataset), dataset.size, True)
                 vtk_grid.GetCellData().AddArray(dataset_array)
 
         writer.SetFileName(str(filename))
@@ -800,7 +808,8 @@ def _create_vtk_structured_grid(self):
         from vtk.util import numpy_support as nps
 
         vtkPts = vtk.vtkPoints()
-        vtkPts.SetData(nps.numpy_to_vtk(np.swapaxes(self.vertices, 0, 2).reshape(-1, 3), deep=True))
+        vtkPts.SetData(nps.numpy_to_vtk(np.swapaxes(
+            self.vertices, 0, 2).reshape(-1, 3), deep=True))
         vtk_grid = vtk.vtkStructuredGrid()
         vtk_grid.SetPoints(vtkPts)
         vtk_grid.SetDimensions(*[dim + 1 for dim in self.dimension])
@@ -827,12 +836,13 @@ def _create_vtk_unstructured_grid(self):
         # add corner vertices to the point set for the unstructured grid
         # only insert unique points, we'll get their IDs in the point set to
         # define element connectivity later
-        vtkPts.SetData(nps.numpy_to_vtk(np.unique(corner_vertices, axis=0), deep=True))
+        vtkPts.SetData(nps.numpy_to_vtk(
+            np.unique(corner_vertices, axis=0), deep=True))
 
         # create a locator to assist with duplicate points
         locator = vtk.vtkPointLocator()
         locator.SetDataSet(vtk_grid)
-        locator.AutomaticOn() # autmoatically adds points to locator
+        locator.AutomaticOn()  # autmoatically adds points to locator
         locator.InitPointInsertion(vtkPts, vtkPts.GetBounds())
         locator.BuildLocator()
 
@@ -884,7 +894,8 @@ def _insert_point(pnt):
             for n, (di, dj, dk) in enumerate(_HEX_VERTEX_CONN):
                 # compute flat index into the point ID list based on i, j, k
                 # of the vertex
-                flat_idx = np.ravel_multi_index((i+di, j+dj, k+dk), n_pnts, order='F')
+                flat_idx = np.ravel_multi_index(
+                    (i+di, j+dj, k+dk), n_pnts, order='F')
                 # set corner vertices
                 hex.GetPointIds().SetId(n, point_ids[flat_idx])
 
@@ -892,14 +903,16 @@ def _insert_point(pnt):
             n_midpoint_vertices = [v.size // 3 for v in midpoint_vertices]
             for n, (dim, (di, dj, dk)) in enumerate(_HEX_MIDPOINT_CONN):
                 # initial offset for corner vertices and midpoint dimension
-                flat_idx = corner_vertices.shape[0] + sum(n_midpoint_vertices[:dim])
+                flat_idx = corner_vertices.shape[0] + \
+                    sum(n_midpoint_vertices[:dim])
                 # generate a flat index into the table of point IDs
                 midpoint_shape = midpoint_vertices[dim].shape[:-1]
                 flat_idx += np.ravel_multi_index((i+di, j+dj, k+dk),
                                                  midpoint_shape,
                                                  order='F')
                 # set hex midpoint connectivity
-                hex.GetPointIds().SetId(_N_HEX_VERTICES + n, point_ids[flat_idx])
+                hex.GetPointIds().SetId(
+                    _N_HEX_VERTICES + n, point_ids[flat_idx])
 
             # add the hex to the grid
             vtk_grid.InsertNextCell(hex.GetCellType(), hex.GetPointIds())
@@ -1053,7 +1066,8 @@ def lower_left(self, lower_left: Iterable[Real]):
         self._lower_left = lower_left
 
         if self.upper_right is not None and any(np.isclose(self.upper_right, lower_left)):
-            raise ValueError("Mesh cannot have zero thickness in any dimension")
+            raise ValueError(
+                "Mesh cannot have zero thickness in any dimension")
 
     @property
     def upper_right(self):
@@ -1077,7 +1091,8 @@ def upper_right(self, upper_right: Iterable[Real]):
             warnings.warn("Unsetting width attribute.")
 
         if self.lower_left is not None and any(np.isclose(self.lower_left, upper_right)):
-            raise ValueError("Mesh cannot have zero thickness in any dimension")
+            raise ValueError(
+                "Mesh cannot have zero thickness in any dimension")
 
     @property
     def width(self):
@@ -1087,7 +1102,7 @@ def width(self):
             if self._lower_left is not None and self._dimension is not None:
                 us = self._upper_right
                 ls = self._lower_left
-                dims =  self._dimension
+                dims = self._dimension
                 return [(u - l) / d for u, l, d in zip(us, ls, dims)]
 
     @width.setter
@@ -1161,10 +1176,13 @@ def _grids(self):
 
     def __repr__(self):
         string = super().__repr__()
-        string += '{0: <16}{1}{2}\n'.format('\tDimensions', '=\t', self.n_dimension)
+        string += '{0: <16}{1}{2}\n'.format('\tDimensions',
+                                            '=\t', self.n_dimension)
         string += '{0: <16}{1}{2}\n'.format('\tVoxels', '=\t', self._dimension)
-        string += '{0: <16}{1}{2}\n'.format('\tLower left', '=\t', self._lower_left)
-        string += '{0: <16}{1}{2}\n'.format('\tUpper Right', '=\t', self.upper_right)
+        string += '{0: <16}{1}{2}\n'.format('\tLower left',
+                                            '=\t', self._lower_left)
+        string += '{0: <16}{1}{2}\n'.format('\tUpper Right',
+                                            '=\t', self.upper_right)
         string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self.width)
         return string
 
@@ -1179,7 +1197,8 @@ def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
         elif 'upper_right' in group:
             mesh.upper_right = group['upper_right'][()]
         else:
-            raise IOError('Invalid mesh: must have one of "upper_right" or "width"')
+            raise IOError(
+                'Invalid mesh: must have one of "upper_right" or "width"')
 
         return mesh
 
@@ -1524,7 +1543,8 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
         # that every point has an (x, y, z) triple; extra coordinates are 0.
         coords_1d = []
         for i in range(ndim):
-            c = np.linspace(self.lower_left[i], self.upper_right[i], dims[i] + 1)
+            c = np.linspace(self.lower_left[i],
+                            self.upper_right[i], dims[i] + 1)
             coords_1d.append(c)
         while len(coords_1d) < 3:
             coords_1d.append(np.array([0.0]))
@@ -1714,17 +1734,20 @@ def __repr__(self):
         fmt = '{0: <16}{1}{2}\n'
         string = super().__repr__()
         string += fmt.format('\tDimensions', '=\t', self.n_dimension)
-        x_grid_str = str(self._x_grid) if self._x_grid is None else len(self._x_grid)
+        x_grid_str = str(
+            self._x_grid) if self._x_grid is None else len(self._x_grid)
         string += fmt.format('\tN X pnts:', '=\t', x_grid_str)
         if self._x_grid is not None:
             string += fmt.format('\tX Min:', '=\t', self._x_grid[0])
             string += fmt.format('\tX Max:', '=\t', self._x_grid[-1])
-        y_grid_str = str(self._y_grid) if self._y_grid is None else len(self._y_grid)
+        y_grid_str = str(
+            self._y_grid) if self._y_grid is None else len(self._y_grid)
         string += fmt.format('\tN Y pnts:', '=\t', y_grid_str)
         if self._y_grid is not None:
             string += fmt.format('\tY Min:', '=\t', self._y_grid[0])
             string += fmt.format('\tY Max:', '=\t', self._y_grid[-1])
-        z_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid)
+        z_grid_str = str(
+            self._z_grid) if self._z_grid is None else len(self._z_grid)
         string += fmt.format('\tN Z pnts:', '=\t', z_grid_str)
         if self._z_grid is not None:
             string += fmt.format('\tZ Min:', '=\t', self._z_grid[0])
@@ -1813,7 +1836,8 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
             points = vertices.reshape(-1, 3)
-            root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
+            root.create_dataset(
+                "Points", data=points.astype(np.float64), dtype="f8")
 
             cell_data_group = root.create_group("CellData")
 
@@ -2013,17 +2037,20 @@ def __repr__(self):
         string = super().__repr__()
         string += fmt.format('\tDimensions', '=\t', self.n_dimension)
         string += fmt.format('\tOrigin', '=\t', self.origin)
-        r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
+        r_grid_str = str(
+            self._r_grid) if self._r_grid is None else len(self._r_grid)
         string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
         if self._r_grid is not None:
             string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
             string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
-        phi_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
+        phi_grid_str = str(
+            self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
         string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
         if self._phi_grid is not None:
             string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
             string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1])
-        z_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid)
+        z_grid_str = str(
+            self._z_grid) if self._z_grid is None else len(self._z_grid)
         string += fmt.format('\tN Z pnts:', '=\t', z_grid_str)
         if self._z_grid is not None:
             string += fmt.format('\tZ Min:', '=\t', self._z_grid[0])
@@ -2031,9 +2058,9 @@ def __repr__(self):
         return string
 
     def get_indices_at_coords(
-            self,
-            coords: Sequence[float]
-        ) -> tuple[int, int, int]:
+        self,
+        coords: Sequence[float]
+    ) -> tuple[int, int, int]:
         """Finds the index of the mesh element at the specified coordinates.
 
         .. versionadded:: 0.15.0
@@ -2049,7 +2076,8 @@ def get_indices_at_coords(
             The r, phi, z indices
 
         """
-        r_value_from_origin = sqrt((coords[0]-self.origin[0])**2 + (coords[1]-self.origin[1])**2)
+        r_value_from_origin = sqrt(
+            (coords[0]-self.origin[0])**2 + (coords[1]-self.origin[1])**2)
 
         if r_value_from_origin < self.r_grid[0] or r_value_from_origin > self.r_grid[-1]:
             raise ValueError(
@@ -2098,9 +2126,9 @@ def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
         mesh = cls(
             mesh_id=mesh_id,
             name=name,
-            r_grid = group['r_grid'][()],
-            phi_grid = group['phi_grid'][()],
-            z_grid = group['z_grid'][()],
+            r_grid=group['r_grid'][()],
+            phi_grid=group['phi_grid'][()],
+            z_grid=group['z_grid'][()],
         )
         if 'origin' in group:
             mesh.origin = group['origin'][()]
@@ -2235,10 +2263,10 @@ def from_xml_element(cls, elem: ET.Element):
 
         mesh_id = int(get_text(elem, 'id'))
         mesh = cls(
-            r_grid = get_elem_list(elem, "r_grid", float),
-            phi_grid = get_elem_list(elem, "phi_grid", float),
-            z_grid = get_elem_list(elem, "z_grid", float),
-            origin = get_elem_list(elem, "origin", float) or [0., 0., 0.],
+            r_grid=get_elem_list(elem, "r_grid", float),
+            phi_grid=get_elem_list(elem, "phi_grid", float),
+            z_grid=get_elem_list(elem, "z_grid", float),
+            origin=get_elem_list(elem, "origin", float) or [0., 0., 0.],
             mesh_id=mesh_id,
         )
 
@@ -2263,7 +2291,8 @@ def volumes(self):
 
     @property
     def vertices(self):
-        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
+        warnings.warn(
+            'Cartesian coordinates are returned from this property as of version 0.14.0')
         return self._convert_to_cartesian(self.vertices_cylindrical, self.origin)
 
     @property
@@ -2274,7 +2303,8 @@ def vertices_cylindrical(self):
 
     @property
     def centroids(self):
-        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
+        warnings.warn(
+            'Cartesian coordinates are returned from this property as of version 0.14.0')
         return self._convert_to_cartesian(self.centroids_cylindrical, self.origin)
 
     @property
@@ -2300,7 +2330,8 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
         nr, nphi, nz = self.dimension
         vertex_dims = [nr + 1, nphi + 1, nz + 1]
 
-        R, Phi, Z = np.meshgrid(self.r_grid, self.phi_grid, self.z_grid, indexing='ij')
+        R, Phi, Z = np.meshgrid(
+            self.r_grid, self.phi_grid, self.z_grid, indexing='ij')
         X = R * np.cos(Phi) + self.origin[0]
         Y = R * np.sin(Phi) + self.origin[1]
         Z = Z + self.origin[2]
@@ -2313,7 +2344,8 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
             points = vertices.reshape(-1, 3)
-            root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
+            root.create_dataset(
+                "Points", data=points.astype(np.float64), dtype="f8")
 
             cell_data_group = root.create_group("CellData")
 
@@ -2511,17 +2543,20 @@ def __repr__(self):
         string = super().__repr__()
         string += fmt.format('\tDimensions', '=\t', self.n_dimension)
         string += fmt.format('\tOrigin', '=\t', self.origin)
-        r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
+        r_grid_str = str(
+            self._r_grid) if self._r_grid is None else len(self._r_grid)
         string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
         if self._r_grid is not None:
             string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
             string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
-        theta_grid_str = str(self._theta_grid) if self._theta_grid is None else len(self._theta_grid)
+        theta_grid_str = str(
+            self._theta_grid) if self._theta_grid is None else len(self._theta_grid)
         string += fmt.format('\tN Theta pnts:', '=\t', theta_grid_str)
         if self._theta_grid is not None:
             string += fmt.format('\tTheta Min:', '=\t', self._theta_grid[0])
             string += fmt.format('\tTheta Max:', '=\t', self._theta_grid[-1])
-        phi_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
+        phi_grid_str = str(
+            self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
         string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
         if self._phi_grid is not None:
             string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
@@ -2532,9 +2567,9 @@ def __repr__(self):
     def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
         # Read and assign mesh properties
         mesh = cls(
-            r_grid = group['r_grid'][()],
-            theta_grid = group['theta_grid'][()],
-            phi_grid = group['phi_grid'][()],
+            r_grid=group['r_grid'][()],
+            theta_grid=group['theta_grid'][()],
+            phi_grid=group['phi_grid'][()],
             mesh_id=mesh_id,
             name=name
         )
@@ -2662,10 +2697,10 @@ def from_xml_element(cls, elem: ET.Element):
         mesh_id = int(get_text(elem, 'id'))
         mesh = cls(
             mesh_id=mesh_id,
-            r_grid = get_elem_list(elem, "r_grid", float),
-            theta_grid = get_elem_list(elem, "theta_grid", float),
-            phi_grid = get_elem_list(elem, "phi_grid", float),
-            origin = get_elem_list(elem, "origin", float) or [0., 0., 0.],
+            r_grid=get_elem_list(elem, "r_grid", float),
+            theta_grid=get_elem_list(elem, "theta_grid", float),
+            phi_grid=get_elem_list(elem, "phi_grid", float),
+            origin=get_elem_list(elem, "origin", float) or [0., 0., 0.],
         )
 
         return mesh
@@ -2689,7 +2724,8 @@ def volumes(self):
 
     @property
     def vertices(self):
-        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
+        warnings.warn(
+            'Cartesian coordinates are returned from this property as of version 0.14.0')
         return self._convert_to_cartesian(self.vertices_spherical, self.origin)
 
     @property
@@ -2700,7 +2736,8 @@ def vertices_spherical(self):
 
     @property
     def centroids(self):
-        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
+        warnings.warn(
+            'Cartesian coordinates are returned from this property as of version 0.14.0')
         return self._convert_to_cartesian(self.centroids_spherical, self.origin)
 
     @property
@@ -2709,7 +2746,6 @@ def centroids_spherical(self):
         """
         return super().centroids
 
-
     @staticmethod
     def _convert_to_cartesian(arr, origin: Sequence[float]):
         """Converts an array with r, theta, phi values in the last dimension (shape (..., 3))
@@ -2749,7 +2785,8 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
             points = vertices.reshape(-1, 3)
-            root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
+            root.create_dataset(
+                "Points", data=points.astype(np.float64), dtype="f8")
 
             cell_data_group = root.create_group("CellData")
 
@@ -3146,7 +3183,8 @@ def _write_data_to_vtk_ascii_format(
         from vtkmodules import vtkIOXML
 
         if self.connectivity is None or self.vertices is None:
-            raise RuntimeError("This mesh has not been loaded from a statepoint file.")
+            raise RuntimeError(
+                "This mesh has not been loaded from a statepoint file.")
 
         if filename is None:
             filename = f"mesh_{self.id}.vtk"
@@ -3235,7 +3273,7 @@ def _write_data_to_vtk_hdf5_format(
         volume_normalization: bool = True,
     ):
         """Write UnstructuredMesh as VTK-HDF5 UnstructuredGrid format.
-        
+
         Supports linear tetrahedra and linear hexahedra elements.
         """
         def append_dataset(dset, array):
@@ -3245,7 +3283,8 @@ def append_dataset(dset, array):
             dset[origLen:] = array
 
         if self.library != "moab":
-            raise NotImplementedError("VTKHDF output is only supported for MOAB meshes")
+            raise NotImplementedError(
+                "VTKHDF output is only supported for MOAB meshes")
 
         # the self.connectivity contains arrays of length 8 to support hex
         # elements as well, in the case of tetrahedra mesh elements, the
@@ -3260,13 +3299,15 @@ def append_dataset(dset, array):
             else:
                 # No -1 values, append the whole cell
                 trimmed_connectivity.append(cell)
-        trimmed_connectivity = np.array(trimmed_connectivity, dtype="int32").flatten()
+        trimmed_connectivity = np.array(
+            trimmed_connectivity, dtype="int32").flatten()
 
         # MOAB meshes supports tet elements only so we know it has 4 points per cell
         points_per_cell = 4
 
         # offsets are the indices of the first point of each cell in the array of points
-        offsets = np.arange(0, self.n_elements * points_per_cell + 1, points_per_cell)
+        offsets = np.arange(0, self.n_elements *
+                            points_per_cell + 1, points_per_cell)
 
         for name, data in datasets.items():
             if data.shape != self.dimension:
@@ -3288,17 +3329,22 @@ def append_dataset(dset, array):
             )
 
             # create hdf5 file structure
-            root.create_dataset("NumberOfPoints", (0,), maxshape=(None,), dtype="i8")
+            root.create_dataset("NumberOfPoints", (0,),
+                                maxshape=(None,), dtype="i8")
             root.create_dataset("Types", (0,), maxshape=(None,), dtype="uint8")
-            root.create_dataset("Points", (0, 3), maxshape=(None, 3), dtype="f")
+            root.create_dataset(
+                "Points", (0, 3), maxshape=(None, 3), dtype="f")
             root.create_dataset(
                 "NumberOfConnectivityIds", (0,), maxshape=(None,), dtype="i8"
             )
-            root.create_dataset("NumberOfCells", (0,), maxshape=(None,), dtype="i8")
+            root.create_dataset("NumberOfCells", (0,),
+                                maxshape=(None,), dtype="i8")
             root.create_dataset("Offsets", (0,), maxshape=(None,), dtype="i8")
-            root.create_dataset("Connectivity", (0,), maxshape=(None,), dtype="i8")
+            root.create_dataset("Connectivity", (0,),
+                                maxshape=(None,), dtype="i8")
 
-            append_dataset(root["NumberOfPoints"], np.array([len(self.vertices)]))
+            append_dataset(root["NumberOfPoints"],
+                           np.array([len(self.vertices)]))
             append_dataset(root["Points"], self.vertices)
             append_dataset(
                 root["NumberOfConnectivityIds"],
@@ -3309,13 +3355,14 @@ def append_dataset(dset, array):
             append_dataset(root["Offsets"], offsets)
 
             append_dataset(
-                root["Types"], np.full(self.n_elements, self._VTK_TETRA, dtype="uint8")
+                root["Types"], np.full(
+                    self.n_elements, self._VTK_TETRA, dtype="uint8")
             )
 
             cell_data_group = root.create_group("CellData")
 
             for name, data in datasets.items():
-                
+
                 cell_data_group.create_dataset(
                     name, (0,), maxshape=(None,), dtype="float64", chunks=True
                 )
@@ -3453,4 +3500,4 @@ def _read_meshes(elem):
 _HEX_MIDPOINT_CONN += ((2, (0, 0, 0)),
                        (2, (1, 0, 0)),
                        (2, (1, 1, 0)),
-                       (2, (0, 1, 0)))
\ No newline at end of file
+                       (2, (0, 1, 0)))
diff --git a/tests/unit_tests/conftest.py b/tests/unit_tests/conftest.py
index 7844b82580e..165fc599120 100644
--- a/tests/unit_tests/conftest.py
+++ b/tests/unit_tests/conftest.py
@@ -59,7 +59,8 @@ def sphere_model():
     model.settings.particles = 100
     model.settings.batches = 10
     model.settings.run_mode = 'fixed source'
-    model.settings.source = openmc.IndependentSource(space=openmc.stats.Point())
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point())
     return model
 
 
@@ -83,6 +84,7 @@ def cell_with_lattice():
             [m_inside[0], m_inside[1], m_inside[3], m_outside],
             univ, lattice)
 
+
 @pytest.fixture
 def mixed_lattice_model(uo2, water):
     cyl = openmc.ZCylinder(r=0.4)
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index c44cf8db5b1..a617b5789c7 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -51,7 +51,7 @@ def test_regular_mesh_bounding_box():
     mesh.upper_right = [2, 3, 5]
     bb = mesh.bounding_box
     assert isinstance(bb, openmc.BoundingBox)
-    np.testing.assert_array_equal(bb.lower_left, (-2, -3 ,-5))
+    np.testing.assert_array_equal(bb.lower_left, (-2, -3, -5))
     np.testing.assert_array_equal(bb.upper_right, (2, 3, 5))
 
 
@@ -62,7 +62,7 @@ def test_rectilinear_mesh_bounding_box():
     mesh.z_grid = [-100., 0., 100.]
     bb = mesh.bounding_box
     assert isinstance(bb, openmc.BoundingBox)
-    np.testing.assert_array_equal(bb.lower_left, (0., -10. ,-100.))
+    np.testing.assert_array_equal(bb.lower_left, (0., -10., -100.))
     np.testing.assert_array_equal(bb.upper_right, (10., 0., 100.))
 
 
@@ -202,14 +202,16 @@ def test_invalid_cylindrical_mesh_errors():
         openmc.CylindricalMesh(r_grid=[5, 1], phi_grid=[0, pi], z_grid=[0, 10])
 
     with pytest.raises(ValueError):
-        openmc.CylindricalMesh(r_grid=[1, 2, 4, 3], phi_grid=[0, pi], z_grid=[0, 10])
+        openmc.CylindricalMesh(r_grid=[1, 2, 4, 3], phi_grid=[
+                               0, pi], z_grid=[0, 10])
 
     with pytest.raises(ValueError):
         openmc.CylindricalMesh(r_grid=[1], phi_grid=[0, pi], z_grid=[0, 10])
 
     # Test invalid phi_grid values
     with pytest.raises(ValueError):
-        openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[-1, 3], z_grid=[0, 10])
+        openmc.CylindricalMesh(
+            r_grid=[0, 1, 2], phi_grid=[-1, 3], z_grid=[0, 10])
 
     with pytest.raises(ValueError):
         openmc.CylindricalMesh(
@@ -226,10 +228,12 @@ def test_invalid_cylindrical_mesh_errors():
         openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[0, pi], z_grid=[5])
 
     with pytest.raises(ValueError):
-        openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[0, pi], z_grid=[5, 1])
+        openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[
+                               0, pi], z_grid=[5, 1])
 
     with pytest.raises(ValueError):
-        openmc.CylindricalMesh(r_grid=[1, 2, 4, 3], phi_grid=[0, pi], z_grid=[0, 10, 5])
+        openmc.CylindricalMesh(r_grid=[1, 2, 4, 3], phi_grid=[
+                               0, pi], z_grid=[0, 10, 5])
 
 
 def test_centroids():
@@ -248,21 +252,26 @@ def test_centroids():
     np.testing.assert_array_almost_equal(mesh.centroids[0, 0, 0], [6., 7., 8.])
 
     # cylindrical mesh
-    mesh = openmc.CylindricalMesh(r_grid=(0, 10), z_grid=(0, 10), phi_grid=(0, np.pi))
-    np.testing.assert_array_almost_equal(mesh.centroids[0, 0, 0], [0.0, 5.0, 5.0])
+    mesh = openmc.CylindricalMesh(
+        r_grid=(0, 10), z_grid=(0, 10), phi_grid=(0, np.pi))
+    np.testing.assert_array_almost_equal(
+        mesh.centroids[0, 0, 0], [0.0, 5.0, 5.0])
     # ensure that setting an origin is handled correctly
     mesh.origin = (5.0, 0, -10)
-    np.testing.assert_array_almost_equal(mesh.centroids[0, 0, 0], [5.0, 5.0, -5.0])
+    np.testing.assert_array_almost_equal(
+        mesh.centroids[0, 0, 0], [5.0, 5.0, -5.0])
 
     # spherical mesh, single element xyz-positive octant
-    mesh = openmc.SphericalMesh(r_grid=[0, 10], theta_grid=[0, 0.5*np.pi], phi_grid=[0, np.pi])
+    mesh = openmc.SphericalMesh(r_grid=[0, 10], theta_grid=[
+                                0, 0.5*np.pi], phi_grid=[0, np.pi])
     x = 5.*np.cos(0.5*np.pi)*np.sin(0.25*np.pi)
     y = 5.*np.sin(0.5*np.pi)*np.sin(0.25*np.pi)
     z = 5.*np.sin(0.25*np.pi)
     np.testing.assert_array_almost_equal(mesh.centroids[0, 0, 0], [x, y, z])
 
     mesh.origin = (-5.0, -5.0, 5.0)
-    np.testing.assert_array_almost_equal(mesh.centroids[0, 0, 0], [x-5.0, y-5.0, z+5.0])
+    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'))
@@ -285,7 +294,7 @@ def test_mesh_vertices(mesh_type):
         np.testing.assert_equal(mesh.vertices[ijk], exp_i_j_k)
 
         # shift the mesh using the llc
-        shift  = np.asarray((3.0, 6.0, 10.0))
+        shift = np.asarray((3.0, 6.0, 10.0))
         mesh.lower_left += shift
         np.testing.assert_equal(mesh.vertices[ijk], exp_i_j_k+shift)
     elif mesh_type == 'rectilinear':
@@ -302,15 +311,19 @@ def test_mesh_vertices(mesh_type):
         r_grid = np.linspace(0, 5, 10)
         z_grid = np.linspace(-10, 10, 20)
         phi_grid = np.linspace(0, 2*np.pi, 8)
-        mesh = openmc.CylindricalMesh(r_grid=r_grid, z_grid=z_grid, phi_grid=phi_grid)
-        exp_vert = np.asarray((mesh.r_grid[2], mesh.phi_grid[3], mesh.z_grid[2]))
+        mesh = openmc.CylindricalMesh(
+            r_grid=r_grid, z_grid=z_grid, phi_grid=phi_grid)
+        exp_vert = np.asarray(
+            (mesh.r_grid[2], mesh.phi_grid[3], mesh.z_grid[2]))
         np.testing.assert_equal(mesh.vertices_cylindrical[ijk], exp_vert)
     elif mesh_type == 'spherical':
         r_grid = np.linspace(0, 13, 14)
         theta_grid = np.linspace(0, np.pi, 11)
         phi_grid = np.linspace(0, 2*np.pi, 7)
-        mesh = openmc.SphericalMesh(r_grid=r_grid, theta_grid=theta_grid, phi_grid=phi_grid)
-        exp_vert = np.asarray((mesh.r_grid[2], mesh.theta_grid[3], mesh.phi_grid[2]))
+        mesh = openmc.SphericalMesh(
+            r_grid=r_grid, theta_grid=theta_grid, phi_grid=phi_grid)
+        exp_vert = np.asarray(
+            (mesh.r_grid[2], mesh.theta_grid[3], mesh.phi_grid[2]))
         np.testing.assert_equal(mesh.vertices_spherical[ijk], exp_vert)
 
 
@@ -325,9 +338,11 @@ def test_CylindricalMesh_get_indices_at_coords():
     assert mesh.get_indices_at_coords([-2, -2, 9]) == (0, 0, 1)
 
     with pytest.raises(ValueError):
-        assert mesh.get_indices_at_coords([8, 8, 1])  # resulting r value to large
+        assert mesh.get_indices_at_coords(
+            [8, 8, 1])  # resulting r value to large
     with pytest.raises(ValueError):
-        assert mesh.get_indices_at_coords([-8, -8, 1])  # resulting r value to large
+        assert mesh.get_indices_at_coords(
+            [-8, -8, 1])  # resulting r value to large
     with pytest.raises(ValueError):
         assert mesh.get_indices_at_coords([1, 0, -1])  # z value below range
     with pytest.raises(ValueError):
@@ -341,11 +356,16 @@ def test_CylindricalMesh_get_indices_at_coords():
         phi_grid=(0, 0.5 * pi, pi, 1.5 * pi, 1.9 * pi),
         z_grid=(-5, 0, 5, 10),
     )
-    assert mesh.get_indices_at_coords([1, 1, 1]) == (0, 0, 1)  # first angle quadrant
-    assert mesh.get_indices_at_coords([2, 2, 6]) == (0, 0, 2)  # first angle quadrant
-    assert mesh.get_indices_at_coords([-2, 0.1, -1]) == (0, 1, 0)  # second angle quadrant
-    assert mesh.get_indices_at_coords([-2, -0.1, -1]) == (0, 2, 0)  # third angle quadrant
-    assert mesh.get_indices_at_coords([2, -0.9, -1]) == (0, 3, 0)  # forth angle quadrant
+    assert mesh.get_indices_at_coords([1, 1, 1]) == (
+        0, 0, 1)  # first angle quadrant
+    assert mesh.get_indices_at_coords([2, 2, 6]) == (
+        0, 0, 2)  # first angle quadrant
+    assert mesh.get_indices_at_coords(
+        [-2, 0.1, -1]) == (0, 1, 0)  # second angle quadrant
+    assert mesh.get_indices_at_coords(
+        [-2, -0.1, -1]) == (0, 2, 0)  # third angle quadrant
+    assert mesh.get_indices_at_coords(
+        [2, -0.9, -1]) == (0, 3, 0)  # forth angle quadrant
 
     with pytest.raises(ValueError):
         assert mesh.get_indices_at_coords([2, -0.1, 1])  # outside of phi range
@@ -357,11 +377,16 @@ def test_CylindricalMesh_get_indices_at_coords():
         z_grid=(-5, 0, 5, 10),
         origin=(100, 200, 300),
     )
-    assert mesh.get_indices_at_coords([101, 201, 301]) == (0, 0, 1)  # first angle quadrant
-    assert mesh.get_indices_at_coords([102, 202, 306]) == (0, 0, 2)  # first angle quadrant
-    assert mesh.get_indices_at_coords([98, 200.1, 299]) == (0, 1, 0)  # second angle quadrant
-    assert mesh.get_indices_at_coords([98, 199.9, 299]) == (0, 2, 0)  # third angle quadrant
-    assert mesh.get_indices_at_coords([102, 199.1, 299]) == (0, 3, 0)  # forth angle quadrant
+    assert mesh.get_indices_at_coords([101, 201, 301]) == (
+        0, 0, 1)  # first angle quadrant
+    assert mesh.get_indices_at_coords([102, 202, 306]) == (
+        0, 0, 2)  # first angle quadrant
+    assert mesh.get_indices_at_coords([98, 200.1, 299]) == (
+        0, 1, 0)  # second angle quadrant
+    assert mesh.get_indices_at_coords([98, 199.9, 299]) == (
+        0, 2, 0)  # third angle quadrant
+    assert mesh.get_indices_at_coords([102, 199.1, 299]) == (
+        0, 3, 0)  # forth angle quadrant
 
 
 def test_mesh_name_roundtrip(run_in_tmpdir):
@@ -386,7 +411,8 @@ def test_mesh_name_roundtrip(run_in_tmpdir):
 
 
 def test_umesh_roundtrip(run_in_tmpdir, request):
-    umesh = openmc.UnstructuredMesh(request.path.parent / 'test_mesh_tets.e', 'moab')
+    umesh = openmc.UnstructuredMesh(
+        request.path.parent / 'test_mesh_tets.e', 'moab')
     umesh.output = True
 
     # create a tally using this mesh
@@ -418,10 +444,10 @@ def simple_umesh(request):
     geometry = openmc.Geometry([cell1])
 
     umesh = openmc.UnstructuredMesh(
-       filename=request.path.parent.parent
+        filename=request.path.parent.parent
         / "regression_tests/external_moab/test_mesh_tets.h5m",
-       library="moab",
-       mesh_id=1
+        library="moab",
+        mesh_id=1
     )
     # setting ID to make it easier to get the mesh from the statepoint later
     mesh_filter = openmc.MeshFilter(umesh)
@@ -467,7 +493,8 @@ def test_umesh(run_in_tmpdir, simple_umesh, export_type):
     rng = np.random.default_rng()
     ref_data = rng.random(simple_umesh.dimension)
     filename = f"test_mesh{export_type}"
-    simple_umesh.write_data_to_vtk(datasets={"mean": ref_data}, filename=filename)
+    simple_umesh.write_data_to_vtk(
+        datasets={"mean": ref_data}, filename=filename)
 
     assert Path(filename).exists()
 
@@ -485,7 +512,8 @@ def test_umesh(run_in_tmpdir, simple_umesh, export_type):
 
     # attempt to apply a dataset with an improper size to a VTK write
     with pytest.raises(ValueError, match='Cannot apply dataset "mean"'):
-        simple_umesh.write_data_to_vtk(datasets={'mean': ref_data[:-2]}, filename=filename)
+        simple_umesh.write_data_to_vtk(
+            datasets={'mean': ref_data[:-2]}, filename=filename)
 
 
 @pytest.mark.skipif(not openmc.lib._dagmc_enabled(), reason="DAGMC not enabled.")
@@ -504,7 +532,7 @@ def test_write_vtkhdf(request, run_in_tmpdir):
     umesh = openmc.UnstructuredMesh(
         request.path.parent / "test_mesh_dagmc_tets.vtk",
         "moab",
-        mesh_id = 1
+        mesh_id=1
     )
     mesh_filter = openmc.MeshFilter(umesh)
 
@@ -526,13 +554,15 @@ def test_write_vtkhdf(request, run_in_tmpdir):
 
     umesh_from_sp = statepoint.meshes[umesh.id]
 
-    datasets={
+    datasets = {
         "mean": my_tally.mean.flatten(),
         "std_dev": my_tally.std_dev.flatten()
     }
 
-    umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtkhdf")
-    umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtk")
+    umesh_from_sp.write_data_to_vtk(
+        datasets=datasets, filename="test_mesh.vtkhdf")
+    umesh_from_sp.write_data_to_vtk(
+        datasets=datasets, filename="test_mesh.vtk")
 
     with pytest.raises(ValueError, match="Unsupported file extension"):
         # Supported file extensions are vtk or vtkhdf, not hdf5, so this should raise an error
@@ -543,7 +573,7 @@ def test_write_vtkhdf(request, run_in_tmpdir):
     with pytest.raises(ValueError, match="Cannot apply dataset"):
         # The shape of the data should match the shape of the mesh, so this should raise an error
         umesh_from_sp.write_data_to_vtk(
-            datasets={'incorrectly_shaped_data': np.array(([1,2,3]))},
+            datasets={'incorrectly_shaped_data': np.array(([1, 2, 3]))},
             filename="test_mesh_incorrect_shape.vtkhdf",
         )
 
@@ -554,6 +584,7 @@ def test_write_vtkhdf(request, run_in_tmpdir):
     with h5py.File("test_mesh.vtkhdf", "r"):
         ...
 
+
 def test_mesh_get_homogenized_materials():
     """Test the get_homogenized_materials method"""
     # Simple model with 1 cm of Fe56 next to 1 cm of H1
@@ -638,9 +669,12 @@ def test_material_volumes_regular_mesh(sphere_model, n_rays):
     mesh.dimension = (2, 2, 2)
     volumes = mesh.material_volumes(sphere_model, n_rays)
     mats = sphere_model.materials
-    np.testing.assert_almost_equal(volumes[mats[0].id], [0., 0., 0., 0., 0., 1., 0., 1.])
-    np.testing.assert_almost_equal(volumes[mats[1].id], [0., 0., 0., 0., 1., 0., 1., 0.])
-    np.testing.assert_almost_equal(volumes[mats[2].id], [1., 1., 1., 1., 0., 0., 0., 0.])
+    np.testing.assert_almost_equal(
+        volumes[mats[0].id], [0., 0., 0., 0., 0., 1., 0., 1.])
+    np.testing.assert_almost_equal(
+        volumes[mats[1].id], [0., 0., 0., 0., 1., 0., 1., 0.])
+    np.testing.assert_almost_equal(
+        volumes[mats[2].id], [1., 1., 1., 1., 0., 0., 0., 0.])
     assert volumes.by_element(4) == [(mats[1].id, 1.)]
     assert volumes.by_element(0) == [(mats[2].id, 1.)]
 
@@ -694,13 +728,14 @@ def test_mesh_material_volumes_boundary_conditions(sphere_model):
     """Test the material volumes method using a regular mesh
     that overlaps with a vacuum boundary condition."""
 
-    mesh = openmc.SphericalMesh.from_domain(sphere_model.geometry, dimension=(1, 1, 1))
+    mesh = openmc.SphericalMesh.from_domain(
+        sphere_model.geometry, dimension=(1, 1, 1))
     # extend mesh beyond the outer sphere surface to test rays crossing the boundary condition
     mesh.r_grid[-1] += 5.0
 
     # add a new cell to the modelthat occupies the outside of the sphere
     sphere_surfaces = list(filter(lambda s: isinstance(s, openmc.Sphere),
-                            sphere_model.geometry.get_all_surfaces().values()))
+                                  sphere_model.geometry.get_all_surfaces().values()))
     outer_cell = openmc.Cell(region=+sphere_surfaces[0])
     sphere_model.geometry.root_universe.add_cell(outer_cell)
 
@@ -744,7 +779,7 @@ def test_raytrace_mesh_infinite_loop(run_in_tmpdir):
     )
     model.settings.run_mode = 'fixed source'
     model.settings.particles = 10
-    model.settings.batches =  1
+    model.settings.batches = 1
 
     # Run the model; this should not cause an infinite loop
     model.run()
@@ -905,6 +940,7 @@ def test_regular_mesh_get_indices_at_coords():
     assert len(result_1d) == 1
     assert result_1d == (5,)
 
+
 def test_write_vtkhdf_regular_mesh(run_in_tmpdir):
     """Test writing a regular mesh to VTKHDF format."""
     mesh = openmc.RegularMesh()
@@ -1193,4 +1229,4 @@ def test_write_ascii_vtk_unchanged(run_in_tmpdir):
     # Get data from file
     arr = reader.GetOutput().GetCellData().GetArray("data")
     read_data = np.array([arr.GetTuple1(i) for i in range(ref_data.size)])
-    np.testing.assert_almost_equal(read_data, ref_data.ravel())
\ No newline at end of file
+    np.testing.assert_almost_equal(read_data, ref_data.ravel())

From fa0f1264f881e3380d9152a42cc97334ab14f4b6 Mon Sep 17 00:00:00 2001
From: Jarvis2001 <moremakarand2001@gmail.com>
Date: Sat, 14 Mar 2026 00:29:00 +0530
Subject: [PATCH 03/20] Some reverting and small fix in filter.py

---
 openmc/filter.py              | 25 ++++++----
 openmc/mesh.py                | 37 +++++++++-----
 tests/unit_tests/test_mesh.py | 94 ++++++++++++++++++++++++++++++++++-
 3 files changed, 133 insertions(+), 23 deletions(-)

diff --git a/openmc/filter.py b/openmc/filter.py
index 87aeb70c3a7..f70c536f3a9 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -37,7 +37,6 @@
 )
 
 
-
 class FilterMeta(ABCMeta):
     """Metaclass for filters that ensures class names are appropriate."""
 
@@ -1113,12 +1112,20 @@ def from_volumes(cls, mesh: openmc.MeshBase, volumes: openmc.MeshMaterialVolumes
             A new MeshMaterialFilter instance
 
         """
-        # Get flat arrays of material IDs and element indices
-        mat_ids = volumes._materials[volumes._materials > -1]
-        elems, _ = np.where(volumes._materials > -1)
-
-        # Stack them into a 2D array of (element, material) pairs
-        bins = np.column_stack((elems, mat_ids))
+        # Build bins sorted by volume ascending (mat ID ascending as tiebreaker)
+        # so the bin ordering is deterministic regardless of the ray-traversal
+        # order stored in the raw _materials array.
+        bins = []
+        for i in range(volumes.num_elements):
+            entries = volumes.by_element(i)
+            # sort by (volume asc, mat_id asc), excluding void (None)
+            entries = sorted(
+                ((mat_id, vol)
+                 for mat_id, vol in entries if mat_id is not None),
+                key=lambda t: (t[1], t[0])
+            )
+            for mat_id, _ in entries:
+                bins.append((i, mat_id))
         return cls(mesh, bins)
 
     def __hash__(self):
@@ -1861,7 +1868,7 @@ def __init__(self, particles, energies=None, filter_id=None):
     def __repr__(self):
         string = type(self).__name__ + '\n'
         string += '{: <16}=\t{}\n'.format('\tParticles',
-            [str(p) for p in self.particles])
+                                          [str(p) for p in self.particles])
         if self.energies is not None:
             string += '{: <16}=\t{}\n'.format('\tEnergies', self.energies)
         string += '{: <16}=\t{}\n'.format('\tID', self.id)
@@ -2171,7 +2178,7 @@ def paths(self):
         if self._paths is None:
             if not hasattr(self, '_geometry'):
                 raise ValueError(
-                    "Model must be exported before the 'paths' attribute is" \
+                    "Model must be exported before the 'paths' attribute is"
                     "available for a DistribcellFilter.")
 
             # Determine paths for cell instances
diff --git a/openmc/mesh.py b/openmc/mesh.py
index 3c94784caad..216cdf21c44 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -693,7 +693,7 @@ def num_mesh_cells(self):
     def write_data_to_vtk(self,
                           filename: PathLike,
                           datasets: dict | None = None,
-                          volume_normalization: bool = False,
+                          volume_normalization: bool | None = None,
                           curvilinear: bool = False):
         """Creates a VTK object of the mesh
 
@@ -709,9 +709,10 @@ def write_data_to_vtk(self,
             with structured indexing in "C" ordering. See the "expand_dims" flag
             of :meth:`~openmc.Tally.get_reshaped_data` on reshaping tally data when using
             :class:`~openmc.MeshFilter`'s.
-        volume_normalization : bool, optional
+        volume_normalization : bool or None, optional
             Whether or not to normalize the data by the volume of the mesh
-            elements.
+            elements. When None (default), the format-appropriate default is
+            used: True for legacy ASCII .vtk files, False for .vtkhdf files.
         curvilinear : bool
             Whether or not to write curvilinear elements. Only applies to
             ``SphericalMesh`` and ``CylindricalMesh``.
@@ -748,10 +749,14 @@ def write_data_to_vtk(self,
             if write_impl is None:
                 raise NotImplementedError(
                     f"VTKHDF output not implemented for {type(self).__name__}")
-            # write_impl is a bound method – do NOT pass self again
-            write_impl(filename, datasets, volume_normalization)
+            # vtkhdf default is False; explicit value is respected
+            norm = volume_normalization if volume_normalization is not None else False
+            write_impl(filename, datasets, norm)
             return None
 
+        # legacy ASCII .vtk default is True
+        norm = volume_normalization if volume_normalization is not None else True
+
         # vtk is an optional dependency only needed for the legacy ASCII path
         import vtk
         from vtk.util import numpy_support as nps
@@ -769,20 +774,26 @@ def write_data_to_vtk(self,
             # maintain a list of the datasets as added to the VTK arrays to
             # ensure they persist in memory until the file is written
             datasets_out = []
+            # Regular/Rectilinear meshes store data in C (ijk) order which
+            # matches VTK's expected ordering directly — no transpose needed.
+            # Curvilinear meshes (Cylindrical, Spherical) require a transpose
+            # to convert from C ordering to the Fortran (kji) ordering that VTK
+            # expects.
+            # TODO: update to "C" ordering throughout
+            needs_transpose = not isinstance(
+                self, (RegularMesh, RectilinearMesh))
             for label, dataset in datasets.items():
                 dataset = self._reshape_vtk_dataset(dataset)
                 self._check_vtk_dataset(label, dataset)
-                # If the array data is 3D, assume is in C ordering and transpose
-                # before flattening to match the ordering expected by the VTK
-                # array based on the way mesh indices are ordered in the Python
-                # API
-                # TODO: update to "C" ordering throughout
                 if dataset.ndim == 3:
-                    dataset = dataset.ravel()
+                    dataset = dataset.T.ravel() if needs_transpose else dataset.ravel()
                 datasets_out.append(dataset)
 
-                if volume_normalization:
-                    dataset /= self.volumes.ravel()
+                if norm:
+                    if needs_transpose:
+                        dataset /= self.volumes.T.ravel()
+                    else:
+                        dataset /= self.volumes.ravel()
 
                 dataset_array = vtk.vtkDoubleArray()
                 dataset_array.SetName(label)
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index a617b5789c7..1a91bf96d6e 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -7,6 +7,8 @@
 from scipy.stats import chi2
 import pytest
 import openmc
+import random
+import itertools
 import openmc.lib
 from openmc.utility_funcs import change_directory
 from uncertainties.unumpy import uarray, nominal_values, std_devs
@@ -723,6 +725,48 @@ def test_mesh_material_volumes_serialize():
     assert new_volumes.by_element(2) == [(2, 0.5), (1, 0.5)]
     assert new_volumes.by_element(3) == [(2, 1.0)]
 
+def test_mesh_material_volumes_serialize_with_bboxes():
+    materials = np.array([
+        [1, -1, -2],
+        [-1, -2, -2],
+        [2, 1, -2],
+        [2, -2, -2]
+    ])
+    volumes = np.array([
+        [0.5, 0.5, 0.0],
+        [1.0, 0.0, 0.0],
+        [0.5, 0.5, 0.0],
+        [1.0, 0.0, 0.0]
+    ])
+
+    # (xmin, ymin, zmin, xmax, ymax, zmax)
+    bboxes = np.empty((4, 3, 6))
+    bboxes[..., 0:3] = np.inf
+    bboxes[..., 3:6] = -np.inf
+    bboxes[0, 0] = [-1.0, -2.0, -3.0, 1.0, 2.0, 3.0]    # material 1
+    bboxes[0, 1] = [-5.0, -6.0, -7.0, 5.0, 6.0, 7.0]    # void
+    bboxes[1, 0] = [0.0, 0.0, 0.0, 10.0, 1.0, 2.0]      # void
+    bboxes[2, 0] = [-1.0, -1.0, -1.0, 0.0, 0.0, 0.0]    # material 2
+    bboxes[2, 1] = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]       # material 1
+    bboxes[3, 0] = [-2.0, -2.0, -2.0, 2.0, 2.0, 2.0]    # material 2
+
+    mmv = openmc.MeshMaterialVolumes(materials, volumes, bboxes)
+    with TemporaryDirectory() as tmpdir:
+        path = f'{tmpdir}/volumes_bboxes.npz'
+        mmv.save(path)
+        loaded = openmc.MeshMaterialVolumes.from_npz(path)
+
+    assert loaded.has_bounding_boxes
+    first = loaded.by_element(0, include_bboxes=True)[0][2]
+    assert isinstance(first, openmc.BoundingBox)
+    np.testing.assert_array_equal(first.lower_left, (-1.0, -2.0, -3.0))
+    np.testing.assert_array_equal(first.upper_right, (1.0, 2.0, 3.0))
+
+    second = loaded.by_element(0, include_bboxes=True)[1][2]
+    assert isinstance(second, openmc.BoundingBox)
+    np.testing.assert_array_equal(second.lower_left, (-5.0, -6.0, -7.0))
+    np.testing.assert_array_equal(second.upper_right, (5.0, 6.0, 7.0))
+
 
 def test_mesh_material_volumes_boundary_conditions(sphere_model):
     """Test the material volumes method using a regular mesh
@@ -751,6 +795,51 @@ def test_mesh_material_volumes_boundary_conditions(sphere_model):
         assert evaluated[0] == expected[0]
         assert evaluated[1] == pytest.approx(expected[1], rel=1e-2)
 
+def test_mesh_material_volumes_bounding_boxes():
+    # Create a model with 8 spherical cells at known locations with random radii
+    box = openmc.model.RectangularParallelepiped(
+        -10, 10, -10, 10, -10, 10, boundary_type='vacuum')
+
+    mat = openmc.Material()
+    mat.add_nuclide('H1', 1.0)
+
+    sph_cells = []
+    for x, y, z in itertools.product((-5., 5.), repeat=3):
+        mat_i = mat.clone()
+        sph = openmc.Sphere(x, y, z, r=random.uniform(0.5, 1.5))
+        sph_cells.append(openmc.Cell(region=-sph, fill=mat_i))
+    background = openmc.Cell(region=-box & openmc.Intersection([~c.region for c in sph_cells]))
+
+    model = openmc.Model()
+    model.geometry = openmc.Geometry(sph_cells + [background])
+    model.settings.particles = 1000
+    model.settings.batches = 10
+
+    # Create a one-element mesh that encompasses the entire geometry
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = (-10., -10., -10.)
+    mesh.upper_right = (10., 10., 10.)
+    mesh.dimension = (1, 1, 1)
+
+    # Run material volume calculation with bounding boxes
+    n_samples = (400, 400, 400)
+    mmv = mesh.material_volumes(model, n_samples, max_materials=10, bounding_boxes=True)
+    assert mmv.has_bounding_boxes
+
+    # Create a mapping of material ID to bounding box
+    bbox_by_mat = {
+        mat_id: bbox
+        for mat_id, vol, bbox in mmv.by_element(0, include_bboxes=True)
+        if mat_id is not None and vol > 0.0
+    }
+
+    # Match the mesh ray spacing used for the bounding box estimator.
+    tol = 0.5 * mesh.bounding_box.width[0] / n_samples[0]
+    for cell in sph_cells:
+        bbox = bbox_by_mat[cell.fill.id]
+        cell_bbox = cell.bounding_box
+        np.testing.assert_allclose(bbox.lower_left, cell_bbox.lower_left, atol=tol)
+        np.testing.assert_allclose(bbox.upper_right, cell_bbox.upper_right, atol=tol)
 
 def test_raytrace_mesh_infinite_loop(run_in_tmpdir):
     # Create a model with one large spherical cell
@@ -1217,7 +1306,10 @@ def test_write_ascii_vtk_unchanged(run_in_tmpdir):
 
     # This should work without requiring vtk module changes
     vtkIOLegacy = pytest.importorskip("vtkmodules.vtkIOLegacy")
-    mesh.write_data_to_vtk(datasets={"data": ref_data}, filename=filename)
+    mesh.write_data_to_vtk(datasets={"data": ref_data},
+                           filename=filename,
+                           volume_normalization=False
+                           )
 
     assert Path(filename).exists()
 

From 357bf6286bde1294b5bab1ec7b6782cec1403c56 Mon Sep 17 00:00:00 2001
From: Makarand More <40858007+Jarvis2001@users.noreply.github.com>
Date: Sat, 14 Mar 2026 22:44:41 +0530
Subject: [PATCH 04/20] Update test_mesh.py

---
 tests/unit_tests/test_mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index b69f5d2e7f5..c4fbc2de656 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1353,4 +1353,4 @@ def test_rectilinear_mesh_get_indices_at_coords():
     with pytest.raises(ValueError):
         mesh.get_indices_at_coords([0.5, -0.5, 110.])
     with pytest.raises(ValueError):
-        mesh.get_indices_at_coords([0.5, -20., 110
+        mesh.get_indices_at_coords([0.5, -20., 110.])

From 7017c1dbd37d8bb9f3032134e56be71267302058 Mon Sep 17 00:00:00 2001
From: Jarvis2001 <moremakarand2001@gmail.com>
Date: Sat, 21 Mar 2026 05:23:56 +0530
Subject: [PATCH 05/20] Reverted Formatting

---
 openmc/filter.py              |  20 +--
 openmc/mesh.py                | 218 +++++++++++----------------------
 tests/unit_tests/conftest.py  |  15 +--
 tests/unit_tests/test_mesh.py | 222 +++++++++++++++++++++-------------
 4 files changed, 218 insertions(+), 257 deletions(-)

diff --git a/openmc/filter.py b/openmc/filter.py
index f70c536f3a9..9472ca7ef93 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -1112,20 +1112,12 @@ def from_volumes(cls, mesh: openmc.MeshBase, volumes: openmc.MeshMaterialVolumes
             A new MeshMaterialFilter instance
 
         """
-        # Build bins sorted by volume ascending (mat ID ascending as tiebreaker)
-        # so the bin ordering is deterministic regardless of the ray-traversal
-        # order stored in the raw _materials array.
-        bins = []
-        for i in range(volumes.num_elements):
-            entries = volumes.by_element(i)
-            # sort by (volume asc, mat_id asc), excluding void (None)
-            entries = sorted(
-                ((mat_id, vol)
-                 for mat_id, vol in entries if mat_id is not None),
-                key=lambda t: (t[1], t[0])
-            )
-            for mat_id, _ in entries:
-                bins.append((i, mat_id))
+        # Get flat arrays of material IDs and element indices
+        mat_ids = volumes._materials[volumes._materials > -1]
+        elems, _ = np.where(volumes._materials > -1)
+
+        # Stack them into a 2D array of (element, material) pairs
+        bins = np.column_stack((elems, mat_ids))
         return cls(mesh, bins)
 
     def __hash__(self):
diff --git a/openmc/mesh.py b/openmc/mesh.py
index a9728307a37..50fbc4203fa 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -70,7 +70,6 @@ class MeshMaterialVolumes(Mapping):
     [(2, 31.87963824195591), (1, 6.129949130817542)]
 
     """
-
     def __init__(
         self,
         materials: np.ndarray,
@@ -142,8 +141,7 @@ def by_element(
         """
         table_size = self._volumes.shape[1]
         if include_bboxes and self._bboxes is None:
-            raise ValueError(
-                'Bounding boxes were not computed for this object.')
+            raise ValueError('Bounding boxes were not computed for this object.')
 
         results = []
         for i in range(table_size):
@@ -290,8 +288,7 @@ def from_hdf5(cls, group: h5py.Group):
             Instance of a MeshBase subclass
 
         """
-        mesh_type = 'regular' if 'type' not in group.keys(
-        ) else group['type'][()].decode()
+        mesh_type = 'regular' if 'type' not in group.keys() else group['type'][()].decode()
         mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
         mesh_name = '' if not 'name' in group else group['name'][()].decode()
 
@@ -399,8 +396,7 @@ def get_homogenized_materials(
             vols = self.material_volumes(model, n_samples, **kwargs)
         else:
             vols = material_volumes
-        mat_volume_by_element = [vols.by_element(
-            i) for i in range(vols.num_elements)]
+        mat_volume_by_element = [vols.by_element(i) for i in range(vols.num_elements)]
 
         # Get dictionary of all materials
         materials = model._get_all_materials()
@@ -506,33 +502,7 @@ def material_volumes(
 
         # Restore original tallies
         model.tallies = original_tallies
-
-        # Sort each element's materials by material ID (ascending), with
-        # None (-1) after positive IDs and empty slots (-2) last.  This
-        # gives a deterministic ordering that is independent of the ray
-        # traversal order used by the C library.
-        mat_arr = volumes._materials.copy()
-        vol_arr = volumes._volumes.copy()
-        bbox_arr = volumes._bboxes.copy() if volumes._bboxes is not None else None
-        n_elements, table_size = mat_arr.shape
-
-        def _sort_key(m):
-            m = int(m)
-            if m == -2:
-                return (2, 0)   # empty slot – always last
-            if m == -1:
-                return (1, 0)   # vacuum / None – just before empty
-            return (0, m)       # real material – ascending ID
-
-        for i in range(n_elements):
-            indices = sorted(range(table_size),
-                             key=lambda j: _sort_key(mat_arr[i, j]))
-            mat_arr[i] = mat_arr[i, indices]
-            vol_arr[i] = vol_arr[i, indices]
-            if bbox_arr is not None:
-                bbox_arr[i] = bbox_arr[i, indices]
-
-        return MeshMaterialVolumes(mat_arr, vol_arr, bbox_arr)
+        return volumes
 
 
 class StructuredMesh(MeshBase):
@@ -638,8 +608,7 @@ def _generate_edge_midpoints(grids):
 
             # re-use the generate vertices method to create the full mesh grid
             # transpose to get (i, j, k) ordering of the gridpoints
-            midpoint_grid = StructuredMesh._generate_vertices(
-                i_grid, j_grid, k_grid)
+            midpoint_grid = StructuredMesh._generate_vertices(i_grid, j_grid, k_grid)
             midpoint_grids.append(midpoint_grid)
 
         return midpoint_grids
@@ -747,11 +716,9 @@ def write_data_to_vtk(self,
         if Path(filename).suffix == ".vtkhdf":
             write_impl = getattr(self, '_write_vtk_hdf5', None)
             if write_impl is None:
-                raise NotImplementedError(
-                    f"VTKHDF output not implemented for {type(self).__name__}")
-            # vtkhdf default is False; explicit value is respected
-            norm = volume_normalization if volume_normalization is not None else False
-            write_impl(filename, datasets, norm)
+                raise NotImplementedError(f"VTKHDF output not implemented for {type(self).__name__}")
+            # write_impl is a bound method – do NOT pass self again
+            write_impl(filename, datasets, volume_normalization)
             return None
 
         # legacy ASCII .vtk default is True
@@ -797,8 +764,7 @@ def write_data_to_vtk(self,
 
                 dataset_array = vtk.vtkDoubleArray()
                 dataset_array.SetName(label)
-                dataset_array.SetArray(
-                    nps.numpy_to_vtk(dataset), dataset.size, True)
+                dataset_array.SetArray(nps.numpy_to_vtk(dataset), dataset.size, True)
                 vtk_grid.GetCellData().AddArray(dataset_array)
 
         writer.SetFileName(str(filename))
@@ -819,8 +785,7 @@ def _create_vtk_structured_grid(self):
         from vtk.util import numpy_support as nps
 
         vtkPts = vtk.vtkPoints()
-        vtkPts.SetData(nps.numpy_to_vtk(np.swapaxes(
-            self.vertices, 0, 2).reshape(-1, 3), deep=True))
+        vtkPts.SetData(nps.numpy_to_vtk(np.swapaxes(self.vertices, 0, 2).reshape(-1, 3), deep=True))
         vtk_grid = vtk.vtkStructuredGrid()
         vtk_grid.SetPoints(vtkPts)
         vtk_grid.SetDimensions(*[dim + 1 for dim in self.dimension])
@@ -847,13 +812,12 @@ def _create_vtk_unstructured_grid(self):
         # add corner vertices to the point set for the unstructured grid
         # only insert unique points, we'll get their IDs in the point set to
         # define element connectivity later
-        vtkPts.SetData(nps.numpy_to_vtk(
-            np.unique(corner_vertices, axis=0), deep=True))
+        vtkPts.SetData(nps.numpy_to_vtk(np.unique(corner_vertices, axis=0), deep=True))
 
         # create a locator to assist with duplicate points
         locator = vtk.vtkPointLocator()
         locator.SetDataSet(vtk_grid)
-        locator.AutomaticOn()  # autmoatically adds points to locator
+        locator.AutomaticOn() # autmoatically adds points to locator
         locator.InitPointInsertion(vtkPts, vtkPts.GetBounds())
         locator.BuildLocator()
 
@@ -905,8 +869,7 @@ def _insert_point(pnt):
             for n, (di, dj, dk) in enumerate(_HEX_VERTEX_CONN):
                 # compute flat index into the point ID list based on i, j, k
                 # of the vertex
-                flat_idx = np.ravel_multi_index(
-                    (i+di, j+dj, k+dk), n_pnts, order='F')
+                flat_idx = np.ravel_multi_index((i+di, j+dj, k+dk), n_pnts, order='F')
                 # set corner vertices
                 hex.GetPointIds().SetId(n, point_ids[flat_idx])
 
@@ -914,16 +877,14 @@ def _insert_point(pnt):
             n_midpoint_vertices = [v.size // 3 for v in midpoint_vertices]
             for n, (dim, (di, dj, dk)) in enumerate(_HEX_MIDPOINT_CONN):
                 # initial offset for corner vertices and midpoint dimension
-                flat_idx = corner_vertices.shape[0] + \
-                    sum(n_midpoint_vertices[:dim])
+                flat_idx = corner_vertices.shape[0] + sum(n_midpoint_vertices[:dim])
                 # generate a flat index into the table of point IDs
                 midpoint_shape = midpoint_vertices[dim].shape[:-1]
                 flat_idx += np.ravel_multi_index((i+di, j+dj, k+dk),
                                                  midpoint_shape,
                                                  order='F')
                 # set hex midpoint connectivity
-                hex.GetPointIds().SetId(
-                    _N_HEX_VERTICES + n, point_ids[flat_idx])
+                hex.GetPointIds().SetId(_N_HEX_VERTICES + n, point_ids[flat_idx])
 
             # add the hex to the grid
             vtk_grid.InsertNextCell(hex.GetCellType(), hex.GetPointIds())
@@ -1077,8 +1038,7 @@ def lower_left(self, lower_left: Iterable[Real]):
         self._lower_left = lower_left
 
         if self.upper_right is not None and any(np.isclose(self.upper_right, lower_left)):
-            raise ValueError(
-                "Mesh cannot have zero thickness in any dimension")
+            raise ValueError("Mesh cannot have zero thickness in any dimension")
 
     @property
     def upper_right(self):
@@ -1102,8 +1062,7 @@ def upper_right(self, upper_right: Iterable[Real]):
             warnings.warn("Unsetting width attribute.")
 
         if self.lower_left is not None and any(np.isclose(self.lower_left, upper_right)):
-            raise ValueError(
-                "Mesh cannot have zero thickness in any dimension")
+            raise ValueError("Mesh cannot have zero thickness in any dimension")
 
     @property
     def width(self):
@@ -1113,7 +1072,7 @@ def width(self):
             if self._lower_left is not None and self._dimension is not None:
                 us = self._upper_right
                 ls = self._lower_left
-                dims = self._dimension
+                dims =  self._dimension
                 return [(u - l) / d for u, l, d in zip(us, ls, dims)]
 
     @width.setter
@@ -1187,13 +1146,10 @@ def _grids(self):
 
     def __repr__(self):
         string = super().__repr__()
-        string += '{0: <16}{1}{2}\n'.format('\tDimensions',
-                                            '=\t', self.n_dimension)
+        string += '{0: <16}{1}{2}\n'.format('\tDimensions', '=\t', self.n_dimension)
         string += '{0: <16}{1}{2}\n'.format('\tVoxels', '=\t', self._dimension)
-        string += '{0: <16}{1}{2}\n'.format('\tLower left',
-                                            '=\t', self._lower_left)
-        string += '{0: <16}{1}{2}\n'.format('\tUpper Right',
-                                            '=\t', self.upper_right)
+        string += '{0: <16}{1}{2}\n'.format('\tLower left', '=\t', self._lower_left)
+        string += '{0: <16}{1}{2}\n'.format('\tUpper Right', '=\t', self.upper_right)
         string += '{0: <16}{1}{2}\n'.format('\tWidth', '=\t', self.width)
         return string
 
@@ -1208,8 +1164,7 @@ def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
         elif 'upper_right' in group:
             mesh.upper_right = group['upper_right'][()]
         else:
-            raise IOError(
-                'Invalid mesh: must have one of "upper_right" or "width"')
+            raise IOError('Invalid mesh: must have one of "upper_right" or "width"')
 
         return mesh
 
@@ -1554,8 +1509,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
         # that every point has an (x, y, z) triple; extra coordinates are 0.
         coords_1d = []
         for i in range(ndim):
-            c = np.linspace(self.lower_left[i],
-                            self.upper_right[i], dims[i] + 1)
+            c = np.linspace(self.lower_left[i], self.upper_right[i], dims[i] + 1)
             coords_1d.append(c)
         while len(coords_1d) < 3:
             coords_1d.append(np.array([0.0]))
@@ -1745,20 +1699,17 @@ def __repr__(self):
         fmt = '{0: <16}{1}{2}\n'
         string = super().__repr__()
         string += fmt.format('\tDimensions', '=\t', self.n_dimension)
-        x_grid_str = str(
-            self._x_grid) if self._x_grid is None else len(self._x_grid)
+        x_grid_str = str(self._x_grid) if self._x_grid is None else len(self._x_grid)
         string += fmt.format('\tN X pnts:', '=\t', x_grid_str)
         if self._x_grid is not None:
             string += fmt.format('\tX Min:', '=\t', self._x_grid[0])
             string += fmt.format('\tX Max:', '=\t', self._x_grid[-1])
-        y_grid_str = str(
-            self._y_grid) if self._y_grid is None else len(self._y_grid)
+        y_grid_str = str(self._y_grid) if self._y_grid is None else len(self._y_grid)
         string += fmt.format('\tN Y pnts:', '=\t', y_grid_str)
         if self._y_grid is not None:
             string += fmt.format('\tY Min:', '=\t', self._y_grid[0])
             string += fmt.format('\tY Max:', '=\t', self._y_grid[-1])
-        z_grid_str = str(
-            self._z_grid) if self._z_grid is None else len(self._z_grid)
+        z_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid)
         string += fmt.format('\tN Z pnts:', '=\t', z_grid_str)
         if self._z_grid is not None:
             string += fmt.format('\tZ Min:', '=\t', self._z_grid[0])
@@ -1884,8 +1835,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
             points = vertices.reshape(-1, 3)
-            root.create_dataset(
-                "Points", data=points.astype(np.float64), dtype="f8")
+            root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
 
             cell_data_group = root.create_group("CellData")
 
@@ -2085,20 +2035,17 @@ def __repr__(self):
         string = super().__repr__()
         string += fmt.format('\tDimensions', '=\t', self.n_dimension)
         string += fmt.format('\tOrigin', '=\t', self.origin)
-        r_grid_str = str(
-            self._r_grid) if self._r_grid is None else len(self._r_grid)
+        r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
         string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
         if self._r_grid is not None:
             string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
             string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
-        phi_grid_str = str(
-            self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
+        phi_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
         string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
         if self._phi_grid is not None:
             string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
             string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1])
-        z_grid_str = str(
-            self._z_grid) if self._z_grid is None else len(self._z_grid)
+        z_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid)
         string += fmt.format('\tN Z pnts:', '=\t', z_grid_str)
         if self._z_grid is not None:
             string += fmt.format('\tZ Min:', '=\t', self._z_grid[0])
@@ -2106,9 +2053,9 @@ def __repr__(self):
         return string
 
     def get_indices_at_coords(
-        self,
-        coords: Sequence[float]
-    ) -> tuple[int, int, int]:
+            self,
+            coords: Sequence[float]
+        ) -> tuple[int, int, int]:
         """Finds the index of the mesh element at the specified coordinates.
 
         .. versionadded:: 0.15.0
@@ -2124,8 +2071,7 @@ def get_indices_at_coords(
             The r, phi, z indices
 
         """
-        r_value_from_origin = sqrt(
-            (coords[0]-self.origin[0])**2 + (coords[1]-self.origin[1])**2)
+        r_value_from_origin = sqrt((coords[0]-self.origin[0])**2 + (coords[1]-self.origin[1])**2)
 
         if r_value_from_origin < self.r_grid[0] or r_value_from_origin > self.r_grid[-1]:
             raise ValueError(
@@ -2174,9 +2120,9 @@ def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
         mesh = cls(
             mesh_id=mesh_id,
             name=name,
-            r_grid=group['r_grid'][()],
-            phi_grid=group['phi_grid'][()],
-            z_grid=group['z_grid'][()],
+            r_grid = group['r_grid'][()],
+            phi_grid = group['phi_grid'][()],
+            z_grid = group['z_grid'][()],
         )
         if 'origin' in group:
             mesh.origin = group['origin'][()]
@@ -2311,10 +2257,10 @@ def from_xml_element(cls, elem: ET.Element):
 
         mesh_id = int(get_text(elem, 'id'))
         mesh = cls(
-            r_grid=get_elem_list(elem, "r_grid", float),
-            phi_grid=get_elem_list(elem, "phi_grid", float),
-            z_grid=get_elem_list(elem, "z_grid", float),
-            origin=get_elem_list(elem, "origin", float) or [0., 0., 0.],
+            r_grid = get_elem_list(elem, "r_grid", float),
+            phi_grid = get_elem_list(elem, "phi_grid", float),
+            z_grid = get_elem_list(elem, "z_grid", float),
+            origin = get_elem_list(elem, "origin", float) or [0., 0., 0.],
             mesh_id=mesh_id,
         )
 
@@ -2339,8 +2285,7 @@ def volumes(self):
 
     @property
     def vertices(self):
-        warnings.warn(
-            'Cartesian coordinates are returned from this property as of version 0.14.0')
+        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
         return self._convert_to_cartesian(self.vertices_cylindrical, self.origin)
 
     @property
@@ -2351,8 +2296,7 @@ def vertices_cylindrical(self):
 
     @property
     def centroids(self):
-        warnings.warn(
-            'Cartesian coordinates are returned from this property as of version 0.14.0')
+        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
         return self._convert_to_cartesian(self.centroids_cylindrical, self.origin)
 
     @property
@@ -2378,8 +2322,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
         nr, nphi, nz = self.dimension
         vertex_dims = [nr + 1, nphi + 1, nz + 1]
 
-        R, Phi, Z = np.meshgrid(
-            self.r_grid, self.phi_grid, self.z_grid, indexing='ij')
+        R, Phi, Z = np.meshgrid(self.r_grid, self.phi_grid, self.z_grid, indexing='ij')
         X = R * np.cos(Phi) + self.origin[0]
         Y = R * np.sin(Phi) + self.origin[1]
         Z = Z + self.origin[2]
@@ -2392,8 +2335,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
             points = vertices.reshape(-1, 3)
-            root.create_dataset(
-                "Points", data=points.astype(np.float64), dtype="f8")
+            root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
 
             cell_data_group = root.create_group("CellData")
 
@@ -2591,20 +2533,17 @@ def __repr__(self):
         string = super().__repr__()
         string += fmt.format('\tDimensions', '=\t', self.n_dimension)
         string += fmt.format('\tOrigin', '=\t', self.origin)
-        r_grid_str = str(
-            self._r_grid) if self._r_grid is None else len(self._r_grid)
+        r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
         string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
         if self._r_grid is not None:
             string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
             string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
-        theta_grid_str = str(
-            self._theta_grid) if self._theta_grid is None else len(self._theta_grid)
+        theta_grid_str = str(self._theta_grid) if self._theta_grid is None else len(self._theta_grid)
         string += fmt.format('\tN Theta pnts:', '=\t', theta_grid_str)
         if self._theta_grid is not None:
             string += fmt.format('\tTheta Min:', '=\t', self._theta_grid[0])
             string += fmt.format('\tTheta Max:', '=\t', self._theta_grid[-1])
-        phi_grid_str = str(
-            self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
+        phi_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
         string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
         if self._phi_grid is not None:
             string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
@@ -2615,9 +2554,9 @@ def __repr__(self):
     def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
         # Read and assign mesh properties
         mesh = cls(
-            r_grid=group['r_grid'][()],
-            theta_grid=group['theta_grid'][()],
-            phi_grid=group['phi_grid'][()],
+            r_grid = group['r_grid'][()],
+            theta_grid = group['theta_grid'][()],
+            phi_grid = group['phi_grid'][()],
             mesh_id=mesh_id,
             name=name
         )
@@ -2745,10 +2684,10 @@ def from_xml_element(cls, elem: ET.Element):
         mesh_id = int(get_text(elem, 'id'))
         mesh = cls(
             mesh_id=mesh_id,
-            r_grid=get_elem_list(elem, "r_grid", float),
-            theta_grid=get_elem_list(elem, "theta_grid", float),
-            phi_grid=get_elem_list(elem, "phi_grid", float),
-            origin=get_elem_list(elem, "origin", float) or [0., 0., 0.],
+            r_grid = get_elem_list(elem, "r_grid", float),
+            theta_grid = get_elem_list(elem, "theta_grid", float),
+            phi_grid = get_elem_list(elem, "phi_grid", float),
+            origin = get_elem_list(elem, "origin", float) or [0., 0., 0.],
         )
 
         return mesh
@@ -2772,8 +2711,7 @@ def volumes(self):
 
     @property
     def vertices(self):
-        warnings.warn(
-            'Cartesian coordinates are returned from this property as of version 0.14.0')
+        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
         return self._convert_to_cartesian(self.vertices_spherical, self.origin)
 
     @property
@@ -2784,8 +2722,7 @@ def vertices_spherical(self):
 
     @property
     def centroids(self):
-        warnings.warn(
-            'Cartesian coordinates are returned from this property as of version 0.14.0')
+        warnings.warn('Cartesian coordinates are returned from this property as of version 0.14.0')
         return self._convert_to_cartesian(self.centroids_spherical, self.origin)
 
     @property
@@ -2794,6 +2731,7 @@ def centroids_spherical(self):
         """
         return super().centroids
 
+
     @staticmethod
     def _convert_to_cartesian(arr, origin: Sequence[float]):
         """Converts an array with r, theta, phi values in the last dimension (shape (..., 3))
@@ -2833,8 +2771,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
             points = vertices.reshape(-1, 3)
-            root.create_dataset(
-                "Points", data=points.astype(np.float64), dtype="f8")
+            root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
 
             cell_data_group = root.create_group("CellData")
 
@@ -3231,8 +3168,7 @@ def _write_data_to_vtk_ascii_format(
         from vtkmodules import vtkIOXML
 
         if self.connectivity is None or self.vertices is None:
-            raise RuntimeError(
-                "This mesh has not been loaded from a statepoint file.")
+            raise RuntimeError("This mesh has not been loaded from a statepoint file.")
 
         if filename is None:
             filename = f"mesh_{self.id}.vtk"
@@ -3321,7 +3257,7 @@ def _write_data_to_vtk_hdf5_format(
         volume_normalization: bool = True,
     ):
         """Write UnstructuredMesh as VTK-HDF5 UnstructuredGrid format.
-
+        
         Supports linear tetrahedra and linear hexahedra elements.
         """
         def append_dataset(dset, array):
@@ -3331,8 +3267,7 @@ def append_dataset(dset, array):
             dset[origLen:] = array
 
         if self.library != "moab":
-            raise NotImplementedError(
-                "VTKHDF output is only supported for MOAB meshes")
+            raise NotImplementedError("VTKHDF output is only supported for MOAB meshes")
 
         # the self.connectivity contains arrays of length 8 to support hex
         # elements as well, in the case of tetrahedra mesh elements, the
@@ -3347,15 +3282,13 @@ def append_dataset(dset, array):
             else:
                 # No -1 values, append the whole cell
                 trimmed_connectivity.append(cell)
-        trimmed_connectivity = np.array(
-            trimmed_connectivity, dtype="int32").flatten()
+        trimmed_connectivity = np.array(trimmed_connectivity, dtype="int32").flatten()
 
         # MOAB meshes supports tet elements only so we know it has 4 points per cell
         points_per_cell = 4
 
         # offsets are the indices of the first point of each cell in the array of points
-        offsets = np.arange(0, self.n_elements *
-                            points_per_cell + 1, points_per_cell)
+        offsets = np.arange(0, self.n_elements * points_per_cell + 1, points_per_cell)
 
         for name, data in datasets.items():
             if data.shape != self.dimension:
@@ -3366,6 +3299,7 @@ def append_dataset(dset, array):
                 )
 
         with h5py.File(filename, "w") as f:
+
             root = f.create_group("VTKHDF")
             vtk_file_format_version = (2, 1)
             root.attrs["Version"] = vtk_file_format_version
@@ -3377,22 +3311,17 @@ def append_dataset(dset, array):
             )
 
             # create hdf5 file structure
-            root.create_dataset("NumberOfPoints", (0,),
-                                maxshape=(None,), dtype="i8")
+            root.create_dataset("NumberOfPoints", (0,), maxshape=(None,), dtype="i8")
             root.create_dataset("Types", (0,), maxshape=(None,), dtype="uint8")
-            root.create_dataset(
-                "Points", (0, 3), maxshape=(None, 3), dtype="f")
+            root.create_dataset("Points", (0, 3), maxshape=(None, 3), dtype="f")
             root.create_dataset(
                 "NumberOfConnectivityIds", (0,), maxshape=(None,), dtype="i8"
             )
-            root.create_dataset("NumberOfCells", (0,),
-                                maxshape=(None,), dtype="i8")
+            root.create_dataset("NumberOfCells", (0,), maxshape=(None,), dtype="i8")
             root.create_dataset("Offsets", (0,), maxshape=(None,), dtype="i8")
-            root.create_dataset("Connectivity", (0,),
-                                maxshape=(None,), dtype="i8")
+            root.create_dataset("Connectivity", (0,), maxshape=(None,), dtype="i8")
 
-            append_dataset(root["NumberOfPoints"],
-                           np.array([len(self.vertices)]))
+            append_dataset(root["NumberOfPoints"], np.array([len(self.vertices)]))
             append_dataset(root["Points"], self.vertices)
             append_dataset(
                 root["NumberOfConnectivityIds"],
@@ -3403,14 +3332,13 @@ def append_dataset(dset, array):
             append_dataset(root["Offsets"], offsets)
 
             append_dataset(
-                root["Types"], np.full(
-                    self.n_elements, self._VTK_TETRA, dtype="uint8")
+                root["Types"], np.full(self.n_elements, self._VTK_TETRA, dtype="uint8")
             )
 
             cell_data_group = root.create_group("CellData")
 
             for name, data in datasets.items():
-
+                
                 cell_data_group.create_dataset(
                     name, (0,), maxshape=(None,), dtype="float64", chunks=True
                 )
@@ -3548,4 +3476,4 @@ def _read_meshes(elem):
 _HEX_MIDPOINT_CONN += ((2, (0, 0, 0)),
                        (2, (1, 0, 0)),
                        (2, (1, 1, 0)),
-                       (2, (0, 1, 0)))
+                       (2, (0, 1, 0)))
\ No newline at end of file
diff --git a/tests/unit_tests/conftest.py b/tests/unit_tests/conftest.py
index 165fc599120..6041d898201 100644
--- a/tests/unit_tests/conftest.py
+++ b/tests/unit_tests/conftest.py
@@ -1,20 +1,9 @@
-import os
-
 import openmc
 import pytest
 
 from tests.regression_tests import config
 
 
-@pytest.fixture
-def run_in_tmpdir(tmp_path):
-    """Run test in a temporary directory, restoring cwd afterwards."""
-    orig = os.getcwd()
-    os.chdir(tmp_path)
-    yield tmp_path
-    os.chdir(orig)
-
-
 @pytest.fixture(scope='module')
 def mpi_intracomm():
     if config['mpi']:
@@ -59,8 +48,7 @@ def sphere_model():
     model.settings.particles = 100
     model.settings.batches = 10
     model.settings.run_mode = 'fixed source'
-    model.settings.source = openmc.IndependentSource(
-        space=openmc.stats.Point())
+    model.settings.source = openmc.IndependentSource(space=openmc.stats.Point())
     return model
 
 
@@ -84,7 +72,6 @@ def cell_with_lattice():
             [m_inside[0], m_inside[1], m_inside[3], m_outside],
             univ, lattice)
 
-
 @pytest.fixture
 def mixed_lattice_model(uo2, water):
     cyl = openmc.ZCylinder(r=0.4)
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index c4fbc2de656..8a023a319ac 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1,6 +1,8 @@
 from math import pi
 from tempfile import TemporaryDirectory
 from pathlib import Path
+import itertools
+import random
 
 import h5py
 import numpy as np
@@ -53,7 +55,7 @@ def test_regular_mesh_bounding_box():
     mesh.upper_right = [2, 3, 5]
     bb = mesh.bounding_box
     assert isinstance(bb, openmc.BoundingBox)
-    np.testing.assert_array_equal(bb.lower_left, (-2, -3, -5))
+    np.testing.assert_array_equal(bb.lower_left, (-2, -3 ,-5))
     np.testing.assert_array_equal(bb.upper_right, (2, 3, 5))
 
 
@@ -64,7 +66,7 @@ def test_rectilinear_mesh_bounding_box():
     mesh.z_grid = [-100., 0., 100.]
     bb = mesh.bounding_box
     assert isinstance(bb, openmc.BoundingBox)
-    np.testing.assert_array_equal(bb.lower_left, (0., -10., -100.))
+    np.testing.assert_array_equal(bb.lower_left, (0., -10. ,-100.))
     np.testing.assert_array_equal(bb.upper_right, (10., 0., 100.))
 
 
@@ -204,16 +206,14 @@ def test_invalid_cylindrical_mesh_errors():
         openmc.CylindricalMesh(r_grid=[5, 1], phi_grid=[0, pi], z_grid=[0, 10])
 
     with pytest.raises(ValueError):
-        openmc.CylindricalMesh(r_grid=[1, 2, 4, 3], phi_grid=[
-                               0, pi], z_grid=[0, 10])
+        openmc.CylindricalMesh(r_grid=[1, 2, 4, 3], phi_grid=[0, pi], z_grid=[0, 10])
 
     with pytest.raises(ValueError):
         openmc.CylindricalMesh(r_grid=[1], phi_grid=[0, pi], z_grid=[0, 10])
 
     # Test invalid phi_grid values
     with pytest.raises(ValueError):
-        openmc.CylindricalMesh(
-            r_grid=[0, 1, 2], phi_grid=[-1, 3], z_grid=[0, 10])
+        openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[-1, 3], z_grid=[0, 10])
 
     with pytest.raises(ValueError):
         openmc.CylindricalMesh(
@@ -230,12 +230,10 @@ def test_invalid_cylindrical_mesh_errors():
         openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[0, pi], z_grid=[5])
 
     with pytest.raises(ValueError):
-        openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[
-                               0, pi], z_grid=[5, 1])
+        openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[0, pi], z_grid=[5, 1])
 
     with pytest.raises(ValueError):
-        openmc.CylindricalMesh(r_grid=[1, 2, 4, 3], phi_grid=[
-                               0, pi], z_grid=[0, 10, 5])
+        openmc.CylindricalMesh(r_grid=[1, 2, 4, 3], phi_grid=[0, pi], z_grid=[0, 10, 5])
 
 
 def test_centroids():
@@ -254,26 +252,21 @@ def test_centroids():
     np.testing.assert_array_almost_equal(mesh.centroids[0, 0, 0], [6., 7., 8.])
 
     # cylindrical mesh
-    mesh = openmc.CylindricalMesh(
-        r_grid=(0, 10), z_grid=(0, 10), phi_grid=(0, np.pi))
-    np.testing.assert_array_almost_equal(
-        mesh.centroids[0, 0, 0], [0.0, 5.0, 5.0])
+    mesh = openmc.CylindricalMesh(r_grid=(0, 10), z_grid=(0, 10), phi_grid=(0, np.pi))
+    np.testing.assert_array_almost_equal(mesh.centroids[0, 0, 0], [0.0, 5.0, 5.0])
     # ensure that setting an origin is handled correctly
     mesh.origin = (5.0, 0, -10)
-    np.testing.assert_array_almost_equal(
-        mesh.centroids[0, 0, 0], [5.0, 5.0, -5.0])
+    np.testing.assert_array_almost_equal(mesh.centroids[0, 0, 0], [5.0, 5.0, -5.0])
 
     # spherical mesh, single element xyz-positive octant
-    mesh = openmc.SphericalMesh(r_grid=[0, 10], theta_grid=[
-                                0, 0.5*np.pi], phi_grid=[0, np.pi])
+    mesh = openmc.SphericalMesh(r_grid=[0, 10], theta_grid=[0, 0.5*np.pi], phi_grid=[0, np.pi])
     x = 5.*np.cos(0.5*np.pi)*np.sin(0.25*np.pi)
     y = 5.*np.sin(0.5*np.pi)*np.sin(0.25*np.pi)
     z = 5.*np.sin(0.25*np.pi)
     np.testing.assert_array_almost_equal(mesh.centroids[0, 0, 0], [x, y, z])
 
     mesh.origin = (-5.0, -5.0, 5.0)
-    np.testing.assert_array_almost_equal(
-        mesh.centroids[0, 0, 0], [x-5.0, y-5.0, z+5.0])
+    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'))
@@ -296,7 +289,7 @@ def test_mesh_vertices(mesh_type):
         np.testing.assert_equal(mesh.vertices[ijk], exp_i_j_k)
 
         # shift the mesh using the llc
-        shift = np.asarray((3.0, 6.0, 10.0))
+        shift  = np.asarray((3.0, 6.0, 10.0))
         mesh.lower_left += shift
         np.testing.assert_equal(mesh.vertices[ijk], exp_i_j_k+shift)
     elif mesh_type == 'rectilinear':
@@ -313,19 +306,15 @@ def test_mesh_vertices(mesh_type):
         r_grid = np.linspace(0, 5, 10)
         z_grid = np.linspace(-10, 10, 20)
         phi_grid = np.linspace(0, 2*np.pi, 8)
-        mesh = openmc.CylindricalMesh(
-            r_grid=r_grid, z_grid=z_grid, phi_grid=phi_grid)
-        exp_vert = np.asarray(
-            (mesh.r_grid[2], mesh.phi_grid[3], mesh.z_grid[2]))
+        mesh = openmc.CylindricalMesh(r_grid=r_grid, z_grid=z_grid, phi_grid=phi_grid)
+        exp_vert = np.asarray((mesh.r_grid[2], mesh.phi_grid[3], mesh.z_grid[2]))
         np.testing.assert_equal(mesh.vertices_cylindrical[ijk], exp_vert)
     elif mesh_type == 'spherical':
         r_grid = np.linspace(0, 13, 14)
         theta_grid = np.linspace(0, np.pi, 11)
         phi_grid = np.linspace(0, 2*np.pi, 7)
-        mesh = openmc.SphericalMesh(
-            r_grid=r_grid, theta_grid=theta_grid, phi_grid=phi_grid)
-        exp_vert = np.asarray(
-            (mesh.r_grid[2], mesh.theta_grid[3], mesh.phi_grid[2]))
+        mesh = openmc.SphericalMesh(r_grid=r_grid, theta_grid=theta_grid, phi_grid=phi_grid)
+        exp_vert = np.asarray((mesh.r_grid[2], mesh.theta_grid[3], mesh.phi_grid[2]))
         np.testing.assert_equal(mesh.vertices_spherical[ijk], exp_vert)
 
 
@@ -340,11 +329,9 @@ def test_CylindricalMesh_get_indices_at_coords():
     assert mesh.get_indices_at_coords([-2, -2, 9]) == (0, 0, 1)
 
     with pytest.raises(ValueError):
-        assert mesh.get_indices_at_coords(
-            [8, 8, 1])  # resulting r value to large
+        assert mesh.get_indices_at_coords([8, 8, 1])  # resulting r value to large
     with pytest.raises(ValueError):
-        assert mesh.get_indices_at_coords(
-            [-8, -8, 1])  # resulting r value to large
+        assert mesh.get_indices_at_coords([-8, -8, 1])  # resulting r value to large
     with pytest.raises(ValueError):
         assert mesh.get_indices_at_coords([1, 0, -1])  # z value below range
     with pytest.raises(ValueError):
@@ -358,16 +345,11 @@ def test_CylindricalMesh_get_indices_at_coords():
         phi_grid=(0, 0.5 * pi, pi, 1.5 * pi, 1.9 * pi),
         z_grid=(-5, 0, 5, 10),
     )
-    assert mesh.get_indices_at_coords([1, 1, 1]) == (
-        0, 0, 1)  # first angle quadrant
-    assert mesh.get_indices_at_coords([2, 2, 6]) == (
-        0, 0, 2)  # first angle quadrant
-    assert mesh.get_indices_at_coords(
-        [-2, 0.1, -1]) == (0, 1, 0)  # second angle quadrant
-    assert mesh.get_indices_at_coords(
-        [-2, -0.1, -1]) == (0, 2, 0)  # third angle quadrant
-    assert mesh.get_indices_at_coords(
-        [2, -0.9, -1]) == (0, 3, 0)  # forth angle quadrant
+    assert mesh.get_indices_at_coords([1, 1, 1]) == (0, 0, 1)  # first angle quadrant
+    assert mesh.get_indices_at_coords([2, 2, 6]) == (0, 0, 2)  # first angle quadrant
+    assert mesh.get_indices_at_coords([-2, 0.1, -1]) == (0, 1, 0)  # second angle quadrant
+    assert mesh.get_indices_at_coords([-2, -0.1, -1]) == (0, 2, 0)  # third angle quadrant
+    assert mesh.get_indices_at_coords([2, -0.9, -1]) == (0, 3, 0)  # forth angle quadrant
 
     with pytest.raises(ValueError):
         assert mesh.get_indices_at_coords([2, -0.1, 1])  # outside of phi range
@@ -379,16 +361,11 @@ def test_CylindricalMesh_get_indices_at_coords():
         z_grid=(-5, 0, 5, 10),
         origin=(100, 200, 300),
     )
-    assert mesh.get_indices_at_coords([101, 201, 301]) == (
-        0, 0, 1)  # first angle quadrant
-    assert mesh.get_indices_at_coords([102, 202, 306]) == (
-        0, 0, 2)  # first angle quadrant
-    assert mesh.get_indices_at_coords([98, 200.1, 299]) == (
-        0, 1, 0)  # second angle quadrant
-    assert mesh.get_indices_at_coords([98, 199.9, 299]) == (
-        0, 2, 0)  # third angle quadrant
-    assert mesh.get_indices_at_coords([102, 199.1, 299]) == (
-        0, 3, 0)  # forth angle quadrant
+    assert mesh.get_indices_at_coords([101, 201, 301]) == (0, 0, 1)  # first angle quadrant
+    assert mesh.get_indices_at_coords([102, 202, 306]) == (0, 0, 2)  # first angle quadrant
+    assert mesh.get_indices_at_coords([98, 200.1, 299]) == (0, 1, 0)  # second angle quadrant
+    assert mesh.get_indices_at_coords([98, 199.9, 299]) == (0, 2, 0)  # third angle quadrant
+    assert mesh.get_indices_at_coords([102, 199.1, 299]) == (0, 3, 0)  # forth angle quadrant
 
 
 def test_mesh_name_roundtrip(run_in_tmpdir):
@@ -413,8 +390,7 @@ def test_mesh_name_roundtrip(run_in_tmpdir):
 
 
 def test_umesh_roundtrip(run_in_tmpdir, request):
-    umesh = openmc.UnstructuredMesh(
-        request.path.parent / 'test_mesh_tets.e', 'moab')
+    umesh = openmc.UnstructuredMesh(request.path.parent / 'test_mesh_tets.e', 'moab')
     umesh.output = True
 
     # create a tally using this mesh
@@ -446,10 +422,10 @@ def simple_umesh(request):
     geometry = openmc.Geometry([cell1])
 
     umesh = openmc.UnstructuredMesh(
-        filename=request.path.parent.parent
+       filename=request.path.parent.parent
         / "regression_tests/external_moab/test_mesh_tets.h5m",
-        library="moab",
-        mesh_id=1
+       library="moab",
+       mesh_id=1
     )
     # setting ID to make it easier to get the mesh from the statepoint later
     mesh_filter = openmc.MeshFilter(umesh)
@@ -495,8 +471,7 @@ def test_umesh(run_in_tmpdir, simple_umesh, export_type):
     rng = np.random.default_rng()
     ref_data = rng.random(simple_umesh.dimension)
     filename = f"test_mesh{export_type}"
-    simple_umesh.write_data_to_vtk(
-        datasets={"mean": ref_data}, filename=filename)
+    simple_umesh.write_data_to_vtk(datasets={"mean": ref_data}, filename=filename)
 
     assert Path(filename).exists()
 
@@ -514,8 +489,7 @@ def test_umesh(run_in_tmpdir, simple_umesh, export_type):
 
     # attempt to apply a dataset with an improper size to a VTK write
     with pytest.raises(ValueError, match='Cannot apply dataset "mean"'):
-        simple_umesh.write_data_to_vtk(
-            datasets={'mean': ref_data[:-2]}, filename=filename)
+        simple_umesh.write_data_to_vtk(datasets={'mean': ref_data[:-2]}, filename=filename)
 
 
 @pytest.mark.skipif(not openmc.lib._dagmc_enabled(), reason="DAGMC not enabled.")
@@ -534,7 +508,7 @@ def test_write_vtkhdf(request, run_in_tmpdir):
     umesh = openmc.UnstructuredMesh(
         request.path.parent / "test_mesh_dagmc_tets.vtk",
         "moab",
-        mesh_id=1
+        mesh_id = 1
     )
     mesh_filter = openmc.MeshFilter(umesh)
 
@@ -556,15 +530,13 @@ def test_write_vtkhdf(request, run_in_tmpdir):
 
     umesh_from_sp = statepoint.meshes[umesh.id]
 
-    datasets = {
+    datasets={
         "mean": my_tally.mean.flatten(),
         "std_dev": my_tally.std_dev.flatten()
     }
 
-    umesh_from_sp.write_data_to_vtk(
-        datasets=datasets, filename="test_mesh.vtkhdf")
-    umesh_from_sp.write_data_to_vtk(
-        datasets=datasets, filename="test_mesh.vtk")
+    umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtkhdf")
+    umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtk")
 
     with pytest.raises(ValueError, match="Unsupported file extension"):
         # Supported file extensions are vtk or vtkhdf, not hdf5, so this should raise an error
@@ -575,7 +547,7 @@ def test_write_vtkhdf(request, run_in_tmpdir):
     with pytest.raises(ValueError, match="Cannot apply dataset"):
         # The shape of the data should match the shape of the mesh, so this should raise an error
         umesh_from_sp.write_data_to_vtk(
-            datasets={'incorrectly_shaped_data': np.array(([1, 2, 3]))},
+            datasets={'incorrectly_shaped_data': np.array(([1,2,3]))},
             filename="test_mesh_incorrect_shape.vtkhdf",
         )
 
@@ -586,7 +558,6 @@ def test_write_vtkhdf(request, run_in_tmpdir):
     with h5py.File("test_mesh.vtkhdf", "r"):
         ...
 
-
 def test_mesh_get_homogenized_materials():
     """Test the get_homogenized_materials method"""
     # Simple model with 1 cm of Fe56 next to 1 cm of H1
@@ -642,6 +613,7 @@ def test_mesh_get_homogenized_materials():
 
 @pytest.fixture
 def sphere_model():
+    openmc.reset_auto_ids()
     # Model with three materials separated by planes x=0 and z=0
     mats = []
     for i in range(3):
@@ -671,12 +643,9 @@ def test_material_volumes_regular_mesh(sphere_model, n_rays):
     mesh.dimension = (2, 2, 2)
     volumes = mesh.material_volumes(sphere_model, n_rays)
     mats = sphere_model.materials
-    np.testing.assert_almost_equal(
-        volumes[mats[0].id], [0., 0., 0., 0., 0., 1., 0., 1.])
-    np.testing.assert_almost_equal(
-        volumes[mats[1].id], [0., 0., 0., 0., 1., 0., 1., 0.])
-    np.testing.assert_almost_equal(
-        volumes[mats[2].id], [1., 1., 1., 1., 0., 0., 0., 0.])
+    np.testing.assert_almost_equal(volumes[mats[0].id], [0., 0., 0., 0., 0., 1., 0., 1.])
+    np.testing.assert_almost_equal(volumes[mats[1].id], [0., 0., 0., 0., 1., 0., 1., 0.])
+    np.testing.assert_almost_equal(volumes[mats[2].id], [1., 1., 1., 1., 0., 0., 0., 0.])
     assert volumes.by_element(4) == [(mats[1].id, 1.)]
     assert volumes.by_element(0) == [(mats[2].id, 1.)]
 
@@ -768,18 +737,60 @@ def test_mesh_material_volumes_serialize_with_bboxes():
     np.testing.assert_array_equal(second.upper_right, (5.0, 6.0, 7.0))
 
 
+def test_mesh_material_volumes_serialize_with_bboxes():
+    materials = np.array([
+        [1, -1, -2],
+        [-1, -2, -2],
+        [2, 1, -2],
+        [2, -2, -2]
+    ])
+    volumes = np.array([
+        [0.5, 0.5, 0.0],
+        [1.0, 0.0, 0.0],
+        [0.5, 0.5, 0.0],
+        [1.0, 0.0, 0.0]
+    ])
+
+    # (xmin, ymin, zmin, xmax, ymax, zmax)
+    bboxes = np.empty((4, 3, 6))
+    bboxes[..., 0:3] = np.inf
+    bboxes[..., 3:6] = -np.inf
+    bboxes[0, 0] = [-1.0, -2.0, -3.0, 1.0, 2.0, 3.0]    # material 1
+    bboxes[0, 1] = [-5.0, -6.0, -7.0, 5.0, 6.0, 7.0]    # void
+    bboxes[1, 0] = [0.0, 0.0, 0.0, 10.0, 1.0, 2.0]      # void
+    bboxes[2, 0] = [-1.0, -1.0, -1.0, 0.0, 0.0, 0.0]    # material 2
+    bboxes[2, 1] = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]       # material 1
+    bboxes[3, 0] = [-2.0, -2.0, -2.0, 2.0, 2.0, 2.0]    # material 2
+
+    mmv = openmc.MeshMaterialVolumes(materials, volumes, bboxes)
+    with TemporaryDirectory() as tmpdir:
+        path = f'{tmpdir}/volumes_bboxes.npz'
+        mmv.save(path)
+        loaded = openmc.MeshMaterialVolumes.from_npz(path)
+
+    assert loaded.has_bounding_boxes
+    first = loaded.by_element(0, include_bboxes=True)[0][2]
+    assert isinstance(first, openmc.BoundingBox)
+    np.testing.assert_array_equal(first.lower_left, (-1.0, -2.0, -3.0))
+    np.testing.assert_array_equal(first.upper_right, (1.0, 2.0, 3.0))
+
+    second = loaded.by_element(0, include_bboxes=True)[1][2]
+    assert isinstance(second, openmc.BoundingBox)
+    np.testing.assert_array_equal(second.lower_left, (-5.0, -6.0, -7.0))
+    np.testing.assert_array_equal(second.upper_right, (5.0, 6.0, 7.0))
+
+
 def test_mesh_material_volumes_boundary_conditions(sphere_model):
     """Test the material volumes method using a regular mesh
     that overlaps with a vacuum boundary condition."""
 
-    mesh = openmc.SphericalMesh.from_domain(
-        sphere_model.geometry, dimension=(1, 1, 1))
+    mesh = openmc.SphericalMesh.from_domain(sphere_model.geometry, dimension=(1, 1, 1))
     # extend mesh beyond the outer sphere surface to test rays crossing the boundary condition
     mesh.r_grid[-1] += 5.0
 
     # add a new cell to the modelthat occupies the outside of the sphere
     sphere_surfaces = list(filter(lambda s: isinstance(s, openmc.Sphere),
-                                  sphere_model.geometry.get_all_surfaces().values()))
+                            sphere_model.geometry.get_all_surfaces().values()))
     outer_cell = openmc.Cell(region=+sphere_surfaces[0])
     sphere_model.geometry.root_universe.add_cell(outer_cell)
 
@@ -841,6 +852,53 @@ def test_mesh_material_volumes_bounding_boxes():
         np.testing.assert_allclose(bbox.lower_left, cell_bbox.lower_left, atol=tol)
         np.testing.assert_allclose(bbox.upper_right, cell_bbox.upper_right, atol=tol)
 
+def test_mesh_material_volumes_bounding_boxes():
+    # Create a model with 8 spherical cells at known locations with random radii
+    box = openmc.model.RectangularParallelepiped(
+        -10, 10, -10, 10, -10, 10, boundary_type='vacuum')
+
+    mat = openmc.Material()
+    mat.add_nuclide('H1', 1.0)
+
+    sph_cells = []
+    for x, y, z in itertools.product((-5., 5.), repeat=3):
+        mat_i = mat.clone()
+        sph = openmc.Sphere(x, y, z, r=random.uniform(0.5, 1.5))
+        sph_cells.append(openmc.Cell(region=-sph, fill=mat_i))
+    background = openmc.Cell(region=-box & openmc.Intersection([~c.region for c in sph_cells]))
+
+    model = openmc.Model()
+    model.geometry = openmc.Geometry(sph_cells + [background])
+    model.settings.particles = 1000
+    model.settings.batches = 10
+
+    # Create a one-element mesh that encompasses the entire geometry
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = (-10., -10., -10.)
+    mesh.upper_right = (10., 10., 10.)
+    mesh.dimension = (1, 1, 1)
+
+    # Run material volume calculation with bounding boxes
+    n_samples = (400, 400, 400)
+    mmv = mesh.material_volumes(model, n_samples, max_materials=10, bounding_boxes=True)
+    assert mmv.has_bounding_boxes
+
+    # Create a mapping of material ID to bounding box
+    bbox_by_mat = {
+        mat_id: bbox
+        for mat_id, vol, bbox in mmv.by_element(0, include_bboxes=True)
+        if mat_id is not None and vol > 0.0
+    }
+
+    # Match the mesh ray spacing used for the bounding box estimator.
+    tol = 0.5 * mesh.bounding_box.width[0] / n_samples[0]
+    for cell in sph_cells:
+        bbox = bbox_by_mat[cell.fill.id]
+        cell_bbox = cell.bounding_box
+        np.testing.assert_allclose(bbox.lower_left, cell_bbox.lower_left, atol=tol)
+        np.testing.assert_allclose(bbox.upper_right, cell_bbox.upper_right, atol=tol)
+
+
 def test_raytrace_mesh_infinite_loop(run_in_tmpdir):
     # Create a model with one large spherical cell
     sphere = openmc.Sphere(r=100, boundary_type='vacuum')
@@ -868,7 +926,7 @@ def test_raytrace_mesh_infinite_loop(run_in_tmpdir):
     )
     model.settings.run_mode = 'fixed source'
     model.settings.particles = 10
-    model.settings.batches = 1
+    model.settings.batches =  1
 
     # Run the model; this should not cause an infinite loop
     model.run()
@@ -953,10 +1011,6 @@ def test_filter_time_mesh(run_in_tmpdir):
     )
 
 
-# =============================================================================
-# VTKHDF Format Tests for StructuredMeshes
-# =============================================================================
-
 def test_regular_mesh_get_indices_at_coords():
     """Test get_indices_at_coords method for RegularMesh"""
     # Create a 10x10x10 mesh from (0,0,0) to (1,1,1)

From 15b0c5eea9503d001393f723ccd72cd58725c957 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Thu, 2 Apr 2026 11:26:23 +0200
Subject: [PATCH 06/20] removed unrelated changes

---
 openmc/filter.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/openmc/filter.py b/openmc/filter.py
index 9472ca7ef93..87aeb70c3a7 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -37,6 +37,7 @@
 )
 
 
+
 class FilterMeta(ABCMeta):
     """Metaclass for filters that ensures class names are appropriate."""
 
@@ -1860,7 +1861,7 @@ def __init__(self, particles, energies=None, filter_id=None):
     def __repr__(self):
         string = type(self).__name__ + '\n'
         string += '{: <16}=\t{}\n'.format('\tParticles',
-                                          [str(p) for p in self.particles])
+            [str(p) for p in self.particles])
         if self.energies is not None:
             string += '{: <16}=\t{}\n'.format('\tEnergies', self.energies)
         string += '{: <16}=\t{}\n'.format('\tID', self.id)
@@ -2170,7 +2171,7 @@ def paths(self):
         if self._paths is None:
             if not hasattr(self, '_geometry'):
                 raise ValueError(
-                    "Model must be exported before the 'paths' attribute is"
+                    "Model must be exported before the 'paths' attribute is" \
                     "available for a DistribcellFilter.")
 
             # Determine paths for cell instances

From ac916d3c0020d09e47f423ff8768a70cfc205ac6 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Thu, 2 Apr 2026 11:26:33 +0200
Subject: [PATCH 07/20] removed unrelated changes

---
 openmc/mesh.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 50fbc4203fa..c9ebfbfcf8a 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -3338,7 +3338,7 @@ def append_dataset(dset, array):
             cell_data_group = root.create_group("CellData")
 
             for name, data in datasets.items():
-                
+
                 cell_data_group.create_dataset(
                     name, (0,), maxshape=(None,), dtype="float64", chunks=True
                 )
@@ -3476,4 +3476,4 @@ def _read_meshes(elem):
 _HEX_MIDPOINT_CONN += ((2, (0, 0, 0)),
                        (2, (1, 0, 0)),
                        (2, (1, 1, 0)),
-                       (2, (0, 1, 0)))
\ No newline at end of file
+                       (2, (0, 1, 0)))

From 014dcc0e0339ae1b3609e34d60b578aceb4184ab Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Thu, 2 Apr 2026 11:27:29 +0200
Subject: [PATCH 08/20] removed duplication of existing test

---
 tests/unit_tests/test_mesh.py | 92 +----------------------------------
 1 file changed, 1 insertion(+), 91 deletions(-)

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 8a023a319ac..6eb76b5a8f2 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -9,8 +9,6 @@
 from scipy.stats import chi2
 import pytest
 import openmc
-import random
-import itertools
 import openmc.lib
 from openmc.utility_funcs import change_directory
 from uncertainties.unumpy import uarray, nominal_values, std_devs
@@ -558,6 +556,7 @@ def test_write_vtkhdf(request, run_in_tmpdir):
     with h5py.File("test_mesh.vtkhdf", "r"):
         ...
 
+
 def test_mesh_get_homogenized_materials():
     """Test the get_homogenized_materials method"""
     # Simple model with 1 cm of Fe56 next to 1 cm of H1
@@ -737,49 +736,6 @@ def test_mesh_material_volumes_serialize_with_bboxes():
     np.testing.assert_array_equal(second.upper_right, (5.0, 6.0, 7.0))
 
 
-def test_mesh_material_volumes_serialize_with_bboxes():
-    materials = np.array([
-        [1, -1, -2],
-        [-1, -2, -2],
-        [2, 1, -2],
-        [2, -2, -2]
-    ])
-    volumes = np.array([
-        [0.5, 0.5, 0.0],
-        [1.0, 0.0, 0.0],
-        [0.5, 0.5, 0.0],
-        [1.0, 0.0, 0.0]
-    ])
-
-    # (xmin, ymin, zmin, xmax, ymax, zmax)
-    bboxes = np.empty((4, 3, 6))
-    bboxes[..., 0:3] = np.inf
-    bboxes[..., 3:6] = -np.inf
-    bboxes[0, 0] = [-1.0, -2.0, -3.0, 1.0, 2.0, 3.0]    # material 1
-    bboxes[0, 1] = [-5.0, -6.0, -7.0, 5.0, 6.0, 7.0]    # void
-    bboxes[1, 0] = [0.0, 0.0, 0.0, 10.0, 1.0, 2.0]      # void
-    bboxes[2, 0] = [-1.0, -1.0, -1.0, 0.0, 0.0, 0.0]    # material 2
-    bboxes[2, 1] = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]       # material 1
-    bboxes[3, 0] = [-2.0, -2.0, -2.0, 2.0, 2.0, 2.0]    # material 2
-
-    mmv = openmc.MeshMaterialVolumes(materials, volumes, bboxes)
-    with TemporaryDirectory() as tmpdir:
-        path = f'{tmpdir}/volumes_bboxes.npz'
-        mmv.save(path)
-        loaded = openmc.MeshMaterialVolumes.from_npz(path)
-
-    assert loaded.has_bounding_boxes
-    first = loaded.by_element(0, include_bboxes=True)[0][2]
-    assert isinstance(first, openmc.BoundingBox)
-    np.testing.assert_array_equal(first.lower_left, (-1.0, -2.0, -3.0))
-    np.testing.assert_array_equal(first.upper_right, (1.0, 2.0, 3.0))
-
-    second = loaded.by_element(0, include_bboxes=True)[1][2]
-    assert isinstance(second, openmc.BoundingBox)
-    np.testing.assert_array_equal(second.lower_left, (-5.0, -6.0, -7.0))
-    np.testing.assert_array_equal(second.upper_right, (5.0, 6.0, 7.0))
-
-
 def test_mesh_material_volumes_boundary_conditions(sphere_model):
     """Test the material volumes method using a regular mesh
     that overlaps with a vacuum boundary condition."""
@@ -852,52 +808,6 @@ def test_mesh_material_volumes_bounding_boxes():
         np.testing.assert_allclose(bbox.lower_left, cell_bbox.lower_left, atol=tol)
         np.testing.assert_allclose(bbox.upper_right, cell_bbox.upper_right, atol=tol)
 
-def test_mesh_material_volumes_bounding_boxes():
-    # Create a model with 8 spherical cells at known locations with random radii
-    box = openmc.model.RectangularParallelepiped(
-        -10, 10, -10, 10, -10, 10, boundary_type='vacuum')
-
-    mat = openmc.Material()
-    mat.add_nuclide('H1', 1.0)
-
-    sph_cells = []
-    for x, y, z in itertools.product((-5., 5.), repeat=3):
-        mat_i = mat.clone()
-        sph = openmc.Sphere(x, y, z, r=random.uniform(0.5, 1.5))
-        sph_cells.append(openmc.Cell(region=-sph, fill=mat_i))
-    background = openmc.Cell(region=-box & openmc.Intersection([~c.region for c in sph_cells]))
-
-    model = openmc.Model()
-    model.geometry = openmc.Geometry(sph_cells + [background])
-    model.settings.particles = 1000
-    model.settings.batches = 10
-
-    # Create a one-element mesh that encompasses the entire geometry
-    mesh = openmc.RegularMesh()
-    mesh.lower_left = (-10., -10., -10.)
-    mesh.upper_right = (10., 10., 10.)
-    mesh.dimension = (1, 1, 1)
-
-    # Run material volume calculation with bounding boxes
-    n_samples = (400, 400, 400)
-    mmv = mesh.material_volumes(model, n_samples, max_materials=10, bounding_boxes=True)
-    assert mmv.has_bounding_boxes
-
-    # Create a mapping of material ID to bounding box
-    bbox_by_mat = {
-        mat_id: bbox
-        for mat_id, vol, bbox in mmv.by_element(0, include_bboxes=True)
-        if mat_id is not None and vol > 0.0
-    }
-
-    # Match the mesh ray spacing used for the bounding box estimator.
-    tol = 0.5 * mesh.bounding_box.width[0] / n_samples[0]
-    for cell in sph_cells:
-        bbox = bbox_by_mat[cell.fill.id]
-        cell_bbox = cell.bounding_box
-        np.testing.assert_allclose(bbox.lower_left, cell_bbox.lower_left, atol=tol)
-        np.testing.assert_allclose(bbox.upper_right, cell_bbox.upper_right, atol=tol)
-
 
 def test_raytrace_mesh_infinite_loop(run_in_tmpdir):
     # Create a model with one large spherical cell

From 73844f8ef7c7cf305992e397b7a707fa7b5c74f1 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Thu, 2 Apr 2026 11:48:05 +0200
Subject: [PATCH 09/20] Fix VTKHDF Type attribute

---
 openmc/mesh.py | 21 ++++++++++++++++++---
 1 file changed, 18 insertions(+), 3 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index c9ebfbfcf8a..828b58c89b2 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -1831,7 +1831,12 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
         with h5py.File(filename, "w") as f:
             root = f.create_group("VTKHDF")
             root.attrs["Version"] = (2, 1)
-            root.attrs["Type"] = b"StructuredGrid"
+            _type = "StructuredGrid".encode("ascii")
+            root.attrs.create(
+                "Type",
+                _type,
+                dtype=h5py.string_dtype("ascii", len(_type)),
+            )
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
             points = vertices.reshape(-1, 3)
@@ -2331,7 +2336,12 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
         with h5py.File(filename, "w") as f:
             root = f.create_group("VTKHDF")
             root.attrs["Version"] = (2, 1)
-            root.attrs["Type"] = b"StructuredGrid"
+            _type = "StructuredGrid".encode("ascii")
+            root.attrs.create(
+                "Type",
+                _type,
+                dtype=h5py.string_dtype("ascii", len(_type)),
+            )
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
             points = vertices.reshape(-1, 3)
@@ -2767,7 +2777,12 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
         with h5py.File(filename, "w") as f:
             root = f.create_group("VTKHDF")
             root.attrs["Version"] = (2, 1)
-            root.attrs["Type"] = b"StructuredGrid"
+            _type = "StructuredGrid".encode("ascii")
+            root.attrs.create(
+                "Type",
+                _type,
+                dtype=h5py.string_dtype("ascii", len(_type)),
+            )
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
             points = vertices.reshape(-1, 3)

From 799098bc9a856348c9c2294a9a5c0b74d56ebf2d Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Thu, 2 Apr 2026 11:55:39 +0200
Subject: [PATCH 10/20] swapped order or elements

---
 openmc/mesh.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 828b58c89b2..b76aa6a4c70 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -1839,7 +1839,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
             )
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
-            points = vertices.reshape(-1, 3)
+            points = np.swapaxes(vertices, 0, 2).reshape(-1, 3)
             root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
 
             cell_data_group = root.create_group("CellData")
@@ -2344,7 +2344,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
             )
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
-            points = vertices.reshape(-1, 3)
+            points = np.swapaxes(vertices, 0, 2).reshape(-1, 3)
             root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
 
             cell_data_group = root.create_group("CellData")
@@ -2785,7 +2785,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
             )
             root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
 
-            points = vertices.reshape(-1, 3)
+            points = np.swapaxes(vertices, 0, 2).reshape(-1, 3)
             root.create_dataset("Points", data=points.astype(np.float64), dtype="f8")
 
             cell_data_group = root.create_group("CellData")

From 51575354bac6fc9d379d118cbf158e03fa0452e1 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Thu, 2 Apr 2026 12:05:09 +0200
Subject: [PATCH 11/20] reverted line removal

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

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 6eb76b5a8f2..747afd053b2 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -693,6 +693,7 @@ def test_mesh_material_volumes_serialize():
     assert new_volumes.by_element(2) == [(2, 0.5), (1, 0.5)]
     assert new_volumes.by_element(3) == [(2, 1.0)]
 
+
 def test_mesh_material_volumes_serialize_with_bboxes():
     materials = np.array([
         [1, -1, -2],
@@ -762,6 +763,7 @@ def test_mesh_material_volumes_boundary_conditions(sphere_model):
         assert evaluated[0] == expected[0]
         assert evaluated[1] == pytest.approx(expected[1], rel=1e-2)
 
+
 def test_mesh_material_volumes_bounding_boxes():
     # Create a model with 8 spherical cells at known locations with random radii
     box = openmc.model.RectangularParallelepiped(

From 95f3ba2357d70e97ef3491900ed87b1f48d1704c Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Thu, 2 Apr 2026 12:10:34 +0200
Subject: [PATCH 12/20] make use of check_vtk_dataset instead of making new
 method

---
 openmc/mesh.py | 48 ++++--------------------------------------------
 1 file changed, 4 insertions(+), 44 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index b76aa6a4c70..6967c2b0566 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -1542,17 +1542,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
 
             for name, data in datasets.items():
                 data = self._reshape_vtk_dataset(data)
-                if data.ndim > 1 and data.shape != self.dimension:
-                    raise ValueError(
-                        f'Cannot apply multidimensional dataset "{name}" with '
-                        f"shape {data.shape} to mesh {self.id} "
-                        f"with dimensions {self.dimension}"
-                    )
-                if data.size != self.n_elements:
-                    raise ValueError(
-                        f"The size of the dataset '{name}' ({data.size}) should be"
-                        f" equal to the number of mesh cells ({self.n_elements})"
-                    )
+                self._check_vtk_dataset(name, data)
 
                 if volume_normalization:
                     data = data / self.volumes
@@ -1849,17 +1839,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
 
             for name, data in datasets.items():
                 data = self._reshape_vtk_dataset(data)
-                if data.ndim > 1 and data.shape != self.dimension:
-                    raise ValueError(
-                        f'Cannot apply dataset "{name}" with '
-                        f"shape {data.shape} to mesh {self.id} "
-                        f"with dimensions {self.dimension}"
-                    )
-                if data.size != self.n_elements:
-                    raise ValueError(
-                        f"The size of the dataset '{name}' ({data.size}) should be"
-                        f" equal to the number of mesh cells ({self.n_elements})"
-                    )
+                self._check_vtk_dataset(name, data)
 
                 if volume_normalization:
                     data = data / self.volumes
@@ -2354,17 +2334,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
 
             for name, data in datasets.items():
                 data = self._reshape_vtk_dataset(data)
-                if data.ndim > 1 and data.shape != self.dimension:
-                    raise ValueError(
-                        f'Cannot apply dataset "{name}" with '
-                        f"shape {data.shape} to mesh {self.id} "
-                        f"with dimensions {self.dimension}"
-                    )
-                if data.size != self.n_elements:
-                    raise ValueError(
-                        f"The size of the dataset '{name}' ({data.size}) should be"
-                        f" equal to the number of mesh cells ({self.n_elements})"
-                    )
+                self._check_vtk_dataset(name, data)
 
                 if volume_normalization:
                     data = data / self.volumes
@@ -2795,17 +2765,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
 
             for name, data in datasets.items():
                 data = self._reshape_vtk_dataset(data)
-                if data.ndim > 1 and data.shape != self.dimension:
-                    raise ValueError(
-                        f'Cannot apply dataset "{name}" with '
-                        f"shape {data.shape} to mesh {self.id} "
-                        f"with dimensions {self.dimension}"
-                    )
-                if data.size != self.n_elements:
-                    raise ValueError(
-                        f"The size of the dataset '{name}' ({data.size}) should be"
-                        f" equal to the number of mesh cells ({self.n_elements})"
-                    )
+                self._check_vtk_dataset(name, data)
 
                 if volume_normalization:
                     data = data / self.volumes

From 8e616826d780445147cb68c254f7d139531f6ce8 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Thu, 2 Apr 2026 12:13:24 +0200
Subject: [PATCH 13/20] test shape instead of size

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

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 747afd053b2..88dbe9cb934 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1203,8 +1203,8 @@ def test_write_vtkhdf_invalid_data_shape(run_in_tmpdir):
     mesh.upper_right = [1., 1., 1.]
     mesh.dimension = [2, 2, 2]
 
-    # Create data with wrong shape
-    wrong_data = np.ones((3, 3, 3))
+    # Create data with wrong shape but correct total size
+    wrong_data = np.ones((4, 2, 1))
     filename = "test_invalid_shape.vtkhdf"
 
     with pytest.raises(ValueError, match="Cannot apply multidimensional dataset"):

From 776be71825b482e8225d9aea8a81bd195bc10392 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Thu, 2 Apr 2026 12:16:34 +0200
Subject: [PATCH 14/20] added vtkhdf to docstring

---
 openmc/mesh.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 6967c2b0566..c13b0301f9a 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -664,12 +664,16 @@ def write_data_to_vtk(self,
                           datasets: dict | None = None,
                           volume_normalization: bool | None = None,
                           curvilinear: bool = False):
-        """Creates a VTK object of the mesh
+        """Creates a VTK object of the mesh and writes it to a file.
+
+        Supported formats are legacy ASCII ``.vtk`` (requires the ``vtk``
+        package) and ``.vtkhdf`` (requires only ``h5py``).
 
         Parameters
         ----------
-        filename : str
-            Name of the VTK file to write.
+        filename : str or PathLike
+            Name of the VTK file to write. Use a ``.vtkhdf`` extension for
+            the VTKHDF format or ``.vtk`` for the legacy ASCII format.
         datasets : dict
             Dictionary whose keys are the data labels and values are the data
             sets. 1D datasets are expected to be extracted directly from

From ab03026c575f2f8f1c2ebedf084041c850ccbd02 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Thu, 2 Apr 2026 12:16:50 +0200
Subject: [PATCH 15/20] tidy up space

---
 tests/unit_tests/test_mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 88dbe9cb934..2e59e9f336d 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1287,7 +1287,7 @@ def test_write_ascii_vtk_unchanged(run_in_tmpdir):
     arr = reader.GetOutput().GetCellData().GetArray("data")
     read_data = np.array([arr.GetTuple1(i) for i in range(ref_data.size)])
     np.testing.assert_almost_equal(read_data, ref_data.ravel())
-    
+
 def test_rectilinear_mesh_get_indices_at_coords():
     """Test get_indices_at_coords method for RectilinearMesh"""
     # Create a 3x2x2 rectilinear mesh with non-uniform spacing

From e24d852f640dcb169a61f6f5480a019e226063ae Mon Sep 17 00:00:00 2001
From: Makarand More <40858007+Jarvis2001@users.noreply.github.com>
Date: Thu, 2 Apr 2026 13:50:38 +0000
Subject: [PATCH 16/20] fixed structured mesh output, expanded rectilinear
 coverage and added some new tests

---
 openmc/mesh.py                |  21 +++++--
 tests/unit_tests/test_mesh.py | 115 ++++++++++++++++++++++++++++++++++
 2 files changed, 130 insertions(+), 6 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index c13b0301f9a..c5f642c6526 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -952,12 +952,21 @@ def _check_vtk_dataset(self, label: str, dataset: np.ndarray):
         if dataset.ndim == 1:
             return
 
-        if dataset.shape != self.dimension:
-            raise ValueError(
-                f'Cannot apply multidimensional dataset "{label}" with '
-                f"shape {dataset.shape} to mesh {self.id} "
-                f"with dimensions {self.dimension}"
-            )
+        if dataset.shape == self.dimension:
+            return
+
+        # allow datasets in lower dimension when the missing dimensions are
+        # singleton (e.g., 2D data on a 3D mesh with one depth layer)
+        if dataset.ndim < len(self.dimension):
+            padded_shape = tuple(dataset.shape) + (1,) * (len(self.dimension) - dataset.ndim)
+            if padded_shape == tuple(self.dimension):
+                return
+
+        raise ValueError(
+            f'Cannot apply multidimensional dataset "{label}" with '
+            f"shape {dataset.shape} to mesh {self.id} "
+            f"with dimensions {self.dimension}"
+        )
 
 
 class HasBoundingBox(Protocol):
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 2e59e9f336d..7e62dc4a574 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1113,6 +1113,49 @@ def test_write_vtkhdf_spherical_mesh(run_in_tmpdir):
         assert "density" in root["CellData"]
 
 
+@pytest.mark.parametrize(
+    "mesh",
+    [
+        openmc.RectilinearMesh(),
+        openmc.CylindricalMesh(
+            r_grid=[0.0, 1.0, 2.0],
+            z_grid=[0.0, 1.0, 2.0],
+            phi_grid=[0.0, np.pi / 2, np.pi],
+        ),
+        openmc.SphericalMesh(
+            r_grid=[0.0, 1.0, 2.0],
+            theta_grid=[0.0, np.pi / 2, np.pi],
+            phi_grid=[0.0, np.pi / 2, np.pi],
+        ),
+    ],
+    ids=["rectilinear", "cylindrical", "spherical"],
+)
+def test_write_vtkhdf_structuredgrid_points_order_and_type(run_in_tmpdir, mesh):
+    """Test VTKHDF _write_vtk_hdf5 point ordering and Type attribute."""
+    if isinstance(mesh, openmc.RectilinearMesh):
+        mesh.x_grid = [0.0, 1.0, 2.0]
+        mesh.y_grid = [0.0, 1.0, 2.0]
+        mesh.z_grid = [0.0, 1.0, 2.0]
+    data = np.arange(mesh.n_elements).reshape(mesh.dimension)
+    mesh.write_data_to_vtk(
+        datasets={"values": data},
+        filename="test_structured.vtkhdf",
+        volume_normalization=False,
+    )
+
+    with h5py.File("test_structured.vtkhdf", "r") as f:
+        root = f["VTKHDF"]
+
+        points = root["Points"][()]
+        expected_points = np.swapaxes(mesh.vertices, 0, 2).reshape(-1, 3)
+        np.testing.assert_allclose(points, expected_points)
+
+        type_attr = root.attrs["Type"]
+        if isinstance(type_attr, bytes):
+            type_attr = type_attr.decode("ascii")
+        assert type_attr == "StructuredGrid"
+
+
 def test_write_vtkhdf_volume_normalization(run_in_tmpdir):
     """Test volume normalization in VTKHDF format."""
     mesh = openmc.RegularMesh()
@@ -1153,6 +1196,78 @@ def test_write_vtkhdf_volume_normalization(run_in_tmpdir):
     np.testing.assert_allclose(unnormalized_data, ref_data.ravel())
 
 
+def test_write_vtkhdf_default_volume_normalization_for_vtkhdf(run_in_tmpdir):
+    """Ensure VTKHDF default keeps data un-normalized by volume."""
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = [0.0, 0.0, 0.0]
+    mesh.upper_right = [10.0, 10.0, 10.0]
+    mesh.dimension = [2, 2, 2]
+
+    data = np.arange(mesh.n_elements).reshape(mesh.dimension)
+    mesh.write_data_to_vtk(
+        datasets={"flux": data},
+        filename="test_default_norm.vtkhdf",
+    )
+
+    with h5py.File("test_default_norm.vtkhdf", "r") as f:
+        saved = f["VTKHDF"]["CellData"]["flux"][()]
+
+    np.testing.assert_allclose(saved, data.T.ravel())
+
+
+@pytest.mark.parametrize(
+    "x_grid,y_grid,z_grid,data,expected", [
+        (
+            [0.0, 1.0, 2.0],
+            [0.0, 1.0],
+            [0.0, 1.0],
+            np.arange(2),
+            np.arange(2),
+        ),
+        (
+            [0.0, 1.0, 2.0],
+            [0.0, 1.0, 2.0],
+            [0.0, 1.0],
+            np.arange(4).reshape(2,2),
+            np.arange(4),
+        ),
+        (
+            [0.0, 1.0, 2.0],
+            [0.0, 1.0, 2.0],
+            [0.0, 1.0, 2.0],
+            np.arange(8).reshape(2,2,2),
+            np.arange(8),
+        ),
+    ],
+    ids=["1d", "2d", "3d"],
+)
+def test_write_vtkhdf_rectilinear_multi_dimensional_data(run_in_tmpdir, x_grid, y_grid, z_grid, data, expected):
+    """Test RectilinearMesh VTKHDF accepts 1D/2D/3D data with singleton padding."""
+    mesh = openmc.RectilinearMesh()
+    mesh.x_grid = x_grid
+    mesh.y_grid = y_grid
+    mesh.z_grid = z_grid
+
+    mesh.write_data_to_vtk(
+        datasets={"flux": data},
+        filename="test_rectilinear_dims.vtkhdf",
+        volume_normalization=False,
+    )
+
+    with h5py.File("test_rectilinear_dims.vtkhdf", "r") as f:
+        saved = f["VTKHDF"]["CellData"]["flux"][()]
+
+    expected_written = np.swapaxes(data, 0, data.ndim-1).ravel() if data.ndim > 1 else data.ravel()
+    if data.ndim == 1:
+        expected_written = expected
+    elif data.ndim == 2:
+        expected_written = data.T.ravel()
+    elif data.ndim == 3:
+        expected_written = data.T.ravel()
+
+    np.testing.assert_allclose(saved, expected_written)
+
+
 def test_write_vtkhdf_multiple_datasets(run_in_tmpdir):
     """Test writing multiple datasets to VTKHDF format."""
     mesh = openmc.RegularMesh()

From 3a92e49e041c7b4a594a7129b6b19524fc9de463 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Sun, 5 Apr 2026 19:53:12 +0200
Subject: [PATCH 17/20] volume_normalization back to True by default

---
 openmc/mesh.py | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 695da41a5a7..fecb8850f02 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -662,7 +662,7 @@ def num_mesh_cells(self):
     def write_data_to_vtk(self,
                           filename: PathLike,
                           datasets: dict | None = None,
-                          volume_normalization: bool | None = None,
+                          volume_normalization: bool = True,
                           curvilinear: bool = False):
         """Creates a VTK object of the mesh and writes it to a file.
 
@@ -682,10 +682,9 @@ def write_data_to_vtk(self,
             with structured indexing in "C" ordering. See the "expand_dims" flag
             of :meth:`~openmc.Tally.get_reshaped_data` on reshaping tally data when using
             :class:`~openmc.MeshFilter`'s.
-        volume_normalization : bool or None, optional
+        volume_normalization : bool, optional
             Whether or not to normalize the data by the volume of the mesh
-            elements. When None (default), the format-appropriate default is
-            used: True for legacy ASCII .vtk files, False for .vtkhdf files.
+            elements. Defaults to True.
         curvilinear : bool
             Whether or not to write curvilinear elements. Only applies to
             ``SphericalMesh`` and ``CylindricalMesh``.
@@ -725,9 +724,6 @@ def write_data_to_vtk(self,
             write_impl(filename, datasets, volume_normalization)
             return None
 
-        # legacy ASCII .vtk default is True
-        norm = volume_normalization if volume_normalization is not None else True
-
         # vtk is an optional dependency only needed for the legacy ASCII path
         import vtk
         from vtk.util import numpy_support as nps
@@ -760,7 +756,7 @@ def write_data_to_vtk(self,
                     dataset = dataset.T.ravel() if needs_transpose else dataset.ravel()
                 datasets_out.append(dataset)
 
-                if norm:
+                if volume_normalization:
                     if needs_transpose:
                         dataset /= self.volumes.T.ravel()
                     else:

From ab1e59d3fef49a84027e333f791a882ec9096b77 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Sun, 5 Apr 2026 20:01:09 +0200
Subject: [PATCH 18/20] update return in doc string

---
 openmc/mesh.py                | 22 +++++++++++-----------
 tests/unit_tests/test_mesh.py |  9 +--------
 2 files changed, 12 insertions(+), 19 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index fecb8850f02..517b4e0c2f0 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -696,8 +696,9 @@ def write_data_to_vtk(self,
 
         Returns
         -------
-        vtk.StructuredGrid or vtk.UnstructuredGrid
-            a VTK grid object representing the mesh
+        vtk.StructuredGrid or vtk.UnstructuredGrid or None
+            A VTK grid object representing the mesh for the legacy ASCII
+            format, or None for the VTKHDF format.
 
         Examples
         --------
@@ -1590,13 +1591,12 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
             coords_1d.append(np.array([0.0]))
 
         # np.meshgrid with indexing='ij' → axis 0 = x, axis 1 = y, axis 2 = z
-        xx, yy, zz = np.meshgrid(*coords_1d, indexing='ij')
-        # Flatten in Fortran (x-fastest) order for VTK point ordering
-        points = np.column_stack([
-            xx.ravel(order='F'),
-            yy.ravel(order='F'),
-            zz.ravel(order='F'),
-        ]).astype(np.float64)   # shape (n_points, 3)
+        vertices = np.stack(
+            np.meshgrid(*coords_1d, indexing='ij'), axis=-1
+        )
+        # Swap first and last spatial axes then flatten, matching the
+        # approach used by RectilinearMesh/CylindricalMesh/SphericalMesh.
+        points = np.swapaxes(vertices, 0, 2).reshape(-1, 3).astype(np.float64)
 
         with h5py.File(filename, "w") as f:
             root = f.create_group("VTKHDF")
@@ -1923,7 +1923,7 @@ def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
                 cell_data_group.create_dataset(
                     name, data=flat_data, dtype="float64", chunks=True
                 )
-    
+
     @classmethod
     def from_bounding_box(
         cls,
@@ -3320,7 +3320,7 @@ def _write_data_to_vtk_hdf5_format(
         volume_normalization: bool = True,
     ):
         """Write UnstructuredMesh as VTK-HDF5 UnstructuredGrid format.
-        
+
         Supports linear tetrahedra and linear hexahedra elements.
         """
         def append_dataset(dset, array):
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 7e62dc4a574..3eddbbd34b2 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1257,14 +1257,7 @@ def test_write_vtkhdf_rectilinear_multi_dimensional_data(run_in_tmpdir, x_grid,
     with h5py.File("test_rectilinear_dims.vtkhdf", "r") as f:
         saved = f["VTKHDF"]["CellData"]["flux"][()]
 
-    expected_written = np.swapaxes(data, 0, data.ndim-1).ravel() if data.ndim > 1 else data.ravel()
-    if data.ndim == 1:
-        expected_written = expected
-    elif data.ndim == 2:
-        expected_written = data.T.ravel()
-    elif data.ndim == 3:
-        expected_written = data.T.ravel()
-
+    expected_written = data.T.ravel() if data.ndim > 1 else data.ravel()
     np.testing.assert_allclose(saved, expected_written)
 
 

From dda304dd4e7c6ac5109b13d2cd404fdccdaebbe5 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Sun, 5 Apr 2026 20:06:34 +0200
Subject: [PATCH 19/20] update tests for defaut true vol norm

---
 tests/unit_tests/test_mesh.py | 23 ++++++++++++++++-------
 1 file changed, 16 insertions(+), 7 deletions(-)

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 3eddbbd34b2..8d9568218b0 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1116,6 +1116,10 @@ def test_write_vtkhdf_spherical_mesh(run_in_tmpdir):
 @pytest.mark.parametrize(
     "mesh",
     [
+        openmc.RegularMesh.from_domain(
+            openmc.BoundingBox([0.0, 0.0, 0.0], [2.0, 2.0, 2.0]),
+            dimension=[2, 2, 2],
+        ),
         openmc.RectilinearMesh(),
         openmc.CylindricalMesh(
             r_grid=[0.0, 1.0, 2.0],
@@ -1128,7 +1132,7 @@ def test_write_vtkhdf_spherical_mesh(run_in_tmpdir):
             phi_grid=[0.0, np.pi / 2, np.pi],
         ),
     ],
-    ids=["rectilinear", "cylindrical", "spherical"],
+    ids=["regular", "rectilinear", "cylindrical", "spherical"],
 )
 def test_write_vtkhdf_structuredgrid_points_order_and_type(run_in_tmpdir, mesh):
     """Test VTKHDF _write_vtk_hdf5 point ordering and Type attribute."""
@@ -1197,13 +1201,13 @@ def test_write_vtkhdf_volume_normalization(run_in_tmpdir):
 
 
 def test_write_vtkhdf_default_volume_normalization_for_vtkhdf(run_in_tmpdir):
-    """Ensure VTKHDF default keeps data un-normalized by volume."""
+    """Ensure VTKHDF default normalizes data by volume (same as legacy)."""
     mesh = openmc.RegularMesh()
     mesh.lower_left = [0.0, 0.0, 0.0]
     mesh.upper_right = [10.0, 10.0, 10.0]
     mesh.dimension = [2, 2, 2]
 
-    data = np.arange(mesh.n_elements).reshape(mesh.dimension)
+    data = np.ones(mesh.dimension) * 100.0
     mesh.write_data_to_vtk(
         datasets={"flux": data},
         filename="test_default_norm.vtkhdf",
@@ -1212,7 +1216,9 @@ def test_write_vtkhdf_default_volume_normalization_for_vtkhdf(run_in_tmpdir):
     with h5py.File("test_default_norm.vtkhdf", "r") as f:
         saved = f["VTKHDF"]["CellData"]["flux"][()]
 
-    np.testing.assert_allclose(saved, data.T.ravel())
+    # Each cell volume is 5*5*5 = 125, default normalization divides by it
+    cell_volume = 125.0
+    np.testing.assert_allclose(saved, (data / cell_volume).T.ravel())
 
 
 @pytest.mark.parametrize(
@@ -1277,7 +1283,8 @@ def test_write_vtkhdf_multiple_datasets(run_in_tmpdir):
     filename = "test_multiple_datasets.vtkhdf"
     mesh.write_data_to_vtk(
         datasets={"flux": data1, "power": data2, "heating": data3},
-        filename=filename
+        filename=filename,
+        volume_normalization=False
     )
 
     assert Path(filename).exists()
@@ -1332,7 +1339,8 @@ def test_write_vtkhdf_1d_mesh(run_in_tmpdir):
     rng = np.random.default_rng(42)
     ref_data = rng.random(mesh.dimension)
     filename = "test_1d_mesh.vtkhdf"
-    mesh.write_data_to_vtk(datasets={"data": ref_data}, filename=filename)
+    mesh.write_data_to_vtk(datasets={"data": ref_data}, filename=filename,
+                           volume_normalization=False)
 
     assert Path(filename).exists()
 
@@ -1354,7 +1362,8 @@ def test_write_vtkhdf_2d_mesh(run_in_tmpdir):
     rng = np.random.default_rng(42)
     ref_data = rng.random(mesh.dimension)
     filename = "test_2d_mesh.vtkhdf"
-    mesh.write_data_to_vtk(datasets={"data": ref_data}, filename=filename)
+    mesh.write_data_to_vtk(datasets={"data": ref_data}, filename=filename,
+                           volume_normalization=False)
 
     assert Path(filename).exists()
 

From e3e4684f5a5fbe6c2cb2edfba0f9c8b7302b7157 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Sun, 5 Apr 2026 20:15:17 +0200
Subject: [PATCH 20/20] added doc strings

---
 openmc/mesh.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 517b4e0c2f0..3d8da0a39ec 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -1572,7 +1572,7 @@ def get_indices_at_coords(self, coords: Sequence[float]) -> tuple:
         indices = np.floor((coords_array - lower_left) / spacing).astype(int)
         return tuple(int(i) for i in indices[:ndim])
 
-    def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
+    def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None, volume_normalization: bool) -> None:
         """Write RegularMesh as VTKHDF StructuredGrid format."""
         dims = self.dimension
         ndim = len(dims)
@@ -1880,7 +1880,7 @@ def get_indices_at_coords(self, coords: Sequence[float]) -> tuple[int, int, int]
 
         return tuple(indices)
 
-    def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
+    def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None, volume_normalization: bool) -> None:
         """Write RectilinearMesh as VTK-HDF5 StructuredGrid format.
 
         Note: vtkRectilinearGrid is not part of the VTKHDF spec yet, so
@@ -2402,7 +2402,7 @@ def _convert_to_cartesian(arr, origin: Sequence[float]):
         arr[..., 2] += origin[2]
         return arr
 
-    def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
+    def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None, volume_normalization: bool) -> None:
         """Write CylindricalMesh as VTK-HDF5 StructuredGrid format."""
         nr, nphi, nz = self.dimension
         vertex_dims = [nr + 1, nphi + 1, nz + 1]
@@ -2819,7 +2819,7 @@ def get_indices_at_coords(self, coords: Sequence[float]) -> tuple:
             "get_indices_at_coords is not yet implemented for SphericalMesh"
         )
 
-    def _write_vtk_hdf5(self, filename, datasets, volume_normalization):
+    def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None, volume_normalization: bool) -> None:
         """Write SphericalMesh as VTK-HDF5 StructuredGrid format."""
         nr, ntheta, nphi = self.dimension
         vertex_dims = [nr + 1, ntheta + 1, nphi + 1]