From b295f4a0952b5cbdcf700a73cd992e9be98f27c6 Mon Sep 17 00:00:00 2001 From: Jon Shimwell Date: Wed, 8 Jan 2025 17:30:18 +0100 Subject: [PATCH 01/22] added option to write vtkhdf files for umesh --- openmc/mesh.py | 225 ++++++++++++---- .../test_mesh_dagmc_tets.vtk | 253 ++++++++++++++++++ tests/unit_tests/test_mesh.py | 68 +++++ 3 files changed, 489 insertions(+), 57 deletions(-) create mode 100644 tests/regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk diff --git a/openmc/mesh.py b/openmc/mesh.py index 0b6f6b84e4b..019c44c8b46 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -5,6 +5,7 @@ from functools import wraps from math import pi, sqrt, atan2 from numbers import Integral, Real +from pathlib import Path import h5py import lxml.etree as ET @@ -2368,7 +2369,9 @@ def write_data_to_vtk( Parameters ---------- filename : str or pathlib.Path - Name of the VTK file to write + Name of the VTK file to write. If the filename ends in '.hdf' then a VTKHDF + format file will be written and can be opened with Paraview versions 5.13.0 + and above, if the filename ends in '.vtk' then a .vtk file will be written. datasets : dict Dictionary whose keys are the data labels and values are numpy appropriately sized arrays @@ -2377,78 +2380,186 @@ def write_data_to_vtk( Whether or not to normalize the data by the volume of the mesh elements """ - import vtk - from vtk.util import numpy_support as nps - if self.connectivity is None or self.vertices is None: - raise RuntimeError('This mesh has not been ' - 'loaded from a statepoint file.') + if Path(filename).suffix == ".hdf": + + def append_dataset(dset, array): + """Convenience function to append data to an HDF5 dataset""" + origLen = dset.shape[0] + dset.resize(origLen + array.shape[0], axis=0) + dset[origLen:] = array + + if self.library != "moab": + raise NotImplemented("VTKHDF output is only supported for MOAB meshes") + + # the self.connectivity contains an arrays of length 8, in the case of + # DAGMC tetrahedra mesh elements, the last 4 values are -1 and can be removed + trimmed_connectivity = [] + for cell in self.connectivity: + # Find the index of the first -1 value, if any + first_negative_index = np.where(cell == -1)[0] + if first_negative_index.size > 0: + # Slice the array up to the first -1 value + trimmed_connectivity.append(cell[: first_negative_index[0]]) + else: + # No -1 values, append the whole cell + trimmed_connectivity.append(cell) + trimmed_connectivity = np.array( + trimmed_connectivity, dtype="int32" + ).flatten() + + # DAGMC supports tet meshes only so we know it has 4 points per cell + points_per_cell = 4 + + # offsets are the indices of the first point of each cell in the array of points + offsets = np.arange( + 0, self.n_elements * points_per_cell + 1, points_per_cell + ) - if filename is None: - filename = f'mesh_{self.id}.vtk' + with h5py.File(filename, "w") as f: + + root = f.create_group("VTKHDF") + root.attrs["Version"] = (2, 1) + ascii_type = "UnstructuredGrid".encode("ascii") + root.attrs.create( + "Type", + ascii_type, + dtype=h5py.string_dtype("ascii", len(ascii_type)), + ) + + # create hdf5 file structure + root.create_dataset( + "NumberOfPoints", (0,), maxshape=(None,), dtype="i8" + ) + root.create_dataset("Types", (0,), maxshape=(None,), dtype="uint8") + root.create_dataset("Points", (0, 3), maxshape=(None, 3), dtype="f") + root.create_dataset( + "NumberOfConnectivityIds", (0,), maxshape=(None,), dtype="i8" + ) + root.create_dataset("NumberOfCells", (0,), maxshape=(None,), dtype="i8") + root.create_dataset("Offsets", (0,), maxshape=(None,), dtype="i8") + root.create_dataset("Connectivity", (0,), maxshape=(None,), dtype="i8") + + append_dataset(root["NumberOfPoints"], np.array([len(self.vertices)])) + append_dataset(root["Points"], self.vertices) + append_dataset( + root["NumberOfConnectivityIds"], + np.array([len(trimmed_connectivity)]), + ) + append_dataset(root["Connectivity"], trimmed_connectivity) + append_dataset(root["NumberOfCells"], np.array([self.n_elements])) + append_dataset(root["Offsets"], offsets) + + # VTK_TETRA type is known as DAGMC only supports tet meshes + append_dataset( + root["Types"], np.full(self.n_elements, 10, dtype="uint8") + ) + + cell_data_group = root.create_group("CellData") - writer = vtk.vtkUnstructuredGridWriter() + for name, data in datasets.items(): - writer.SetFileName(str(filename)) + if data.shape != self.dimension: + raise ValueError( + f'Cannot apply dataset "{name}" with ' + f"shape {data.shape} to mesh {self.id} " + f"with dimensions {self.dimension}" + ) - grid = vtk.vtkUnstructuredGrid() + cell_data_group.create_dataset( + name, (0,), maxshape=(None,), dtype="float64", chunks=True + ) - vtk_pnts = vtk.vtkPoints() - vtk_pnts.SetData(nps.numpy_to_vtk(self.vertices)) - grid.SetPoints(vtk_pnts) + if volume_normalization: + data /= self.volumes + append_dataset(cell_data_group[name], data) - n_skipped = 0 - for elem_type, conn in zip(self.element_types, self.connectivity): - if elem_type == self._LINEAR_TET: - elem = vtk.vtkTetra() - elif elem_type == self._LINEAR_HEX: - elem = vtk.vtkHexahedron() - elif elem_type == self._UNSUPPORTED_ELEM: - n_skipped += 1 - else: - raise RuntimeError(f'Invalid element type {elem_type} found') - for i, c in enumerate(conn): - if c == -1: - break - elem.GetPointIds().SetId(i, c) + elif Path(filename).suffix == ".vtk": - grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds()) + import vtk + from vtk.util import numpy_support as nps - if n_skipped > 0: - warnings.warn(f'{n_skipped} elements were not written because ' - 'they are not of type linear tet/hex') + if self.connectivity is None or self.vertices is None: + raise RuntimeError( + "This mesh has not been " "loaded from a statepoint file." + ) - # check that datasets are the correct size - datasets_out = [] - if datasets is not None: - for name, data in datasets.items(): - if data.shape != self.dimension: - raise ValueError(f'Cannot apply dataset "{name}" with ' - f'shape {data.shape} to mesh {self.id} ' - f'with dimensions {self.dimension}') + if filename is None: + filename = f"mesh_{self.id}.vtk" + + writer = vtk.vtkUnstructuredGridWriter() - if volume_normalization: + writer.SetFileName(str(filename)) + + grid = vtk.vtkUnstructuredGrid() + + vtk_pnts = vtk.vtkPoints() + vtk_pnts.SetData(nps.numpy_to_vtk(self.vertices)) + grid.SetPoints(vtk_pnts) + + n_skipped = 0 + for elem_type, conn in zip(self.element_types, self.connectivity): + if elem_type == self._LINEAR_TET: + elem = vtk.vtkTetra() + elif elem_type == self._LINEAR_HEX: + elem = vtk.vtkHexahedron() + elif elem_type == self._UNSUPPORTED_ELEM: + n_skipped += 1 + else: + raise RuntimeError(f"Invalid element type {elem_type} found") + for i, c in enumerate(conn): + if c == -1: + break + elem.GetPointIds().SetId(i, c) + + grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds()) + + if n_skipped > 0: + warnings.warn( + f"{n_skipped} elements were not written because " + "they are not of type linear tet/hex" + ) + + # check that datasets are the correct size + datasets_out = [] + if datasets is not None: for name, data in datasets.items(): - if np.issubdtype(data.dtype, np.integer): - warnings.warn(f'Integer data set "{name}" will ' - 'not be volume-normalized.') - continue - data /= self.volumes + if data.shape != self.dimension: + raise ValueError( + f'Cannot apply dataset "{name}" with ' + f"shape {data.shape} to mesh {self.id} " + f"with dimensions {self.dimension}" + ) - # add data to the mesh - for name, data in datasets.items(): - datasets_out.append(data) - arr = vtk.vtkDoubleArray() - arr.SetName(name) - arr.SetNumberOfTuples(data.size) + if volume_normalization: + for name, data in datasets.items(): + if np.issubdtype(data.dtype, np.integer): + warnings.warn( + f'Integer data set "{name}" will ' + "not be volume-normalized." + ) + continue + data /= self.volumes + + # add data to the mesh + for name, data in datasets.items(): + datasets_out.append(data) + arr = vtk.vtkDoubleArray() + arr.SetName(name) + arr.SetNumberOfTuples(data.size) - for i in range(data.size): - arr.SetTuple1(i, data.flat[i]) - grid.GetCellData().AddArray(arr) + for i in range(data.size): + arr.SetTuple1(i, data.flat[i]) + grid.GetCellData().AddArray(arr) - writer.SetInputData(grid) + writer.SetInputData(grid) - writer.Write() + writer.Write() + + else: + raise ValueError( + f"Unsupported file extension for '{filename}'. Extension must be '.hdf' or '.vtk'." + ) @classmethod def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): diff --git a/tests/regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk b/tests/regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk new file mode 100644 index 00000000000..daad9e15a75 --- /dev/null +++ b/tests/regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk @@ -0,0 +1,253 @@ +# vtk DataFile Version 3.0 +MOAB 5.5.1 +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 84 double +-0.5 -0.5 0.5 +-0.5 -0.5 -0.5 +-0.5 0.5 0.5 +-0.5 0.5 -0.5 +0.5 -0.5 0.5 +0.5 -0.5 -0.5 +0.5 0.5 0.5 +0.5 0.5 -0.5 +-0.5 0 0 +0.5 0 0 +0 -0.5 0 +0 0.5 0 +0 0 -0.5 +0 0 0.5 +-0.5 -0.5 0.5 +-0.5 -0.5 -0.5 +-0.5 0.5 0.5 +-0.5 0.5 -0.5 +0.5 -0.5 0.5 +0.5 -0.5 -0.5 +0.5 0.5 0.5 +0.5 0.5 -0.5 +-0.5 0 0 +0.5 0 0 +0 -0.5 0 +0 0.5 0 +0 0 -0.5 +0 0 0.5 +-0.5 -0.5 0.5 +-0.5 -0.5 -0.5 +-0.5 0.5 0.5 +-0.5 0.5 -0.5 +0.5 -0.5 0.5 +0.5 -0.5 -0.5 +0.5 0.5 0.5 +0.5 0.5 -0.5 +-0.5 0 0 +0.5 0 0 +0 -0.5 0 +0 0.5 0 +0 0 -0.5 +0 0 0.5 +-0.5 -0.5 0.5 +-0.5 -0.5 -0.5 +-0.5 0.5 0.5 +-0.5 0.5 -0.5 +0.5 -0.5 0.5 +0.5 -0.5 -0.5 +0.5 0.5 0.5 +0.5 0.5 -0.5 +-0.5 0 0 +0.5 0 0 +0 -0.5 0 +0 0.5 0 +0 0 -0.5 +0 0 0.5 +-0.5 -0.5 0.5 +-0.5 -0.5 -0.5 +-0.5 0.5 0.5 +-0.5 0.5 -0.5 +0.5 -0.5 0.5 +0.5 -0.5 -0.5 +0.5 0.5 0.5 +0.5 0.5 -0.5 +-0.5 0 0 +0.5 0 0 +0 -0.5 0 +0 0.5 0 +0 0 -0.5 +0 0 0.5 +-0.5 -0.5 0.5 +-0.5 -0.5 -0.5 +-0.5 0.5 0.5 +-0.5 0.5 -0.5 +0.5 -0.5 0.5 +0.5 -0.5 -0.5 +0.5 0.5 0.5 +0.5 0.5 -0.5 +-0.5 0 0 +0.5 0 0 +0 -0.5 0 +0 0.5 0 +0 0 -0.5 +0 0 0.5 +CELLS 24 96 +3 1 0 8 +3 2 3 8 +3 8 3 1 +3 8 0 2 +3 19 23 18 +3 20 23 21 +3 23 19 21 +3 23 20 18 +3 28 29 38 +3 33 32 38 +3 38 32 28 +3 38 29 33 +3 44 53 45 +3 49 53 48 +3 53 44 48 +3 53 49 45 +3 68 57 59 +3 68 63 61 +3 59 63 68 +3 61 57 68 +3 83 72 70 +3 83 74 76 +3 72 83 76 +3 74 83 70 +CELL_TYPES 24 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +POINT_DATA 84 +SCALARS GLOBAL_ID int 1 +LOOKUP_TABLE default +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +CELL_DATA 24 +SCALARS GLOBAL_ID int 1 +LOOKUP_TABLE default +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 +-1 diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 27322165f64..1cbb1c4a8ea 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -3,6 +3,7 @@ import numpy as np import pytest import openmc +from pathlib import Path @pytest.mark.parametrize("val_left,val_right", [(0, 0), (-1., -1.), (2.0, 2)]) @@ -449,3 +450,70 @@ 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_umesh(run_in_tmpdir, request): + """Performs a minimal UnstructuredMesh simulation, reads in the resulting statepoint + file and writes the mesh data to vtk and hdf files. It is necessary to read in the + unstructured mesh from a statepoint file to ensure it has all the required attributes + """ + + surf1 = openmc.Sphere(R=1000.0, boundary_type="vacuum") + cell1 = openmc.Cell(region=-surf1) + my_geometry = openmc.Geometry([cell1]) + + umesh = openmc.UnstructuredMesh( + request.path.parent.parent + / "regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk", + "moab", + ) + umesh.id = ( + 1 # setting ID to make it easier to get the mesh from the statepoint later + ) + mesh_filter = openmc.MeshFilter(umesh) + + # Create flux mesh tally to score alpha production + mesh_tally = openmc.Tally(name="test_tally") + mesh_tally.filters = [mesh_filter] + mesh_tally.scores = ["flux"] + + tallies = openmc.Tallies([mesh_tally]) + + my_source = openmc.IndependentSource() + my_source.space = openmc.stats.Point((0.4, 0, 0.4)) + + settings = openmc.Settings() + settings.run_mode = "fixed source" + settings.batches = 1 + settings.particles = 100 + settings.source = my_source + + my_model = openmc.Model( + materials=None, geometry=my_geometry, settings=settings, tallies=tallies + ) + + statepoint_file = my_model.run() + + statepoint = openmc.StatePoint(statepoint_file) + + my_tally = statepoint.get_tally(name="test_tally") + + umesh_from_sp = statepoint.meshes[1] + + umesh_from_sp.write_data_to_vtk( + datasets={"mean": my_tally.mean.flatten()}, + filename="test_mesh.vtk", + ) + umesh_from_sp.write_data_to_vtk( + datasets={"mean": my_tally.mean.flatten()}, + filename="test_mesh.hdf", + ) + with pytest.raises(ValueError): + # Supported file extensions are vtk or hdf, not hdf5, so this should raise an error + umesh_from_sp.write_data_to_vtk( + datasets={"mean": my_tally.mean.flatten()}, + filename="test_mesh.hdf5", + ) + + assert Path("test_mesh.vtk").exists() + assert Path("test_mesh.hdf").exists() From f1ee07da387a9b21768e11e65eb4b765427f42f3 Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 8 Jan 2025 17:41:49 +0100 Subject: [PATCH 02/22] testing data shape --- tests/unit_tests/test_mesh.py | 41 ++++++++++++++++++++--------------- 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 1cbb1c4a8ea..9c652fcb97a 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -453,9 +453,10 @@ def test_mesh_get_homogenized_materials(): def test_umesh(run_in_tmpdir, request): - """Performs a minimal UnstructuredMesh simulation, reads in the resulting statepoint - file and writes the mesh data to vtk and hdf files. It is necessary to read in the - unstructured mesh from a statepoint file to ensure it has all the required attributes + """Performs a minimal UnstructuredMesh simulation, reads in the resulting + statepoint file and writes the mesh data to vtk and hdf files. It is + necessary to read in the unstructured mesh from a statepoint file to ensure + it has all the required attributes """ surf1 = openmc.Sphere(R=1000.0, boundary_type="vacuum") @@ -467,9 +468,8 @@ def test_umesh(run_in_tmpdir, request): / "regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk", "moab", ) - umesh.id = ( - 1 # setting ID to make it easier to get the mesh from the statepoint later - ) + # setting ID to make it easier to get the mesh from the statepoint later + umesh.id = 1 mesh_filter = openmc.MeshFilter(umesh) # Create flux mesh tally to score alpha production @@ -480,12 +480,11 @@ def test_umesh(run_in_tmpdir, request): tallies = openmc.Tallies([mesh_tally]) my_source = openmc.IndependentSource() - my_source.space = openmc.stats.Point((0.4, 0, 0.4)) settings = openmc.Settings() settings.run_mode = "fixed source" - settings.batches = 1 - settings.particles = 100 + settings.batches = 2 + settings.particles = 10 settings.source = my_source my_model = openmc.Model( @@ -500,20 +499,26 @@ def test_umesh(run_in_tmpdir, request): umesh_from_sp = statepoint.meshes[1] - umesh_from_sp.write_data_to_vtk( - datasets={"mean": my_tally.mean.flatten()}, - filename="test_mesh.vtk", - ) - umesh_from_sp.write_data_to_vtk( - datasets={"mean": my_tally.mean.flatten()}, - filename="test_mesh.hdf", - ) + datasets={ + "mean": my_tally.mean.flatten(), + "std_dev": my_tally.std_dev.flatten() + } + + umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtk") + umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.hdf") + with pytest.raises(ValueError): # Supported file extensions are vtk or hdf, not hdf5, so this should raise an error umesh_from_sp.write_data_to_vtk( - datasets={"mean": my_tally.mean.flatten()}, + datasets=datasets, filename="test_mesh.hdf5", ) + with pytest.raises(ValueError): + # Supported file extensions are vtk or hdf, not hdf5, so this should raise an error + umesh_from_sp.write_data_to_vtk( + datasets={'incorrectly_shaped_data': np.array(([1,2,3]))}, + filename="test_mesh.hdf", + ) assert Path("test_mesh.vtk").exists() assert Path("test_mesh.hdf").exists() From 52e6b8a9c07fe28b4f098e711b8d0382a4e0c61c Mon Sep 17 00:00:00 2001 From: shimwell Date: Wed, 8 Jan 2025 17:46:52 +0100 Subject: [PATCH 03/22] improved tests, added match string --- tests/unit_tests/test_mesh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 9c652fcb97a..6876ad89563 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -507,14 +507,14 @@ def test_umesh(run_in_tmpdir, request): umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtk") umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.hdf") - with pytest.raises(ValueError): + with pytest.raises(ValueError, match="Unsupported file extension"): # Supported file extensions are vtk or hdf, not hdf5, so this should raise an error umesh_from_sp.write_data_to_vtk( datasets=datasets, filename="test_mesh.hdf5", ) - with pytest.raises(ValueError): - # Supported file extensions are vtk or hdf, not hdf5, so this should raise an error + with pytest.raises(ValueError, match="Cannot apply dataset"): + # The shape of the data should match the shape of the mesh, so this should raise an error umesh_from_sp.write_data_to_vtk( datasets={'incorrectly_shaped_data': np.array(([1,2,3]))}, filename="test_mesh.hdf", From b2c775c55ab13737c94bf0aa70107d00a7edffba Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 9 Jan 2025 09:32:10 +0100 Subject: [PATCH 04/22] skip single test if no dagmc --- tests/unit_tests/test_mesh.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 6876ad89563..c87b193b9b2 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -451,7 +451,11 @@ def test_mesh_get_homogenized_materials(): model, n_samples=1000, include_void=False) assert m5.get_mass_density('H1') == pytest.approx(1.0) +skip_if_no_dagmc = pytest.mark.skipif( + not openmc.lib._dagmc_enabled(), + reason="DAGMC CAD geometry is not enabled.") +@skip_if_no_dagmc def test_umesh(run_in_tmpdir, request): """Performs a minimal UnstructuredMesh simulation, reads in the resulting statepoint file and writes the mesh data to vtk and hdf files. It is From aa1383571ef7b0136c5da8fccae6d4e606235442 Mon Sep 17 00:00:00 2001 From: Jon Shimwell Date: Thu, 9 Jan 2025 10:28:18 +0100 Subject: [PATCH 05/22] review improvements by mwestphal --- openmc/mesh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 019c44c8b46..f681ac79c54 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2369,7 +2369,7 @@ def write_data_to_vtk( Parameters ---------- filename : str or pathlib.Path - Name of the VTK file to write. If the filename ends in '.hdf' then a VTKHDF + Name of the VTK file to write. If the filename ends in '.vtkhdf' then a VTKHDF format file will be written and can be opened with Paraview versions 5.13.0 and above, if the filename ends in '.vtk' then a .vtk file will be written. datasets : dict @@ -2381,7 +2381,7 @@ def write_data_to_vtk( volume of the mesh elements """ - if Path(filename).suffix == ".hdf": + if Path(filename).suffix == ".vtkhdf": def append_dataset(dset, array): """Convenience function to append data to an HDF5 dataset""" @@ -2558,7 +2558,7 @@ def append_dataset(dset, array): else: raise ValueError( - f"Unsupported file extension for '{filename}'. Extension must be '.hdf' or '.vtk'." + f"Unsupported file extension for '{filename}'. Extension must be '.vtkhdf' or '.vtk'." ) @classmethod From 177e8dc8b9ede14a6a4aa3bad5cb3e07a010db93 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 9 Jan 2025 14:45:24 +0100 Subject: [PATCH 06/22] changed .hdf to .vtkhdf --- tests/unit_tests/test_mesh.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index c87b193b9b2..8b4ad89faa1 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -509,10 +509,10 @@ def test_umesh(run_in_tmpdir, request): } umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtk") - umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.hdf") + umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtkhdf") with pytest.raises(ValueError, match="Unsupported file extension"): - # Supported file extensions are vtk or hdf, not hdf5, so this should raise an error + # Supported file extensions are vtk or vtkhdf, not hdf5, so this should raise an error umesh_from_sp.write_data_to_vtk( datasets=datasets, filename="test_mesh.hdf5", @@ -521,8 +521,8 @@ def test_umesh(run_in_tmpdir, request): # The shape of the data should match the shape of the mesh, so this should raise an error umesh_from_sp.write_data_to_vtk( datasets={'incorrectly_shaped_data': np.array(([1,2,3]))}, - filename="test_mesh.hdf", + filename="test_mesh.vtkhdf", ) assert Path("test_mesh.vtk").exists() - assert Path("test_mesh.hdf").exists() + assert Path("test_mesh.vtkhdf").exists() From ba4c097db9cdedc11b1f36db6a1dfedc7fd53b7c Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 9 Jan 2025 14:46:20 +0100 Subject: [PATCH 07/22] changed .hdf to .vtkhdf --- tests/unit_tests/test_mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 8b4ad89faa1..9ae5ae6f649 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -458,7 +458,7 @@ def test_mesh_get_homogenized_materials(): @skip_if_no_dagmc def test_umesh(run_in_tmpdir, request): """Performs a minimal UnstructuredMesh simulation, reads in the resulting - statepoint file and writes the mesh data to vtk and hdf files. It is + statepoint file and writes the mesh data to vtk and vtkhdf files. It is necessary to read in the unstructured mesh from a statepoint file to ensure it has all the required attributes """ From 6f397fdd55da155c6426f844e9e972134ffe0aa3 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 4 Feb 2025 12:13:00 +0100 Subject: [PATCH 08/22] using vtk mob mesh not vtk output mesh Co-authored-by: rherrero-pf <156206440+rherrero-pf@users.noreply.github.com> --- .../test_mesh_dagmc_tets.vtk | 376 +++++++----------- 1 file changed, 141 insertions(+), 235 deletions(-) diff --git a/tests/regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk b/tests/regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk index daad9e15a75..2ddd228cf13 100644 --- a/tests/regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk +++ b/tests/regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk @@ -1,8 +1,8 @@ -# vtk DataFile Version 3.0 -MOAB 5.5.1 +# vtk DataFile Version 2.0 +made_with_cad_to_dagmc_package, Created by Gmsh 4.12.1 ASCII DATASET UNSTRUCTURED_GRID -POINTS 84 double +POINTS 14 double -0.5 -0.5 0.5 -0.5 -0.5 -0.5 -0.5 0.5 0.5 @@ -17,237 +17,143 @@ POINTS 84 double 0 0.5 0 0 0 -0.5 0 0 0.5 --0.5 -0.5 0.5 --0.5 -0.5 -0.5 --0.5 0.5 0.5 --0.5 0.5 -0.5 -0.5 -0.5 0.5 -0.5 -0.5 -0.5 -0.5 0.5 0.5 -0.5 0.5 -0.5 --0.5 0 0 -0.5 0 0 -0 -0.5 0 -0 0.5 0 -0 0 -0.5 -0 0 0.5 --0.5 -0.5 0.5 --0.5 -0.5 -0.5 --0.5 0.5 0.5 --0.5 0.5 -0.5 -0.5 -0.5 0.5 -0.5 -0.5 -0.5 -0.5 0.5 0.5 -0.5 0.5 -0.5 --0.5 0 0 -0.5 0 0 -0 -0.5 0 -0 0.5 0 -0 0 -0.5 -0 0 0.5 --0.5 -0.5 0.5 --0.5 -0.5 -0.5 --0.5 0.5 0.5 --0.5 0.5 -0.5 -0.5 -0.5 0.5 -0.5 -0.5 -0.5 -0.5 0.5 0.5 -0.5 0.5 -0.5 --0.5 0 0 -0.5 0 0 -0 -0.5 0 -0 0.5 0 -0 0 -0.5 -0 0 0.5 --0.5 -0.5 0.5 --0.5 -0.5 -0.5 --0.5 0.5 0.5 --0.5 0.5 -0.5 -0.5 -0.5 0.5 -0.5 -0.5 -0.5 -0.5 0.5 0.5 -0.5 0.5 -0.5 --0.5 0 0 -0.5 0 0 -0 -0.5 0 -0 0.5 0 -0 0 -0.5 -0 0 0.5 --0.5 -0.5 0.5 --0.5 -0.5 -0.5 --0.5 0.5 0.5 --0.5 0.5 -0.5 -0.5 -0.5 0.5 -0.5 -0.5 -0.5 -0.5 0.5 0.5 -0.5 0.5 -0.5 --0.5 0 0 -0.5 0 0 -0 -0.5 0 -0 0.5 0 -0 0 -0.5 -0 0 0.5 -CELLS 24 96 + +CELLS 68 268 +1 0 +1 1 +1 2 +1 3 +1 4 +1 5 +1 6 +1 7 +2 1 0 +2 0 2 +2 3 2 +2 1 3 +2 5 4 +2 4 6 +2 7 6 +2 5 7 +2 1 5 +2 0 4 +2 3 7 +2 2 6 3 1 0 8 +3 0 2 8 +3 3 1 8 3 2 3 8 -3 8 3 1 -3 8 0 2 -3 19 23 18 -3 20 23 21 -3 23 19 21 -3 23 20 18 -3 28 29 38 -3 33 32 38 -3 38 32 28 -3 38 29 33 -3 44 53 45 -3 49 53 48 -3 53 44 48 -3 53 49 45 -3 68 57 59 -3 68 63 61 -3 59 63 68 -3 61 57 68 -3 83 72 70 -3 83 74 76 -3 72 83 76 -3 74 83 70 -CELL_TYPES 24 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5 -POINT_DATA 84 -SCALARS GLOBAL_ID int 1 -LOOKUP_TABLE default --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 -CELL_DATA 24 -SCALARS GLOBAL_ID int 1 -LOOKUP_TABLE default --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 --1 +3 5 9 4 +3 4 9 6 +3 7 9 5 +3 6 9 7 +3 0 1 10 +3 4 0 10 +3 1 5 10 +3 5 4 10 +3 2 11 3 +3 6 11 2 +3 3 11 7 +3 7 11 6 +3 1 3 12 +3 5 1 12 +3 3 7 12 +3 7 5 12 +3 0 13 2 +3 4 13 0 +3 2 13 6 +3 6 13 4 +4 13 8 12 10 +4 11 8 12 13 +4 12 11 13 9 +4 10 12 13 9 +4 12 3 11 7 +4 13 2 8 0 +4 8 12 1 3 +4 11 3 8 2 +4 10 8 1 0 +4 0 10 13 4 +4 1 12 10 5 +4 11 2 13 6 +4 4 9 13 6 +4 6 9 11 7 +4 10 9 4 5 +4 7 9 12 5 +4 3 8 12 11 +4 1 12 8 10 +4 8 11 2 13 +4 10 13 8 0 +4 13 10 9 4 +4 11 13 9 6 +4 12 11 9 7 +4 9 10 12 5 + +CELL_TYPES 68 +1 +1 +1 +1 +1 +1 +1 +1 +3 +3 +3 +3 +3 +3 +3 +3 +3 +3 +3 +3 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 +10 From 2426da89eb6df91f3dc9b97a02233b16eff1895c Mon Sep 17 00:00:00 2001 From: Jon Shimwell Date: Tue, 4 Feb 2025 12:20:55 +0100 Subject: [PATCH 09/22] moved data check earlier to avoid empty file creation --- openmc/mesh.py | 18 +++++++++++------- tests/unit_tests/test_mesh.py | 6 +++--- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index f681ac79c54..88e69573ac7 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2383,6 +2383,8 @@ def write_data_to_vtk( if Path(filename).suffix == ".vtkhdf": + print('Writing VTKHDF file') + def append_dataset(dset, array): """Convenience function to append data to an HDF5 dataset""" origLen = dset.shape[0] @@ -2416,6 +2418,15 @@ def append_dataset(dset, array): 0, self.n_elements * points_per_cell + 1, points_per_cell ) + for name, data in datasets.items(): + + if data.shape != self.dimension: + raise ValueError( + f'Cannot apply dataset "{name}" with ' + f"shape {data.shape} to mesh {self.id} " + f"with dimensions {self.dimension}" + ) + with h5py.File(filename, "w") as f: root = f.create_group("VTKHDF") @@ -2459,13 +2470,6 @@ def append_dataset(dset, array): for name, data in datasets.items(): - if data.shape != self.dimension: - raise ValueError( - f'Cannot apply dataset "{name}" with ' - f"shape {data.shape} to mesh {self.id} " - f"with dimensions {self.dimension}" - ) - cell_data_group.create_dataset( name, (0,), maxshape=(None,), dtype="float64", chunks=True ) diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 9ae5ae6f649..6d06c4ba0f6 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -456,7 +456,7 @@ def test_mesh_get_homogenized_materials(): reason="DAGMC CAD geometry is not enabled.") @skip_if_no_dagmc -def test_umesh(run_in_tmpdir, request): +def test_umesh(request): """Performs a minimal UnstructuredMesh simulation, reads in the resulting statepoint file and writes the mesh data to vtk and vtkhdf files. It is necessary to read in the unstructured mesh from a statepoint file to ensure @@ -508,8 +508,8 @@ def test_umesh(run_in_tmpdir, request): "std_dev": my_tally.std_dev.flatten() } - umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtk") umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtkhdf") + umesh_from_sp.write_data_to_vtk(datasets=datasets, filename="test_mesh.vtk") with pytest.raises(ValueError, match="Unsupported file extension"): # Supported file extensions are vtk or vtkhdf, not hdf5, so this should raise an error @@ -521,7 +521,7 @@ def test_umesh(run_in_tmpdir, request): # The shape of the data should match the shape of the mesh, so this should raise an error umesh_from_sp.write_data_to_vtk( datasets={'incorrectly_shaped_data': np.array(([1,2,3]))}, - filename="test_mesh.vtkhdf", + filename="test_mesh_incorrect_shape.vtkhdf", ) assert Path("test_mesh.vtk").exists() From 1bb33e619e2b32ee702a084d4fade9be21e72068 Mon Sep 17 00:00:00 2001 From: shimwell Date: Sun, 9 Mar 2025 02:08:52 +0100 Subject: [PATCH 10/22] minimal git diff --- openmc/mesh.py | 1 - 1 file changed, 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 1b41a5828ed..2b5aa62cdf4 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2576,7 +2576,6 @@ def append_dataset(dset, array): ) for name, data in datasets.items(): - if data.shape != self.dimension: raise ValueError( f'Cannot apply dataset "{name}" with ' From d746783c623644c2e8002f4a0dbec07a6fe789de Mon Sep 17 00:00:00 2001 From: shimwell Date: Sun, 9 Mar 2025 12:31:27 +0100 Subject: [PATCH 11/22] removed print line --- openmc/mesh.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 2b5aa62cdf4..0a9cdac207e 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2540,8 +2540,6 @@ def write_data_to_vtk( if Path(filename).suffix == ".vtkhdf": - print('Writing VTKHDF file') - def append_dataset(dset, array): """Convenience function to append data to an HDF5 dataset""" origLen = dset.shape[0] From bdc3a6ead7831917d0931a9c8feba822bbb71d13 Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 10 Mar 2025 12:12:20 +0100 Subject: [PATCH 12/22] added value error for unsupported file suffix --- openmc/mesh.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 0a9cdac207e..006e66656be 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2632,7 +2632,7 @@ def append_dataset(dset, array): data /= self.volumes append_dataset(cell_data_group[name], data) - else: + elif Path(filename).suffix == ".vtk" or Path(filename).suffix == ".vtu": from vtkmodules.util import numpy_support from vtkmodules import vtkCommonCore @@ -2724,6 +2724,12 @@ def append_dataset(dset, array): writer.SetInputData(grid) writer.Write() + + else: + raise ValueError( + "Unsupported file extension, The filename must end with " + "'.vtkhdf', '.vtu' or '.vtk'" + ) @classmethod def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): From 70203b7e02d002a0009e4d36e93132331ee9b306 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 24 Oct 2025 12:06:29 +0200 Subject: [PATCH 13/22] [skip ci] improved comments from code review Co-authored-by: Patrick Shriwise --- openmc/mesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 5590608fcb3..2c5de7fca37 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2698,7 +2698,7 @@ def append_dataset(dset, array): raise NotImplemented("VTKHDF output is only supported for MOAB meshes") # the self.connectivity contains an arrays of length 8, in the case of - # DAGMC tetrahedra mesh elements, the last 4 values are -1 and can be removed + # MOAB tetrahedra mesh elements, the last 4 values are -1 and can be removed trimmed_connectivity = [] for cell in self.connectivity: # Find the index of the first -1 value, if any @@ -2713,7 +2713,7 @@ def append_dataset(dset, array): trimmed_connectivity, dtype="int32" ).flatten() - # DAGMC supports tet meshes only so we know it has 4 points per cell + # MOAB meshes supports tet elements only so we know it has 4 points per cell points_per_cell = 4 # offsets are the indices of the first point of each cell in the array of points From 098d79003bf218d2162230791b3716c338b38c69 Mon Sep 17 00:00:00 2001 From: Jon Shimwell Date: Sun, 26 Oct 2025 17:26:22 +0000 Subject: [PATCH 14/22] write_data_to_vtk split into 2 functions --- openmc/mesh.py | 377 ++++++++++++++++++++++++++----------------------- 1 file changed, 201 insertions(+), 176 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 2c5de7fca37..e2414e5b291 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2656,7 +2656,8 @@ def write_vtk_mesh(self, **kwargs): warnings.warn( "The 'UnstructuredMesh.write_vtk_mesh' method has been renamed " "to 'write_data_to_vtk' and will be removed in a future version " - " of OpenMC.", FutureWarning + " of OpenMC.", + FutureWarning, ) self.write_data_to_vtk(**kwargs) @@ -2688,220 +2689,244 @@ def write_data_to_vtk( if Path(filename).suffix == ".vtkhdf": - def append_dataset(dset, array): - """Convenience function to append data to an HDF5 dataset""" - origLen = dset.shape[0] - dset.resize(origLen + array.shape[0], axis=0) - dset[origLen:] = array - - if self.library != "moab": - raise NotImplemented("VTKHDF output is only supported for MOAB meshes") - - # the self.connectivity contains an arrays of length 8, in the case of - # MOAB tetrahedra mesh elements, the last 4 values are -1 and can be removed - trimmed_connectivity = [] - for cell in self.connectivity: - # Find the index of the first -1 value, if any - first_negative_index = np.where(cell == -1)[0] - if first_negative_index.size > 0: - # Slice the array up to the first -1 value - trimmed_connectivity.append(cell[: first_negative_index[0]]) - else: - # No -1 values, append the whole cell - trimmed_connectivity.append(cell) - trimmed_connectivity = np.array( - trimmed_connectivity, dtype="int32" - ).flatten() - - # MOAB meshes supports tet elements only so we know it has 4 points per cell - points_per_cell = 4 - - # offsets are the indices of the first point of each cell in the array of points - offsets = np.arange( - 0, self.n_elements * points_per_cell + 1, points_per_cell + self._write_data_to_vtk_vtk_format( + filename=filename, + datasets=datasets, + volume_normalization=volume_normalization, ) - for name, data in datasets.items(): - if data.shape != self.dimension: - raise ValueError( - f'Cannot apply dataset "{name}" with ' - f"shape {data.shape} to mesh {self.id} " - f"with dimensions {self.dimension}" - ) + elif Path(filename).suffix == ".vtk" or Path(filename).suffix == ".vtu": - with h5py.File(filename, "w") as f: + self._write_data_to_vtk_hdf5_format( + filename=filename, + datasets=datasets, + volume_normalization=volume_normalization, + ) - root = f.create_group("VTKHDF") - root.attrs["Version"] = (2, 1) - ascii_type = "UnstructuredGrid".encode("ascii") - root.attrs.create( - "Type", - ascii_type, - dtype=h5py.string_dtype("ascii", len(ascii_type)), - ) + else: + raise ValueError( + "Unsupported file extension, The filename must end with " + "'.vtkhdf', '.vtu' or '.vtk'" + ) - # create hdf5 file structure - root.create_dataset( - "NumberOfPoints", (0,), maxshape=(None,), dtype="i8" - ) - root.create_dataset("Types", (0,), maxshape=(None,), dtype="uint8") - root.create_dataset("Points", (0, 3), maxshape=(None, 3), dtype="f") - root.create_dataset( - "NumberOfConnectivityIds", (0,), maxshape=(None,), dtype="i8" - ) - root.create_dataset("NumberOfCells", (0,), maxshape=(None,), dtype="i8") - root.create_dataset("Offsets", (0,), maxshape=(None,), dtype="i8") - root.create_dataset("Connectivity", (0,), maxshape=(None,), dtype="i8") - - append_dataset(root["NumberOfPoints"], np.array([len(self.vertices)])) - append_dataset(root["Points"], self.vertices) - append_dataset( - root["NumberOfConnectivityIds"], - np.array([len(trimmed_connectivity)]), - ) - append_dataset(root["Connectivity"], trimmed_connectivity) - append_dataset(root["NumberOfCells"], np.array([self.n_elements])) - append_dataset(root["Offsets"], offsets) + def _write_data_to_vtk_vtk_format( + self, + filename: PathLike | None = None, + datasets: dict | None = None, + volume_normalization: bool = True, + ): + from vtkmodules.util import numpy_support + from vtkmodules import vtkCommonCore + from vtkmodules import vtkCommonDataModel + from vtkmodules import vtkIOLegacy + from vtkmodules import vtkIOXML - # VTK_TETRA type is known as DAGMC only supports tet meshes - append_dataset( - root["Types"], np.full(self.n_elements, 10, dtype="uint8") - ) + if self.connectivity is None or self.vertices is None: + raise RuntimeError("This mesh has not been loaded from a statepoint file.") - cell_data_group = root.create_group("CellData") + if filename is None: + filename = f"mesh_{self.id}.vtk" - for name, data in datasets.items(): + if Path(filename).suffix == ".vtk": + writer = vtkIOLegacy.vtkUnstructuredGridWriter() - cell_data_group.create_dataset( - name, (0,), maxshape=(None,), dtype="float64", chunks=True - ) + elif Path(filename).suffix == ".vtu": + writer = vtkIOXML.vtkXMLUnstructuredGridWriter() + writer.SetCompressorTypeToZLib() + writer.SetDataModeToBinary() - if volume_normalization: - data /= self.volumes - append_dataset(cell_data_group[name], data) + writer.SetFileName(str(filename)) - elif Path(filename).suffix == ".vtk" or Path(filename).suffix == ".vtu": - - from vtkmodules.util import numpy_support - from vtkmodules import vtkCommonCore - from vtkmodules import vtkCommonDataModel - from vtkmodules import vtkIOLegacy - from vtkmodules import vtkIOXML - - if self.connectivity is None or self.vertices is None: - raise RuntimeError( - "This mesh has not been loaded from a statepoint file." - ) + grid = vtkCommonDataModel.vtkUnstructuredGrid() + + points = vtkCommonCore.vtkPoints() + points.SetData(numpy_support.numpy_to_vtk(self.vertices)) + grid.SetPoints(points) + + n_skipped = 0 + for elem_type, conn in zip(self.element_types, self.connectivity): + if elem_type == self._LINEAR_TET: + elem = vtkCommonDataModel.vtkTetra() + elif elem_type == self._LINEAR_HEX: + elem = vtkCommonDataModel.vtkHexahedron() + elif elem_type == self._UNSUPPORTED_ELEM: + n_skipped += 1 + continue + else: + raise RuntimeError(f"Invalid element type {elem_type} found") - if filename is None: - filename = f"mesh_{self.id}.vtk" + for i, c in enumerate(conn): + if c == -1: + break + elem.GetPointIds().SetId(i, c) - if Path(filename).suffix == ".vtk": - writer = vtkIOLegacy.vtkUnstructuredGridWriter() + grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds()) - elif Path(filename).suffix == ".vtu": - writer = vtkIOXML.vtkXMLUnstructuredGridWriter() - writer.SetCompressorTypeToZLib() - writer.SetDataModeToBinary() + if n_skipped > 0: + warnings.warn( + f"{n_skipped} elements were not written because " + "they are not of type linear tet/hex" + ) - writer.SetFileName(str(filename)) + # check that datasets are the correct size + datasets_out = [] + if datasets is not None: + for name, data in datasets.items(): + if data.shape != self.dimension: + raise ValueError( + f'Cannot apply dataset "{name}" with ' + f"shape {data.shape} to mesh {self.id} " + f"with dimensions {self.dimension}" + ) - grid = vtkCommonDataModel.vtkUnstructuredGrid() + if volume_normalization: + for name, data in datasets.items(): + if np.issubdtype(data.dtype, np.integer): + warnings.warn( + f'Integer data set "{name}" will ' + "not be volume-normalized." + ) + continue + data /= self.volumes - points = vtkCommonCore.vtkPoints() - points.SetData(numpy_support.numpy_to_vtk(self.vertices)) - grid.SetPoints(points) + # add data to the mesh + for name, data in datasets.items(): + datasets_out.append(data) + arr = vtkCommonCore.vtkDoubleArray() + arr.SetName(name) + arr.SetNumberOfTuples(data.size) - n_skipped = 0 - for elem_type, conn in zip(self.element_types, self.connectivity): - if elem_type == self._LINEAR_TET: - elem = vtkCommonDataModel.vtkTetra() - elif elem_type == self._LINEAR_HEX: - elem = vtkCommonDataModel.vtkHexahedron() - elif elem_type == self._UNSUPPORTED_ELEM: - n_skipped += 1 - continue - else: - raise RuntimeError(f"Invalid element type {elem_type} found") + for i in range(data.size): + arr.SetTuple1(i, data.flat[i]) + grid.GetCellData().AddArray(arr) - for i, c in enumerate(conn): - if c == -1: - break - elem.GetPointIds().SetId(i, c) + writer.SetInputData(grid) - grid.InsertNextCell(elem.GetCellType(), elem.GetPointIds()) + writer.Write() - if n_skipped > 0: - warnings.warn( - f"{n_skipped} elements were not written because " - "they are not of type linear tet/hex" + def _write_data_to_vtk_hdf_format( + self, + filename: PathLike | None = None, + datasets: dict | None = None, + volume_normalization: bool = True, + ): + def append_dataset(dset, array): + """Convenience function to append data to an HDF5 dataset""" + origLen = dset.shape[0] + dset.resize(origLen + array.shape[0], axis=0) + dset[origLen:] = array + + if self.library != "moab": + raise NotImplemented("VTKHDF output is only supported for MOAB meshes") + + # the self.connectivity contains an arrays of length 8, in the case of + # MOAB tetrahedra mesh elements, the last 4 values are -1 and can be removed + trimmed_connectivity = [] + for cell in self.connectivity: + # Find the index of the first -1 value, if any + first_negative_index = np.where(cell == -1)[0] + if first_negative_index.size > 0: + # Slice the array up to the first -1 value + trimmed_connectivity.append(cell[: first_negative_index[0]]) + else: + # No -1 values, append the whole cell + trimmed_connectivity.append(cell) + trimmed_connectivity = np.array(trimmed_connectivity, dtype="int32").flatten() + + # MOAB meshes supports tet elements only so we know it has 4 points per cell + points_per_cell = 4 + + # offsets are the indices of the first point of each cell in the array of points + offsets = np.arange(0, self.n_elements * points_per_cell + 1, points_per_cell) + + for name, data in datasets.items(): + if data.shape != self.dimension: + raise ValueError( + f'Cannot apply dataset "{name}" with ' + f"shape {data.shape} to mesh {self.id} " + f"with dimensions {self.dimension}" ) - # check that datasets are the correct size - datasets_out = [] - if datasets is not None: - for name, data in datasets.items(): - if data.shape != self.dimension: - raise ValueError( - f'Cannot apply dataset "{name}" with ' - f"shape {data.shape} to mesh {self.id} " - f"with dimensions {self.dimension}" - ) + with h5py.File(filename, "w") as f: - if volume_normalization: - for name, data in datasets.items(): - if np.issubdtype(data.dtype, np.integer): - warnings.warn( - f'Integer data set "{name}" will ' - "not be volume-normalized." - ) - continue - data /= self.volumes - - # add data to the mesh - for name, data in datasets.items(): - datasets_out.append(data) - arr = vtkCommonCore.vtkDoubleArray() - arr.SetName(name) - arr.SetNumberOfTuples(data.size) + root = f.create_group("VTKHDF") + vtk_file_format_version = (2, 1) + root.attrs["Version"] = vtk_file_format_version + ascii_type = "UnstructuredGrid".encode("ascii") + root.attrs.create( + "Type", + ascii_type, + dtype=h5py.string_dtype("ascii", len(ascii_type)), + ) - for i in range(data.size): - arr.SetTuple1(i, data.flat[i]) - grid.GetCellData().AddArray(arr) + # create hdf5 file structure + root.create_dataset("NumberOfPoints", (0,), maxshape=(None,), dtype="i8") + root.create_dataset("Types", (0,), maxshape=(None,), dtype="uint8") + root.create_dataset("Points", (0, 3), maxshape=(None, 3), dtype="f") + root.create_dataset( + "NumberOfConnectivityIds", (0,), maxshape=(None,), dtype="i8" + ) + root.create_dataset("NumberOfCells", (0,), maxshape=(None,), dtype="i8") + root.create_dataset("Offsets", (0,), maxshape=(None,), dtype="i8") + root.create_dataset("Connectivity", (0,), maxshape=(None,), dtype="i8") + + append_dataset(root["NumberOfPoints"], np.array([len(self.vertices)])) + append_dataset(root["Points"], self.vertices) + append_dataset( + root["NumberOfConnectivityIds"], + np.array([len(trimmed_connectivity)]), + ) + append_dataset(root["Connectivity"], trimmed_connectivity) + append_dataset(root["NumberOfCells"], np.array([self.n_elements])) + append_dataset(root["Offsets"], offsets) + + # VTK_TETRA type is known as DAGMC only supports tet meshes + element_type = 10 # VTK_TETRA + append_dataset( + root["Types"], np.full(self.n_elements, element_type, dtype="uint8") + ) - writer.SetInputData(grid) + cell_data_group = root.create_group("CellData") - writer.Write() - - else: - raise ValueError( - "Unsupported file extension, The filename must end with " - "'.vtkhdf', '.vtu' or '.vtk'" - ) + for name, data in datasets.items(): + + cell_data_group.create_dataset( + name, (0,), maxshape=(None,), dtype="float64", chunks=True + ) + + if volume_normalization: + data /= self.volumes + append_dataset(cell_data_group[name], data) @classmethod def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): - filename = group['filename'][()].decode() - library = group['library'][()].decode() - if 'options' in group.attrs: - options = group.attrs['options'].decode() + filename = group["filename"][()].decode() + library = group["library"][()].decode() + if "options" in group.attrs: + options = group.attrs["options"].decode() else: options = None - mesh = cls(filename=filename, library=library, mesh_id=mesh_id, name=name, options=options) + mesh = cls( + filename=filename, + library=library, + mesh_id=mesh_id, + name=name, + options=options, + ) mesh._has_statepoint_data = True - vol_data = group['volumes'][()] + vol_data = group["volumes"][()] mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],)) mesh.n_elements = mesh.volumes.size - vertices = group['vertices'][()] + vertices = group["vertices"][()] mesh._vertices = vertices.reshape((-1, 3)) - connectivity = group['connectivity'][()] + connectivity = group["connectivity"][()] mesh._connectivity = connectivity.reshape((-1, 8)) - mesh._element_types = group['element_types'][()] + mesh._element_types = group["element_types"][()] - if 'length_multiplier' in group: - mesh.length_multiplier = group['length_multiplier'][()] + if "length_multiplier" in group: + mesh.length_multiplier = group["length_multiplier"][()] return mesh @@ -2920,7 +2945,7 @@ def to_xml_element(self): element.set("library", self._library) if self.options is not None: - element.set('options', self.options) + element.set("options", self.options) subelement = ET.SubElement(element, "filename") subelement.text = str(self.filename) From 698d1c69a6203b71c9b2620d1789bc9ea364e4ac Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 30 Oct 2025 17:57:33 +0100 Subject: [PATCH 15/22] [skip ci] Apply suggestions from code review Co-authored-by: Patrick Shriwise --- openmc/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index e2414e5b291..71f03f08e84 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2709,7 +2709,7 @@ def write_data_to_vtk( "'.vtkhdf', '.vtu' or '.vtk'" ) - def _write_data_to_vtk_vtk_format( + def _write_data_to_vtk_ascii_format( self, filename: PathLike | None = None, datasets: dict | None = None, From da337a05a1f2692e4e84c5d1174d4d4cba000a85 Mon Sep 17 00:00:00 2001 From: shimwell Date: Thu, 30 Oct 2025 18:06:36 +0100 Subject: [PATCH 16/22] review comments on vtkhdf --- openmc/mesh.py | 9 ++++----- tests/unit_tests/test_mesh.py | 5 ++--- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 71f03f08e84..021de58559b 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2447,6 +2447,7 @@ class UnstructuredMesh(MeshBase): _UNSUPPORTED_ELEM = -1 _LINEAR_TET = 0 _LINEAR_HEX = 1 + _VTK_TETRA = 10 # VTK_TETRA type is known as DAGMC only supports tet meshes def __init__(self, filename: PathLike, library: str, mesh_id: int | None = None, name: str = '', length_multiplier: float = 1.0, @@ -2689,7 +2690,7 @@ def write_data_to_vtk( if Path(filename).suffix == ".vtkhdf": - self._write_data_to_vtk_vtk_format( + self._write_data_to_vtk_ascii_format( filename=filename, datasets=datasets, volume_normalization=volume_normalization, @@ -2804,7 +2805,7 @@ def _write_data_to_vtk_ascii_format( writer.Write() - def _write_data_to_vtk_hdf_format( + def _write_data_to_vtk_hdf5_format( self, filename: PathLike | None = None, datasets: dict | None = None, @@ -2880,10 +2881,8 @@ def append_dataset(dset, array): append_dataset(root["NumberOfCells"], np.array([self.n_elements])) append_dataset(root["Offsets"], offsets) - # VTK_TETRA type is known as DAGMC only supports tet meshes - element_type = 10 # VTK_TETRA append_dataset( - root["Types"], np.full(self.n_elements, element_type, dtype="uint8") + root["Types"], np.full(self.n_elements, _VTK_TETRA, dtype="uint8") ) cell_data_group = root.create_group("CellData") diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 841d9159a2a..20872f98c5f 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -526,9 +526,8 @@ def test_umesh(request): statepoint_file = my_model.run() - statepoint = openmc.StatePoint(statepoint_file) - - my_tally = statepoint.get_tally(name="test_tally") + with openmc.StatePoint(statepoint_file) as statepoint: + my_tally = statepoint.get_tally(name="test_tally") umesh_from_sp = statepoint.meshes[1] From 0cd21ed2ed5072444906ae3979690337fe60651f Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 31 Oct 2025 09:32:51 +0100 Subject: [PATCH 17/22] corrected function name to vtkhdf --- openmc/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 021de58559b..b14ed2a3d19 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2690,7 +2690,7 @@ def write_data_to_vtk( if Path(filename).suffix == ".vtkhdf": - self._write_data_to_vtk_ascii_format( + self._write_data_to_vtk_hdf5_format( filename=filename, datasets=datasets, volume_normalization=volume_normalization, From 27157b613aad2d86ae9721bb63f423e14c8f858d Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 3 Nov 2025 15:20:11 +0100 Subject: [PATCH 18/22] added missing self --- openmc/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index b14ed2a3d19..6d022fe0734 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2882,7 +2882,7 @@ def append_dataset(dset, array): append_dataset(root["Offsets"], offsets) append_dataset( - root["Types"], np.full(self.n_elements, _VTK_TETRA, dtype="uint8") + root["Types"], np.full(self.n_elements, self._VTK_TETRA, dtype="uint8") ) cell_data_group = root.create_group("CellData") From 226d1812003514800abd48378640a6abc9ab6d32 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 3 Nov 2025 21:55:50 -0600 Subject: [PATCH 19/22] Updates to test_mesh. Adding symbolic link to file for less verbose filepath. --- tests/unit_tests/test_mesh.py | 33 +++++++++-------------- tests/unit_tests/test_mesh_dagmc_tets.vtk | 1 + 2 files changed, 13 insertions(+), 21 deletions(-) create mode 120000 tests/unit_tests/test_mesh_dagmc_tets.vtk diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 20872f98c5f..84c4a8570bd 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -485,24 +485,23 @@ def test_umesh(run_in_tmpdir, simple_umesh, export_type): simple_umesh.write_data_to_vtk(datasets={'mean': ref_data[:-2]}, filename=filename) @pytest.mark.skipif(not openmc.lib._dagmc_enabled(), reason="DAGMC not enabled.") -def test_umesh(request): +def test_write_vtkhdf(request, run_in_tmpdir): """Performs a minimal UnstructuredMesh simulation, reads in the resulting statepoint file and writes the mesh data to vtk and vtkhdf files. It is necessary to read in the unstructured mesh from a statepoint file to ensure it has all the required attributes """ + model = openmc.Model() surf1 = openmc.Sphere(R=1000.0, boundary_type="vacuum") cell1 = openmc.Cell(region=-surf1) - my_geometry = openmc.Geometry([cell1]) + model.geometry = openmc.Geometry([cell1]) umesh = openmc.UnstructuredMesh( - request.path.parent.parent - / "regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk", + request.path.parent / "test_mesh_dagmc_tets.vtk", "moab", + mesh_id = 1 ) - # setting ID to make it easier to get the mesh from the statepoint later - umesh.id = 1 mesh_filter = openmc.MeshFilter(umesh) # Create flux mesh tally to score alpha production @@ -510,26 +509,18 @@ def test_umesh(request): mesh_tally.filters = [mesh_filter] mesh_tally.scores = ["flux"] - tallies = openmc.Tallies([mesh_tally]) - - my_source = openmc.IndependentSource() + model.tallies = [mesh_tally] - settings = openmc.Settings() - settings.run_mode = "fixed source" - settings.batches = 2 - settings.particles = 10 - settings.source = my_source - - my_model = openmc.Model( - materials=None, geometry=my_geometry, settings=settings, tallies=tallies - ) + model.settings.run_mode = "fixed source" + model.settings.batches = 2 + model.settings.particles = 10 - statepoint_file = my_model.run() + statepoint_file = model.run() with openmc.StatePoint(statepoint_file) as statepoint: my_tally = statepoint.get_tally(name="test_tally") - umesh_from_sp = statepoint.meshes[1] + umesh_from_sp = statepoint.meshes[umesh.id] datasets={ "mean": my_tally.mean.flatten(), @@ -554,7 +545,7 @@ def test_umesh(request): assert Path("test_mesh.vtk").exists() assert Path("test_mesh.vtkhdf").exists() - + def test_mesh_get_homogenized_materials(): """Test the get_homogenized_materials method""" diff --git a/tests/unit_tests/test_mesh_dagmc_tets.vtk b/tests/unit_tests/test_mesh_dagmc_tets.vtk new file mode 120000 index 00000000000..9f7000175d7 --- /dev/null +++ b/tests/unit_tests/test_mesh_dagmc_tets.vtk @@ -0,0 +1 @@ +../regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk \ No newline at end of file From b8a72730fc0fbb9cb0f42f0a81ec9bc623b440e0 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 3 Nov 2025 23:26:09 -0600 Subject: [PATCH 20/22] Correct call to internal method for non-hdf file formats --- openmc/mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 6d022fe0734..b87cd1a0948 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2698,7 +2698,7 @@ def write_data_to_vtk( elif Path(filename).suffix == ".vtk" or Path(filename).suffix == ".vtu": - self._write_data_to_vtk_hdf5_format( + self._write_data_to_vtk_ascii_format( filename=filename, datasets=datasets, volume_normalization=volume_normalization, From b4a815914ca43248aa3f1f749e9f964e3fcb9725 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 4 Nov 2025 10:51:28 -0600 Subject: [PATCH 21/22] Comment updates and simplifying PR diff --- openmc/mesh.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index b87cd1a0948..54d796f59c2 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2447,7 +2447,7 @@ class UnstructuredMesh(MeshBase): _UNSUPPORTED_ELEM = -1 _LINEAR_TET = 0 _LINEAR_HEX = 1 - _VTK_TETRA = 10 # VTK_TETRA type is known as DAGMC only supports tet meshes + _VTK_TETRA = 10 def __init__(self, filename: PathLike, library: str, mesh_id: int | None = None, name: str = '', length_multiplier: float = 1.0, @@ -2820,8 +2820,9 @@ def append_dataset(dset, array): if self.library != "moab": raise NotImplemented("VTKHDF output is only supported for MOAB meshes") - # the self.connectivity contains an arrays of length 8, in the case of - # MOAB tetrahedra mesh elements, the last 4 values are -1 and can be removed + # the self.connectivity contains arrays of length 8 to support hex + # elements as well, in the case of tetrahedra mesh elements, the + # last 4 values are -1 and are removed trimmed_connectivity = [] for cell in self.connectivity: # Find the index of the first -1 value, if any @@ -2902,7 +2903,7 @@ def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): filename = group["filename"][()].decode() library = group["library"][()].decode() if "options" in group.attrs: - options = group.attrs["options"].decode() + options = group.attrs['options'].decode() else: options = None From 3c71decac2a17498f13bbe181e19d552e2dd8a68 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Nov 2025 07:01:15 -0600 Subject: [PATCH 22/22] Small fixes --- openmc/mesh.py | 4 ++-- tests/unit_tests/test_mesh.py | 8 +++++++- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 54d796f59c2..bc47b2b8f0e 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -2678,7 +2678,7 @@ def write_data_to_vtk( filename : str or pathlib.Path Name of the VTK file to write. If the filename ends in '.vtkhdf' then a VTKHDF format file will be written. If the filename ends in - '.vtu' then a binary VTU format file will be written, if the + '.vtu' then a binary VTU format file will be written. If the filename ends in '.vtk' then a legacy VTK file will be written. datasets : dict Dictionary whose keys are the data labels and values are numpy @@ -2818,7 +2818,7 @@ def append_dataset(dset, array): dset[origLen:] = array if self.library != "moab": - raise NotImplemented("VTKHDF output is only supported for MOAB meshes") + raise NotImplementedError("VTKHDF output is only supported for MOAB meshes") # the self.connectivity contains arrays of length 8 to support hex # elements as well, in the case of tetrahedra mesh elements, the diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 84c4a8570bd..2bbd447fa65 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -2,6 +2,7 @@ from tempfile import TemporaryDirectory from pathlib import Path +import h5py import numpy as np import pytest import openmc @@ -484,6 +485,7 @@ def test_umesh(run_in_tmpdir, simple_umesh, export_type): with pytest.raises(ValueError, match='Cannot apply dataset "mean"') as e: simple_umesh.write_data_to_vtk(datasets={'mean': ref_data[:-2]}, filename=filename) + @pytest.mark.skipif(not openmc.lib._dagmc_enabled(), reason="DAGMC not enabled.") def test_write_vtkhdf(request, run_in_tmpdir): """Performs a minimal UnstructuredMesh simulation, reads in the resulting @@ -493,7 +495,7 @@ def test_write_vtkhdf(request, run_in_tmpdir): """ model = openmc.Model() - surf1 = openmc.Sphere(R=1000.0, boundary_type="vacuum") + surf1 = openmc.Sphere(r=1000.0, boundary_type="vacuum") cell1 = openmc.Cell(region=-surf1) model.geometry = openmc.Geometry([cell1]) @@ -546,6 +548,10 @@ def test_write_vtkhdf(request, run_in_tmpdir): assert Path("test_mesh.vtk").exists() assert Path("test_mesh.vtkhdf").exists() + # just ensure we can open the file without error + with h5py.File("test_mesh.vtkhdf", "r"): + ... + def test_mesh_get_homogenized_materials(): """Test the get_homogenized_materials method"""