8#ifndef GRIDFORMAT_VTK_HDF_UNSTRUCTURED_GRID_WRITER_HPP_
9#define GRIDFORMAT_VTK_HDF_UNSTRUCTURED_GRID_WRITER_HPP_
10#if GRIDFORMAT_HAVE_HIGH_FIVE
15#include <gridformat/common/exceptions.hpp>
16#include <gridformat/common/md_layout.hpp>
17#include <gridformat/common/concepts.hpp>
18#include <gridformat/common/lvalue_reference.hpp>
20#include <gridformat/parallel/communication.hpp>
21#include <gridformat/parallel/concepts.hpp>
31template<
bool is_transient,
32 Concepts::UnstructuredGrid G,
33 Concepts::Communicator Communicator = NullCommunicator>
35 static constexpr int root_rank = 0;
36 static constexpr std::size_t vtk_space_dim = 3;
38 using CT = CoordinateType<G>;
39 template<
typename T>
using Vector = std::array<T, vtk_space_dim>;
40 template<
typename T>
using Tensor = std::array<T, vtk_space_dim*vtk_space_dim>;
43 using HDF5File = HDF5::File<Communicator>;
47 .append_null_terminator_to_strings =
true
50 struct TimeSeriesOffsets {
51 std::size_t cell_offset;
52 std::size_t connectivity_offset;
53 std::size_t point_offset;
60 requires(!is_transient && std::is_same_v<Communicator, NullCommunicator>)
65 requires(!is_transient && std::is_copy_constructible_v<Communicator>)
71 std::string filename_without_extension,
73 requires(is_transient && std::is_same_v<Communicator, NullCommunicator>)
76 , _timeseries_filename{std::move(filename_without_extension) +
".hdf"}
77 , _transient_opts{std::move(opts)}
81 const Communicator& comm,
82 std::string filename_without_extension,
84 requires(is_transient && std::is_copy_constructible_v<Communicator>)
87 , _timeseries_filename{std::move(filename_without_extension) +
".hdf"}
88 , _transient_opts{std::move(opts)}
91 const Communicator& communicator()
const {
96 void _write(std::ostream&)
const {
97 throw InvalidState(
"VTKHDFUnstructuredGridWriter does not support export into stream");
100 void _write(
const std::string& filename_with_ext)
const {
101 if constexpr (is_transient)
102 throw InvalidState(
"This overload only works for non-transient output");
103 HDF5File file{filename_with_ext, _comm, HDF5File::overwrite};
107 std::string _write(
double t) {
108 if constexpr (!is_transient)
109 throw InvalidState(
"This overload only works for transient output");
111 if (this->_step_count == 0)
112 HDF5File::clear(_timeseries_filename, _comm);
114 HDF5File file{_timeseries_filename, _comm, HDF5File::append};
115 const auto offsets = _write_to(file);
117 file.write_attribute(this->_step_count+1,
"/VTKHDF/Steps/NSteps");
118 file.write(std::array{t},
"VTKHDF/Steps/Values");
119 file.write(std::vector{offsets.point_offset},
"VTKHDF/Steps/PointOffsets");
120 file.write(std::vector{std::array{offsets.cell_offset}},
"/VTKHDF/Steps/CellOffsets");
121 file.write(std::vector{std::array{offsets.connectivity_offset}},
"VTKHDF/Steps/ConnectivityIdOffsets");
123 file.write(std::vector{Parallel::size(_comm)},
"/VTKHDF/Steps/NumberOfParts");
124 if (this->_step_count > 0 && _transient_opts.
static_grid) {
126 std::vector{_get_last_step_data(file,
"PartOffsets")},
127 "/VTKHDF/Steps/PartOffsets"
130 const std::size_t offset = this->_step_count == 0
132 : _get_last_step_data(file,
"PartOffsets") + Parallel::size(_comm);
133 file.write(std::vector{offset},
"/VTKHDF/Steps/PartOffsets");
136 return _timeseries_filename;
139 TimeSeriesOffsets _write_to(HDF5File& file)
const {
140 file.write_attribute(std::array<std::size_t, 2>{(is_transient ? 2 : 1), 0},
"/VTKHDF/Version");
141 file.write_attribute(
"UnstructuredGrid",
"/VTKHDF/Type");
143 TimeSeriesOffsets offsets;
144 const auto context = IOContext::from(this->grid(), _comm, root_rank);
145 _write_num_cells_and_points(file, context);
146 offsets.point_offset = _write_coordinates(file, context);
147 offsets.connectivity_offset = _write_connectivity(file, context);
148 offsets.cell_offset = _write_types(file, context);
149 _write_offsets(file, context);
150 _write_meta_data(file);
151 _write_point_fields(file, context);
152 _write_cell_fields(file, context);
157 void _write_num_cells_and_points(HDF5File& file,
const IOContext& context)
const {
158 _write_values(file,
"/VTKHDF/NumberOfPoints", std::vector{number_of_points(this->grid())}, context);
159 _write_values(file,
"/VTKHDF/NumberOfCells", std::vector{number_of_cells(this->grid())}, context);
162 std::size_t _write_coordinates(HDF5File& file,
const IOContext& context)
const {
163 if constexpr (is_transient) {
164 if (this->_step_count > 0 && _transient_opts.
static_grid)
165 return _get_last_step_data(file,
"PointOffsets");
167 const auto coords_field = VTK::make_coordinates_field<CT>(this->grid(),
false);
168 const auto offset = _get_current_offset(file,
"/VTKHDF/Points");
169 _write_point_field(file,
"/VTKHDF/Points", *coords_field, context);
173 std::size_t _write_connectivity(HDF5File& file,
const IOContext& context)
const {
174 if constexpr (is_transient) {
175 if (this->_step_count > 0 && _transient_opts.
static_grid)
176 return _get_last_step_data(file,
"ConnectivityIdOffsets");
178 const auto point_id_map = make_point_id_map(this->grid());
179 const auto connectivity_field = VTK::make_connectivity_field(this->grid(), point_id_map);
180 const auto num_entries = connectivity_field->layout().number_of_entries();
181 const auto my_num_ids = std::vector{
static_cast<long>(num_entries)};
182 std::vector<long> connectivity(num_entries);
183 connectivity_field->export_to(connectivity);
184 const auto offset = _get_current_offset(file,
"/VTKHDF/Connectivity");
185 _write_values(file,
"/VTKHDF/Connectivity", connectivity, context);
186 _write_values(file,
"/VTKHDF/NumberOfConnectivityIds", my_num_ids, context);
190 std::size_t _write_types(HDF5File& file,
const IOContext& context)
const {
191 if constexpr (is_transient) {
192 if (this->_step_count > 0 && _transient_opts.
static_grid)
193 return _get_last_step_data(file,
"CellOffsets");
195 const auto types_field = VTK::make_cell_types_field(this->grid());
196 std::vector<std::uint8_t> types(types_field->layout().number_of_entries());
197 types_field->export_to(types);
198 const auto offset = _get_current_offset(file,
"VTKHDF/Types");
199 _write_values(file,
"VTKHDF/Types", types, context);
203 std::size_t _write_offsets(HDF5File& file,
const IOContext& context)
const {
204 if constexpr (is_transient) {
205 if (this->_step_count > 0 && _transient_opts.
static_grid)
206 return _get_last_step_data(file,
"CellOffsets");
208 const auto offsets_field = VTK::make_offsets_field(this->grid());
209 const auto num_offset_entries = offsets_field->layout().number_of_entries() + 1;
210 std::vector<long> offsets(num_offset_entries);
211 offsets_field->export_to(std::ranges::subrange(std::next(offsets.begin()), offsets.end()));
212 offsets[0] =
long{0};
213 const auto offset = _get_current_offset(file,
"/VTKHDF/Offsets");
214 _write_values(file,
"/VTKHDF/Offsets", offsets, context);
218 void _write_meta_data(HDF5File& file)
const {
219 std::ranges::for_each(this->_meta_data_field_names(), [&] (
const std::string& name) {
220 if constexpr (is_transient) {
222 file.write(std::array{0},
"/VTKHDF/Steps/FieldDataOffsets/" + name);
225 file.write(std::array{this->_step_count},
"/VTKHDF/Steps/FieldDataOffsets/" + name);
228 TransformedField sub{this->_get_meta_data_field_ptr(name), FieldTransformation::as_sub_field};
229 file.write(sub,
"/VTKHDF/FieldData/" + name);
231 file.write(*this->_get_meta_data_field_ptr(name),
"/VTKHDF/FieldData/" + name);
236 void _write_point_fields(HDF5File& file,
const IOContext& context)
const {
237 std::ranges::for_each(this->_point_field_names(), [&] (
const std::string& name) {
238 if constexpr (is_transient)
241 _get_current_offset(file,
"/VTKHDF/PointData/" + name),
242 "/VTKHDF/Steps/PointDataOffsets/" + name
244 auto reshaped = _reshape(VTK::make_vtk_field(this->_get_point_field_ptr(name)));
245 _write_point_field(file,
"/VTKHDF/PointData/" + name, *reshaped, context);
249 void _write_cell_fields(HDF5File& file,
const IOContext& context)
const {
250 std::ranges::for_each(this->_cell_field_names(), [&] (
const std::string& name) {
251 if constexpr (is_transient)
254 _get_current_offset(file,
"VTKHDF/CellData/" + name),
255 "/VTKHDF/Steps/CellDataOffsets/" + name
257 auto reshaped = _reshape(VTK::make_vtk_field(this->_get_cell_field_ptr(name)));
258 _write_cell_field(file,
"/VTKHDF/CellData/" + name, *reshaped, context);
264 const auto layout = field_ptr->layout();
265 if (layout.dimension() > 2)
268 MDLayout{{layout.extent(0), layout.number_of_entries(1)}}
273 template<Concepts::Scalar T>
274 void _write_values(HDF5File& file,
275 const std::string& path,
276 const std::vector<T>& values,
278 if (context.is_parallel) {
279 const auto num_values = values.size();
280 const auto total_num_values = Parallel::sum(_comm, num_values, root_rank);
281 const auto my_total_num_values = Parallel::broadcast(_comm, total_num_values, root_rank);
282 const auto my_offset = _accumulate_rank_offset(num_values);
284 .offset = std::vector{my_offset},
285 .count = std::vector{num_values},
286 .total_size = std::vector{my_total_num_values}
288 file.write(values, path, slice);
290 file.write(values, path);
294 void _write_point_field(HDF5File& file,
295 const std::string& path,
298 _write_field(file, path, field, context.is_parallel, context.my_point_offset, context.num_points_total);
301 void _write_cell_field(HDF5File& file,
302 const std::string& path,
305 _write_field(file, path, field, context.is_parallel, context.my_cell_offset, context.num_cells_total);
308 void _write_field(HDF5File& file,
309 const std::string& path,
312 std::size_t main_offset,
313 std::size_t main_size)
const {
315 const auto layout = field.
layout();
316 std::vector<std::size_t> count(layout.dimension());
317 layout.export_to(count);
319 std::vector<std::size_t> size = count;
320 std::vector<std::size_t> offset(layout.dimension(), 0);
321 offset.at(0) = main_offset;
322 size.at(0) = main_size;
324 file.write(field, path, HDF5::Slice{
325 .offset = std::move(offset),
326 .count = std::move(count),
327 .total_size = std::move(size)
330 file.write(field, path);
334 std::size_t _accumulate_rank_offset(std::size_t my_size)
const {
335 auto all_offsets = Parallel::gather(_comm, my_size, root_rank);
336 if (Parallel::rank(_comm) == root_rank) {
337 std::partial_sum(all_offsets.begin(), std::prev(all_offsets.end()), all_offsets.begin());
338 std::copy(std::next(all_offsets.rbegin()), all_offsets.rend(), all_offsets.rbegin());
341 const auto my_offset = Parallel::scatter(_comm, all_offsets, root_rank);
342 if (my_offset.size() != 1)
343 throw ValueError(
"Unexpected scatter result");
347 void _write_step_offset(HDF5File& file,
348 const std::size_t offset,
349 const std::string& path)
const {
350 file.write(std::array{offset}, path);
353 std::size_t _get_current_offset(HDF5File& file,
const std::string& path)
const {
354 const auto cur_dimensions = file.get_dimensions(path);
355 return cur_dimensions ? (*cur_dimensions)[0] : 0;
358 std::size_t _get_last_step_data(
const HDF5File& file,
const std::string& sub_path)
const {
359 if (this->_step_count == 0)
360 throw ValueError(
"Last step data can only be read after at least one write");
361 const std::string path =
"VTKHDF/Steps/" + sub_path;
363 auto step_dimensions = file.get_dimensions(path).value();
364 std::ranges::fill(step_dimensions, std::size_t{1});
366 auto access_offset = step_dimensions;
367 std::ranges::fill(access_offset, std::size_t{0});
368 access_offset.at(0) = this->_step_count - 1;
370 return file.template read_dataset_to<std::size_t>(path, HDF5::Slice{
371 .offset = std::move(access_offset),
372 .count = std::move(step_dimensions)
377 std::string _timeseries_filename =
"";
385template<Concepts::UnstructuredGr
id G, Concepts::Communicator C = NullCommunicator>
389 using ParentType::ParentType;
396template<Concepts::UnstructuredGr
id G, Concepts::Communicator C = NullCommunicator>
400 using ParentType::ParentType;
403template<Concepts::UnstructuredGr
id G>
405template<Concepts::UnstructuredGr
id G, Concepts::Communicator C>
408template<Concepts::UnstructuredGr
id G>
410template<Concepts::UnstructuredGr
id G, Concepts::Communicator C>
411VTKHDFUnstructuredTimeSeriesWriter(
const G&,
const C&, std::string, VTK::HDFTransientOptions = {}) -> VTKHDFUnstructuredTimeSeriesWriter<G, C>;
Base classes for grid data writers.
std::shared_ptr< const Field > FieldPtr
Pointer type used by writers/readers for fields.
Definition: field.hpp:186
FieldPtr make_field_ptr(F &&f)
Factory function for field pointers.
Definition: field.hpp:192
Common functionality for writing VTK HDF files.
Helper class to store processor offsets in parallel writes.
Definition: hdf_common.hpp:37
Common functionality for VTK writers.