From f10d9ac7b629e984fb6fbc9adfde0217a1dacb50 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Fri, 25 Feb 2022 16:56:56 -0600
Subject: [PATCH] Making volume retrieval consistent in Python API

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

diff --git a/openmc/mesh.py b/openmc/mesh.py
index 6c210ae596..9a217497f1 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -203,6 +203,22 @@ class RegularMesh(MeshBase):
     def num_mesh_cells(self):
         return np.prod(self._dimension)
 
+    @property
+    def volumes(self):
+        """Return Volumes for every mesh cell
+
+        Returns
+        -------
+        volumes : Iterable of float
+            Volumes
+
+        """
+        return np.full(self.dimension, np.prod(self.width))
+
+    @property
+    def total_volume(self):
+        return np.prod(self.dimension) * np.prod(self.width)
+
     @property
     def indices(self):
         ndim = len(self._dimension)
@@ -575,6 +591,26 @@ class RectilinearMesh(MeshBase):
     def z_grid(self):
         return self._z_grid
 
+    @property
+    def volumes(self):
+        """Return Volumes for every mesh cell
+
+        Returns
+        -------
+        volumes : Iterable of float
+            Volumes
+
+        """
+        V_x = np.diff(self.x_grid)
+        V_y = np.diff(self.y_grid)
+        V_z = np.diff(self.z_grid)
+
+        return np.multiply.outer(np.outer(V_x, V_y), V_z)
+
+    @property
+    def total_volume(self):
+        return np.sum(self.volumes)
+
     @property
     def indices(self):
         nx = len(self.x_grid) - 1
@@ -681,6 +717,7 @@ class RectilinearMesh(MeshBase):
 
         return element
 
+
 class CylindricalMesh(MeshBase):
     """A 3D cylindrical mesh
 
@@ -850,7 +887,8 @@ class CylindricalMesh(MeshBase):
         mesh.z_grid = [float(x) for x in get_text(elem, "z_grid").split()]
         return mesh
 
-    def calc_mesh_volumes(self):
+    @property
+    def volumes(self):
         """Return Volumes for every mesh cell
 
         Returns
@@ -1037,7 +1075,8 @@ class SphericalMesh(MeshBase):
         mesh.phi_grid = [float(x) for x in get_text(elem, "phi_grid").split()]
         return mesh
 
-    def calc_mesh_volumes(self):
+    @property
+    def volumes(self):
         """Return Volumes for every mesh cell
 
         Returns
@@ -1150,6 +1189,15 @@ class UnstructuredMesh(MeshBase):
 
     @property
     def volumes(self):
+        """Return Volumes for every mesh cell if
+        populated by a StatePoint file
+
+        Returns
+        -------
+        volumes : Iterable of float
+            Volumes
+
+        """
         return self._volumes
 
     @volumes.setter