Skip to content

Commit e5c7d0c

Browse files
jon-proximafusionshimwellrherrero-pfpshriwisepaulromano
authored
Adding vtkhdf option to write vtk data (#3252)
Co-authored-by: shimwell <mail@jshimwell.com> Co-authored-by: Jonathan Shimwell <drshimwell@gmail.com> Co-authored-by: rherrero-pf <156206440+rherrero-pf@users.noreply.github.com> Co-authored-by: Patrick Shriwise <pshriwise@gmail.com> Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
1 parent bd76fc0 commit e5c7d0c

File tree

4 files changed

+377
-18
lines changed

4 files changed

+377
-18
lines changed

‎openmc/mesh.py‎

Lines changed: 148 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
from functools import wraps
66
from math import pi, sqrt, atan2
77
from numbers import Integral, Real
8+
from pathlib import Path
89
from typing import Protocol
910

1011
import h5py
@@ -2443,6 +2444,7 @@ class UnstructuredMesh(MeshBase):
24432444
_UNSUPPORTED_ELEM = -1
24442445
_LINEAR_TET = 0
24452446
_LINEAR_HEX = 1
2447+
_VTK_TETRA = 10
24462448

24472449
def __init__(self, filename: PathLike, library: str, mesh_id: int | None = None,
24482450
name: str = '', length_multiplier: float = 1.0,
@@ -2652,7 +2654,8 @@ def write_vtk_mesh(self, **kwargs):
26522654
warnings.warn(
26532655
"The 'UnstructuredMesh.write_vtk_mesh' method has been renamed "
26542656
"to 'write_data_to_vtk' and will be removed in a future version "
2655-
" of OpenMC.", FutureWarning
2657+
" of OpenMC.",
2658+
FutureWarning,
26562659
)
26572660
self.write_data_to_vtk(**kwargs)
26582661

@@ -2670,26 +2673,54 @@ def write_data_to_vtk(
26702673
Parameters
26712674
----------
26722675
filename : str or pathlib.Path
2673-
Name of the VTK file to write. If the filename ends in '.vtu' then a
2674-
binary VTU format file will be written, if the filename ends in
2675-
'.vtk' then a legacy VTK file will be written.
2676+
Name of the VTK file to write. If the filename ends in '.vtkhdf'
2677+
then a VTKHDF format file will be written. If the filename ends in
2678+
'.vtu' then a binary VTU format file will be written. If the
2679+
filename ends in '.vtk' then a legacy VTK file will be written.
26762680
datasets : dict
26772681
Dictionary whose keys are the data labels and values are numpy
26782682
appropriately sized arrays of the data
26792683
volume_normalization : bool
26802684
Whether or not to normalize the data by the volume of the mesh
26812685
elements
26822686
"""
2687+
2688+
if Path(filename).suffix == ".vtkhdf":
2689+
2690+
self._write_data_to_vtk_hdf5_format(
2691+
filename=filename,
2692+
datasets=datasets,
2693+
volume_normalization=volume_normalization,
2694+
)
2695+
2696+
elif Path(filename).suffix == ".vtk" or Path(filename).suffix == ".vtu":
2697+
2698+
self._write_data_to_vtk_ascii_format(
2699+
filename=filename,
2700+
datasets=datasets,
2701+
volume_normalization=volume_normalization,
2702+
)
2703+
2704+
else:
2705+
raise ValueError(
2706+
"Unsupported file extension, The filename must end with "
2707+
"'.vtkhdf', '.vtu' or '.vtk'"
2708+
)
2709+
2710+
def _write_data_to_vtk_ascii_format(
2711+
self,
2712+
filename: PathLike | None = None,
2713+
datasets: dict | None = None,
2714+
volume_normalization: bool = True,
2715+
):
26832716
from vtkmodules.util import numpy_support
26842717
from vtkmodules import vtkCommonCore
26852718
from vtkmodules import vtkCommonDataModel
26862719
from vtkmodules import vtkIOLegacy
26872720
from vtkmodules import vtkIOXML
26882721

26892722
if self.connectivity is None or self.vertices is None:
2690-
raise RuntimeError(
2691-
"This mesh has not been loaded from a statepoint file."
2692-
)
2723+
raise RuntimeError("This mesh has not been loaded from a statepoint file.")
26932724

26942725
if filename is None:
26952726
filename = f"mesh_{self.id}.vtk"
@@ -2771,29 +2802,128 @@ def write_data_to_vtk(
27712802

27722803
writer.Write()
27732804

2805+
def _write_data_to_vtk_hdf5_format(
2806+
self,
2807+
filename: PathLike | None = None,
2808+
datasets: dict | None = None,
2809+
volume_normalization: bool = True,
2810+
):
2811+
def append_dataset(dset, array):
2812+
"""Convenience function to append data to an HDF5 dataset"""
2813+
origLen = dset.shape[0]
2814+
dset.resize(origLen + array.shape[0], axis=0)
2815+
dset[origLen:] = array
2816+
2817+
if self.library != "moab":
2818+
raise NotImplementedError("VTKHDF output is only supported for MOAB meshes")
2819+
2820+
# the self.connectivity contains arrays of length 8 to support hex
2821+
# elements as well, in the case of tetrahedra mesh elements, the
2822+
# last 4 values are -1 and are removed
2823+
trimmed_connectivity = []
2824+
for cell in self.connectivity:
2825+
# Find the index of the first -1 value, if any
2826+
first_negative_index = np.where(cell == -1)[0]
2827+
if first_negative_index.size > 0:
2828+
# Slice the array up to the first -1 value
2829+
trimmed_connectivity.append(cell[: first_negative_index[0]])
2830+
else:
2831+
# No -1 values, append the whole cell
2832+
trimmed_connectivity.append(cell)
2833+
trimmed_connectivity = np.array(trimmed_connectivity, dtype="int32").flatten()
2834+
2835+
# MOAB meshes supports tet elements only so we know it has 4 points per cell
2836+
points_per_cell = 4
2837+
2838+
# offsets are the indices of the first point of each cell in the array of points
2839+
offsets = np.arange(0, self.n_elements * points_per_cell + 1, points_per_cell)
2840+
2841+
for name, data in datasets.items():
2842+
if data.shape != self.dimension:
2843+
raise ValueError(
2844+
f'Cannot apply dataset "{name}" with '
2845+
f"shape {data.shape} to mesh {self.id} "
2846+
f"with dimensions {self.dimension}"
2847+
)
2848+
2849+
with h5py.File(filename, "w") as f:
2850+
2851+
root = f.create_group("VTKHDF")
2852+
vtk_file_format_version = (2, 1)
2853+
root.attrs["Version"] = vtk_file_format_version
2854+
ascii_type = "UnstructuredGrid".encode("ascii")
2855+
root.attrs.create(
2856+
"Type",
2857+
ascii_type,
2858+
dtype=h5py.string_dtype("ascii", len(ascii_type)),
2859+
)
2860+
2861+
# create hdf5 file structure
2862+
root.create_dataset("NumberOfPoints", (0,), maxshape=(None,), dtype="i8")
2863+
root.create_dataset("Types", (0,), maxshape=(None,), dtype="uint8")
2864+
root.create_dataset("Points", (0, 3), maxshape=(None, 3), dtype="f")
2865+
root.create_dataset(
2866+
"NumberOfConnectivityIds", (0,), maxshape=(None,), dtype="i8"
2867+
)
2868+
root.create_dataset("NumberOfCells", (0,), maxshape=(None,), dtype="i8")
2869+
root.create_dataset("Offsets", (0,), maxshape=(None,), dtype="i8")
2870+
root.create_dataset("Connectivity", (0,), maxshape=(None,), dtype="i8")
2871+
2872+
append_dataset(root["NumberOfPoints"], np.array([len(self.vertices)]))
2873+
append_dataset(root["Points"], self.vertices)
2874+
append_dataset(
2875+
root["NumberOfConnectivityIds"],
2876+
np.array([len(trimmed_connectivity)]),
2877+
)
2878+
append_dataset(root["Connectivity"], trimmed_connectivity)
2879+
append_dataset(root["NumberOfCells"], np.array([self.n_elements]))
2880+
append_dataset(root["Offsets"], offsets)
2881+
2882+
append_dataset(
2883+
root["Types"], np.full(self.n_elements, self._VTK_TETRA, dtype="uint8")
2884+
)
2885+
2886+
cell_data_group = root.create_group("CellData")
2887+
2888+
for name, data in datasets.items():
2889+
2890+
cell_data_group.create_dataset(
2891+
name, (0,), maxshape=(None,), dtype="float64", chunks=True
2892+
)
2893+
2894+
if volume_normalization:
2895+
data /= self.volumes
2896+
append_dataset(cell_data_group[name], data)
2897+
27742898
@classmethod
27752899
def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str):
2776-
filename = group['filename'][()].decode()
2777-
library = group['library'][()].decode()
2778-
if 'options' in group.attrs:
2900+
filename = group["filename"][()].decode()
2901+
library = group["library"][()].decode()
2902+
if "options" in group.attrs:
27792903
options = group.attrs['options'].decode()
27802904
else:
27812905
options = None
27822906

2783-
mesh = cls(filename=filename, library=library, mesh_id=mesh_id, name=name, options=options)
2907+
mesh = cls(
2908+
filename=filename,
2909+
library=library,
2910+
mesh_id=mesh_id,
2911+
name=name,
2912+
options=options,
2913+
)
27842914
mesh._has_statepoint_data = True
2785-
vol_data = group['volumes'][()]
2915+
vol_data = group["volumes"][()]
27862916
mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],))
27872917
mesh.n_elements = mesh.volumes.size
27882918

2789-
vertices = group['vertices'][()]
2919+
vertices = group["vertices"][()]
27902920
mesh._vertices = vertices.reshape((-1, 3))
2791-
connectivity = group['connectivity'][()]
2921+
connectivity = group["connectivity"][()]
27922922
mesh._connectivity = connectivity.reshape((-1, 8))
2793-
mesh._element_types = group['element_types'][()]
2923+
mesh._element_types = group["element_types"][()]
27942924

2795-
if 'length_multiplier' in group:
2796-
mesh.length_multiplier = group['length_multiplier'][()]
2925+
if "length_multiplier" in group:
2926+
mesh.length_multiplier = group["length_multiplier"][()]
27972927

27982928
return mesh
27992929

@@ -2812,7 +2942,7 @@ def to_xml_element(self):
28122942

28132943
element.set("library", self._library)
28142944
if self.options is not None:
2815-
element.set('options', self.options)
2945+
element.set("options", self.options)
28162946
subelement = ET.SubElement(element, "filename")
28172947
subelement.text = str(self.filename)
28182948

‎tests/regression_tests/unstructured_mesh/test_mesh_dagmc_tets.vtk‎

Lines changed: 159 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,159 @@
1+
# vtk DataFile Version 2.0
2+
made_with_cad_to_dagmc_package, Created by Gmsh 4.12.1
3+
ASCII
4+
DATASET UNSTRUCTURED_GRID
5+
POINTS 14 double
6+
-0.5 -0.5 0.5
7+
-0.5 -0.5 -0.5
8+
-0.5 0.5 0.5
9+
-0.5 0.5 -0.5
10+
0.5 -0.5 0.5
11+
0.5 -0.5 -0.5
12+
0.5 0.5 0.5
13+
0.5 0.5 -0.5
14+
-0.5 0 0
15+
0.5 0 0
16+
0 -0.5 0
17+
0 0.5 0
18+
0 0 -0.5
19+
0 0 0.5
20+
21+
CELLS 68 268
22+
1 0
23+
1 1
24+
1 2
25+
1 3
26+
1 4
27+
1 5
28+
1 6
29+
1 7
30+
2 1 0
31+
2 0 2
32+
2 3 2
33+
2 1 3
34+
2 5 4
35+
2 4 6
36+
2 7 6
37+
2 5 7
38+
2 1 5
39+
2 0 4
40+
2 3 7
41+
2 2 6
42+
3 1 0 8
43+
3 0 2 8
44+
3 3 1 8
45+
3 2 3 8
46+
3 5 9 4
47+
3 4 9 6
48+
3 7 9 5
49+
3 6 9 7
50+
3 0 1 10
51+
3 4 0 10
52+
3 1 5 10
53+
3 5 4 10
54+
3 2 11 3
55+
3 6 11 2
56+
3 3 11 7
57+
3 7 11 6
58+
3 1 3 12
59+
3 5 1 12
60+
3 3 7 12
61+
3 7 5 12
62+
3 0 13 2
63+
3 4 13 0
64+
3 2 13 6
65+
3 6 13 4
66+
4 13 8 12 10
67+
4 11 8 12 13
68+
4 12 11 13 9
69+
4 10 12 13 9
70+
4 12 3 11 7
71+
4 13 2 8 0
72+
4 8 12 1 3
73+
4 11 3 8 2
74+
4 10 8 1 0
75+
4 0 10 13 4
76+
4 1 12 10 5
77+
4 11 2 13 6
78+
4 4 9 13 6
79+
4 6 9 11 7
80+
4 10 9 4 5
81+
4 7 9 12 5
82+
4 3 8 12 11
83+
4 1 12 8 10
84+
4 8 11 2 13
85+
4 10 13 8 0
86+
4 13 10 9 4
87+
4 11 13 9 6
88+
4 12 11 9 7
89+
4 9 10 12 5
90+
91+
CELL_TYPES 68
92+
1
93+
1
94+
1
95+
1
96+
1
97+
1
98+
1
99+
1
100+
3
101+
3
102+
3
103+
3
104+
3
105+
3
106+
3
107+
3
108+
3
109+
3
110+
3
111+
3
112+
5
113+
5
114+
5
115+
5
116+
5
117+
5
118+
5
119+
5
120+
5
121+
5
122+
5
123+
5
124+
5
125+
5
126+
5
127+
5
128+
5
129+
5
130+
5
131+
5
132+
5
133+
5
134+
5
135+
5
136+
10
137+
10
138+
10
139+
10
140+
10
141+
10
142+
10
143+
10
144+
10
145+
10
146+
10
147+
10
148+
10
149+
10
150+
10
151+
10
152+
10
153+
10
154+
10
155+
10
156+
10
157+
10
158+
10
159+
10

0 commit comments

Comments
 (0)