Add openmc_mesh_get_volumes C API function (#2869)

include/openmc/capi.h

@@ -106,6 +106,7 @@ int openmc_mesh_filter_set_translation(int32_t index, double translation[3]);
 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_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);

openmc/lib/mesh.py

@@ -39,6 +39,9 @@ class _MaterialVolume(Structure):
 _dll.openmc_mesh_get_n_elements.argtypes = [c_int32, POINTER(c_size_t)]
 _dll.openmc_mesh_get_n_elements.restype = c_int
 _dll.openmc_mesh_get_n_elements.errcheck = _error_handler
+_dll.openmc_mesh_get_volumes.argtypes = [c_int32, POINTER(c_double)]
+_dll.openmc_mesh_get_volumes.restype = c_int
+_dll.openmc_mesh_get_volumes.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)]
@@ -149,11 +152,18 @@ def id(self, mesh_id):
         _dll.openmc_mesh_set_id(self._index, mesh_id)
 
     @property
-    def n_elements(self):
+    def n_elements(self) -> int:
         n = c_size_t()
         _dll.openmc_mesh_get_n_elements(self._index, n)
         return n.value
 
+    @property
+    def volumes(self) -> np.ndarray:
+        volumes = np.empty((self.n_elements,))
+        _dll.openmc_mesh_get_volumes(
+            self._index, volumes.ctypes.data_as(POINTER(c_double)))
+        return volumes
+
     def material_volumes(
             self,
             n_samples: int = 10_000,
@@ -276,6 +286,10 @@ class RegularMesh(Mesh):
         are given, it is assumed that the mesh is an x-y mesh.
     width : numpy.ndarray
         The width of mesh cells in each direction.
+    n_elements : int
+        Total number of mesh elements.
+    volumes : numpy.ndarray
+        Volume of each mesh element in [cm^3]
 
     """
     mesh_type = 'regular'
@@ -358,6 +372,10 @@ class RectilinearMesh(Mesh):
         The upper-right corner of the structrued mesh.
     width : numpy.ndarray
         The width of mesh cells in each direction.
+    n_elements : int
+        Total number of mesh elements.
+    volumes : numpy.ndarray
+        Volume of each mesh element in [cm^3]
 
     """
     mesh_type = 'rectilinear'
@@ -457,6 +475,10 @@ class CylindricalMesh(Mesh):
         The upper-right corner of the structrued mesh.
     width : numpy.ndarray
         The width of mesh cells in each direction.
+    n_elements : int
+        Total number of mesh elements.
+    volumes : numpy.ndarray
+        Volume of each mesh element in [cm^3]
 
     """
     mesh_type = 'cylindrical'
@@ -555,6 +577,10 @@ class SphericalMesh(Mesh):
         The upper-right corner of the structrued mesh.
     width : numpy.ndarray
         The width of mesh cells in each direction.
+    n_elements : int
+        Total number of mesh elements.
+    volumes : numpy.ndarray
+        Volume of each mesh element in [cm^3]
 
     """
     mesh_type = 'spherical'

openmc/material.py

@@ -571,7 +571,7 @@ def add_components(self, components: dict, percent_type: str = 'ao'):
 
         for component, params in components.items():
             cv.check_type('component', component, str)
-            if isinstance(params, float):
+            if isinstance(params, Real):
                 params = {'percent': params}
 
             else:

src/mesh.cpp

@@ -1881,6 +1881,17 @@ extern "C" int openmc_mesh_get_n_elements(int32_t index, size_t* n)
   return 0;
 }
 
+//! Get the volume of each element in the mesh
+extern "C" int openmc_mesh_get_volumes(int32_t index, double* volumes)
+{
+  if (int err = check_mesh(index))
+    return err;
+  for (int i = 0; i < model::meshes[index]->n_bins(); ++i) {
+    volumes[i] = model::meshes[index]->volume(i);
+  }
+  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)
 {

tests/unit_tests/test_lib.py

@@ -566,6 +566,8 @@ def test_regular_mesh(lib_init):
     assert mesh.upper_right == pytest.approx(ur)
     assert mesh.width == pytest.approx(width)
 
+    np.testing.assert_allclose(mesh.volumes, 1.0)
+
     meshes = openmc.lib.meshes
     assert isinstance(meshes, Mapping)
     assert len(meshes) == 1
@@ -644,6 +646,8 @@ def test_rectilinear_mesh(lib_init):
             for k, diff_z in enumerate(np.diff(z_grid)):
                 assert np.all(mesh.width[i, j, k, :] == (10, 10, 10))
 
+    np.testing.assert_allclose(mesh.volumes, 1000.0)
+
     with pytest.raises(exc.AllocationError):
         mesh2 = openmc.lib.RectilinearMesh(mesh.id)
 
@@ -688,6 +692,9 @@ def test_cylindrical_mesh(lib_init):
             for k, _ in enumerate(np.diff(z_grid)):
                 assert np.allclose(mesh.width[i, j, k, :], (5, deg2rad(10), 10))
 
+    np.testing.assert_allclose(mesh.volumes[::2], 10/360 * pi * 5**2 * 10)
+    np.testing.assert_allclose(mesh.volumes[1::2], 10/360 * pi * (10**2 - 5**2) * 10)
+
     with pytest.raises(exc.AllocationError):
         mesh2 = openmc.lib.CylindricalMesh(mesh.id)
 
@@ -734,6 +741,13 @@ def test_spherical_mesh(lib_init):
             for k, _ in enumerate(np.diff(phi_grid)):
                 assert np.allclose(mesh.width[i, j, k, :], (5, deg2rad(10), deg2rad(10)))
 
+    dtheta = lambda d1, d2: np.cos(deg2rad(d1)) - np.cos(deg2rad(d2))
+    f = 1/3 * deg2rad(10.)
+    np.testing.assert_allclose(mesh.volumes[::4],  f * 5**3 * dtheta(0., 10.))
+    np.testing.assert_allclose(mesh.volumes[1::4], f * (10**3 - 5**3) * dtheta(0., 10.))
+    np.testing.assert_allclose(mesh.volumes[2::4], f * 5**3 * dtheta(10., 20.))
+    np.testing.assert_allclose(mesh.volumes[3::4], f * (10**3 - 5**3) * dtheta(10., 20.))
+
     with pytest.raises(exc.AllocationError):
         mesh2 = openmc.lib.SphericalMesh(mesh.id)