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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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/32] 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]

From 3d29e6eb2a4000be475cb0d55c58528bc7d38c69 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Tue, 14 Apr 2026 15:45:37 +0200
Subject: [PATCH 21/32] revert ravel code changes in vtk branch

---
 openmc/mesh.py                | 20 +++++++-------------
 tests/unit_tests/test_mesh.py |  2 +-
 2 files changed, 8 insertions(+), 14 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 3d8da0a39ec..0da943755c1 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -742,26 +742,20 @@ 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.T.ravel() if needs_transpose else dataset.ravel()
+                    dataset = dataset.T.ravel()
                 datasets_out.append(dataset)
 
                 if volume_normalization:
-                    if needs_transpose:
-                        dataset /= self.volumes.T.ravel()
-                    else:
-                        dataset /= self.volumes.ravel()
+                    dataset /= self.volumes.T.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 8d9568218b0..1990b24c027 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1403,7 +1403,7 @@ 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())
+    np.testing.assert_almost_equal(read_data, ref_data.T.ravel())
 
 def test_rectilinear_mesh_get_indices_at_coords():
     """Test get_indices_at_coords method for RectilinearMesh"""

From 17c9b2ab51945dc37fd2facd3857e8da2a6a2fba Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Tue, 14 Apr 2026 15:50:40 +0200
Subject: [PATCH 22/32] same VTKHDF name in all docstrings

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

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 0da943755c1..a8f2ac92681 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -1875,7 +1875,7 @@ def get_indices_at_coords(self, coords: Sequence[float]) -> tuple[int, int, int]
         return tuple(indices)
 
     def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None, volume_normalization: bool) -> None:
-        """Write RectilinearMesh as VTK-HDF5 StructuredGrid format.
+        """Write RectilinearMesh as VTKHDF StructuredGrid format.
 
         Note: vtkRectilinearGrid is not part of the VTKHDF spec yet, so
         StructuredGrid with explicit point coordinates is used instead.
@@ -2397,7 +2397,7 @@ def _convert_to_cartesian(arr, origin: Sequence[float]):
         return arr
 
     def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None, volume_normalization: bool) -> None:
-        """Write CylindricalMesh as VTK-HDF5 StructuredGrid format."""
+        """Write CylindricalMesh as VTKHDF StructuredGrid format."""
         nr, nphi, nz = self.dimension
         vertex_dims = [nr + 1, nphi + 1, nz + 1]
 
@@ -2814,7 +2814,7 @@ def get_indices_at_coords(self, coords: Sequence[float]) -> tuple:
         )
 
     def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None, volume_normalization: bool) -> None:
-        """Write SphericalMesh as VTK-HDF5 StructuredGrid format."""
+        """Write SphericalMesh as VTKHDF StructuredGrid format."""
         nr, ntheta, nphi = self.dimension
         vertex_dims = [nr + 1, ntheta + 1, nphi + 1]
 

From e046963ca9c0975d0f3efe2708bf5eede874e4ac Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Tue, 14 Apr 2026 15:53:20 +0200
Subject: [PATCH 23/32] test avoids uniform data values

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

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 1990b24c027..e69c491c247 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1167,8 +1167,9 @@ def test_write_vtkhdf_volume_normalization(run_in_tmpdir):
     mesh.upper_right = [10., 10., 10.]
     mesh.dimension = [2, 2, 2]
 
-    # Create data with known values
-    ref_data = np.ones(mesh.dimension) * 100.0
+    # Use non-uniform data so that .T.ravel() and .ravel() differ,
+    # catching any ordering mistakes in the implementation.
+    ref_data = np.arange(8, dtype=float).reshape(mesh.dimension)
     filename_with_norm = "test_normalized.vtkhdf"
     filename_without_norm = "test_unnormalized.vtkhdf"
 
@@ -1195,9 +1196,9 @@ def test_write_vtkhdf_volume_normalization(run_in_tmpdir):
 
     # Volume for each cell is 5 x 5 x 5 = 125
     cell_volume = 125.0
-    expected_normalized = ref_data.ravel() / cell_volume
+    expected_normalized = ref_data.T.ravel() / cell_volume
     np.testing.assert_allclose(normalized_data, expected_normalized)
-    np.testing.assert_allclose(unnormalized_data, ref_data.ravel())
+    np.testing.assert_allclose(unnormalized_data, ref_data.T.ravel())
 
 
 def test_write_vtkhdf_default_volume_normalization_for_vtkhdf(run_in_tmpdir):

From f5e7d583696d20c112578d1d66f8ea01b3989f0e Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Tue, 14 Apr 2026 15:59:50 +0200
Subject: [PATCH 24/32] all classes have method so no guard needed

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

diff --git a/openmc/mesh.py b/openmc/mesh.py
index a8f2ac92681..953e3258140 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -718,11 +718,7 @@ def write_data_to_vtk(self,
            >>> 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)
+            self._write_vtk_hdf5(filename, datasets, volume_normalization)
             return None
 
         # vtk is an optional dependency only needed for the legacy ASCII path

From 774c057fcfe6e58162368204476509c9fae93fc6 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Wed, 6 May 2026 21:38:10 +0200
Subject: [PATCH 25/32] dim check now passing

---
 tests/unit_tests/test_mesh.py | 27 ++++++++++++++-------------
 1 file changed, 14 insertions(+), 13 deletions(-)

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index e69c491c247..cd94d7d75b2 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1337,20 +1337,20 @@ def test_write_vtkhdf_1d_mesh(run_in_tmpdir):
     mesh.upper_right = [10.]
     mesh.dimension = [5]
 
-    rng = np.random.default_rng(42)
-    ref_data = rng.random(mesh.dimension)
+    ref_data = np.arange(5, dtype=float)
     filename = "test_1d_mesh.vtkhdf"
-    mesh.write_data_to_vtk(datasets={"data": ref_data}, filename=filename,
-                           volume_normalization=False)
+    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
+        assert len(dims) == 1
+        saved = root["CellData"]["data"][()]
+
+    # Default volume_normalization=True; cell width = 2.0 → volume = 2.0
+    np.testing.assert_allclose(saved, ref_data / 2.0)
 
 
 def test_write_vtkhdf_2d_mesh(run_in_tmpdir):
@@ -1360,20 +1360,21 @@ def test_write_vtkhdf_2d_mesh(run_in_tmpdir):
     mesh.upper_right = [10., 10.]
     mesh.dimension = [5, 3]
 
-    rng = np.random.default_rng(42)
-    ref_data = rng.random(mesh.dimension)
+    ref_data = np.arange(15, dtype=float).reshape(mesh.dimension)
     filename = "test_2d_mesh.vtkhdf"
-    mesh.write_data_to_vtk(datasets={"data": ref_data}, filename=filename,
-                           volume_normalization=False)
+    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
+        saved = root["CellData"]["data"][()]
+
+    # Default volume_normalization=True; cell volume = 2.0 * (10/3)
+    cell_volume = 2.0 * (10.0 / 3.0)
+    np.testing.assert_allclose(saved, (ref_data / cell_volume).T.ravel())
 
 
 def test_write_ascii_vtk_unchanged(run_in_tmpdir):

From 6cc9f54e01aab0ce8d10148fec1cf0eda7982094 Mon Sep 17 00:00:00 2001
From: Makarand More <40858007+Jarvis2001@users.noreply.github.com>
Date: Tue, 26 May 2026 23:30:49 +0000
Subject: [PATCH 26/32] Resolved merge conflicts and fixed UnstructuredMesh
 VTKHDF Offsets array length

The VTKHDF UnstructuredGrid writer produced Offsets arrays with only
n_cells entries. vtkHDFReader requires n_cells + 1 entries per the
VTKHDF spec, so it loaded zero cells despite valid connectivity data.
---
 openmc/mesh.py | 251 ++++++++++++++++++++++++++++++-------------------
 1 file changed, 152 insertions(+), 99 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 0adac2ba196..e467da4a531 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):
@@ -289,7 +291,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()
 
@@ -397,7 +400,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()
@@ -609,7 +613,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
@@ -663,7 +668,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 | None = None,
                           curvilinear: bool = False):
         """Creates a VTK object of the mesh and writes it to a file.
 
@@ -683,7 +688,7 @@ 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. Defaults to True.
         curvilinear : bool
@@ -739,24 +744,31 @@ 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.T.ravel()
+                    dataset = dataset.T.ravel() if needs_transpose else dataset.ravel()
                 datasets_out.append(dataset)
 
                 if volume_normalization:
-                    dataset /= self.volumes.T.ravel()
+                    if needs_transpose:
+                        dataset = dataset / self.volumes.T.ravel()
+                    else:
+                        dataset = dataset / self.volumes.ravel()
 
                 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))
@@ -777,7 +789,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])
@@ -804,12 +817,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()
 
@@ -861,7 +875,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])
 
@@ -869,14 +884,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())
@@ -930,32 +947,29 @@ def _check_vtk_dataset(self, label: str, dataset: np.ndarray):
         """
         cv.check_type('data label', label, str)
 
+        # If this is a multidimensional dataset, validate its shape first so
+        # that the user receives the more specific error message.
+        if dataset.ndim > 1:
+            if dataset.shape == self.dimension:
+                return
+
+            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}"
+            )
+
         if dataset.size != self.n_elements:
             raise ValueError(
                 f"The size of the dataset '{label}' ({dataset.size}) should be"
                 f" equal to the number of mesh cells ({self.n_elements})"
             )
 
-        # accept a flat array as-is, assuming it is in the correct order
-        if dataset.ndim == 1:
-            return
-
-        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}"
-        )
-
     @classmethod
     def from_domain(
         cls,
@@ -1120,7 +1134,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):
@@ -1144,7 +1159,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):
@@ -1154,7 +1170,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
@@ -1228,10 +1244,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
 
@@ -1246,7 +1265,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
 
@@ -1755,17 +1775,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])
@@ -2128,17 +2151,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])
@@ -2146,9 +2172,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
@@ -2164,7 +2190,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(
@@ -2213,9 +2240,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'][()]
@@ -2333,10 +2360,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,
         )
 
@@ -2361,7 +2388,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
@@ -2372,7 +2400,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
@@ -2604,17 +2633,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])
@@ -2625,9 +2657,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
         )
@@ -2743,10 +2775,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
@@ -2770,7 +2802,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
@@ -2781,7 +2814,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
@@ -2790,7 +2824,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))
@@ -3223,7 +3256,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"
@@ -3326,6 +3360,8 @@ def append_dataset(dset, array):
         if self.library != "moab":
             raise NotImplementedError("VTKHDF output is only supported for MOAB meshes")
 
+        datasets = {} if datasets is None else datasets
+
         # This writer supports linear tetrahedra and linear hexahedra elements
         conn_list = []           # flattened connectivity ids
         cell_sizes = []          # number of points per cell
@@ -3356,12 +3392,15 @@ def append_dataset(dset, array):
             )
 
         connectivity = np.asarray(conn_list, dtype=np.int64)
+        n_cells = len(vtk_types)
 
-        # Offsets must be length (numCells + 1) with a leading 0 and
-        # cumulative end-indices thereafter, per VTK's layout
-        cell_sizes_arr = np.asarray(cell_sizes, dtype=np.int64)
-        offsets = np.zeros(cell_sizes_arr.size + 1, dtype=np.int64)
-        np.cumsum(cell_sizes_arr, out=offsets[1:])
+        # Offsets are the index of the first point of each cell in the
+        # flattened connectivity array.
+        if cell_sizes:
+            # VTKHDF requires n_cells + 1 offset entries per partition
+            offsets = np.concatenate(([0], np.cumsum(cell_sizes, dtype=np.int64)))
+        else:
+            offsets = np.array([0], dtype=np.int64)
 
         for name, data in datasets.items():
             if data.shape != self.dimension:
@@ -3372,7 +3411,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
@@ -3383,24 +3421,39 @@ def append_dataset(dset, array):
                 dtype=h5py.string_dtype("ascii", len(ascii_type)),
             )
 
-            # Create HDF5 file structure compliant with VTKHDF UnstructuredGrid
-            n_points = int(len(self.vertices))
-            n_cells = int(len(cell_sizes))
-            n_conn_ids = int(len(connectivity))
-
-            root.create_dataset("NumberOfPoints", data=(n_points,), dtype="i8")
-            root.create_dataset("NumberOfCells", data=(n_cells,), dtype="i8")
-            root.create_dataset("NumberOfConnectivityIds", data=(n_conn_ids,), dtype="i8")
-            root.create_dataset("Points", data=self.vertices.astype(np.float64, copy=False), dtype="f8")
-            root.create_dataset("Types", data=np.asarray(vtk_types, dtype=np.uint8), dtype="uint8")
-            root.create_dataset("Offsets", data=offsets.astype("i8"), dtype="i8")
-            root.create_dataset("Connectivity", data=connectivity.astype("i8"), dtype="i8")
+            # create hdf5 file structure
+            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(
+                "NumberOfConnectivityIds", (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")
+
+            append_dataset(root["NumberOfPoints"],
+                           np.array([len(self.vertices)]))
+            append_dataset(root["Points"], self.vertices)
+            append_dataset(
+                root["NumberOfConnectivityIds"],
+                np.array([len(connectivity)]),
+            )
+            append_dataset(root["Connectivity"], connectivity)
+            append_dataset(root["NumberOfCells"], np.array([n_cells]))
+            append_dataset(root["Offsets"], offsets)
+            append_dataset(root["Types"], np.asarray(vtk_types, dtype="uint8"))
 
             cell_data_group = root.create_group("CellData")
 
             for name, data in datasets.items():
+                data = np.asarray(data, dtype="float64")
                 if volume_normalization:
-                    data /= self.volumes
+                    data = data / self.volumes
                 cell_data_group.create_dataset(
                     name, data=data, dtype="float64", chunks=True
                 )

From cc947e8903bdd856e6b2ae820bcb2349aa24970e Mon Sep 17 00:00:00 2001
From: Makarand More <40858007+Jarvis2001@users.noreply.github.com>
Date: Wed, 27 May 2026 08:18:08 +0000
Subject: [PATCH 27/32] test should pass now

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

diff --git a/openmc/mesh.py b/openmc/mesh.py
index e467da4a531..1cbbf7ad847 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -723,6 +723,9 @@ def write_data_to_vtk(self,
            >>> heating = tally.get_reshaped_data(expand_dims=True)
            >>> mesh.write_data_to_vtk({'heating': heating})
         """
+        if volume_normalization is None:
+            volume_normalization = True
+
         if Path(filename).suffix == ".vtkhdf":
             self._write_vtk_hdf5(filename, datasets, volume_normalization)
             return None
@@ -744,26 +747,20 @@ 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.T.ravel() if needs_transpose else dataset.ravel()
+                    dataset = dataset.T.ravel()
                 datasets_out.append(dataset)
 
                 if volume_normalization:
-                    if needs_transpose:
-                        dataset = dataset / self.volumes.T.ravel()
-                    else:
-                        dataset = dataset / self.volumes.ravel()
+                    dataset /= self.volumes.T.ravel()
 
                 dataset_array = vtk.vtkDoubleArray()
                 dataset_array.SetName(label)

From 7a0982d48c947d1f5f0350980d64eb9086bb339f Mon Sep 17 00:00:00 2001
From: Makarand More <40858007+Jarvis2001@users.noreply.github.com>
Date: Wed, 27 May 2026 08:18:08 +0000
Subject: [PATCH 28/32] test should pass now

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

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 1cbbf7ad847..b40105aa732 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):
@@ -291,8 +289,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()
 
@@ -400,8 +397,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()
@@ -613,8 +609,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
@@ -764,8 +759,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))
@@ -786,8 +780,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])
@@ -814,13 +807,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()
 
@@ -872,8 +864,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])
 
@@ -881,16 +872,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())
@@ -1131,8 +1120,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):
@@ -1156,8 +1144,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):
@@ -1241,13 +1228,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
 
@@ -1262,8 +1246,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
 
@@ -1772,20 +1755,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])
@@ -2148,27 +2128,24 @@ 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])
             string += fmt.format('\tZ Max:', '=\t', self._z_grid[-1])
         return string
 
-    def get_indices_at_coords(
+    def get_indices_at_coords(        
         self,
         coords: Sequence[float]
     ) -> tuple[int, int, int]:
@@ -2187,8 +2164,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(
@@ -2237,9 +2213,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'][()]
@@ -2357,10 +2333,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,
         )
 
@@ -2385,8 +2361,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
@@ -2397,8 +2372,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
@@ -2630,20 +2604,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])
@@ -2654,9 +2625,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
         )
@@ -2772,10 +2743,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
@@ -2799,8 +2770,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
@@ -2811,8 +2781,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
@@ -3253,8 +3222,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"
@@ -3354,9 +3322,6 @@ def append_dataset(dset, array):
             dset.resize(origLen + array.shape[0], axis=0)
             dset[origLen:] = array
 
-        if self.library != "moab":
-            raise NotImplementedError("VTKHDF output is only supported for MOAB meshes")
-
         datasets = {} if datasets is None else datasets
 
         # This writer supports linear tetrahedra and linear hexahedra elements
@@ -3450,7 +3415,7 @@ def append_dataset(dset, array):
             for name, data in datasets.items():
                 data = np.asarray(data, dtype="float64")
                 if volume_normalization:
-                    data = data / self.volumes
+                    data /= self.volumes
                 cell_data_group.create_dataset(
                     name, data=data, dtype="float64", chunks=True
                 )

From dd2481809da543d98377999cddfb61427bf77bcb Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Mon, 3 Aug 2026 14:39:09 +0200
Subject: [PATCH 29/32] Revert volume_normalization default to True

The None default was reintroduced in 6cc9f54e0 when resolving merge
conflicts, undoing the earlier revert in 3a92e49e0. Restore the plain
bool default so the public signature and behaviour match develop.
---
 openmc/mesh.py | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index e4fb3710e3b..3745cf7cd23 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -661,7 +661,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.
 
@@ -681,9 +681,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. Defaults to True.
+            elements.
         curvilinear : bool
             Whether or not to write curvilinear elements. Only applies to
             ``SphericalMesh`` and ``CylindricalMesh``.
@@ -716,9 +716,6 @@ def write_data_to_vtk(self,
            >>> heating = tally.get_reshaped_data(expand_dims=True)
            >>> mesh.write_data_to_vtk({'heating': heating})
         """
-        if volume_normalization is None:
-            volume_normalization = True
-
         if Path(filename).suffix == ".vtkhdf":
             self._write_vtk_hdf5(filename, datasets, volume_normalization)
             return None

From 05ecdbc9cea25483584240d699fd1ff6e14195b9 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Mon, 3 Aug 2026 14:52:10 +0200
Subject: [PATCH 30/32] Write structured meshes as VTKHDF types that VTK
 supports

VTKHDF defines no StructuredGrid or RectilinearGrid dataset type, so the
files written for the four structured mesh types were rejected by every
VTK reader with "Unknown data set type: StructuredGrid" and loaded in
ParaView as an empty object.

A regular mesh is a uniform grid, so it is now written as ImageData,
described by an extent, origin and spacing. This needs no point
coordinates at all and is roughly four times smaller than before. The
extent is padded to six entries so 1D and 2D meshes are valid.

Rectilinear, cylindrical and spherical meshes are written as an
UnstructuredGrid of linear hexahedra, matching what UnstructuredMesh
already does. Points come from the existing vertices property, which
reproduces the coordinates the four separate implementations were
computing by hand, so those collapse into one shared method.

Element ordering is unchanged, with the first mesh index varying fastest.
Volume normalization now divides a flat dataset by volumes in the same
order as the data rather than relying on broadcasting.

The new tests read every file back with vtkHDFReader and check both the
element data and the element geometry against the mesh vertices. All
eight fail on the previous output and pass on this.
---
 openmc/mesh.py                | 336 +++++++++++++++++-----------------
 tests/unit_tests/test_mesh.py | 231 ++++++++++++++---------
 2 files changed, 305 insertions(+), 262 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index daf75431929..1be270c8098 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -22,6 +22,24 @@
 from .utility_funcs import input_path
 
 
+# VTKHDF format version and element type used for structured mesh output
+_VTKHDF_VERSION = (2, 1)
+_VTK_HEXAHEDRON = 12
+
+
+def _write_vtkhdf_type(group, name: str):
+    """Write the VTKHDF dataset type as a fixed length ASCII string
+
+    A variable length string is read back as ``str`` rather than ``bytes``,
+    which does not match what VTK readers expect.
+
+    """
+    ascii_name = name.encode("ascii")
+    group.attrs.create(
+        "Type", ascii_name, dtype=h5py.string_dtype("ascii", len(ascii_name))
+    )
+
+
 class MeshMaterialVolumes(Mapping):
     """Results from a material volume in mesh calculation.
 
@@ -958,6 +976,112 @@ def _check_vtk_dataset(self, label: str, dataset: np.ndarray):
                 f" equal to the number of mesh cells ({self.n_elements})"
             )
 
+    def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None,
+                        volume_normalization: bool) -> None:
+        """Write the mesh and data in the VTKHDF UnstructuredGrid format
+
+        VTKHDF defines no StructuredGrid or RectilinearGrid dataset type, so
+        meshes with non-uniform or curvilinear elements are written as an
+        unstructured grid of linear hexahedra. :class:`RegularMesh` overrides
+        this with the more compact ImageData format.
+
+        Parameters
+        ----------
+        filename : str or PathLike
+            Name of the VTKHDF file to write.
+        datasets : dict or None
+            Dictionary whose keys are the data labels and values are the data
+            sets.
+        volume_normalization : bool
+            Whether or not to normalize the data by the volume of the mesh
+            elements.
+
+        """
+        n_i, n_j, n_k = self.dimension
+        vertex_shape = (n_i + 1, n_j + 1, n_k + 1)
+        points = np.asarray(self.vertices).reshape(-1, 3)
+
+        # Index of the (i, j, k) corner of each element, ordered with i varying
+        # fastest to match the cell ordering used elsewhere for VTK output
+        i, j, k = np.meshgrid(np.arange(n_i), np.arange(n_j), np.arange(n_k),
+                              indexing='ij')
+        corner = np.ravel_multi_index((i, j, k), vertex_shape)
+        corner = corner.transpose(2, 1, 0).ravel()
+
+        # Offset from the (i, j, k) corner to each vertex of a VTK hexahedron
+        stride_i = (n_j + 1) * (n_k + 1)
+        stride_j = n_k + 1
+        offsets = np.array([
+            0,
+            stride_i,
+            stride_i + stride_j,
+            stride_j,
+            1,
+            stride_i + 1,
+            stride_i + stride_j + 1,
+            stride_j + 1,
+        ])
+        connectivity = (corner[:, None] + offsets[None, :]).ravel()
+
+        n_cells = corner.size
+        with h5py.File(filename, "w") as f:
+            root = f.create_group("VTKHDF")
+            root.attrs["Version"] = _VTKHDF_VERSION
+            _write_vtkhdf_type(root, "UnstructuredGrid")
+            root.create_dataset("NumberOfPoints", data=[len(points)], dtype="i8")
+            root.create_dataset("Points", data=points, dtype="f8")
+            root.create_dataset("NumberOfCells", data=[n_cells], dtype="i8")
+            root.create_dataset("NumberOfConnectivityIds",
+                                data=[connectivity.size], dtype="i8")
+            root.create_dataset("Connectivity", data=connectivity, dtype="i8")
+            root.create_dataset(
+                "Offsets",
+                data=np.arange(0, connectivity.size + 1, _N_HEX_VERTICES),
+                dtype="i8",
+            )
+            root.create_dataset("Types",
+                                data=np.full(n_cells, _VTK_HEXAHEDRON),
+                                dtype="uint8")
+
+            cell_data_group = root.create_group("CellData")
+            for label, dataset in self._vtk_hdf5_datasets(datasets,
+                                                          volume_normalization):
+                cell_data_group.create_dataset(label, data=dataset,
+                                               dtype="float64", chunks=True)
+
+    def _vtk_hdf5_datasets(self, datasets: dict | None,
+                           volume_normalization: bool):
+        """Validate and flatten datasets for VTKHDF output
+
+        Yields
+        ------
+        tuple of (str, numpy.ndarray)
+            Data label and the corresponding flat array of element values,
+            ordered with the first mesh index varying fastest.
+
+        """
+        if not datasets:
+            return
+
+        for label, dataset in datasets.items():
+            dataset = self._reshape_vtk_dataset(dataset)
+            self._check_vtk_dataset(label, dataset)
+
+            if volume_normalization:
+                # a flat dataset is already ordered with the first index
+                # varying fastest, so the volumes have to match that ordering
+                if dataset.ndim > 1:
+                    dataset = dataset / self.volumes
+                else:
+                    dataset = dataset / self.volumes.T.ravel()
+
+            if dataset.ndim > 1:
+                dataset = dataset.T.ravel()
+            else:
+                dataset = dataset.ravel()
+
+            yield label, dataset
+
     @classmethod
     def from_domain(
         cls,
@@ -1568,59 +1692,54 @@ 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: PathLike, datasets: dict | None, volume_normalization: bool) -> None:
-        """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
-        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)
+    def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None,
+                        volume_normalization: bool) -> None:
+        """Write the mesh and data in the VTKHDF ImageData format
+
+        A regular mesh is a uniform grid, so it is fully described by an origin
+        and a spacing and no point coordinates need to be stored. Meshes with
+        fewer than three dimensions are padded to three, as required by the
+        VTKHDF ImageData specification.
+
+        Parameters
+        ----------
+        filename : str or PathLike
+            Name of the VTKHDF file to write.
+        datasets : dict or None
+            Dictionary whose keys are the data labels and values are the data
+            sets.
+        volume_normalization : bool
+            Whether or not to normalize the data by the volume of the mesh
+            elements.
+
+        """
+        ndim = len(self.dimension)
+        shape = list(self.dimension) + [1] * (3 - ndim)
+        origin = list(self.lower_left) + [0.0] * (3 - ndim)
+        spacing = list(self.width) + [1.0] * (3 - ndim)
 
         with h5py.File(filename, "w") as f:
             root = f.create_group("VTKHDF")
-            root.attrs["Version"] = (2, 1)
-            _type = "StructuredGrid".encode("ascii")
+            root.attrs["Version"] = _VTKHDF_VERSION
+            _write_vtkhdf_type(root, "ImageData")
             root.attrs.create(
-                "Type",
-                _type,
-                dtype=h5py.string_dtype("ascii", len(_type)),
+                "WholeExtent",
+                np.array([0, shape[0], 0, shape[1], 0, shape[2]], dtype="i8"),
             )
-            root.create_dataset("Dimensions", data=vertex_dims, dtype="i8")
-            root.create_dataset("Points", data=points, dtype="f8")
+            root.attrs.create("Origin", np.array(origin, dtype="f8"))
+            root.attrs.create("Spacing", np.array(spacing, dtype="f8"))
+            root.attrs.create("Direction", np.eye(3).ravel())
 
             cell_data_group = root.create_group("CellData")
-
-            if not datasets:
-                return
-
-            for name, data in datasets.items():
-                data = self._reshape_vtk_dataset(data)
-                self._check_vtk_dataset(name, data)
-
-                if volume_normalization:
-                    data = data / self.volumes
-
-                flat_data = data.T.ravel() if data.ndim > 1 else data.ravel()
+            for label, dataset in self._vtk_hdf5_datasets(datasets,
+                                                          volume_normalization):
+                # ImageData element data is stored with the last index varying
+                # fastest, i.e. with shape (nz, ny, nx)
                 cell_data_group.create_dataset(
-                    name, data=flat_data, dtype="float64", chunks=True
+                    label,
+                    data=dataset.reshape(shape[2], shape[1], shape[0]),
+                    dtype="float64",
+                    chunks=True,
                 )
 
 
@@ -1878,49 +1997,6 @@ def get_indices_at_coords(self, coords: Sequence[float]) -> tuple[int, int, int]
 
         return tuple(indices)
 
-    def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None, volume_normalization: bool) -> None:
-        """Write RectilinearMesh as VTKHDF 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)
-            _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 = 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")
-
-            if not datasets:
-                return
-
-            for name, data in datasets.items():
-                data = self._reshape_vtk_dataset(data)
-                self._check_vtk_dataset(name, data)
-
-                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
-                )
 
     @classmethod
     def from_bounding_box(
@@ -2402,47 +2478,6 @@ def _convert_to_cartesian(arr, origin: Sequence[float]):
         arr[..., 2] += origin[2]
         return arr
 
-    def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None, volume_normalization: bool) -> None:
-        """Write CylindricalMesh as VTKHDF 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)
-            _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 = 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")
-
-            if not datasets:
-                return
-
-            for name, data in datasets.items():
-                data = self._reshape_vtk_dataset(data)
-                self._check_vtk_dataset(name, data)
-
-                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):
@@ -2892,49 +2927,6 @@ def get_indices_at_coords(
 
         return (r_index, theta_index, phi_index)
 
-    def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None, volume_normalization: bool) -> None:
-        """Write SphericalMesh as VTKHDF 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)
-            _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 = 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")
-
-            if not datasets:
-                return
-
-            for name, data in datasets.items():
-                data = self._reshape_vtk_dataset(data)
-                self._check_vtk_dataset(name, data)
-
-                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):
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 4efb12d429c..d7dea66b4ee 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1032,14 +1032,94 @@ def test_regular_mesh_get_indices_at_coords():
     assert len(result_1d) == 1
     assert result_1d == (5,)
 
+def _read_vtkhdf(filename):
+    """Read a VTKHDF file with VTK, skipping the test if VTK is unavailable"""
+    vtk = pytest.importorskip("vtk")
+    reader = vtk.vtkHDFReader()
+    reader.SetFileName(str(filename))
+    reader.Update()
+    return reader.GetOutput()
+
+
+VTKHDF_MESHES = {
+    "regular": lambda: openmc.RegularMesh.from_domain(
+        openmc.BoundingBox([0.0, 0.0, 0.0], [2.0, 3.0, 4.0]),
+        dimension=[2, 3, 4],
+    ),
+    "rectilinear": lambda: openmc.RectilinearMesh.from_domain(
+        openmc.BoundingBox([0.0, 0.0, 0.0], [2.0, 3.0, 4.0]),
+        dimension=[2, 3, 4],
+    ),
+    "cylindrical": lambda: openmc.CylindricalMesh(
+        r_grid=[0.0, 1.0, 2.0],
+        phi_grid=np.linspace(0.0, 2 * pi, 4),
+        z_grid=[0.0, 1.0, 2.0, 3.0, 4.0],
+    ),
+    "spherical": lambda: openmc.SphericalMesh(
+        r_grid=[0.0, 1.0, 2.0],
+        theta_grid=np.linspace(0.0, pi, 4),
+        phi_grid=np.linspace(0.0, 2 * pi, 6),
+    ),
+}
+
+
+@pytest.mark.parametrize("kind", list(VTKHDF_MESHES))
+def test_write_vtkhdf_readable_by_vtk(run_in_tmpdir, kind):
+    """VTKHDF output must be readable by VTK, not merely well formed HDF5.
+
+    VTKHDF has no StructuredGrid or RectilinearGrid dataset type, so a file
+    claiming to be one is silently rejected by every VTK reader. Assert that
+    the file loads and that the element data survives the round trip.
+    """
+    mesh = VTKHDF_MESHES[kind]()
+    # distinct values per element so any reordering is detected
+    data = np.arange(float(mesh.n_elements)).reshape(mesh.dimension)
+
+    filename = f"test_readable_{kind}.vtkhdf"
+    mesh.write_data_to_vtk(datasets={"values": data}, filename=filename,
+                           volume_normalization=False)
+
+    grid = _read_vtkhdf(filename)
+    assert grid.GetNumberOfCells() == mesh.n_elements
+    assert grid.GetNumberOfPoints() == np.prod([d + 1 for d in mesh.dimension])
+
+    array = grid.GetCellData().GetArray("values")
+    assert array is not None
+    read_data = np.array([array.GetTuple1(i) for i in range(mesh.n_elements)])
+    np.testing.assert_allclose(read_data, data.T.ravel())
+
+
+@pytest.mark.parametrize("kind", list(VTKHDF_MESHES))
+def test_write_vtkhdf_element_geometry(run_in_tmpdir, kind):
+    """Each element read back by VTK must span the expected mesh vertices"""
+    mesh = VTKHDF_MESHES[kind]()
+    mesh.write_data_to_vtk(datasets=None, filename=f"test_geom_{kind}.vtkhdf")
+
+    grid = _read_vtkhdf(f"test_geom_{kind}.vtkhdf")
+    n_i, n_j, _ = mesh.dimension
+    vertices = np.asarray(mesh.vertices)
+
+    for cell_id in range(grid.GetNumberOfCells()):
+        i = cell_id % n_i
+        j = (cell_id // n_i) % n_j
+        k = cell_id // (n_i * n_j)
+        cell = grid.GetCell(cell_id)
+        found = np.array([cell.GetPoints().GetPoint(n) for n in range(8)])
+        expected = vertices[i:i+2, j:j+2, k:k+2].reshape(-1, 3)
+        np.testing.assert_allclose(
+            np.array(sorted(map(tuple, np.round(found, 8)))),
+            np.array(sorted(map(tuple, np.round(expected, 8)))),
+            atol=1e-8,
+        )
+
+
 def test_write_vtkhdf_regular_mesh(run_in_tmpdir):
-    """Test writing a regular mesh to VTKHDF format."""
+    """A regular mesh is uniform, so it is written as VTKHDF ImageData."""
     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"
@@ -1047,34 +1127,33 @@ def test_write_vtkhdf_regular_mesh(run_in_tmpdir):
 
     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 root.attrs["Type"] == b"ImageData"
         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
+        # ImageData is described by an extent, origin and spacing rather than
+        # explicit point coordinates
+        assert "Points" not in root
+        np.testing.assert_array_equal(root.attrs["WholeExtent"],
+                                      [0, 2, 0, 3, 0, 4])
+        np.testing.assert_allclose(root.attrs["Origin"], [-5., -5., -5.])
+        np.testing.assert_allclose(root.attrs["Spacing"], [5., 10/3, 2.5])
 
-        # Check data shape
-        cell_data = root["CellData"]["data"][()]
-        assert cell_data.shape == (ref_data.size,)
+        # element data is stored with the last index varying fastest
+        assert root["CellData"]["data"].shape == (4, 3, 2)
 
 
 def test_write_vtkhdf_rectilinear_mesh(run_in_tmpdir):
-    """Test writing a rectilinear mesh to VTKHDF format."""
+    """A rectilinear mesh is written as a VTKHDF UnstructuredGrid."""
     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"
@@ -1082,16 +1161,18 @@ def test_write_vtkhdf_rectilinear_mesh(run_in_tmpdir):
 
     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"UnstructuredGrid"
         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
+        assert root["CellData"]["flux"].size == ref_data.size
+        assert root["NumberOfCells"][0] == mesh.n_elements
+        # all elements are linear hexahedra
+        assert np.all(root["Types"][()] == 12)
+        assert root["Offsets"].size == mesh.n_elements + 1
 
 
 def test_write_vtkhdf_cylindrical_mesh(run_in_tmpdir):
@@ -1103,7 +1184,6 @@ def test_write_vtkhdf_cylindrical_mesh(run_in_tmpdir):
         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"
@@ -1111,17 +1191,16 @@ def test_write_vtkhdf_cylindrical_mesh(run_in_tmpdir):
 
     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"UnstructuredGrid"
         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)
+        # one point per mesh vertex, r x phi x z
+        assert root["NumberOfPoints"][0] == (3+1) * (4+1) * (2+1)
+        assert root["NumberOfCells"][0] == 3 * 4 * 2
 
 
 def test_write_vtkhdf_spherical_mesh(run_in_tmpdir):
@@ -1133,7 +1212,6 @@ def test_write_vtkhdf_spherical_mesh(run_in_tmpdir):
         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"
@@ -1141,60 +1219,34 @@ def test_write_vtkhdf_spherical_mesh(run_in_tmpdir):
 
     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"UnstructuredGrid"
         assert "Points" in root
         assert "CellData" in root
         assert "density" in root["CellData"]
 
 
 @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],
-            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=["regular", "rectilinear", "cylindrical", "spherical"],
+    "kind", ["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)
+def test_write_vtkhdf_points_match_vertices(run_in_tmpdir, kind):
+    """Points written for an unstructured grid must be the mesh vertices."""
+    mesh = VTKHDF_MESHES[kind]()
+    data = np.arange(float(mesh.n_elements)).reshape(mesh.dimension)
     mesh.write_data_to_vtk(
         datasets={"values": data},
-        filename="test_structured.vtkhdf",
+        filename="test_unstructured.vtkhdf",
         volume_normalization=False,
     )
 
-    with h5py.File("test_structured.vtkhdf", "r") as f:
+    with h5py.File("test_unstructured.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"
+        np.testing.assert_allclose(
+            root["Points"][()], np.asarray(mesh.vertices).reshape(-1, 3)
+        )
+        assert root.attrs["Type"] == b"UnstructuredGrid"
 
 
 def test_write_vtkhdf_volume_normalization(run_in_tmpdir):
@@ -1210,21 +1262,17 @@ def test_write_vtkhdf_volume_normalization(run_in_tmpdir):
     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"][()]
 
@@ -1233,9 +1281,11 @@ def test_write_vtkhdf_volume_normalization(run_in_tmpdir):
 
     # Volume for each cell is 5 x 5 x 5 = 125
     cell_volume = 125.0
-    expected_normalized = ref_data.T.ravel() / cell_volume
-    np.testing.assert_allclose(normalized_data, expected_normalized)
-    np.testing.assert_allclose(unnormalized_data, ref_data.T.ravel())
+    shape = (2, 2, 2)
+    np.testing.assert_allclose(
+        normalized_data, (ref_data.T.ravel() / cell_volume).reshape(shape))
+    np.testing.assert_allclose(
+        unnormalized_data, ref_data.T.ravel().reshape(shape))
 
 
 def test_write_vtkhdf_default_volume_normalization_for_vtkhdf(run_in_tmpdir):
@@ -1256,7 +1306,8 @@ def test_write_vtkhdf_default_volume_normalization_for_vtkhdf(run_in_tmpdir):
 
     # 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())
+    np.testing.assert_allclose(
+        saved, (data / cell_volume).T.ravel().reshape(2, 2, 2))
 
 
 @pytest.mark.parametrize(
@@ -1312,7 +1363,6 @@ def test_write_vtkhdf_multiple_datasets(run_in_tmpdir):
     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)
@@ -1327,26 +1377,19 @@ def test_write_vtkhdf_multiple_datasets(run_in_tmpdir):
 
     assert Path(filename).exists()
 
-    # Verify all datasets are present
+    shape = (2, 2, 2)
     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()
-        )
+            root["CellData"]["flux"][()], data1.T.ravel().reshape(shape))
         np.testing.assert_allclose(
-            root["CellData"]["power"][()],
-            data2.T.ravel()
-        )
+            root["CellData"]["power"][()], data2.T.ravel().reshape(shape))
         np.testing.assert_allclose(
-            root["CellData"]["heating"][()],
-            data3.T.ravel()
-        )
+            root["CellData"]["heating"][()], data3.T.ravel().reshape(shape))
 
 
 def test_write_vtkhdf_invalid_data_shape(run_in_tmpdir):
@@ -1382,12 +1425,16 @@ def test_write_vtkhdf_1d_mesh(run_in_tmpdir):
 
     with h5py.File(filename, "r") as f:
         root = f["VTKHDF"]
-        dims = root["Dimensions"][()]
-        assert len(dims) == 1
+        # a VTKHDF extent always has six entries, flat axes are padded to one
+        np.testing.assert_array_equal(root.attrs["WholeExtent"],
+                                      [0, 5, 0, 1, 0, 1])
         saved = root["CellData"]["data"][()]
 
-    # Default volume_normalization=True; cell width = 2.0 → volume = 2.0
-    np.testing.assert_allclose(saved, ref_data / 2.0)
+    # Default volume_normalization=True; cell width = 2.0 -> volume = 2.0
+    np.testing.assert_allclose(saved, (ref_data / 2.0).reshape(1, 1, 5))
+
+    grid = _read_vtkhdf(filename)
+    assert grid.GetNumberOfCells() == 5
 
 
 def test_write_vtkhdf_2d_mesh(run_in_tmpdir):
@@ -1405,13 +1452,17 @@ def test_write_vtkhdf_2d_mesh(run_in_tmpdir):
 
     with h5py.File(filename, "r") as f:
         root = f["VTKHDF"]
-        dims = root["Dimensions"][()]
-        assert len(dims) == 2
+        np.testing.assert_array_equal(root.attrs["WholeExtent"],
+                                      [0, 5, 0, 3, 0, 1])
         saved = root["CellData"]["data"][()]
 
     # Default volume_normalization=True; cell volume = 2.0 * (10/3)
     cell_volume = 2.0 * (10.0 / 3.0)
-    np.testing.assert_allclose(saved, (ref_data / cell_volume).T.ravel())
+    np.testing.assert_allclose(
+        saved, (ref_data / cell_volume).T.ravel().reshape(1, 3, 5))
+
+    grid = _read_vtkhdf(filename)
+    assert grid.GetNumberOfCells() == 15
 
 
 def test_write_ascii_vtk_unchanged(run_in_tmpdir):

From 2c7c65054491f5b4fc2e78bbdc8a94bcd3df3f9a Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Mon, 3 Aug 2026 14:58:33 +0200
Subject: [PATCH 31/32] Inline the VTK hexahedron cell type

Used in one place, so a named constant does not earn its keep.
---
 openmc/mesh.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 1be270c8098..7afa2240fe3 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -22,9 +22,8 @@
 from .utility_funcs import input_path
 
 
-# VTKHDF format version and element type used for structured mesh output
+# VTKHDF format version used for structured mesh output
 _VTKHDF_VERSION = (2, 1)
-_VTK_HEXAHEDRON = 12
 
 
 def _write_vtkhdf_type(group, name: str):
@@ -1039,8 +1038,8 @@ def _write_vtk_hdf5(self, filename: PathLike, datasets: dict | None,
                 data=np.arange(0, connectivity.size + 1, _N_HEX_VERTICES),
                 dtype="i8",
             )
-            root.create_dataset("Types",
-                                data=np.full(n_cells, _VTK_HEXAHEDRON),
+            # 12 is the VTK cell type for a linear hexahedron
+            root.create_dataset("Types", data=np.full(n_cells, 12),
                                 dtype="uint8")
 
             cell_data_group = root.create_group("CellData")

From aee4de8306e700e60845c303d3b1fb88a7ee7c62 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Sat, 10 Oct 2026 16:27:40 +0200
Subject: [PATCH 32/32] Reshape padded datasets before VTKHDF volume
 normalization

A multidimensional dataset with trailing axes of size one omitted, such
as (2, 2) for a mesh of dimension (2, 2, 1), broadcast against the
volumes and wrote more cell values than there are cells.
---
 openmc/mesh.py                |  4 ++++
 tests/unit_tests/test_mesh.py | 13 ++++++++++---
 2 files changed, 14 insertions(+), 3 deletions(-)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 2b6469adb5c..da37b090d80 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -1068,6 +1068,10 @@ def _vtk_hdf5_datasets(self, datasets: dict | None,
         for label, dataset in datasets.items():
             dataset = self._reshape_vtk_dataset(dataset)
             self._check_vtk_dataset(label, dataset)
+            # restore any trailing axes of size one so that the dataset
+            # aligns with the volumes instead of broadcasting against them
+            if dataset.ndim > 1:
+                dataset = dataset.reshape(self.dimension)
 
             if volume_normalization:
                 # a flat dataset is already ordered with the first index
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 2cb36551ee3..b1e6b2b2474 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1506,8 +1506,15 @@ def test_write_vtkhdf_default_volume_normalization_for_vtkhdf(run_in_tmpdir):
     ],
     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."""
+@pytest.mark.parametrize("volume_normalization", [False, True])
+def test_write_vtkhdf_rectilinear_multi_dimensional_data(
+    run_in_tmpdir, x_grid, y_grid, z_grid, data, expected, volume_normalization
+):
+    """Test RectilinearMesh VTKHDF accepts 1D/2D/3D data with singleton padding.
+
+    The grids have unit spacing, so every element volume is one and the
+    written values are the same with and without volume normalization.
+    """
     mesh = openmc.RectilinearMesh()
     mesh.x_grid = x_grid
     mesh.y_grid = y_grid
@@ -1516,7 +1523,7 @@ def test_write_vtkhdf_rectilinear_multi_dimensional_data(run_in_tmpdir, x_grid,
     mesh.write_data_to_vtk(
         datasets={"flux": data},
         filename="test_rectilinear_dims.vtkhdf",
-        volume_normalization=False,
+        volume_normalization=volume_normalization,
     )
 
     with h5py.File("test_rectilinear_dims.vtkhdf", "r") as f: