Some reverting and small fix in filter.py

Author: Jarvis2001

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

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)

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()