From d1a4815de1f6159a0098a7a5867f13539af8b059 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sat, 6 Jul 2024 17:23:52 -0500
Subject: [PATCH 01/23] Initial version of mesh material volumes based on
 raytracing

---
 include/openmc/mesh.h |  60 ++++++++++++++++++
 openmc/lib/mesh.py    |  20 ++++++
 src/mesh.cpp          | 137 ++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 217 insertions(+)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 6f3f2b0a480..38b2c4dadeb 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -77,6 +77,55 @@ class Mesh {
     double volume;    //!< volume in [cm^3]
   };
 
+  using shape_type = std::vector<size_t>;
+
+  class MaterialVolumeRT {
+  public:
+    MaterialVolumeRT(int32_t* mats, double* vols, int max_materials)
+      : materials_(mats), volumes_(vols), n_mats_(max_materials)
+    {}
+
+    void add_volume(int index_mesh, int index_material, double volume)
+    {
+      int i;
+      for (i = 0; i < n_mats_; ++i) {
+        if (materials(index_mesh, i) == index_material)
+          break;
+
+        if (materials(index_mesh, i) == -2) {
+          materials(index_mesh, i) = index_material;
+          break;
+        }
+      }
+
+      if (i >= n_mats_) {
+        fatal_error(
+          "Number of materials for mesh matvol raytrace insufficient.");
+      }
+
+      // Accumulate volume
+#pragma omp atomic
+      volumes(index_mesh, i) += volume;
+    }
+
+    int32_t& materials(int i, int j) { return materials_[i * n_mats_ + j]; }
+    const int32_t& materials(int i, int j) const
+    {
+      return materials_[i * n_mats_ + j];
+    }
+
+    double& volumes(int i, int j) { return volumes_[i * n_mats_ + j]; }
+    const double& volumes(int i, int j) const
+    {
+      return volumes_[i * n_mats_ + j];
+    }
+
+  private:
+    int32_t* materials_; //!< material index (bins, max_mats)
+    double* volumes_;    //!< volume in [cm^3] (bins, max_mats)
+    int n_mats_;
+  };
+
   // Constructors and destructor
   Mesh() = default;
   Mesh(pugi::xml_node node);
@@ -186,6 +235,17 @@ class Mesh {
   vector<MaterialVolume> material_volumes(
     int n_sample, int bin, uint64_t* seed) const;
 
+  //! Determine volume of materials within each mesh elemenet by raytracing
+  //
+  //! \param[in] ny Number of samples in y direction
+  //! \param[in] nz Number of samples in z direction
+  //! \param[in] max_materials Maximum number of materials in a single mesh
+  //!                          element
+  //! \param[inout] materials Array storing material indices
+  //! \param[inout] volumes Array storing volumes
+  void material_volumes_raytrace(int ny, int nz, int max_materials,
+    int32_t* materials, double* volumes) const;
+
   //! Determine bounding box of mesh
   //
   //! \return Bounding box of mesh
diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 16bec019863..27f77dcb637 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -29,6 +29,9 @@ class _MaterialVolume(Structure):
     ]
 
 
+arr_2d_int32 = np.ctypeslib.ndpointer(dtype=np.int32, ndim=2, flags='CONTIGUOUS')
+arr_2d_double = np.ctypeslib.ndpointer(dtype=np.double, ndim=2, flags='CONTIGUOUS')
+
 # Mesh functions
 _dll.openmc_extend_meshes.argtypes = [c_int32, c_char_p, POINTER(c_int32),
                                       POINTER(c_int32)]
@@ -55,6 +58,10 @@ class _MaterialVolume(Structure):
     POINTER(c_int), POINTER(c_uint64)]
 _dll.openmc_mesh_material_volumes.restype = c_int
 _dll.openmc_mesh_material_volumes.errcheck = _error_handler
+_dll.openmc_mesh_material_volumes_raytrace.argtypes = [
+    c_int32, c_int, c_int, c_int, arr_2d_int32, arr_2d_double]
+_dll.openmc_mesh_material_volumes_raytrace.restype = c_int
+_dll.openmc_mesh_material_volumes_raytrace.errcheck = _error_handler
 _dll.openmc_mesh_get_plot_bins.argtypes = [
     c_int32, _Position, _Position, c_int, POINTER(c_int), POINTER(c_int32)
 ]
@@ -243,6 +250,19 @@ def material_volumes(
             ])
         return volumes
 
+    def material_volumes_raytrace(
+            self, ny: int = 1000, nz: int = 1000, max_materials: int = 4):
+
+        n = self.n_elements
+        materials = np.full((n, max_materials), -2, dtype=np.int32)
+        volumes = np.zeros((n, max_materials), dtype=np.float64)
+
+        print('Raytracing...')
+        _dll.openmc_mesh_raytrace_material_volumes(
+            self._index, ny, nz, max_materials, materials, volumes)
+
+        return (np.ma.MaskedArray(materials, materials == -2), volumes)
+
     def get_plot_bins(
             self,
             origin: Sequence[float],
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 0f8b4b14e84..d9c02325a4b 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -9,6 +9,7 @@
 #include "mpi.h"
 #endif
 
+#include "xtensor/xadapt.hpp"
 #include "xtensor/xbuilder.hpp"
 #include "xtensor/xeval.hpp"
 #include "xtensor/xmath.hpp"
@@ -236,6 +237,131 @@ vector<Mesh::MaterialVolume> Mesh::material_volumes(
   return result;
 }
 
+void Mesh::material_volumes_raytrace(
+  int ny, int nz, int max_materials, int32_t* materials, double* volumes) const
+{
+  // Create object for keeping track of materials/volumes
+  MaterialVolumeRT result(materials, volumes, max_materials);
+
+  // Determine bounding box
+  auto bbox = this->bounding_box();
+
+  // Loop over rays on left face of bounding box
+  double dy = (bbox.ymax - bbox.ymin) / ny;
+  double dz = (bbox.zmax - bbox.zmin) / nz;
+
+#pragma omp parallel
+  {
+    // Preallocate vector for mesh indices and lenght fractions and p
+    std::vector<int> bins;
+    std::vector<double> length_fractions;
+    Particle p;
+
+    SourceSite site;
+    site.r = {bbox.xmin, 0.0, 0.0};
+    site.u = {1.0, 0.0, 0.0};
+    site.E = 1.0;
+    site.particle = ParticleType::neutron;
+
+#pragma omp for
+    for (int iy = 0; iy < ny; ++iy) {
+      site.r.y = bbox.ymin + (iy + 0.5) * dy;
+      for (int iz = 0; iz < nz; ++iz) {
+        site.r.z = bbox.zmin + (iz + 0.5) * dz;
+
+        p.from_source(&site);
+
+        // Determine particle's location
+        if (!exhaustive_find_cell(p)) {
+          fatal_error("Mesh is not fully contained in geometry.");
+        }
+
+        // Set birth cell attribute
+        if (p.cell_born() == C_NONE)
+          p.cell_born() = p.lowest_coord().cell;
+
+        // Initialize last cells from current cell
+        for (int j = 0; j < p.n_coord(); ++j) {
+          p.cell_last(j) = p.coord(j).cell;
+        }
+        p.n_coord_last() = p.n_coord();
+
+        while (true) {
+          // Ray trace from r_start to r_end
+          Position r0 = p.r();
+          double max_distance = bbox.xmax - r0.x;
+
+          // Find the distance to the nearest boundary
+          BoundaryInfo boundary = distance_to_boundary(p);
+
+          // Advance particle forward
+          double distance = std::min(boundary.distance, max_distance);
+          p.move_distance(distance);
+
+          // Determine what mesh elements were crossed by particle
+          bins.clear();
+          length_fractions.clear();
+          this->bins_crossed(r0, p.r(), p.u(), bins, length_fractions);
+
+          // Add volumes to any mesh elements that were crossed
+          int i_material = p.material();
+          if (i_material != C_NONE) {
+            i_material = model::materials[i_material]->id();
+          }
+          for (int i_bin = 0; i_bin < bins.size(); i_bin++) {
+            int mesh_index = bins[i_bin];
+            double length = distance * length_fractions[i_bin];
+
+            // Add volume to result
+            result.add_volume(mesh_index, i_material, length * dy * dz);
+          }
+
+          if (distance == max_distance)
+            break;
+
+          for (int j = 0; j < p.n_coord(); ++j) {
+            p.cell_last(j) = p.coord(j).cell;
+          }
+          p.n_coord_last() = p.n_coord();
+
+          // Set surface that particle is on and adjust coordinate levels
+          p.surface() = boundary.surface_index;
+          p.n_coord() = boundary.coord_level;
+
+          if (boundary.lattice_translation[0] != 0 ||
+              boundary.lattice_translation[1] != 0 ||
+              boundary.lattice_translation[2] != 0) {
+            // Particle crosses lattice boundary
+
+            cross_lattice(p, boundary);
+          } else {
+            // Particle crosses surface
+            // TODO: off-by-one
+            const auto& surf {model::surfaces[std::abs(p.surface()) - 1].get()};
+            p.cross_surface(*surf);
+          }
+        }
+      }
+    }
+  }
+
+  // for (int i = 0; i < this->n_bins(); ++i) {
+  //   for (int j = 0; j < max_materials; ++j) {
+  //     int i_material = result.materials(i, j);
+  //     if (i_material == -2)
+  //       break;
+
+  //     std::string mat_id =
+  //       (i_material == -1)
+  //         ? "void"
+  //         : fmt::format("{}", model::materials[i_material]->id());
+  //     double volume = result.volumes(i, j);
+  //     fmt::print("Element {}, Material ID={}, volume={}\n", i, mat_id,
+  //     volume);
+  //   }
+  // }
+}
+
 //==============================================================================
 // Structured Mesh implementation
 //==============================================================================
@@ -1962,6 +2088,17 @@ extern "C" int openmc_mesh_material_volumes(int32_t index, int n_sample,
   return (n == -1) ? OPENMC_E_ALLOCATE : 0;
 }
 
+extern "C" int openmc_mesh_material_volumes_raytrace(int32_t index, int ny,
+  int nz, int max_mats, int32_t* materials, double* volumes)
+{
+  if (int err = check_mesh(index))
+    return err;
+
+  model::meshes[index]->material_volumes_raytrace(
+    ny, nz, max_mats, materials, volumes);
+  return 0;
+}
+
 extern "C" int openmc_mesh_get_plot_bins(int32_t index, Position origin,
   Position width, int basis, int* pixels, int32_t* data)
 {

From 97f64481d0f3995933ab9631e355d35efacfc00f Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 18 Aug 2024 16:40:23 -0500
Subject: [PATCH 02/23] Allow Python material_volumes_raytrace to accept n_rays

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

diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 27f77dcb637..96e131bedb9 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -1,6 +1,7 @@
 from collections.abc import Mapping, Sequence
 from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, Structure,
                     create_string_buffer, c_uint64, c_size_t)
+from math import sqrt
 from random import getrandbits
 import sys
 from weakref import WeakValueDictionary
@@ -251,14 +252,28 @@ def material_volumes(
         return volumes
 
     def material_volumes_raytrace(
-            self, ny: int = 1000, nz: int = 1000, max_materials: int = 4):
-
+            self,
+            n_rays: int | tuple[int, int] = 10_000,
+            max_materials: int = 4
+    ):
+        # Determine number of rays in y/z directions
+        if isinstance(n_rays, int):
+            ll, ur = self.bounding_box
+            width_y = ur[1] - ll[1]
+            width_z = ur[2] - ll[2]
+            aspect_ratio = width_y / width_z
+            ny = int(sqrt(n_rays * aspect_ratio))
+            nz = int(sqrt(n_rays / aspect_ratio))
+        else:
+            ny, nz = n_rays
+
+        # Preallocate arrays for material indices and volumes
         n = self.n_elements
         materials = np.full((n, max_materials), -2, dtype=np.int32)
         volumes = np.zeros((n, max_materials), dtype=np.float64)
 
-        print('Raytracing...')
-        _dll.openmc_mesh_raytrace_material_volumes(
+        # Run material volume calculation
+        _dll.openmc_mesh_material_volumes_raytrace(
             self._index, ny, nz, max_materials, materials, volumes)
 
         return (np.ma.MaskedArray(materials, materials == -2), volumes)

From b3c3b54f8b01edbe72478078867b574be8b9d747 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 19 Aug 2024 07:57:09 -0500
Subject: [PATCH 03/23] Implement MeshMaterialVolumes class

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

diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 96e131bedb9..67d4c7ffd2b 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -16,10 +16,11 @@
 from .material import Material
 from .plot import _Position
 from ..bounding_box import BoundingBox
+from ..checkvalue import PathLike
 
 __all__ = [
     'Mesh', 'RegularMesh', 'RectilinearMesh', 'CylindricalMesh',
-    'SphericalMesh', 'UnstructuredMesh', 'meshes'
+    'SphericalMesh', 'UnstructuredMesh', 'meshes', 'MeshMaterialVolumes'
 ]
 
 
@@ -121,6 +122,45 @@ class _MaterialVolume(Structure):
 _dll.openmc_spherical_mesh_set_grid.errcheck = _error_handler
 
 
+class MeshMaterialVolumes(Mapping):
+    def __init__(self, materials, volumes):
+        self._materials = materials
+        self._volumes = volumes
+
+    @property
+    def num_elements(self) -> int:
+        return self._volumes.shape[0]
+
+    def __iter__(self):
+        for mat in np.unique(self._materials):
+            if mat > 0:
+                yield mat
+
+    def __len__(self):
+        return (np.unique(self._materials) > 0).sum()
+
+    def __repr__(self) -> str:
+        ids, counts = np.unique(self._materials, return_counts=True)
+        return '{' + '\n '.join(
+            f'{id}: <{count} nonzero volumes>' for id, count in zip(ids, counts) if id > 0) + '}'
+
+    def __getitem__(self, material_id: int):
+        volumes = np.zeros(self.num_elements)
+        for i in range(self._volumes.shape[1]):
+            indices = (self._materials[:, i] == material_id)
+            volumes[indices] = self._volumes[indices, i]
+        return volumes
+
+    def save(self, filename: PathLike):
+        np.savez_compressed(
+            filename, materials=self._materials, volumes=self._volumes)
+
+    @classmethod
+    def from_npz(cls, filename: PathLike):
+        filedata = np.load(filename)
+        return cls(filedata['materials'], filedata['volumes'])
+
+
 class Mesh(_FortranObjectWithID):
     """Base class to represent mesh objects
 
@@ -255,7 +295,29 @@ def material_volumes_raytrace(
             self,
             n_rays: int | tuple[int, int] = 10_000,
             max_materials: int = 4
-    ):
+    ) -> MeshMaterialVolumes:
+        """Determine volume of materials in each mesh element.
+
+        This method works by raytracing repeatedly through the mesh from one
+        side to count of the estimated volume of each material in all mesh
+        elements.
+
+        Parameters
+        ----------
+        n_rays : int or 2-tuple of int
+            Total number of rays to follow. The rays start on an x plane and are
+            evenly distributed over the y and z dimensions. When specified as a
+            2-tuple, it is interpreted as the number of rays in the y and z
+            dimensions.
+        max_materials : int, optional
+            Maximum number of materials in any given mesh element.
+
+        Returns
+        -------
+        Dictionary-like object that maps material IDs to an array of volumes
+        equal in size to the number of mesh elements.
+
+        """
         # Determine number of rays in y/z directions
         if isinstance(n_rays, int):
             ll, ur = self.bounding_box
@@ -276,7 +338,7 @@ def material_volumes_raytrace(
         _dll.openmc_mesh_material_volumes_raytrace(
             self._index, ny, nz, max_materials, materials, volumes)
 
-        return (np.ma.MaskedArray(materials, materials == -2), volumes)
+        return MeshMaterialVolumes(materials, volumes)
 
     def get_plot_bins(
             self,
@@ -754,7 +816,6 @@ def __getitem__(self, key):
             raise KeyError(str(e))
         return _get_mesh(index.value)
 
-
     def __iter__(self):
         for i in range(len(self)):
             yield _get_mesh(i).id

From d794d6547b820080fe2859dd2a3927e94aaceee0 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 19 Aug 2024 10:53:34 -0500
Subject: [PATCH 04/23] Normalize material volumes based on known element
 volumes

---
 src/mesh.cpp | 29 ++++++++++++++---------------
 1 file changed, 14 insertions(+), 15 deletions(-)

diff --git a/src/mesh.cpp b/src/mesh.cpp
index d9c02325a4b..4e3f9b046ae 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -345,21 +345,20 @@ void Mesh::material_volumes_raytrace(
     }
   }
 
-  // for (int i = 0; i < this->n_bins(); ++i) {
-  //   for (int j = 0; j < max_materials; ++j) {
-  //     int i_material = result.materials(i, j);
-  //     if (i_material == -2)
-  //       break;
-
-  //     std::string mat_id =
-  //       (i_material == -1)
-  //         ? "void"
-  //         : fmt::format("{}", model::materials[i_material]->id());
-  //     double volume = result.volumes(i, j);
-  //     fmt::print("Element {}, Material ID={}, volume={}\n", i, mat_id,
-  //     volume);
-  //   }
-  // }
+  // Normalize based on known volumes of elements
+  for (int i = 0; i < this->n_bins(); ++i) {
+    // Estimated total volume in element i
+    double volume = 0.0;
+    for (int j = 0; j < max_materials; ++j) {
+      volume += result.volumes(i, j);
+    }
+
+    // Renormalize volumes based on known volume of elemnet i
+    double norm = this->volume(i) / volume;
+    for (int j = 0; j < max_materials; ++j) {
+      result.volumes(i, j) *= norm;
+    }
+  }
 }
 
 //==============================================================================

From 56fb913447c11e476e88f9f14a7b12d5c8a5c3a1 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 19 Aug 2024 22:52:11 -0500
Subject: [PATCH 05/23] Fix typos

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

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 38b2c4dadeb..91c43acb11c 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -216,7 +216,7 @@ class Mesh {
 
   virtual std::string get_mesh_type() const = 0;
 
-  //! Determine volume of materials within a single mesh elemenet
+  //! Determine volume of materials within a single mesh element
   //
   //! \param[in] n_sample Number of samples within each element
   //! \param[in] bin Index of mesh element
@@ -226,7 +226,7 @@ class Mesh {
   int material_volumes(int n_sample, int bin, gsl::span<MaterialVolume> volumes,
     uint64_t* seed) const;
 
-  //! Determine volume of materials within a single mesh elemenet
+  //! Determine volume of materials within a single mesh element
   //
   //! \param[in] n_sample Number of samples within each element
   //! \param[in] bin Index of mesh element
@@ -235,7 +235,7 @@ class Mesh {
   vector<MaterialVolume> material_volumes(
     int n_sample, int bin, uint64_t* seed) const;
 
-  //! Determine volume of materials within each mesh elemenet by raytracing
+  //! Determine volume of materials within each mesh element by raytracing
   //
   //! \param[in] ny Number of samples in y direction
   //! \param[in] nz Number of samples in z direction

From 7d8647981f7985e7ce7f9c9d47dfe2a6e98c1cb5 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 25 Aug 2024 16:44:50 -0500
Subject: [PATCH 06/23] Add openmc.Mesh.material_volumes method

---
 openmc/lib/mesh.py |  38 +-------------
 openmc/mesh.py     | 126 ++++++++++++++++++++++++++++++++++++++++++++-
 2 files changed, 125 insertions(+), 39 deletions(-)

diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 67d4c7ffd2b..7b3d398b093 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -17,6 +17,7 @@
 from .plot import _Position
 from ..bounding_box import BoundingBox
 from ..checkvalue import PathLike
+from ..mesh import MeshMaterialVolumes
 
 __all__ = [
     'Mesh', 'RegularMesh', 'RectilinearMesh', 'CylindricalMesh',
@@ -122,43 +123,6 @@ class _MaterialVolume(Structure):
 _dll.openmc_spherical_mesh_set_grid.errcheck = _error_handler
 
 
-class MeshMaterialVolumes(Mapping):
-    def __init__(self, materials, volumes):
-        self._materials = materials
-        self._volumes = volumes
-
-    @property
-    def num_elements(self) -> int:
-        return self._volumes.shape[0]
-
-    def __iter__(self):
-        for mat in np.unique(self._materials):
-            if mat > 0:
-                yield mat
-
-    def __len__(self):
-        return (np.unique(self._materials) > 0).sum()
-
-    def __repr__(self) -> str:
-        ids, counts = np.unique(self._materials, return_counts=True)
-        return '{' + '\n '.join(
-            f'{id}: <{count} nonzero volumes>' for id, count in zip(ids, counts) if id > 0) + '}'
-
-    def __getitem__(self, material_id: int):
-        volumes = np.zeros(self.num_elements)
-        for i in range(self._volumes.shape[1]):
-            indices = (self._materials[:, i] == material_id)
-            volumes[indices] = self._volumes[indices, i]
-        return volumes
-
-    def save(self, filename: PathLike):
-        np.savez_compressed(
-            filename, materials=self._materials, volumes=self._volumes)
-
-    @classmethod
-    def from_npz(cls, filename: PathLike):
-        filedata = np.load(filename)
-        return cls(filedata['materials'], filedata['volumes'])
 
 
 class Mesh(_FortranObjectWithID):
diff --git a/openmc/mesh.py b/openmc/mesh.py
index a706b8fa811..6cf0f50aa65 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -1,12 +1,11 @@
 from __future__ import annotations
 import warnings
 from abc import ABC, abstractmethod
-from collections.abc import Iterable, Sequence
+from collections.abc import Iterable, Sequence, Mapping
 from functools import wraps
 from math import pi, sqrt, atan2
 from numbers import Integral, Real
 from pathlib import Path
-import tempfile
 
 import h5py
 import lxml.etree as ET
@@ -21,6 +20,70 @@
 from .surface import _BOUNDARY_TYPES
 
 
+class MeshMaterialVolumes(Mapping):
+    """Results from a material volume in mesh calculation.
+
+    Parameters
+    ----------
+    materials : np.ndarray
+        Array of shape (elements, max_materials) storing material IDs
+    volumes : np.ndarray
+        Array of shape (elements, max_materials) storing material volumes
+
+    """
+    def __init__(self, materials: np.ndarray, volumes):
+        self._materials = materials
+        self._volumes = volumes
+
+    @property
+    def num_elements(self) -> int:
+        return self._volumes.shape[0]
+
+    def __iter__(self):
+        for mat in np.unique(self._materials):
+            if mat > 0:
+                yield mat
+
+    def __len__(self) -> int:
+        return (np.unique(self._materials) > 0).sum()
+
+    def __repr__(self) -> str:
+        ids, counts = np.unique(self._materials, return_counts=True)
+        return '{' + '\n '.join(
+            f'{id}: <{count} nonzero volumes>' for id, count in zip(ids, counts) if id > 0) + '}'
+
+    def __getitem__(self, material_id: int) -> np.ndarray:
+        volumes = np.zeros(self.num_elements)
+        for i in range(self._volumes.shape[1]):
+            indices = (self._materials[:, i] == material_id)
+            volumes[indices] = self._volumes[indices, i]
+        return volumes
+
+    def save(self, filename: PathLike):
+        """Save material volumes to a .npz file.
+
+        Parameters
+        ----------
+        filename : path-like
+            Filename where data will be saved
+        """
+        np.savez_compressed(
+            filename, materials=self._materials, volumes=self._volumes)
+
+    @classmethod
+    def from_npz(cls, filename: PathLike) -> MeshMaterialVolumes:
+        """Generate material volumes from a .npz file
+
+        Parameters
+        ----------
+        filename : path-like
+            File where data will be read from
+
+        """
+        filedata = np.load(filename)
+        return cls(filedata['materials'], filedata['volumes'])
+
+
 class MeshBase(IDManagerMixin, ABC):
     """A mesh that partitions geometry for tallying purposes.
 
@@ -252,6 +315,65 @@ def get_homogenized_materials(
 
         return homogenized_materials
 
+    def material_volumes(
+            self,
+            model: openmc.Model,
+            n_rays: int | tuple[int, int] = 10_000,
+            max_materials: int = 4,
+            **kwargs
+    ) -> list[openmc.Material]:
+        """Determine volume of materials in each mesh element.
+
+        This method works by raytracing repeatedly through the mesh from one
+        side to count of the estimated volume of each material in all mesh
+        elements.
+
+        Parameters
+        ----------
+        model : openmc.Model
+            Model containing materials.
+        n_rays : int or 2-tuple of int
+            Total number of rays to follow. The rays start on an x plane and are
+            evenly distributed over the y and z dimensions. When specified as a
+            2-tuple, it is interpreted as the number of rays in the y and z
+            dimensions.
+        max_materials : int, optional
+            Maximum number of materials in any given mesh element.
+        **kwargs : dict
+            Keyword arguments passed to :func:`openmc.lib.init`
+
+        Returns
+        -------
+        Dictionary-like object that maps material IDs to an array of volumes
+        equal in size to the number of mesh elements.
+
+        """
+        import openmc.lib
+
+        with change_directory(tmpdir=True):
+            # In order to get mesh into model, we temporarily replace the
+            # tallies with a single mesh tally using the current mesh
+            original_tallies = model.tallies
+            new_tally = openmc.Tally()
+            new_tally.filters = [openmc.MeshFilter(self)]
+            new_tally.scores = ['flux']
+            model.tallies = [new_tally]
+
+            # Export model to XML
+            model.export_to_model_xml()
+
+            # Get material volume fractions
+            kwargs.setdefault('output', False)
+            openmc.lib.init(['-c'], **kwargs)
+            mesh = openmc.lib.tallies[new_tally.id].filters[0].mesh
+            volumes = mesh.material_volumes_raytrace(n_rays, max_materials)
+            openmc.lib.finalize()
+
+            # Restore original tallies
+            model.tallies = original_tallies
+
+        return volumes
+
 
 class StructuredMesh(MeshBase):
     """A base class for structured mesh functionality

From c35460809ec7bf1654b00d4dc3091b5a64ebb432 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 25 Aug 2024 17:13:51 -0500
Subject: [PATCH 07/23] Automatically reallocate if max_materials is
 insufficient

---
 include/openmc/mesh.h | 11 +++++++++--
 openmc/lib/mesh.py    | 14 ++++++++++++--
 src/mesh.cpp          | 16 ++++++++++++++--
 3 files changed, 35 insertions(+), 6 deletions(-)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 91c43acb11c..4649df57f71 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -89,18 +89,22 @@ class Mesh {
     {
       int i;
       for (i = 0; i < n_mats_; ++i) {
+        // Check whether material already is present
         if (materials(index_mesh, i) == index_material)
           break;
 
+        // If not already present, check for unused position (-2)
         if (materials(index_mesh, i) == -2) {
           materials(index_mesh, i) = index_material;
           break;
         }
       }
 
+      // If maximum number of materials exceeded, set a flag that can be checked
+      // later
       if (i >= n_mats_) {
-        fatal_error(
-          "Number of materials for mesh matvol raytrace insufficient.");
+        too_many_mats_ = true;
+        return;
       }
 
       // Accumulate volume
@@ -120,10 +124,13 @@ class Mesh {
       return volumes_[i * n_mats_ + j];
     }
 
+    bool too_many_mats() const { return too_many_mats_; }
+
   private:
     int32_t* materials_; //!< material index (bins, max_mats)
     double* volumes_;    //!< volume in [cm^3] (bins, max_mats)
     int n_mats_;
+    bool too_many_mats_ = false;
   };
 
   // Constructors and destructor
diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 7b3d398b093..48cb8621bf0 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -299,8 +299,18 @@ def material_volumes_raytrace(
         volumes = np.zeros((n, max_materials), dtype=np.float64)
 
         # Run material volume calculation
-        _dll.openmc_mesh_material_volumes_raytrace(
-            self._index, ny, nz, max_materials, materials, volumes)
+        while True:
+            try:
+                _dll.openmc_mesh_material_volumes_raytrace(
+                    self._index, ny, nz, max_materials, materials, volumes)
+            except AllocationError:
+                # Increase size of result array and try again
+                max_materials *= 2
+                materials = np.full((n, max_materials), -2, dtype=np.int32)
+                volumes = np.zeros((n, max_materials), dtype=np.float64)
+            else:
+                # If no error, break out of loop
+                break
 
         return MeshMaterialVolumes(materials, volumes)
 
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 4e3f9b046ae..e57aa516bd6 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -345,6 +345,12 @@ void Mesh::material_volumes_raytrace(
     }
   }
 
+  // Check whether max number of materials was exceeded
+  if (result.too_many_mats()) {
+    throw std::runtime_error("Maximum number of materials for mesh material "
+                             "volume calculation insufficient.");
+  }
+
   // Normalize based on known volumes of elements
   for (int i = 0; i < this->n_bins(); ++i) {
     // Estimated total volume in element i
@@ -2093,8 +2099,14 @@ extern "C" int openmc_mesh_material_volumes_raytrace(int32_t index, int ny,
   if (int err = check_mesh(index))
     return err;
 
-  model::meshes[index]->material_volumes_raytrace(
-    ny, nz, max_mats, materials, volumes);
+  try {
+    model::meshes[index]->material_volumes_raytrace(
+      ny, nz, max_mats, materials, volumes);
+  } catch (const std::exception& e) {
+    set_errmsg(e.what());
+    return OPENMC_E_ALLOCATE;
+  }
+
   return 0;
 }
 

From a288054e8f97b2767d58d018970e9e74f6fae257 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 25 Aug 2024 22:00:45 -0500
Subject: [PATCH 08/23] Avoid potential race condition on material indexing

---
 include/openmc/mesh.h | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 4649df57f71..0090fc61ba6 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -85,17 +85,18 @@ class Mesh {
       : materials_(mats), volumes_(vols), n_mats_(max_materials)
     {}
 
-    void add_volume(int index_mesh, int index_material, double volume)
+    void add_volume(int index_elem, int index_material, double volume)
     {
       int i;
+#pragma omp critical(MeshMatVol)
       for (i = 0; i < n_mats_; ++i) {
         // Check whether material already is present
-        if (materials(index_mesh, i) == index_material)
+        if (materials(index_elem, i) == index_material)
           break;
 
         // If not already present, check for unused position (-2)
-        if (materials(index_mesh, i) == -2) {
-          materials(index_mesh, i) = index_material;
+        if (materials(index_elem, i) == -2) {
+          materials(index_elem, i) = index_material;
           break;
         }
       }
@@ -109,7 +110,7 @@ class Mesh {
 
       // Accumulate volume
 #pragma omp atomic
-      volumes(index_mesh, i) += volume;
+      volumes(index_elem, i) += volume;
     }
 
     int32_t& materials(int i, int j) { return materials_[i * n_mats_ + j]; }

From 5b9041ae992fb26b5c22d68a3dacebccbd3f386c Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 26 Aug 2024 23:03:44 -0500
Subject: [PATCH 09/23] Add by_element method on MeshMaterialVolumes

---
 openmc/mesh.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 6cf0f50aa65..1128ec72c27 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -59,6 +59,26 @@ def __getitem__(self, material_id: int) -> np.ndarray:
             volumes[indices] = self._volumes[indices, i]
         return volumes
 
+    def by_element(self, index_elem: int) -> list[tuple[int | None, float]]:
+        """Get a list of volumes for each material within a specific element.
+
+        Parameters
+        ----------
+        index_elem : int
+            Mesh element index
+
+        Returns
+        -------
+        list of tuple of (material ID, volume)
+
+        """
+        max_mats = self._volumes.shape[1]
+        return [
+            (m if m > -1 else None, self._volumes[index_elem, i])
+            for i in range(max_mats)
+            if (m := self._materials[index_elem, i]) != -2
+        ]
+
     def save(self, filename: PathLike):
         """Save material volumes to a .npz file.
 

From 0afbd9b22456d15c0d7cc2b39524717c5a18d84c Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 26 Aug 2024 23:04:22 -0500
Subject: [PATCH 10/23] Avoid fatal_error in middle of
 material_volumes_raytrace

---
 src/mesh.cpp | 19 +++++++++++++++----
 1 file changed, 15 insertions(+), 4 deletions(-)

diff --git a/src/mesh.cpp b/src/mesh.cpp
index e57aa516bd6..2ad4d804438 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -34,6 +34,7 @@
 #include "openmc/random_dist.h"
 #include "openmc/search.h"
 #include "openmc/settings.h"
+#include "openmc/string_utils.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
 #include "openmc/volume_calc.h"
@@ -250,6 +251,9 @@ void Mesh::material_volumes_raytrace(
   double dy = (bbox.ymax - bbox.ymin) / ny;
   double dz = (bbox.zmax - bbox.zmin) / nz;
 
+  // Set flag for mesh being contained within model
+  bool out_of_model = false;
+
 #pragma omp parallel
   {
     // Preallocate vector for mesh indices and lenght fractions and p
@@ -273,7 +277,8 @@ void Mesh::material_volumes_raytrace(
 
         // Determine particle's location
         if (!exhaustive_find_cell(p)) {
-          fatal_error("Mesh is not fully contained in geometry.");
+          out_of_model = true;
+          continue;
         }
 
         // Set birth cell attribute
@@ -345,8 +350,10 @@ void Mesh::material_volumes_raytrace(
     }
   }
 
-  // Check whether max number of materials was exceeded
-  if (result.too_many_mats()) {
+  // Check for errors
+  if (out_of_model) {
+    throw std::runtime_error("Mesh not fully contained in geometry.");
+  } else if (result.too_many_mats()) {
     throw std::runtime_error("Maximum number of materials for mesh material "
                              "volume calculation insufficient.");
   }
@@ -2104,7 +2111,11 @@ extern "C" int openmc_mesh_material_volumes_raytrace(int32_t index, int ny,
       ny, nz, max_mats, materials, volumes);
   } catch (const std::exception& e) {
     set_errmsg(e.what());
-    return OPENMC_E_ALLOCATE;
+    if (starts_with(e.what(), "Mesh")) {
+      return OPENMC_E_GEOMETRY;
+    } else {
+      return OPENMC_E_ALLOCATE;
+    }
   }
 
   return 0;

From af9dfaa8c4a069f705802ace1db75f212025ae5f Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 26 Aug 2024 23:29:40 -0500
Subject: [PATCH 11/23] Update tests to use material_volumes_raytrace

---
 tests/unit_tests/test_lib.py | 44 ++++++++++++++++++------------------
 1 file changed, 22 insertions(+), 22 deletions(-)

diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py
index 64c16c238ea..b7c8761af34 100644
--- a/tests/unit_tests/test_lib.py
+++ b/tests/unit_tests/test_lib.py
@@ -600,19 +600,19 @@ def test_regular_mesh(lib_init):
     mesh.dimension = (2, 2, 1)
     mesh.set_parameters(lower_left=(-0.63, -0.63, -0.5),
                         upper_right=(0.63, 0.63, 0.5))
-    vols = mesh.material_volumes()
-    assert len(vols) == 4
-    for elem_vols in vols:
+    vols = mesh.material_volumes_raytrace()
+    assert vols.num_elements == 4
+    for i in range(vols.num_elements):
+        elem_vols = vols.by_element(i)
         assert sum(f[1] for f in elem_vols) == pytest.approx(1.26 * 1.26 / 4)
 
-    # If the mesh extends beyond the boundaries of the model, the volumes should
-    # still be reported correctly
+    # If the mesh extends beyond the boundaries of the model, we should get a
+    # GeometryError
     mesh.dimension = (1, 1, 1)
     mesh.set_parameters(lower_left=(-1.0, -1.0, -0.5),
                         upper_right=(1.0, 1.0, 0.5))
-    vols = mesh.material_volumes(100_000)
-    for elem_vols in vols:
-        assert sum(f[1] for f in elem_vols) == pytest.approx(1.26 * 1.26, 1e-2)
+    with pytest.raises(exc.GeometryError, match="not fully contained"):
+        vols = mesh.material_volumes_raytrace()
 
 
 def test_regular_mesh_get_plot_bins(lib_init):
@@ -682,12 +682,12 @@ def test_rectilinear_mesh(lib_init):
     w = 1.26
     mesh.set_grid([-w/2, -w/4, w/2], [-w/2, -w/4, w/2], [-0.5, 0.5])
 
-    vols = mesh.material_volumes()
-    assert len(vols) == 4
-    assert sum(f[1] for f in vols[0]) == pytest.approx(w/4 * w/4)
-    assert sum(f[1] for f in vols[1]) == pytest.approx(w/4 * 3*w/4)
-    assert sum(f[1] for f in vols[2]) == pytest.approx(3*w/4 * w/4)
-    assert sum(f[1] for f in vols[3]) == pytest.approx(3*w/4 * 3*w/4)
+    vols = mesh.material_volumes_raytrace()
+    assert vols.num_elements == 4
+    assert sum(f[1] for f in vols.by_element(0)) == pytest.approx(w/4 * w/4)
+    assert sum(f[1] for f in vols.by_element(1)) == pytest.approx(w/4 * 3*w/4)
+    assert sum(f[1] for f in vols.by_element(2)) == pytest.approx(3*w/4 * w/4)
+    assert sum(f[1] for f in vols.by_element(3)) == pytest.approx(3*w/4 * 3*w/4)
 
 
 def test_cylindrical_mesh(lib_init):
@@ -736,12 +736,12 @@ def test_cylindrical_mesh(lib_init):
     z_grid = (-0.5, 0.5)
     mesh.set_grid(r_grid, phi_grid, z_grid)
 
-    vols = mesh.material_volumes()
-    assert len(vols) == 6
+    vols = mesh.material_volumes_raytrace()
+    assert vols.num_elements == 6
     for i in range(0, 6, 2):
-        assert sum(f[1] for f in vols[i]) == pytest.approx(pi * 0.25**2 / 3)
+        assert sum(f[1] for f in vols.by_element(i)) == pytest.approx(pi * 0.25**2 / 3)
     for i in range(1, 6, 2):
-        assert sum(f[1] for f in vols[i]) == pytest.approx(pi * (0.5**2 - 0.25**2) / 3)
+        assert sum(f[1] for f in vols.by_element(i)) == pytest.approx(pi * (0.5**2 - 0.25**2) / 3)
 
 
 def test_spherical_mesh(lib_init):
@@ -794,15 +794,15 @@ def test_spherical_mesh(lib_init):
     phi_grid = np.linspace(0., 2.0*pi, 4)
     mesh.set_grid(r_grid, theta_grid, phi_grid)
 
-    vols = mesh.material_volumes()
-    assert len(vols) == 12
+    vols = mesh.material_volumes_raytrace()
+    assert vols.num_elements == 12
     d_theta = theta_grid[1] - theta_grid[0]
     d_phi = phi_grid[1] - phi_grid[0]
     for i in range(0, 12, 2):
-        assert sum(f[1] for f in vols[i]) == pytest.approx(
+        assert sum(f[1] for f in vols.by_element(i)) == pytest.approx(
             0.25**3 / 3 * d_theta * d_phi * 2/pi)
     for i in range(1, 12, 2):
-        assert sum(f[1] for f in vols[i]) == pytest.approx(
+        assert sum(f[1] for f in vols.by_element(i)) == pytest.approx(
             (0.5**3 - 0.25**3) / 3 * d_theta * d_phi * 2/pi)
 
 

From 7f0a3df9f7408bfb6fee540613b16a8703c092ad Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Tue, 27 Aug 2024 00:13:19 -0500
Subject: [PATCH 12/23] Use raytraced version in get_homogenized_materials

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

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 1128ec72c27..24ae0b62bbf 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -231,8 +231,7 @@ def from_xml_element(cls, elem: ET.Element):
     def get_homogenized_materials(
             self,
             model: openmc.Model,
-            n_samples: int = 10_000,
-            prn_seed: int | None = None,
+            n_samples: int | tuple[int, int] = 10_000,
             include_void: bool = True,
             **kwargs
     ) -> list[openmc.Material]:
@@ -245,15 +244,15 @@ def get_homogenized_materials(
         model : openmc.Model
             Model containing materials to be homogenized and the associated
             geometry.
-        n_samples : int
-            Number of samples in each mesh element.
-        prn_seed : int, optional
-            Pseudorandom number generator (PRNG) seed; if None, one will be
-            generated randomly.
+        n_samples : int or 2-tuple of int
+            Total number of rays to sample. The rays start on an x plane and are
+            evenly distributed over the y and z dimensions. When specified as a
+            2-tuple, it is interpreted as the number of rays in the y and z
+            dimensions.
         include_void : bool, optional
             Whether homogenization should include voids.
         **kwargs
-            Keyword-arguments passed to :func:`openmc.lib.init`.
+            Keyword-arguments passed to :meth:`MeshBase.material_volumes`.
 
         Returns
         -------
@@ -261,34 +260,8 @@ def get_homogenized_materials(
             Homogenized material in each mesh element
 
         """
-        import openmc.lib
-
-        with change_directory(tmpdir=True):
-            # In order to get mesh into model, we temporarily replace the
-            # tallies with a single mesh tally using the current mesh
-            original_tallies = model.tallies
-            new_tally = openmc.Tally()
-            new_tally.filters = [openmc.MeshFilter(self)]
-            new_tally.scores = ['flux']
-            model.tallies = [new_tally]
-
-            # Export model to XML
-            model.export_to_model_xml()
-
-            # Get material volume fractions
-            openmc.lib.init(**kwargs)
-            mesh = openmc.lib.tallies[new_tally.id].filters[0].mesh
-            mat_volume_by_element = [
-                [
-                    (mat.id if mat is not None else None, volume)
-                    for mat, volume in mat_volume_list
-                ]
-                for mat_volume_list in mesh.material_volumes(n_samples, prn_seed)
-            ]
-            openmc.lib.finalize()
-
-            # Restore original tallies
-            model.tallies = original_tallies
+        vols = self.material_volumes(model, n_samples, **kwargs)
+        mat_volume_by_element = [vols.by_element(i) for i in range(vols.num_elements)]
 
         # Create homogenized material for each element
         materials = model.geometry.get_all_materials()
@@ -341,7 +314,7 @@ def material_volumes(
             n_rays: int | tuple[int, int] = 10_000,
             max_materials: int = 4,
             **kwargs
-    ) -> list[openmc.Material]:
+    ) -> MeshMaterialVolumes:
         """Determine volume of materials in each mesh element.
 
         This method works by raytracing repeatedly through the mesh from one

From 93e517101da213e8c1d07fbe5c9e8bbb80ddd982 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 28 Aug 2024 22:53:44 -0500
Subject: [PATCH 13/23] Rename n_rays to n_samples

---
 openmc/lib/mesh.py | 14 +++++++-------
 openmc/mesh.py     |  8 ++++----
 2 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 48cb8621bf0..0020636a20f 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -257,7 +257,7 @@ def material_volumes(
 
     def material_volumes_raytrace(
             self,
-            n_rays: int | tuple[int, int] = 10_000,
+            n_samples: int | tuple[int, int] = 10_000,
             max_materials: int = 4
     ) -> MeshMaterialVolumes:
         """Determine volume of materials in each mesh element.
@@ -268,8 +268,8 @@ def material_volumes_raytrace(
 
         Parameters
         ----------
-        n_rays : int or 2-tuple of int
-            Total number of rays to follow. The rays start on an x plane and are
+        n_samples : int or 2-tuple of int
+            Total number of rays to sample. The rays start on an x plane and are
             evenly distributed over the y and z dimensions. When specified as a
             2-tuple, it is interpreted as the number of rays in the y and z
             dimensions.
@@ -283,15 +283,15 @@ def material_volumes_raytrace(
 
         """
         # Determine number of rays in y/z directions
-        if isinstance(n_rays, int):
+        if isinstance(n_samples, int):
             ll, ur = self.bounding_box
             width_y = ur[1] - ll[1]
             width_z = ur[2] - ll[2]
             aspect_ratio = width_y / width_z
-            ny = int(sqrt(n_rays * aspect_ratio))
-            nz = int(sqrt(n_rays / aspect_ratio))
+            ny = int(sqrt(n_samples * aspect_ratio))
+            nz = int(sqrt(n_samples / aspect_ratio))
         else:
-            ny, nz = n_rays
+            ny, nz = n_samples
 
         # Preallocate arrays for material indices and volumes
         n = self.n_elements
diff --git a/openmc/mesh.py b/openmc/mesh.py
index 24ae0b62bbf..a6082e60a84 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -311,7 +311,7 @@ def get_homogenized_materials(
     def material_volumes(
             self,
             model: openmc.Model,
-            n_rays: int | tuple[int, int] = 10_000,
+            n_samples: int | tuple[int, int] = 10_000,
             max_materials: int = 4,
             **kwargs
     ) -> MeshMaterialVolumes:
@@ -325,8 +325,8 @@ def material_volumes(
         ----------
         model : openmc.Model
             Model containing materials.
-        n_rays : int or 2-tuple of int
-            Total number of rays to follow. The rays start on an x plane and are
+        n_samples : int or 2-tuple of int
+            Total number of rays to sample. The rays start on an x plane and are
             evenly distributed over the y and z dimensions. When specified as a
             2-tuple, it is interpreted as the number of rays in the y and z
             dimensions.
@@ -359,7 +359,7 @@ def material_volumes(
             kwargs.setdefault('output', False)
             openmc.lib.init(['-c'], **kwargs)
             mesh = openmc.lib.tallies[new_tally.id].filters[0].mesh
-            volumes = mesh.material_volumes_raytrace(n_rays, max_materials)
+            volumes = mesh.material_volumes_raytrace(n_samples, max_materials)
             openmc.lib.finalize()
 
             # Restore original tallies

From 006cb1f691a74cc34cefe2438c65938906c2973a Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 28 Aug 2024 23:10:13 -0500
Subject: [PATCH 14/23] Remove point sampling version of material_volumes

---
 include/openmc/capi.h        |   4 +-
 include/openmc/mesh.h        |  25 +-------
 openmc/lib/mesh.py           |  64 +-------------------
 openmc/mesh.py               |   2 +-
 src/mesh.cpp                 | 112 ++---------------------------------
 tests/unit_tests/test_lib.py |  10 ++--
 6 files changed, 17 insertions(+), 200 deletions(-)

diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index 9401156a64f..223a6b1877e 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -107,8 +107,8 @@ int openmc_mesh_get_id(int32_t index, int32_t* id);
 int openmc_mesh_set_id(int32_t index, int32_t id);
 int openmc_mesh_get_n_elements(int32_t index, size_t* n);
 int openmc_mesh_get_volumes(int32_t index, double* volumes);
-int openmc_mesh_material_volumes(int32_t index, int n_sample, int bin,
-  int result_size, void* result, int* hits, uint64_t* seed);
+int openmc_mesh_material_volumes(int32_t index, int ny, int nz, int max_mats,
+  int32_t* materials, double* volumes);
 int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
 int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
 int openmc_new_filter(const char* type, int32_t* index);
diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 0090fc61ba6..e12edf2d7e3 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -224,26 +224,7 @@ class Mesh {
 
   virtual std::string get_mesh_type() const = 0;
 
-  //! Determine volume of materials within a single mesh element
-  //
-  //! \param[in] n_sample Number of samples within each element
-  //! \param[in] bin Index of mesh element
-  //! \param[out] Array of (material index, volume) for desired element
-  //! \param[inout] seed Pseudorandom number seed
-  //! \return Number of materials within element
-  int material_volumes(int n_sample, int bin, gsl::span<MaterialVolume> volumes,
-    uint64_t* seed) const;
-
-  //! Determine volume of materials within a single mesh element
-  //
-  //! \param[in] n_sample Number of samples within each element
-  //! \param[in] bin Index of mesh element
-  //! \param[inout] seed Pseudorandom number seed
-  //! \return Vector of (material index, volume) for desired element
-  vector<MaterialVolume> material_volumes(
-    int n_sample, int bin, uint64_t* seed) const;
-
-  //! Determine volume of materials within each mesh element by raytracing
+  //! Determine volume of materials within each mesh element
   //
   //! \param[in] ny Number of samples in y direction
   //! \param[in] nz Number of samples in z direction
@@ -251,8 +232,8 @@ class Mesh {
   //!                          element
   //! \param[inout] materials Array storing material indices
   //! \param[inout] volumes Array storing volumes
-  void material_volumes_raytrace(int ny, int nz, int max_materials,
-    int32_t* materials, double* volumes) const;
+  void material_volumes(int ny, int nz, int max_materials, int32_t* materials,
+    double* volumes) const;
 
   //! Determine bounding box of mesh
   //
diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 0020636a20f..4f4232832f8 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -57,14 +57,9 @@ class _MaterialVolume(Structure):
 _dll.openmc_mesh_bounding_box.restype = c_int
 _dll.openmc_mesh_bounding_box.errcheck = _error_handler
 _dll.openmc_mesh_material_volumes.argtypes = [
-    c_int32, c_int, c_int, c_int, POINTER(_MaterialVolume),
-    POINTER(c_int), POINTER(c_uint64)]
+    c_int32, c_int, c_int, c_int, arr_2d_int32, arr_2d_double]
 _dll.openmc_mesh_material_volumes.restype = c_int
 _dll.openmc_mesh_material_volumes.errcheck = _error_handler
-_dll.openmc_mesh_material_volumes_raytrace.argtypes = [
-    c_int32, c_int, c_int, c_int, arr_2d_int32, arr_2d_double]
-_dll.openmc_mesh_material_volumes_raytrace.restype = c_int
-_dll.openmc_mesh_material_volumes_raytrace.errcheck = _error_handler
 _dll.openmc_mesh_get_plot_bins.argtypes = [
     c_int32, _Position, _Position, c_int, POINTER(c_int), POINTER(c_int32)
 ]
@@ -201,61 +196,6 @@ def bounding_box(self) -> BoundingBox:
         return BoundingBox(ll, ur)
 
     def material_volumes(
-            self,
-            n_samples: int = 10_000,
-            prn_seed: int | None = None
-    ) -> list[list[tuple[Material, float]]]:
-        """Determine volume of materials in each mesh element
-
-        .. versionadded:: 0.15.0
-
-        Parameters
-        ----------
-        n_samples : int
-            Number of samples in each mesh element
-        prn_seed : int
-            Pseudorandom number generator (PRNG) seed; if None, one will be
-            generated randomly.
-
-        Returns
-        -------
-        List of tuple of (material, volume) for each mesh element. Void volume
-        is represented by having a value of None in the first element of a
-        tuple.
-
-        """
-        if n_samples <= 0:
-            raise ValueError("Number of samples must be positive")
-        if prn_seed is None:
-            prn_seed = getrandbits(63)
-        prn_seed = c_uint64(prn_seed)
-
-        # Preallocate space for MaterialVolume results
-        size = 16
-        result = (_MaterialVolume * size)()
-
-        hits = c_int()  # Number of materials hit in a given element
-        volumes = []
-        for i_element in range(self.n_elements):
-            while True:
-                try:
-                    _dll.openmc_mesh_material_volumes(
-                        self._index, n_samples, i_element, size, result, hits, prn_seed)
-                except AllocationError:
-                    # Increase size of result array and try again
-                    size *= 2
-                    result = (_MaterialVolume * size)()
-                else:
-                    # If no error, break out of loop
-                    break
-
-            volumes.append([
-                (Material(index=r.material), r.volume)
-                for r in result[:hits.value]
-            ])
-        return volumes
-
-    def material_volumes_raytrace(
             self,
             n_samples: int | tuple[int, int] = 10_000,
             max_materials: int = 4
@@ -301,7 +241,7 @@ def material_volumes_raytrace(
         # Run material volume calculation
         while True:
             try:
-                _dll.openmc_mesh_material_volumes_raytrace(
+                _dll.openmc_mesh_material_volumes(
                     self._index, ny, nz, max_materials, materials, volumes)
             except AllocationError:
                 # Increase size of result array and try again
diff --git a/openmc/mesh.py b/openmc/mesh.py
index a6082e60a84..35c322bb4f7 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -359,7 +359,7 @@ def material_volumes(
             kwargs.setdefault('output', False)
             openmc.lib.init(['-c'], **kwargs)
             mesh = openmc.lib.tallies[new_tally.id].filters[0].mesh
-            volumes = mesh.material_volumes_raytrace(n_samples, max_materials)
+            volumes = mesh.material_volumes(n_samples, max_materials)
             openmc.lib.finalize()
 
             # Restore original tallies
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 2ad4d804438..a8345624683 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -152,93 +152,7 @@ vector<double> Mesh::volumes() const
   return volumes;
 }
 
-int Mesh::material_volumes(
-  int n_sample, int bin, gsl::span<MaterialVolume> result, uint64_t* seed) const
-{
-  vector<int32_t> materials;
-  vector<int64_t> hits;
-
-#pragma omp parallel
-  {
-    vector<int32_t> local_materials;
-    vector<int64_t> local_hits;
-    GeometryState geom;
-
-#pragma omp for
-    for (int i = 0; i < n_sample; ++i) {
-      // Get seed for i-th sample
-      uint64_t seed_i = future_seed(3 * i, *seed);
-
-      // Sample position and set geometry state
-      geom.r() = this->sample_element(bin, &seed_i);
-      geom.u() = {1., 0., 0.};
-      geom.n_coord() = 1;
-
-      // If this location is not in the geometry at all, move on to next block
-      if (!exhaustive_find_cell(geom))
-        continue;
-
-      int i_material = geom.material();
-
-      // Check if this material was previously hit and if so, increment count
-      auto it =
-        std::find(local_materials.begin(), local_materials.end(), i_material);
-      if (it == local_materials.end()) {
-        local_materials.push_back(i_material);
-        local_hits.push_back(1);
-      } else {
-        local_hits[it - local_materials.begin()]++;
-      }
-    } // omp for
-
-    // Reduce index/hits lists from each thread into a single copy
-    reduce_indices_hits(local_materials, local_hits, materials, hits);
-  } // omp parallel
-
-  // Advance RNG seed
-  advance_prn_seed(3 * n_sample, seed);
-
-  // Make sure span passed in is large enough
-  if (hits.size() > result.size()) {
-    return -1;
-  }
-
-  // Convert hits to fractions
-  for (int i_mat = 0; i_mat < hits.size(); ++i_mat) {
-    double fraction = double(hits[i_mat]) / n_sample;
-    result[i_mat].material = materials[i_mat];
-    result[i_mat].volume = fraction * this->volume(bin);
-  }
-  return hits.size();
-}
-
-vector<Mesh::MaterialVolume> Mesh::material_volumes(
-  int n_sample, int bin, uint64_t* seed) const
-{
-  // Create result vector with space for 8 pairs
-  vector<Mesh::MaterialVolume> result;
-  result.reserve(8);
-
-  int size = -1;
-  while (true) {
-    // Get material volumes
-    size = this->material_volumes(
-      n_sample, bin, {result.data(), result.data() + result.capacity()}, seed);
-
-    // If capacity was sufficient, resize the vector and return
-    if (size >= 0) {
-      result.resize(size);
-      break;
-    }
-
-    // Otherwise, increase capacity of the vector
-    result.reserve(2 * result.capacity());
-  }
-
-  return result;
-}
-
-void Mesh::material_volumes_raytrace(
+void Mesh::material_volumes(
   int ny, int nz, int max_materials, int32_t* materials, double* volumes) const
 {
   // Create object for keeping track of materials/volumes
@@ -2082,32 +1996,14 @@ extern "C" int openmc_mesh_bounding_box(int32_t index, double* ll, double* ur)
   return 0;
 }
 
-extern "C" int openmc_mesh_material_volumes(int32_t index, int n_sample,
-  int bin, int result_size, void* result, int* hits, uint64_t* seed)
-{
-  auto result_ = reinterpret_cast<Mesh::MaterialVolume*>(result);
-  if (!result_) {
-    set_errmsg("Invalid result pointer passed to openmc_mesh_material_volumes");
-    return OPENMC_E_INVALID_ARGUMENT;
-  }
-
-  if (int err = check_mesh(index))
-    return err;
-
-  int n = model::meshes[index]->material_volumes(
-    n_sample, bin, {result_, result_ + result_size}, seed);
-  *hits = n;
-  return (n == -1) ? OPENMC_E_ALLOCATE : 0;
-}
-
-extern "C" int openmc_mesh_material_volumes_raytrace(int32_t index, int ny,
-  int nz, int max_mats, int32_t* materials, double* volumes)
+extern "C" int openmc_mesh_material_volumes(int32_t index, int ny, int nz,
+  int max_mats, int32_t* materials, double* volumes)
 {
   if (int err = check_mesh(index))
     return err;
 
   try {
-    model::meshes[index]->material_volumes_raytrace(
+    model::meshes[index]->material_volumes(
       ny, nz, max_mats, materials, volumes);
   } catch (const std::exception& e) {
     set_errmsg(e.what());
diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py
index b7c8761af34..8ab35335fc0 100644
--- a/tests/unit_tests/test_lib.py
+++ b/tests/unit_tests/test_lib.py
@@ -600,7 +600,7 @@ def test_regular_mesh(lib_init):
     mesh.dimension = (2, 2, 1)
     mesh.set_parameters(lower_left=(-0.63, -0.63, -0.5),
                         upper_right=(0.63, 0.63, 0.5))
-    vols = mesh.material_volumes_raytrace()
+    vols = mesh.material_volumes()
     assert vols.num_elements == 4
     for i in range(vols.num_elements):
         elem_vols = vols.by_element(i)
@@ -612,7 +612,7 @@ def test_regular_mesh(lib_init):
     mesh.set_parameters(lower_left=(-1.0, -1.0, -0.5),
                         upper_right=(1.0, 1.0, 0.5))
     with pytest.raises(exc.GeometryError, match="not fully contained"):
-        vols = mesh.material_volumes_raytrace()
+        vols = mesh.material_volumes()
 
 
 def test_regular_mesh_get_plot_bins(lib_init):
@@ -682,7 +682,7 @@ def test_rectilinear_mesh(lib_init):
     w = 1.26
     mesh.set_grid([-w/2, -w/4, w/2], [-w/2, -w/4, w/2], [-0.5, 0.5])
 
-    vols = mesh.material_volumes_raytrace()
+    vols = mesh.material_volumes()
     assert vols.num_elements == 4
     assert sum(f[1] for f in vols.by_element(0)) == pytest.approx(w/4 * w/4)
     assert sum(f[1] for f in vols.by_element(1)) == pytest.approx(w/4 * 3*w/4)
@@ -736,7 +736,7 @@ def test_cylindrical_mesh(lib_init):
     z_grid = (-0.5, 0.5)
     mesh.set_grid(r_grid, phi_grid, z_grid)
 
-    vols = mesh.material_volumes_raytrace()
+    vols = mesh.material_volumes()
     assert vols.num_elements == 6
     for i in range(0, 6, 2):
         assert sum(f[1] for f in vols.by_element(i)) == pytest.approx(pi * 0.25**2 / 3)
@@ -794,7 +794,7 @@ def test_spherical_mesh(lib_init):
     phi_grid = np.linspace(0., 2.0*pi, 4)
     mesh.set_grid(r_grid, theta_grid, phi_grid)
 
-    vols = mesh.material_volumes_raytrace()
+    vols = mesh.material_volumes()
     assert vols.num_elements == 12
     d_theta = theta_grid[1] - theta_grid[0]
     d_phi = phi_grid[1] - phi_grid[0]

From 812d9a8d3313ba51a8c8a03d2e578288eb45f274 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 28 Aug 2024 23:58:29 -0500
Subject: [PATCH 15/23] Clean up, add documentation

---
 docs/source/pythonapi/base.rst |  1 +
 include/openmc/mesh.h          | 99 +++++++++++++++-------------------
 openmc/lib/mesh.py             |  2 -
 openmc/mesh.py                 | 32 ++++++++++-
 src/mesh.cpp                   | 39 +++++++++++++-
 5 files changed, 112 insertions(+), 61 deletions(-)

diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst
index 5ae9f20edf4..4724544491a 100644
--- a/docs/source/pythonapi/base.rst
+++ b/docs/source/pythonapi/base.rst
@@ -150,6 +150,7 @@ Constructing Tallies
    openmc.CylindricalMesh
    openmc.SphericalMesh
    openmc.UnstructuredMesh
+   openmc.MeshMaterialVolumes
    openmc.Trigger
    openmc.TallyDerivative
    openmc.Tally
diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index e12edf2d7e3..d990209a2fb 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -69,70 +69,57 @@ extern const libMesh::Parallel::Communicator* libmesh_comm;
 } // namespace settings
 #endif
 
-class Mesh {
+//==============================================================================
+//! Helper class for keeping track of volume for each material in a mesh element
+//==============================================================================
+
+namespace detail {
+
+class MaterialVolumes {
 public:
-  // Types, aliases
-  struct MaterialVolume {
-    int32_t material; //!< material index
-    double volume;    //!< volume in [cm^3]
-  };
+  MaterialVolumes(int32_t* mats, double* vols, int max_materials)
+    : materials_(mats), volumes_(vols), n_mats_(max_materials)
+  {}
 
-  using shape_type = std::vector<size_t>;
+  //! Add volume for a given material in a mesh element
+  //
+  //! \param[in] index_elem Index of the mesh element
+  //! \param[in] index_material Index of the material
+  //! \param[in] volume Volume to add
+  void add_volume(int index_elem, int index_material, double volume);
+
+  // Accessors
+  int32_t& materials(int i, int j) { return materials_[i * n_mats_ + j]; }
+  const int32_t& materials(int i, int j) const
+  {
+    return materials_[i * n_mats_ + j];
+  }
 
-  class MaterialVolumeRT {
-  public:
-    MaterialVolumeRT(int32_t* mats, double* vols, int max_materials)
-      : materials_(mats), volumes_(vols), n_mats_(max_materials)
-    {}
+  double& volumes(int i, int j) { return volumes_[i * n_mats_ + j]; }
+  const double& volumes(int i, int j) const
+  {
+    return volumes_[i * n_mats_ + j];
+  }
 
-    void add_volume(int index_elem, int index_material, double volume)
-    {
-      int i;
-#pragma omp critical(MeshMatVol)
-      for (i = 0; i < n_mats_; ++i) {
-        // Check whether material already is present
-        if (materials(index_elem, i) == index_material)
-          break;
-
-        // If not already present, check for unused position (-2)
-        if (materials(index_elem, i) == -2) {
-          materials(index_elem, i) = index_material;
-          break;
-        }
-      }
-
-      // If maximum number of materials exceeded, set a flag that can be checked
-      // later
-      if (i >= n_mats_) {
-        too_many_mats_ = true;
-        return;
-      }
-
-      // Accumulate volume
-#pragma omp atomic
-      volumes(index_elem, i) += volume;
-    }
+  bool too_many_mats() const { return too_many_mats_; }
 
-    int32_t& materials(int i, int j) { return materials_[i * n_mats_ + j]; }
-    const int32_t& materials(int i, int j) const
-    {
-      return materials_[i * n_mats_ + j];
-    }
+private:
+  int32_t* materials_; //!< material index (bins, max_mats)
+  double* volumes_;    //!< volume in [cm^3] (bins, max_mats)
+  int n_mats_;         //!< Maximum number of materials in a single mesh element
+  bool too_many_mats_ = false; //!< Whether the maximum number of materials has
+                               //!< been exceeded
+};
 
-    double& volumes(int i, int j) { return volumes_[i * n_mats_ + j]; }
-    const double& volumes(int i, int j) const
-    {
-      return volumes_[i * n_mats_ + j];
-    }
+} // namespace detail
 
-    bool too_many_mats() const { return too_many_mats_; }
+//==============================================================================
+//! Base mesh class
+//==============================================================================
 
-  private:
-    int32_t* materials_; //!< material index (bins, max_mats)
-    double* volumes_;    //!< volume in [cm^3] (bins, max_mats)
-    int n_mats_;
-    bool too_many_mats_ = false;
-  };
+class Mesh {
+public:
+  // Types, aliases
 
   // Constructors and destructor
   Mesh() = default;
diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 4f4232832f8..19dda46be32 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -118,8 +118,6 @@ class _MaterialVolume(Structure):
 _dll.openmc_spherical_mesh_set_grid.errcheck = _error_handler
 
 
-
-
 class Mesh(_FortranObjectWithID):
     """Base class to represent mesh objects
 
diff --git a/openmc/mesh.py b/openmc/mesh.py
index 35c322bb4f7..7c17227f492 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -23,13 +23,41 @@
 class MeshMaterialVolumes(Mapping):
     """Results from a material volume in mesh calculation.
 
+    This class provides multiple ways of accessing information about material
+    volumes in individual mesh elements. First, the class behaves like a
+    dictionary that maps material IDs to an array of volumes equal in size to
+    the number of mesh elements. Second, the class provides a :meth:`by_element`
+    method that gives all the material volumes for a specific mesh element.
+
     Parameters
     ----------
-    materials : np.ndarray
+    materials : numpy.ndarray
         Array of shape (elements, max_materials) storing material IDs
-    volumes : np.ndarray
+    volumes : numpy.ndarray
         Array of shape (elements, max_materials) storing material volumes
 
+    See Also
+    --------
+    openmc.MeshBase.material_volumes
+
+    Examples
+    --------
+    If you want to get the volume of a specific material in every mesh element,
+    index the object with the material ID:
+
+    >>> volumes = mesh.material_volumes(...)
+    >>> volumes
+    {1: <32121 nonzero volumes>
+     2: <338186 nonzero volumes>
+     3: <49120 nonzero volumes>}
+
+    If you want the volume of all materials in a specific mesh element, use the
+    :meth:`by_element` method:
+
+    >>> volumes = mesh.material_volumes(...)
+    >>> volumes.by_element(42)
+    [(2, 31.87963824195591), (1, 6.129949130817542)]
+
     """
     def __init__(self, materials: np.ndarray, volumes):
         self._materials = materials
diff --git a/src/mesh.cpp b/src/mesh.cpp
index a8345624683..7649451cc3f 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -107,6 +107,43 @@ inline bool check_intersection_point(double x1, double x0, double y1, double y0,
 // Mesh implementation
 //==============================================================================
 
+namespace detail {
+
+void MaterialVolumes::add_volume(
+  int index_elem, int index_material, double volume)
+{
+  int i;
+#pragma omp critical(MeshMatVol)
+  for (i = 0; i < n_mats_; ++i) {
+    // Check whether material already is present
+    if (this->materials(index_elem, i) == index_material)
+      break;
+
+    // If not already present, check for unused position (-2)
+    if (this->materials(index_elem, i) == -2) {
+      this->materials(index_elem, i) = index_material;
+      break;
+    }
+  }
+
+  // If maximum number of materials exceeded, set a flag that can be checked
+  // later
+  if (i >= n_mats_) {
+    too_many_mats_ = true;
+    return;
+  }
+
+  // Accumulate volume
+#pragma omp atomic
+  this->volumes(index_elem, i) += volume;
+}
+
+} // namespace detail
+
+//==============================================================================
+// Mesh implementation
+//==============================================================================
+
 Mesh::Mesh(pugi::xml_node node)
 {
   // Read mesh id
@@ -156,7 +193,7 @@ void Mesh::material_volumes(
   int ny, int nz, int max_materials, int32_t* materials, double* volumes) const
 {
   // Create object for keeping track of materials/volumes
-  MaterialVolumeRT result(materials, volumes, max_materials);
+  detail::MaterialVolumes result(materials, volumes, max_materials);
 
   // Determine bounding box
   auto bbox = this->bounding_box();

From 0fcddbffaa5319ab38791d9c824201a81b8fbcaa Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 2 Sep 2024 12:48:11 -0500
Subject: [PATCH 16/23] Add test for material_volumes

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

diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index ed08be816f5..7356ca99a45 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -402,7 +402,7 @@ def test_mesh_get_homogenized_materials():
     mesh.lower_left = (-1., -1., -1.)
     mesh.upper_right = (1., 1., 1.)
     mesh.dimension = (3, 1, 1)
-    m1, m2, m3 = mesh.get_homogenized_materials(model, n_samples=1_000_000)
+    m1, m2, m3 = mesh.get_homogenized_materials(model, n_samples=10_000)
 
     # Left mesh element should be only Fe56
     assert m1.get_mass_density('Fe56') == pytest.approx(5.0)
@@ -418,7 +418,7 @@ def test_mesh_get_homogenized_materials():
     mesh_void.lower_left = (0.5, 0.5, -1.)
     mesh_void.upper_right = (1.5, 1.5, 1.)
     mesh_void.dimension = (1, 1, 1)
-    m4, = mesh_void.get_homogenized_materials(model, n_samples=1_000_000)
+    m4, = mesh_void.get_homogenized_materials(model, n_samples=10_000)
 
     # Mesh element that overlaps void should have half density
     assert m4.get_mass_density('H1') == pytest.approx(0.5, rel=1e-2)
@@ -428,3 +428,51 @@ def test_mesh_get_homogenized_materials():
     m5, = mesh_void.get_homogenized_materials(
         model, n_samples=1000, include_void=False)
     assert m5.get_mass_density('H1') == pytest.approx(1.0)
+
+
+def test_mesh_material_volumes():
+    """Test the material_volumes method"""
+    # Model with three materials separated by planes x=0 and z=0
+    mats = []
+    for i in range(3):
+        mat = openmc.Material()
+        mat.add_nuclide('H1', 1.0)
+        mat.set_density('g/cm3', float(i + 1))
+        mats.append(mat)
+
+    sph = openmc.Sphere(r=25.0, boundary_type='vacuum')
+    x0 = openmc.XPlane(0.0)
+    z0 = openmc.ZPlane(0.0)
+    cell1 = openmc.Cell(fill=mats[0], region=-sph & +x0 & +z0)
+    cell2 = openmc.Cell(fill=mats[1], region=-sph & -x0 & +z0)
+    cell3 = openmc.Cell(fill=mats[2], region=-sph & -z0)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell1, cell2, cell3])
+
+    # Test material volumes on a regular mesh
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = (-1., -1., -1.)
+    mesh.upper_right = (1., 1., 1.)
+    mesh.dimension = (2, 2, 2)
+    volumes = mesh.material_volumes(model)
+    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.)]
+
+    # Test material volumes on a cylindrical mesh
+    cyl_mesh = openmc.CylindricalMesh([0., 1.], [-1., 0., 1.,], [0.0, pi/4, 3*pi/4, 5*pi/4, 7*pi/4, 2*pi])
+    volumes = cyl_mesh.material_volumes(model)
+    np.testing.assert_almost_equal(volumes[mats[0].id], [
+        0., 0., 0., 0., 0.,
+        pi/8, pi/8, 0., pi/8, pi/8
+    ])
+    np.testing.assert_almost_equal(volumes[mats[1].id], [
+        0., 0., 0., 0., 0.,
+        0., pi/8, pi/4, pi/8, 0.
+    ])
+    np.testing.assert_almost_equal(volumes[mats[2].id], [
+        pi/8, pi/4, pi/4, pi/4, pi/8,
+        0., 0., 0., 0., 0.
+    ])

From 1819a089bf2d2b105cc6cfd3d5d6adfe78518a18 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 2 Sep 2024 15:19:56 -0500
Subject: [PATCH 17/23] Add test for save/from_npz methods

---
 openmc/lib/mesh.py            |  8 +++++++-
 openmc/mesh.py                |  8 ++++++--
 tests/unit_tests/test_mesh.py | 29 ++++++++++++++++++++++++++++-
 3 files changed, 41 insertions(+), 4 deletions(-)

diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 19dda46be32..196a1ddf64f 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -204,6 +204,12 @@ def material_volumes(
         side to count of the estimated volume of each material in all mesh
         elements.
 
+        .. versionadded:: 0.15.0
+
+        .. versionchanged:: 0.15.1
+            Material volumes are now determined by raytracing rather than by
+            point sampling.
+
         Parameters
         ----------
         n_samples : int or 2-tuple of int
@@ -212,7 +218,7 @@ def material_volumes(
             2-tuple, it is interpreted as the number of rays in the y and z
             dimensions.
         max_materials : int, optional
-            Maximum number of materials in any given mesh element.
+            Estimated maximum number of materials in any given mesh element.
 
         Returns
         -------
diff --git a/openmc/mesh.py b/openmc/mesh.py
index 7c17227f492..0d09bbf1f86 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -29,6 +29,8 @@ class MeshMaterialVolumes(Mapping):
     the number of mesh elements. Second, the class provides a :meth:`by_element`
     method that gives all the material volumes for a specific mesh element.
 
+    .. versionadded:: 0.15.1
+
     Parameters
     ----------
     materials : numpy.ndarray
@@ -59,7 +61,7 @@ class MeshMaterialVolumes(Mapping):
     [(2, 31.87963824195591), (1, 6.129949130817542)]
 
     """
-    def __init__(self, materials: np.ndarray, volumes):
+    def __init__(self, materials: np.ndarray, volumes: np.ndarray):
         self._materials = materials
         self._volumes = volumes
 
@@ -349,6 +351,8 @@ def material_volumes(
         side to count of the estimated volume of each material in all mesh
         elements.
 
+        .. versionadded:: 0.15.1
+
         Parameters
         ----------
         model : openmc.Model
@@ -359,7 +363,7 @@ def material_volumes(
             2-tuple, it is interpreted as the number of rays in the y and z
             dimensions.
         max_materials : int, optional
-            Maximum number of materials in any given mesh element.
+            Estimated maximum number of materials in any given mesh element.
         **kwargs : dict
             Keyword arguments passed to :func:`openmc.lib.init`
 
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 7356ca99a45..6e95bdc769f 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -1,4 +1,5 @@
 from math import pi
+from tempfile import TemporaryDirectory
 
 import numpy as np
 import pytest
@@ -462,7 +463,8 @@ def test_mesh_material_volumes():
     assert volumes.by_element(0) == [(mats[2].id, 1.)]
 
     # Test material volumes on a cylindrical mesh
-    cyl_mesh = openmc.CylindricalMesh([0., 1.], [-1., 0., 1.,], [0.0, pi/4, 3*pi/4, 5*pi/4, 7*pi/4, 2*pi])
+    cyl_mesh = openmc.CylindricalMesh(
+        [0., 1.], [-1., 0., 1.,], [0.0, pi/4, 3*pi/4, 5*pi/4, 7*pi/4, 2*pi])
     volumes = cyl_mesh.material_volumes(model)
     np.testing.assert_almost_equal(volumes[mats[0].id], [
         0., 0., 0., 0., 0.,
@@ -476,3 +478,28 @@ def test_mesh_material_volumes():
         pi/8, pi/4, pi/4, pi/4, pi/8,
         0., 0., 0., 0., 0.
     ])
+
+
+def test_mesh_material_volumes_serialize():
+    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]
+    ])
+    volumes = openmc.MeshMaterialVolumes(materials, volumes)
+    with TemporaryDirectory() as tmpdir:
+        path = f'{tmpdir}/volumes.npz'
+        volumes.save(path)
+        new_volumes = openmc.MeshMaterialVolumes.from_npz(path)
+
+    assert new_volumes.by_element(0) == [(1, 0.5), (None, 0.5)]
+    assert new_volumes.by_element(1) == [(None, 1.0)]
+    assert new_volumes.by_element(2) == [(2, 0.5), (1, 0.5)]
+    assert new_volumes.by_element(3) == [(2, 1.0)]

From 762d1f88499d5543aa2abbbc4f608135fa98c7bc Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 13 Sep 2024 23:02:49 -0500
Subject: [PATCH 18/23] Extend material volume calculation to use rays in all
 three directions

---
 include/openmc/bounding_box.h |   4 +
 include/openmc/capi.h         |   4 +-
 include/openmc/mesh.h         |   5 +-
 openmc/lib/mesh.py            |  55 ++++-----
 openmc/mesh.py                |  29 ++---
 src/mesh.cpp                  | 205 +++++++++++++++++++---------------
 tests/unit_tests/test_mesh.py |  23 ++--
 7 files changed, 179 insertions(+), 146 deletions(-)

diff --git a/include/openmc/bounding_box.h b/include/openmc/bounding_box.h
index 40f603583b4..d02d92cb41e 100644
--- a/include/openmc/bounding_box.h
+++ b/include/openmc/bounding_box.h
@@ -4,6 +4,7 @@
 #include <algorithm> // for min, max
 
 #include "openmc/constants.h"
+#include "openmc/position.h"
 
 namespace openmc {
 
@@ -54,6 +55,9 @@ struct BoundingBox {
     zmax = std::max(zmax, other.zmax);
     return *this;
   }
+
+  inline Position min() const { return {xmin, ymin, zmin}; }
+  inline Position max() const { return {xmax, ymax, zmax}; }
 };
 
 } // namespace openmc
diff --git a/include/openmc/capi.h b/include/openmc/capi.h
index 223a6b1877e..382718a47e3 100644
--- a/include/openmc/capi.h
+++ b/include/openmc/capi.h
@@ -107,8 +107,8 @@ int openmc_mesh_get_id(int32_t index, int32_t* id);
 int openmc_mesh_set_id(int32_t index, int32_t id);
 int openmc_mesh_get_n_elements(int32_t index, size_t* n);
 int openmc_mesh_get_volumes(int32_t index, double* volumes);
-int openmc_mesh_material_volumes(int32_t index, int ny, int nz, int max_mats,
-  int32_t* materials, double* volumes);
+int openmc_mesh_material_volumes(int32_t index, int nx, int ny, int nz,
+  int max_mats, int32_t* materials, double* volumes);
 int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
 int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
 int openmc_new_filter(const char* type, int32_t* index);
diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index d990209a2fb..4d28555b88e 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -213,14 +213,15 @@ class Mesh {
 
   //! Determine volume of materials within each mesh element
   //
+  //! \param[in] nx Number of samples in x direction
   //! \param[in] ny Number of samples in y direction
   //! \param[in] nz Number of samples in z direction
   //! \param[in] max_materials Maximum number of materials in a single mesh
   //!                          element
   //! \param[inout] materials Array storing material indices
   //! \param[inout] volumes Array storing volumes
-  void material_volumes(int ny, int nz, int max_materials, int32_t* materials,
-    double* volumes) const;
+  void material_volumes(int nx, int ny, int nz, int max_materials,
+    int32_t* materials, double* volumes) const;
 
   //! Determine bounding box of mesh
   //
diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 196a1ddf64f..28339490d55 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -1,8 +1,7 @@
 from collections.abc import Mapping, Sequence
-from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, Structure,
-                    create_string_buffer, c_uint64, c_size_t)
+from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER,
+                    create_string_buffer, c_size_t)
 from math import sqrt
-from random import getrandbits
 import sys
 from weakref import WeakValueDictionary
 
@@ -13,10 +12,8 @@
 from . import _dll
 from .core import _FortranObjectWithID
 from .error import _error_handler
-from .material import Material
 from .plot import _Position
 from ..bounding_box import BoundingBox
-from ..checkvalue import PathLike
 from ..mesh import MeshMaterialVolumes
 
 __all__ = [
@@ -25,13 +22,6 @@
 ]
 
 
-class _MaterialVolume(Structure):
-    _fields_ = [
-        ("material", c_int32),
-        ("volume", c_double)
-    ]
-
-
 arr_2d_int32 = np.ctypeslib.ndpointer(dtype=np.int32, ndim=2, flags='CONTIGUOUS')
 arr_2d_double = np.ctypeslib.ndpointer(dtype=np.double, ndim=2, flags='CONTIGUOUS')
 
@@ -57,7 +47,7 @@ class _MaterialVolume(Structure):
 _dll.openmc_mesh_bounding_box.restype = c_int
 _dll.openmc_mesh_bounding_box.errcheck = _error_handler
 _dll.openmc_mesh_material_volumes.argtypes = [
-    c_int32, c_int, c_int, c_int, arr_2d_int32, arr_2d_double]
+    c_int32, c_int, c_int, c_int, c_int, arr_2d_int32, arr_2d_double]
 _dll.openmc_mesh_material_volumes.restype = c_int
 _dll.openmc_mesh_material_volumes.errcheck = _error_handler
 _dll.openmc_mesh_get_plot_bins.argtypes = [
@@ -195,14 +185,15 @@ def bounding_box(self) -> BoundingBox:
 
     def material_volumes(
             self,
-            n_samples: int | tuple[int, int] = 10_000,
+            n_samples: int | tuple[int, int, int] = 10_000,
             max_materials: int = 4
     ) -> MeshMaterialVolumes:
         """Determine volume of materials in each mesh element.
 
-        This method works by raytracing repeatedly through the mesh from one
-        side to count of the estimated volume of each material in all mesh
-        elements.
+        This method works by raytracing repeatedly through the mesh to count the
+        estimated volume of each material in all mesh elements. Three sets of
+        rays are used: one set parallel to the x-axis, one parallel to the
+        y-axis, and one parallel to the z-axis.
 
         .. versionadded:: 0.15.0
 
@@ -212,11 +203,11 @@ def material_volumes(
 
         Parameters
         ----------
-        n_samples : int or 2-tuple of int
-            Total number of rays to sample. The rays start on an x plane and are
-            evenly distributed over the y and z dimensions. When specified as a
-            2-tuple, it is interpreted as the number of rays in the y and z
-            dimensions.
+        n_samples : int or 3-tuple of int
+            Total number of rays to sample. The number of rays in each direction
+            is determined by the aspect ratio of the mesh bounding box. When
+            specified as a 3-tuple, it is interpreted as the number of rays in
+            the x, y, and z dimensions.
         max_materials : int, optional
             Estimated maximum number of materials in any given mesh element.
 
@@ -226,16 +217,18 @@ def material_volumes(
         equal in size to the number of mesh elements.
 
         """
-        # Determine number of rays in y/z directions
         if isinstance(n_samples, int):
-            ll, ur = self.bounding_box
-            width_y = ur[1] - ll[1]
-            width_z = ur[2] - ll[2]
-            aspect_ratio = width_y / width_z
-            ny = int(sqrt(n_samples * aspect_ratio))
-            nz = int(sqrt(n_samples / aspect_ratio))
+            # Determine number of rays in each direction based on aspect ratios
+            # and using the relation (nx*ny + ny*nz + nx*nz) = n_samples
+            width_x, width_y, width_z = self.bounding_box.width
+            ax = width_x / width_z
+            ay = width_y / width_z
+            f = sqrt(n_samples/(ax*ay + ax + ay))
+            nx = round(f * ax)
+            ny = round(f * ay)
+            nz = round(f)
         else:
-            ny, nz = n_samples
+            nx, ny, nz = n_samples
 
         # Preallocate arrays for material indices and volumes
         n = self.n_elements
@@ -246,7 +239,7 @@ def material_volumes(
         while True:
             try:
                 _dll.openmc_mesh_material_volumes(
-                    self._index, ny, nz, max_materials, materials, volumes)
+                    self._index, nx, ny, nz, max_materials, materials, volumes)
             except AllocationError:
                 # Increase size of result array and try again
                 max_materials *= 2
diff --git a/openmc/mesh.py b/openmc/mesh.py
index 0d09bbf1f86..eda32984451 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -261,7 +261,7 @@ def from_xml_element(cls, elem: ET.Element):
     def get_homogenized_materials(
             self,
             model: openmc.Model,
-            n_samples: int | tuple[int, int] = 10_000,
+            n_samples: int | tuple[int, int, int] = 10_000,
             include_void: bool = True,
             **kwargs
     ) -> list[openmc.Material]:
@@ -275,10 +275,10 @@ def get_homogenized_materials(
             Model containing materials to be homogenized and the associated
             geometry.
         n_samples : int or 2-tuple of int
-            Total number of rays to sample. The rays start on an x plane and are
-            evenly distributed over the y and z dimensions. When specified as a
-            2-tuple, it is interpreted as the number of rays in the y and z
-            dimensions.
+            Total number of rays to sample. The number of rays in each direction
+            is determined by the aspect ratio of the mesh bounding box. When
+            specified as a 3-tuple, it is interpreted as the number of rays in
+            the x, y, and z dimensions.
         include_void : bool, optional
             Whether homogenization should include voids.
         **kwargs
@@ -341,15 +341,16 @@ def get_homogenized_materials(
     def material_volumes(
             self,
             model: openmc.Model,
-            n_samples: int | tuple[int, int] = 10_000,
+            n_samples: int | tuple[int, int, int] = 10_000,
             max_materials: int = 4,
             **kwargs
     ) -> MeshMaterialVolumes:
         """Determine volume of materials in each mesh element.
 
-        This method works by raytracing repeatedly through the mesh from one
-        side to count of the estimated volume of each material in all mesh
-        elements.
+        This method works by raytracing repeatedly through the mesh to count the
+        estimated volume of each material in all mesh elements. Three sets of
+        rays are used: one set parallel to the x-axis, one parallel to the
+        y-axis, and one parallel to the z-axis.
 
         .. versionadded:: 0.15.1
 
@@ -357,11 +358,11 @@ def material_volumes(
         ----------
         model : openmc.Model
             Model containing materials.
-        n_samples : int or 2-tuple of int
-            Total number of rays to sample. The rays start on an x plane and are
-            evenly distributed over the y and z dimensions. When specified as a
-            2-tuple, it is interpreted as the number of rays in the y and z
-            dimensions.
+        n_samples : int or 3-tuple of int
+            Total number of rays to sample. The number of rays in each direction
+            is determined by the aspect ratio of the mesh bounding box. When
+            specified as a 3-tuple, it is interpreted as the number of rays in
+            the x, y, and z dimensions.
         max_materials : int, optional
             Estimated maximum number of materials in any given mesh element.
         **kwargs : dict
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 7649451cc3f..53814c6c19c 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -189,8 +189,8 @@ vector<double> Mesh::volumes() const
   return volumes;
 }
 
-void Mesh::material_volumes(
-  int ny, int nz, int max_materials, int32_t* materials, double* volumes) const
+void Mesh::material_volumes(int nx, int ny, int nz, int max_materials,
+  int32_t* materials, double* volumes) const
 {
   // Create object for keeping track of materials/volumes
   detail::MaterialVolumes result(materials, volumes, max_materials);
@@ -198,9 +198,11 @@ void Mesh::material_volumes(
   // Determine bounding box
   auto bbox = this->bounding_box();
 
+  std::array<int, 3> n_rays = {nx, ny, nz};
+
   // Loop over rays on left face of bounding box
-  double dy = (bbox.ymax - bbox.ymin) / ny;
-  double dz = (bbox.zmax - bbox.zmin) / nz;
+  Position width((bbox.xmax - bbox.xmin) / nx, (bbox.ymax - bbox.ymin) / ny,
+    (bbox.zmax - bbox.zmin) / nz);
 
   // Set flag for mesh being contained within model
   bool out_of_model = false;
@@ -213,88 +215,105 @@ void Mesh::material_volumes(
     Particle p;
 
     SourceSite site;
-    site.r = {bbox.xmin, 0.0, 0.0};
-    site.u = {1.0, 0.0, 0.0};
     site.E = 1.0;
     site.particle = ParticleType::neutron;
 
-#pragma omp for
-    for (int iy = 0; iy < ny; ++iy) {
-      site.r.y = bbox.ymin + (iy + 0.5) * dy;
-      for (int iz = 0; iz < nz; ++iz) {
-        site.r.z = bbox.zmin + (iz + 0.5) * dz;
-
-        p.from_source(&site);
-
-        // Determine particle's location
-        if (!exhaustive_find_cell(p)) {
-          out_of_model = true;
-          continue;
-        }
-
-        // Set birth cell attribute
-        if (p.cell_born() == C_NONE)
-          p.cell_born() = p.lowest_coord().cell;
-
-        // Initialize last cells from current cell
-        for (int j = 0; j < p.n_coord(); ++j) {
-          p.cell_last(j) = p.coord(j).cell;
-        }
-        p.n_coord_last() = p.n_coord();
-
-        while (true) {
-          // Ray trace from r_start to r_end
-          Position r0 = p.r();
-          double max_distance = bbox.xmax - r0.x;
-
-          // Find the distance to the nearest boundary
-          BoundaryInfo boundary = distance_to_boundary(p);
-
-          // Advance particle forward
-          double distance = std::min(boundary.distance, max_distance);
-          p.move_distance(distance);
+    for (int axis = 0; axis < 3; ++axis) {
+      // Set starting position and direction
+      site.r = {0.0, 0.0, 0.0};
+      site.r[axis] = bbox.min()[axis];
+      site.u = {0.0, 0.0, 0.0};
+      site.u[axis] = 1.0;
+
+      // Determine width of rays and number of rays in other directions
+      int ax1 = (axis + 1) % 3;
+      int ax2 = (axis + 2) % 3;
+      double min1 = bbox.min()[ax1];
+      double min2 = bbox.min()[ax2];
+      double d1 = width[ax1];
+      double d2 = width[ax2];
+      int n1 = n_rays[ax1];
+      int n2 = n_rays[ax2];
 
-          // Determine what mesh elements were crossed by particle
-          bins.clear();
-          length_fractions.clear();
-          this->bins_crossed(r0, p.r(), p.u(), bins, length_fractions);
+#pragma omp for
+      for (int i1 = 0; i1 < n1; ++i1) {
+        site.r[ax1] = min1 + (i1 + 0.5) * d1;
+        for (int i2 = 0; i2 < n2; ++i2) {
+          site.r[ax2] = min2 + (i2 + 0.5) * d2;
 
-          // Add volumes to any mesh elements that were crossed
-          int i_material = p.material();
-          if (i_material != C_NONE) {
-            i_material = model::materials[i_material]->id();
-          }
-          for (int i_bin = 0; i_bin < bins.size(); i_bin++) {
-            int mesh_index = bins[i_bin];
-            double length = distance * length_fractions[i_bin];
+          p.from_source(&site);
 
-            // Add volume to result
-            result.add_volume(mesh_index, i_material, length * dy * dz);
+          // Determine particle's location
+          if (!exhaustive_find_cell(p)) {
+            out_of_model = true;
+            continue;
           }
 
-          if (distance == max_distance)
-            break;
+          // Set birth cell attribute
+          if (p.cell_born() == C_NONE)
+            p.cell_born() = p.lowest_coord().cell;
 
+          // Initialize last cells from current cell
           for (int j = 0; j < p.n_coord(); ++j) {
             p.cell_last(j) = p.coord(j).cell;
           }
           p.n_coord_last() = p.n_coord();
 
-          // Set surface that particle is on and adjust coordinate levels
-          p.surface() = boundary.surface_index;
-          p.n_coord() = boundary.coord_level;
-
-          if (boundary.lattice_translation[0] != 0 ||
-              boundary.lattice_translation[1] != 0 ||
-              boundary.lattice_translation[2] != 0) {
-            // Particle crosses lattice boundary
-
-            cross_lattice(p, boundary);
-          } else {
-            // Particle crosses surface
-            // TODO: off-by-one
-            const auto& surf {model::surfaces[std::abs(p.surface()) - 1].get()};
-            p.cross_surface(*surf);
+          while (true) {
+            // Ray trace from r_start to r_end
+            Position r0 = p.r();
+            double max_distance = bbox.max()[axis] - r0[axis];
+
+            // Find the distance to the nearest boundary
+            BoundaryInfo boundary = distance_to_boundary(p);
+
+            // Advance particle forward
+            double distance = std::min(boundary.distance, max_distance);
+            p.move_distance(distance);
+
+            // Determine what mesh elements were crossed by particle
+            bins.clear();
+            length_fractions.clear();
+            this->bins_crossed(r0, p.r(), p.u(), bins, length_fractions);
+
+            // Add volumes to any mesh elements that were crossed
+            int i_material = p.material();
+            if (i_material != C_NONE) {
+              i_material = model::materials[i_material]->id();
+            }
+            for (int i_bin = 0; i_bin < bins.size(); i_bin++) {
+              int mesh_index = bins[i_bin];
+              double length = distance * length_fractions[i_bin];
+
+              // Add volume to result
+              result.add_volume(mesh_index, i_material, length * d1 * d2);
+            }
+
+            if (distance == max_distance)
+              break;
+
+            for (int j = 0; j < p.n_coord(); ++j) {
+              p.cell_last(j) = p.coord(j).cell;
+            }
+            p.n_coord_last() = p.n_coord();
+
+            // Set surface that particle is on and adjust coordinate levels
+            p.surface() = boundary.surface_index;
+            p.n_coord() = boundary.coord_level;
+
+            if (boundary.lattice_translation[0] != 0 ||
+                boundary.lattice_translation[1] != 0 ||
+                boundary.lattice_translation[2] != 0) {
+              // Particle crosses lattice boundary
+
+              cross_lattice(p, boundary);
+            } else {
+              // Particle crosses surface
+              // TODO: off-by-one
+              const auto& surf {
+                model::surfaces[std::abs(p.surface()) - 1].get()};
+              p.cross_surface(*surf);
+            }
           }
         }
       }
@@ -686,8 +705,8 @@ xt::xtensor<double, 1> StructuredMesh::count_sites(
 }
 
 // raytrace through the mesh. The template class T will do the tallying.
-// A modern optimizing compiler can recognize the noop method of T and eleminate
-// that call entirely.
+// A modern optimizing compiler can recognize the noop method of T and
+// eleminate that call entirely.
 template<class T>
 void StructuredMesh::raytrace_mesh(
   Position r0, Position r1, const Direction& u, T tally) const
@@ -755,7 +774,8 @@ void StructuredMesh::raytrace_mesh(
       if (traveled_distance >= total_distance)
         return;
 
-      // If we have not reached r1, we have hit a surface. Tally outward current
+      // If we have not reached r1, we have hit a surface. Tally outward
+      // current
       tally.surface(ijk, k, distances[k].max_surface, false);
 
       // Update cell and calculate distance to next surface in k-direction.
@@ -767,15 +787,16 @@ void StructuredMesh::raytrace_mesh(
       // Check if we have left the interior of the mesh
       in_mesh = ((ijk[k] >= 1) && (ijk[k] <= shape_[k]));
 
-      // If we are still inside the mesh, tally inward current for the next cell
+      // If we are still inside the mesh, tally inward current for the next
+      // cell
       if (in_mesh)
         tally.surface(ijk, k, !distances[k].max_surface, true);
 
     } else { // not inside mesh
 
-      // For all directions outside the mesh, find the distance that we need to
-      // travel to reach the next surface. Use the largest distance, as only
-      // this will cross all outer surfaces.
+      // For all directions outside the mesh, find the distance that we need
+      // to travel to reach the next surface. Use the largest distance, as
+      // only this will cross all outer surfaces.
       int k_max {0};
       for (int k = 0; k < n; ++k) {
         if ((ijk[k] < 1 || ijk[k] > shape_[k]) &&
@@ -789,7 +810,8 @@ void StructuredMesh::raytrace_mesh(
       if (traveled_distance >= total_distance)
         return;
 
-      // Calculate the new cell index and update all distances to next surfaces.
+      // Calculate the new cell index and update all distances to next
+      // surfaces.
       ijk = get_indices(r0 + (traveled_distance + TINY_BIT) * u, in_mesh);
       for (int k = 0; k < n; ++k) {
         distances[k] =
@@ -1666,7 +1688,8 @@ double SphericalMesh::find_theta_crossing(
   const double b = r.dot(u) * cos_t_2 - r.z * u.z;
   const double c = r.dot(r) * cos_t_2 - r.z * r.z;
 
-  // if factor of s^2 is zero, direction of flight is parallel to theta surface
+  // if factor of s^2 is zero, direction of flight is parallel to theta
+  // surface
   if (std::abs(a) < FP_PRECISION) {
     // if b vanishes, direction of flight is within theta surface and crossing
     // is not possible
@@ -1674,7 +1697,8 @@ double SphericalMesh::find_theta_crossing(
       return INFTY;
 
     const double s = -0.5 * c / b;
-    // Check if solution is in positive direction of flight and has correct sign
+    // Check if solution is in positive direction of flight and has correct
+    // sign
     if ((s > l) && (std::signbit(r.z + s * u.z) == sgn))
       return s;
 
@@ -2033,15 +2057,15 @@ extern "C" int openmc_mesh_bounding_box(int32_t index, double* ll, double* ur)
   return 0;
 }
 
-extern "C" int openmc_mesh_material_volumes(int32_t index, int ny, int nz,
-  int max_mats, int32_t* materials, double* volumes)
+extern "C" int openmc_mesh_material_volumes(int32_t index, int nx, int ny,
+  int nz, int max_mats, int32_t* materials, double* volumes)
 {
   if (int err = check_mesh(index))
     return err;
 
   try {
     model::meshes[index]->material_volumes(
-      ny, nz, max_mats, materials, volumes);
+      nx, ny, nz, max_mats, materials, volumes);
   } catch (const std::exception& e) {
     set_errmsg(e.what());
     if (starts_with(e.what(), "Mesh")) {
@@ -2486,7 +2510,8 @@ void MOABMesh::intersect_track(const moab::CartVect& start,
   moab::ErrorCode rval;
   vector<moab::EntityHandle> tris;
   // get all intersections with triangles in the tet mesh
-  // (distances are relative to the start point, not the previous intersection)
+  // (distances are relative to the start point, not the previous
+  // intersection)
   rval = kdtree_->ray_intersect_triangles(kdtree_root_, FP_COINCIDENT,
     dir.array(), start.array(), tris, hits, 0, track_len);
   if (rval != moab::MB_SUCCESS) {
@@ -3050,8 +3075,8 @@ void LibMesh::set_mesh_pointer_from_filename(const std::string& filename)
 void LibMesh::initialize()
 {
   if (!settings::libmesh_comm) {
-    fatal_error(
-      "Attempting to use an unstructured mesh without a libMesh communicator.");
+    fatal_error("Attempting to use an unstructured mesh without a libMesh "
+                "communicator.");
   }
 
   // assuming that unstructured meshes used in OpenMC are 3D
@@ -3304,8 +3329,8 @@ void read_meshes(pugi::xml_node root)
     // Check to make sure multiple meshes in the same file don't share IDs
     int id = std::stoi(get_node_value(node, "id"));
     if (contains(mesh_ids, id)) {
-      fatal_error(fmt::format(
-        "Two or more meshes use the same unique ID '{}' in the same input file",
+      fatal_error(fmt::format("Two or more meshes use the same unique ID "
+                              "'{}' in the same input file",
         id));
     }
     mesh_ids.insert(id);
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 6e95bdc769f..14ad6a0c4c3 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -419,7 +419,7 @@ def test_mesh_get_homogenized_materials():
     mesh_void.lower_left = (0.5, 0.5, -1.)
     mesh_void.upper_right = (1.5, 1.5, 1.)
     mesh_void.dimension = (1, 1, 1)
-    m4, = mesh_void.get_homogenized_materials(model, n_samples=10_000)
+    m4, = mesh_void.get_homogenized_materials(model, n_samples=(100, 100, 0))
 
     # Mesh element that overlaps void should have half density
     assert m4.get_mass_density('H1') == pytest.approx(0.5, rel=1e-2)
@@ -431,8 +431,8 @@ def test_mesh_get_homogenized_materials():
     assert m5.get_mass_density('H1') == pytest.approx(1.0)
 
 
-def test_mesh_material_volumes():
-    """Test the material_volumes method"""
+@pytest.fixture
+def sphere_model():
     # Model with three materials separated by planes x=0 and z=0
     mats = []
     for i in range(3):
@@ -449,23 +449,32 @@ def test_mesh_material_volumes():
     cell3 = openmc.Cell(fill=mats[2], region=-sph & -z0)
     model = openmc.Model()
     model.geometry = openmc.Geometry([cell1, cell2, cell3])
+    model.materials = openmc.Materials(mats)
+    return model
 
-    # Test material volumes on a regular mesh
+
+@pytest.mark.parametrize("n_rays", [1000, (10, 10, 0), (10, 0, 10), (0, 10, 10)])
+def test_material_volumes_regular_mesh(sphere_model, n_rays):
+    """Test the material_volumes method on a regular mesh"""
     mesh = openmc.RegularMesh()
     mesh.lower_left = (-1., -1., -1.)
     mesh.upper_right = (1., 1., 1.)
     mesh.dimension = (2, 2, 2)
-    volumes = mesh.material_volumes(model)
+    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.])
     assert volumes.by_element(4) == [(mats[1].id, 1.)]
     assert volumes.by_element(0) == [(mats[2].id, 1.)]
 
-    # Test material volumes on a cylindrical mesh
+
+def test_material_volumes_cylindrical_mesh(sphere_model):
+    """Test the material_volumes method on a cylindrical mesh"""
     cyl_mesh = openmc.CylindricalMesh(
         [0., 1.], [-1., 0., 1.,], [0.0, pi/4, 3*pi/4, 5*pi/4, 7*pi/4, 2*pi])
-    volumes = cyl_mesh.material_volumes(model)
+    volumes = cyl_mesh.material_volumes(sphere_model, (0, 100, 100))
+    mats = sphere_model.materials
     np.testing.assert_almost_equal(volumes[mats[0].id], [
         0., 0., 0., 0., 0.,
         pi/8, pi/8, 0., pi/8, pi/8

From 190b5463af89031f770c1989c497887fe27a4a75 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 18 Sep 2024 15:54:30 -0500
Subject: [PATCH 19/23] Fix typo in mesh.cpp

Co-authored-by: Olek <45364492+yardasol@users.noreply.github.com>
---
 src/mesh.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mesh.cpp b/src/mesh.cpp
index 53814c6c19c..591f0bb110d 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -706,7 +706,7 @@ xt::xtensor<double, 1> StructuredMesh::count_sites(
 
 // raytrace through the mesh. The template class T will do the tallying.
 // A modern optimizing compiler can recognize the noop method of T and
-// eleminate that call entirely.
+// eliminate that call entirely.
 template<class T>
 void StructuredMesh::raytrace_mesh(
   Position r0, Position r1, const Direction& u, T tally) const

From 095473b44ff46beee70347b3a8319262b0c30f32 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 5 Dec 2024 12:42:20 -0600
Subject: [PATCH 20/23] Fix comments

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

diff --git a/src/mesh.cpp b/src/mesh.cpp
index ca0df92e5ca..ae3ef4f43fa 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -200,7 +200,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int max_materials,
 
   std::array<int, 3> n_rays = {nx, ny, nz};
 
-  // Loop over rays on left face of bounding box
+  // Determine effective width of rays
   Position width((bbox.xmax - bbox.xmin) / nx, (bbox.ymax - bbox.ymin) / ny,
     (bbox.zmax - bbox.zmin) / nz);
 
@@ -235,6 +235,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int max_materials,
       int n1 = n_rays[ax1];
       int n2 = n_rays[ax2];
 
+      // Loop over rays on face of bounding box
 #pragma omp for
       for (int i1 = 0; i1 < n1; ++i1) {
         site.r[ax1] = min1 + (i1 + 0.5) * d1;

From 38b246a494890aa2b5307a0e9b3fc0e2cab44d4b Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sun, 29 Dec 2024 18:07:21 -0500
Subject: [PATCH 21/23] Parallelize mesh material volume calc with MPI

---
 src/mesh.cpp | 85 +++++++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 81 insertions(+), 4 deletions(-)

diff --git a/src/mesh.cpp b/src/mesh.cpp
index ae3ef4f43fa..425cbf3c7fe 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -29,6 +29,7 @@
 #include "openmc/memory.h"
 #include "openmc/message_passing.h"
 #include "openmc/openmp_interface.h"
+#include "openmc/output.h"
 #include "openmc/particle_data.h"
 #include "openmc/plot.h"
 #include "openmc/random_dist.h"
@@ -37,6 +38,7 @@
 #include "openmc/string_utils.h"
 #include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
+#include "openmc/timer.h"
 #include "openmc/volume_calc.h"
 #include "openmc/xml_interface.h"
 
@@ -192,6 +194,20 @@ vector<double> Mesh::volumes() const
 void Mesh::material_volumes(int nx, int ny, int nz, int max_materials,
   int32_t* materials, double* volumes) const
 {
+  if (mpi::master) {
+    header("MESH MATERIAL VOLUMES CALCULATION", 7);
+  }
+  write_message(7, "Number of rays (x) = {}", nx);
+  write_message(7, "Number of rays (y) = {}", ny);
+  write_message(7, "Number of rays (z) = {}", nz);
+  int64_t n_total = nx * ny + ny * nz + nx * nz;
+  write_message(7, "Total number of rays = {}", n_total);
+  write_message(
+    7, "Maximum number of materials per mesh element = {}", max_materials);
+
+  Timer timer;
+  timer.start();
+
   // Create object for keeping track of materials/volumes
   detail::MaterialVolumes result(materials, volumes, max_materials);
 
@@ -235,11 +251,19 @@ void Mesh::material_volumes(int nx, int ny, int nz, int max_materials,
       int n1 = n_rays[ax1];
       int n2 = n_rays[ax2];
 
+      // Divide rays in first direction over MPI processes by computing starting
+      // and ending indices
+      int min_work = n1 / mpi::n_procs;
+      int remainder = n1 % mpi::n_procs;
+      int n1_local = (mpi::rank < remainder) ? min_work + 1 : min_work;
+      int i1_start = mpi::rank * min_work + std::min(mpi::rank, remainder);
+      int i1_end = i1_start + n1_local;
+
       // Loop over rays on face of bounding box
-#pragma omp for
-      for (int i1 = 0; i1 < n1; ++i1) {
-        site.r[ax1] = min1 + (i1 + 0.5) * d1;
+#pragma omp for collapse(2)
+      for (int i1 = i1_start; i1 < i1_end; ++i1) {
         for (int i2 = 0; i2 < n2; ++i2) {
+          site.r[ax1] = min1 + (i1 + 0.5) * d1;
           site.r[ax2] = min2 + (i2 + 0.5) * d2;
 
           p.from_source(&site);
@@ -329,6 +353,46 @@ void Mesh::material_volumes(int nx, int ny, int nz, int max_materials,
                              "volume calculation insufficient.");
   }
 
+  // Compute time for raytracing
+  double t_raytrace = timer.elapsed();
+
+#ifdef OPENMC_MPI
+  // Combine results from multiple MPI processes
+  if (mpi::n_procs > 1) {
+    int total = this->n_bins() * max_materials;
+    if (mpi::master) {
+      // Allocate temporary buffer for receiving data
+      std::vector<int32_t> mats(total);
+      std::vector<double> vols(total);
+
+      for (int i = 1; i < mpi::n_procs; ++i) {
+        // Receive material indices and volumes from process i
+        MPI_Recv(
+          mats.data(), total, MPI_INT, i, i, mpi::intracomm, MPI_STATUS_IGNORE);
+        MPI_Recv(vols.data(), total, MPI_DOUBLE, i, i, mpi::intracomm,
+          MPI_STATUS_IGNORE);
+
+        // Combine with existing results
+        for (int index_elem = 0; index_elem < n_bins(); ++index_elem) {
+          for (int k = 0; k < max_materials; ++k) {
+            int index = index_elem * max_materials + k;
+            result.add_volume(index_elem, mats[index], vols[index]);
+          }
+        }
+      }
+    } else {
+      // Send material indices and volumes to process 0
+      MPI_Send(materials, total, MPI_INT, 0, mpi::rank, mpi::intracomm);
+      MPI_Send(volumes, total, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm);
+    }
+  }
+
+  // Report time for MPI communication
+  double t_mpi = timer.elapsed() - t_raytrace;
+#else
+  double t_mpi = 0.0;
+#endif
+
   // Normalize based on known volumes of elements
   for (int i = 0; i < this->n_bins(); ++i) {
     // Estimated total volume in element i
@@ -337,12 +401,25 @@ void Mesh::material_volumes(int nx, int ny, int nz, int max_materials,
       volume += result.volumes(i, j);
     }
 
-    // Renormalize volumes based on known volume of elemnet i
+    // Renormalize volumes based on known volume of element i
     double norm = this->volume(i) / volume;
     for (int j = 0; j < max_materials; ++j) {
       result.volumes(i, j) *= norm;
     }
   }
+
+  // Show elapsed time
+  timer.stop();
+  double t_total = timer.elapsed();
+  double t_normalize = t_total - t_raytrace - t_mpi;
+  if (mpi::master) {
+    header("Timing Statistics", 7);
+    show_time("Total time elapsed", t_total);
+    show_time("Ray tracing", t_raytrace, 1);
+    show_time("Ray tracing (per ray)", t_raytrace / n_total, 1);
+    show_time("MPI communication", t_mpi, 1);
+    show_time("Normalization", t_normalize, 1);
+  }
 }
 
 //==============================================================================

From 7fc29c9d48c1003405e7b76801060b70e3622ea2 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 3 Jan 2025 12:01:00 -0600
Subject: [PATCH 22/23] Parallelize calls to add_volume when collecting results
 across processes

---
 include/openmc/mesh.h |  1 +
 src/mesh.cpp          | 36 +++++++++++++++++++++++++++++++++---
 2 files changed, 34 insertions(+), 3 deletions(-)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index d03926f0baf..d01c816f8e6 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -87,6 +87,7 @@ class MaterialVolumes {
   //! \param[in] index_material Index of the material
   //! \param[in] volume Volume to add
   void add_volume(int index_elem, int index_material, double volume);
+  void add_volume_unsafe(int index_elem, int index_material, double volume);
 
   // Accessors
   int32_t& materials(int i, int j) { return materials_[i * n_mats_ + j]; }
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 425cbf3c7fe..4ffaeb28932 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -106,7 +106,7 @@ inline bool check_intersection_point(double x1, double x0, double y1, double y0,
 }
 
 //==============================================================================
-// Mesh implementation
+// MaterialVolumes implementation
 //==============================================================================
 
 namespace detail {
@@ -140,6 +140,34 @@ void MaterialVolumes::add_volume(
   this->volumes(index_elem, i) += volume;
 }
 
+// Same as add_volume above, but without the OpenMP critical/atomic section
+void MaterialVolumes::add_volume_unsafe(
+  int index_elem, int index_material, double volume)
+{
+  int i;
+  for (i = 0; i < n_mats_; ++i) {
+    // Check whether material already is present
+    if (this->materials(index_elem, i) == index_material)
+      break;
+
+    // If not already present, check for unused position (-2)
+    if (this->materials(index_elem, i) == -2) {
+      this->materials(index_elem, i) = index_material;
+      break;
+    }
+  }
+
+  // If maximum number of materials exceeded, set a flag that can be checked
+  // later
+  if (i >= n_mats_) {
+    too_many_mats_ = true;
+    return;
+  }
+
+  // Accumulate volume
+  this->volumes(index_elem, i) += volume;
+}
+
 } // namespace detail
 
 //==============================================================================
@@ -372,11 +400,13 @@ void Mesh::material_volumes(int nx, int ny, int nz, int max_materials,
         MPI_Recv(vols.data(), total, MPI_DOUBLE, i, i, mpi::intracomm,
           MPI_STATUS_IGNORE);
 
-        // Combine with existing results
+        // Combine with existing results; we can call thread unsafe version of
+        // add_volume because each thread is operating on a different element
+#pragma omp for
         for (int index_elem = 0; index_elem < n_bins(); ++index_elem) {
           for (int k = 0; k < max_materials; ++k) {
             int index = index_elem * max_materials + k;
-            result.add_volume(index_elem, mats[index], vols[index]);
+            result.add_volume_unsafe(index_elem, mats[index], vols[index]);
           }
         }
       }

From d4ae59bfda7e774fcf92760633d69d2655ef4abe Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 3 Jan 2025 17:04:59 -0600
Subject: [PATCH 23/23] Fix output for Mesh.material_volumes

---
 openmc/lib/mesh.py | 12 ++++++++----
 openmc/mesh.py     |  3 ++-
 src/mesh.cpp       |  1 +
 3 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py
index 28339490d55..0a53fd251e7 100644
--- a/openmc/lib/mesh.py
+++ b/openmc/lib/mesh.py
@@ -10,7 +10,7 @@
 
 from ..exceptions import AllocationError, InvalidIDError
 from . import _dll
-from .core import _FortranObjectWithID
+from .core import _FortranObjectWithID, quiet_dll
 from .error import _error_handler
 from .plot import _Position
 from ..bounding_box import BoundingBox
@@ -186,7 +186,8 @@ def bounding_box(self) -> BoundingBox:
     def material_volumes(
             self,
             n_samples: int | tuple[int, int, int] = 10_000,
-            max_materials: int = 4
+            max_materials: int = 4,
+            output: bool = True,
     ) -> MeshMaterialVolumes:
         """Determine volume of materials in each mesh element.
 
@@ -210,6 +211,8 @@ def material_volumes(
             the x, y, and z dimensions.
         max_materials : int, optional
             Estimated maximum number of materials in any given mesh element.
+        output : bool, optional
+            Whether or not to show output.
 
         Returns
         -------
@@ -238,8 +241,9 @@ def material_volumes(
         # Run material volume calculation
         while True:
             try:
-                _dll.openmc_mesh_material_volumes(
-                    self._index, nx, ny, nz, max_materials, materials, volumes)
+                with quiet_dll(output):
+                    _dll.openmc_mesh_material_volumes(
+                        self._index, nx, ny, nz, max_materials, materials, volumes)
             except AllocationError:
                 # Increase size of result array and try again
                 max_materials *= 2
diff --git a/openmc/mesh.py b/openmc/mesh.py
index 106bde1fadb..4e59c37179f 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -392,7 +392,8 @@ def material_volumes(
             kwargs.setdefault('output', False)
             openmc.lib.init(['-c'], **kwargs)
             mesh = openmc.lib.tallies[new_tally.id].filters[0].mesh
-            volumes = mesh.material_volumes(n_samples, max_materials)
+            volumes = mesh.material_volumes(
+                n_samples, max_materials, output=kwargs['output'])
             openmc.lib.finalize()
 
             # Restore original tallies
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 4ffaeb28932..970edbe595a 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -449,6 +449,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int max_materials,
     show_time("Ray tracing (per ray)", t_raytrace / n_total, 1);
     show_time("MPI communication", t_mpi, 1);
     show_time("Normalization", t_normalize, 1);
+    std::fflush(stdout);
   }
 }