implement to_hdf5 for unstructured meshes using a new c api function

Author: GuySten

openmc/lib/mesh.py

@@ -1,5 +1,5 @@
 from collections.abc import Mapping, Sequence
-from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER, c_void_p,
+from ctypes import (c_int, c_int32, c_int64, c_char_p, c_double, POINTER, c_void_p,
                     create_string_buffer, c_size_t)
 from math import sqrt
 import sys
@@ -18,7 +18,7 @@
 
 __all__ = [
     'Mesh', 'RegularMesh', 'RectilinearMesh', 'CylindricalMesh',
-    'SphericalMesh', 'UnstructuredMesh', 'meshes', 'MeshMaterialVolumes'
+    'SphericalMesh', 'UnstructuredMesh', 'meshes', 'MeshMaterialVolumes', 'export_unstructured_mesh'
 ]
 
 
@@ -108,6 +108,10 @@
 _dll.openmc_spherical_mesh_set_grid.restype = c_int
 _dll.openmc_spherical_mesh_set_grid.errcheck = _error_handler
 
+_dll.openmc_unstructured_mesh_export_hdf5.argtypes = [c_int32, c_int64]
+_dll.openmc_unstructured_mesh_export_hdf5.restype = c_int
+_dll.openmc_unstructured_mesh_export_hdf5.errcheck = _error_handler
+
 
 class Mesh(_FortranObjectWithID):
     """Base class to represent mesh objects
@@ -740,6 +744,12 @@ class UnstructuredMesh(Mesh):
 }
 
 
+def export_unstructured_mesh(mesh, group):
+    index = c_int32()
+    _dll.openmc_get_mesh_index(mesh.id, index)
+    _dll.openmc_unstructured_mesh_export_hdf5(index.value, int(group.id.id))
+
+
 def _get_mesh(index):
     mesh_type = create_string_buffer(20)
     _dll.openmc_mesh_get_type(index, mesh_type)

openmc/mesh.py

@@ -3295,11 +3295,20 @@ def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
         return mesh
 
     def to_hdf5(self, group: h5py.Group):
-        raise NotImplementedError(
-            "UnstructuredMesh.to_hdf5 is not implemented in Python. "
-            "Use openmc.lib.export_weight_windows() to export weight "
-            "windows on unstructured meshes."
+        import openmc.lib
+        mesh_group = super().to_hdf5(group)
+        model = openmc.Model()
+        sph = openmc.Sphere(boundary_type='vacuum')
+        model.geometry = openmc.Geometry([openmc.Cell(region=-sph)])
+        model.settings.particles = 100
+        model.settings.batches = 1
+        wwg = openmc.WeightWindowGenerator(
+              method='magic',
+              mesh=self,
         )
+        model.settings.weight_window_generators = [wwg] 
+        with openmc.lib.TemporarySession(model):
+            openmc.lib.export_unstructured_mesh(self, mesh_group)
 
     def to_xml_element(self):
         """Return XML representation of the mesh

openmc/weight_windows.py

@@ -1081,19 +1081,6 @@ def export_to_hdf5(self, path: PathLike = 'weight_windows.h5', **init_kwargs):
         cv.check_type('path', path, PathLike)
         path = Path(path).resolve()
 
-        # Any unstructured mesh forces the whole list onto the lib fallback.
-        if any(isinstance(ww.mesh, UnstructuredMesh) for ww in self):
-            import openmc.lib
-            model = openmc.Model()
-            sph = openmc.Sphere(boundary_type='vacuum')
-            model.geometry = openmc.Geometry([openmc.Cell(region=-sph)])
-            model.settings.weight_windows = self
-            model.settings.particles = 100
-            model.settings.batches = 1
-            with openmc.lib.TemporarySession(model, **init_kwargs):
-                openmc.lib.export_weight_windows(path)
-            return
-
         with h5py.File(path, 'w') as f:
             f.attrs['filetype'] = np.bytes_('weight_windows')
             f.attrs['version'] = np.asarray([1, 0], dtype=np.int32)

src/mesh.cpp

@@ -2843,6 +2843,27 @@ extern "C" int openmc_spherical_mesh_set_grid(int32_t index,
     index, grid_x, nx, grid_y, ny, grid_z, nz);
 }
 
+extern "C" int openmc_unstructured_mesh_export_hdf5(
+  int32_t index, hid_t mesh_group)
+{
+  if (int err = check_mesh_type<UnstructuredMesh>(index))
+    return err;
+  UnstructuredMesh* m =
+    dynamic_cast<UnstructuredMesh*>(model::meshes[index].get());
+
+  // Write mesh type
+  write_dataset(mesh_group, "type", m->get_mesh_type());
+
+  // Write mesh ID
+  write_attribute(mesh_group, "id", m->id_);
+
+  // Write mesh name
+  write_dataset(mesh_group, "name", m->name_);
+
+  m->to_hdf5_inner(mesh_group);
+  return 0;
+}
+
 #ifdef OPENMC_DAGMC_ENABLED
 
 const std::string MOABMesh::mesh_lib_type = "moab";