diff --git a/openmc/mesh.py b/openmc/mesh.py index e4d83d81d3d..b916757b628 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -9,6 +9,7 @@ import h5py import lxml.etree as ET import numpy as np +from pathlib import Path import openmc import openmc.checkvalue as cv @@ -2358,17 +2359,21 @@ def write_vtk_mesh(self, **kwargs): self.write_data_to_vtk(**kwargs) def write_data_to_vtk( - self, - filename: PathLike | None = None, - datasets: dict | None = None, - volume_normalization: bool = True + self, + filename: PathLike | None = None, + datasets: dict | None = None, + volume_normalization: bool = True, ): """Map data to unstructured VTK mesh elements. + If filename is None, then a filename will be generated based on the mesh ID, + and exported to VTK format. + 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 '.vtu' then a 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 appropriately sized arrays @@ -2377,36 +2382,48 @@ 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 + 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.') + raise RuntimeError( + "This mesh has not been loaded from a statepoint file." + ) if filename is None: - filename = f'mesh_{self.id}.vtk' + filename = f"mesh_{self.id}.vtk" + + if Path(filename).suffix == ".vtk": + writer = vtkIOLegacy.vtkUnstructuredGridWriter() - writer = vtk.vtkUnstructuredGridWriter() + elif Path(filename).suffix == ".vtu": + writer = vtkIOXML.vtkXMLUnstructuredGridWriter() + writer.SetCompressorTypeToZLib() + writer.SetDataModeToBinary() writer.SetFileName(str(filename)) - grid = vtk.vtkUnstructuredGrid() + grid = vtkCommonDataModel.vtkUnstructuredGrid() - vtk_pnts = vtk.vtkPoints() - vtk_pnts.SetData(nps.numpy_to_vtk(self.vertices)) - grid.SetPoints(vtk_pnts) + 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 = vtk.vtkTetra() + elem = vtkCommonDataModel.vtkTetra() elif elem_type == self._LINEAR_HEX: - elem = vtk.vtkHexahedron() + elem = vtkCommonDataModel.vtkHexahedron() elif elem_type == self._UNSUPPORTED_ELEM: n_skipped += 1 + continue else: - raise RuntimeError(f'Invalid element type {elem_type} found') + raise RuntimeError(f"Invalid element type {elem_type} found") + for i, c in enumerate(conn): if c == -1: break @@ -2415,30 +2432,36 @@ def write_data_to_vtk( 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') + 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 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}') + raise ValueError( + f'Cannot apply dataset "{name}" with ' + f"shape {data.shape} to mesh {self.id} " + f"with dimensions {self.dimension}" + ) 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.') + 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 = vtkCommonCore.vtkDoubleArray() arr.SetName(name) arr.SetNumberOfTuples(data.size) 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..2ddd228cf13 --- /dev/null +++ b/tests/regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk @@ -0,0 +1,159 @@ +# vtk DataFile Version 2.0 +made_with_cad_to_dagmc_package, Created by Gmsh 4.12.1 +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 14 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 + +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 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 diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py index 27322165f64..5c6dd517c6d 100644 --- a/tests/unit_tests/test_mesh.py +++ b/tests/unit_tests/test_mesh.py @@ -3,6 +3,9 @@ import numpy as np import pytest import openmc +from pathlib import Path +import tempfile +from vtkmodules import vtkIOLegacy, vtkIOXML @pytest.mark.parametrize("val_left,val_right", [(0, 0), (-1., -1.), (2.0, 2)]) @@ -397,6 +400,85 @@ def test_umesh_roundtrip(run_in_tmpdir, request): assert umesh.id == xml_mesh.id +skip_if_no_dagmc = pytest.mark.skipif( + not openmc.lib._dagmc_enabled(), + reason="DAGMC CAD geometry is not enabled.") + +@skip_if_no_dagmc +@pytest.mark.parametrize('export_type', ('.vtk','.vtu')) +def test_umesh(request, export_type): + """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 + """ + surf1 = openmc.Sphere(r=1.0, boundary_type="vacuum") + material1 = openmc.Material() + material1.add_components(components={"H" : 1.0}) + material1.set_density('g/cm3', 1.0) + + materials = openmc.Materials([material1]) + cell1 = openmc.Cell(region=-surf1, fill=material1) + geometry = openmc.Geometry([cell1]) + + umesh = openmc.UnstructuredMesh( + filename= request.path.parent.parent + / "regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk", + library= "moab", + mesh_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 = ["total"] + + tallies = openmc.Tallies([mesh_tally]) + + source = openmc.IndependentSource() + + settings = openmc.Settings() + settings.run_mode = "fixed source" + settings.batches = 2 + settings.particles = 100 + settings.source = source + + model = openmc.Model( + materials=materials, geometry=geometry, settings=settings, tallies=tallies + ) + + with tempfile.TemporaryDirectory() as tmpdir: + + statepoint_file = model.run(cwd=tmpdir) + + statepoint = openmc.StatePoint(statepoint_file) + tally: openmc.Tally = statepoint.get_tally(name="test_tally") + umesh_from_sp: openmc.UnstructuredMesh = statepoint.meshes[1] + + datasets={ + "mean": tally.mean.squeeze(), + "std_dev": tally.std_dev.squeeze() + } + + filename = Path(tmpdir) / ("test_mesh" + export_type) + umesh_from_sp.write_data_to_vtk(datasets=datasets, filename=filename) + + assert Path(filename).exists() + + if export_type == ".vtk": + reader = vtkIOLegacy.vtkGenericDataObjectReader() + reader.SetFileName(str(filename)) + reader.Update() + assert reader.GetOutput().GetCellData().GetArray("mean") + + elif export_type == ".vtk": + reader = vtkIOXML.vtkGenericDataObjectReader() + reader.SetFileName(str(filename)) + reader.Update() + assert reader.GetOutput().GetCellData().GetArray("mean") + def test_mesh_get_homogenized_materials(): """Test the get_homogenized_materials method"""