8#ifndef GRIDFORMAT_VTK_HDF_IMAGE_GRID_WRITER_HPP_
9#define GRIDFORMAT_VTK_HDF_IMAGE_GRID_WRITER_HPP_
10#if GRIDFORMAT_HAVE_HIGH_FIVE
18#include <gridformat/common/exceptions.hpp>
19#include <gridformat/common/md_layout.hpp>
20#include <gridformat/common/concepts.hpp>
21#include <gridformat/common/matrix.hpp>
22#include <gridformat/common/ranges.hpp>
23#include <gridformat/common/type_traits.hpp>
24#include <gridformat/common/lvalue_reference.hpp>
25#include <gridformat/common/field_transformations.hpp>
27#include <gridformat/parallel/communication.hpp>
28#include <gridformat/parallel/concepts.hpp>
39template<
bool is_transient, Concepts::ImageGr
id Gr
id, Concepts::Communicator Communicator = NullCommunicator>
41 static constexpr int root_rank = 0;
42 static constexpr std::size_t dim = dimension<Grid>;
43 static constexpr std::size_t vtk_space_dim = 3;
44 static constexpr std::array<std::size_t, 2> version{1, 0};
46 using CT = CoordinateType<Grid>;
48 using HDF5File = HDF5::File<Communicator>;
51 const std::array<CT, dim> origin;
52 const std::array<CT, dim> spacing;
53 const std::array<std::size_t, dim> extents;
55 template<
typename O,
typename S,
typename E>
56 ImageSpecs(O&& origin, S&& spacing, E&& extents)
57 : origin{Ranges::to_array<dim, CT>(std::forward<O>(origin))}
58 , spacing{Ranges::to_array<dim, CT>(std::forward<S>(spacing))}
59 , extents{Ranges::to_array<dim, std::size_t>(std::forward<E>(extents))}
65 .append_null_terminator_to_strings =
true
70 requires(std::is_same_v<Communicator, NullCommunicator>)
75 requires(std::is_copy_constructible_v<Communicator>)
81 std::string filename_without_extension,
84 .static_meta_data =
false
86 requires(is_transient && std::is_same_v<Communicator, NullCommunicator>)
91 const Communicator& comm,
92 std::string filename_without_extension,
95 .static_meta_data =
false
97 requires(is_transient && std::is_copy_constructible_v<Communicator>)
100 , _timeseries_filename{std::move(filename_without_extension) +
".hdf"}
101 , _transient_opts{std::move(opts)} {
103 throw ValueError(
"Transient VTK-HDF ImageData files do not support evolving grids");
106 const Communicator& communicator()
const {
111 void _write(std::ostream&)
const {
112 throw InvalidState(
"VTKHDFImageGridWriter does not support export into stream");
115 std::string _write([[maybe_unused]]
double t) {
116 if constexpr (!is_transient)
117 throw InvalidState(
"This overload only works for transient output");
119 if (this->_step_count == 0)
120 HDF5File::clear(_timeseries_filename, _comm);
122 HDF5File file{_timeseries_filename, _comm, HDF5File::Mode::append};
124 file.write_attribute(this->_step_count+1,
"/VTKHDF/Steps/NSteps");
125 file.write(std::array{t},
"/VTKHDF/Steps/Values");
126 return _timeseries_filename;
129 void _write(
const std::string& filename_with_ext)
const {
130 if constexpr (is_transient)
131 throw InvalidState(
"This overload only works for non-transient output");
132 HDF5File file{filename_with_ext, _comm, HDF5File::Mode::overwrite};
136 void _write_to(HDF5File& file)
const {
137 const ImageSpecs my_specs{origin(this->grid()), spacing(this->grid()), extents(this->grid())};
138 const auto [overall_specs, my_offset] = _get_image_specs(my_specs);
139 const auto cell_slice_base = _make_slice(
140 overall_specs.extents,
148 const auto my_point_extents = Ranges::to_array<dim>(
149 std::views::iota(std::size_t{0}, dim)
150 | std::views::transform([&] (std::size_t dir) {
151 const auto overall_end = overall_specs.extents[dir];
152 const auto my_end = my_specs.extents[dir] + my_offset[dir];
153 return my_end < overall_end ? my_specs.extents[dir] : my_specs.extents[dir] + 1;
155 const auto point_slice_base = _make_slice<true>(
156 overall_specs.extents,
161 file.write_attribute(std::array<std::size_t, 2>{(is_transient ? 2 : 1), 0},
"/VTKHDF/Version");
162 file.write_attribute(Ranges::to_array<vtk_space_dim>(overall_specs.origin),
"/VTKHDF/Origin");
163 file.write_attribute(Ranges::to_array<vtk_space_dim>(overall_specs.spacing),
"/VTKHDF/Spacing");
164 file.write_attribute(VTK::CommonDetail::get_extents(overall_specs.extents),
"/VTKHDF/WholeExtent");
165 file.write_attribute(_get_direction(),
"/VTKHDF/Direction");
166 file.write_attribute(
"ImageData",
"/VTKHDF/Type");
168 std::ranges::for_each(this->_meta_data_field_names(), [&] (
const std::string& name) {
169 if constexpr (is_transient) {
171 file.write(std::array{0},
"/VTKHDF/Steps/FieldDataOffsets/" + name);
174 file.write(std::array{this->_step_count},
"/VTKHDF/Steps/FieldDataOffsets/" + name);
176 auto field_ptr = this->_get_meta_data_field_ptr(name);
177 auto sub =
make_field_ptr(TransformedField{field_ptr, FieldTransformation::as_sub_field});
178 _write_field(file, sub,
"/VTKHDF/FieldData/" + name, _slice_from(sub));
180 auto field_ptr = this->_get_meta_data_field_ptr(name);
181 _write_field(file, field_ptr,
"/VTKHDF/FieldData/" + name, _slice_from(field_ptr));
186 std::vector<std::size_t> non_zero_extents;
188 my_specs.extents | std::views::filter([] (
auto e) { return e != 0; }),
189 std::back_inserter(non_zero_extents)
192 std::ranges::for_each(this->_point_field_names(), [&] (
const std::string& name) {
193 auto field_ptr = _reshape(
194 VTK::make_vtk_field(this->_get_point_field_ptr(name)),
195 Ranges::incremented(non_zero_extents, 1) | std::views::reverse,
196 point_slice_base.count
198 _write_field(file, field_ptr,
"/VTKHDF/PointData/" + name, point_slice_base);
201 std::ranges::for_each(this->_cell_field_names(), [&] (
const std::string& name) {
202 auto field_ptr = _reshape(
203 VTK::make_vtk_field(this->_get_cell_field_ptr(name)),
204 non_zero_extents | std::views::reverse,
205 cell_slice_base.count
207 _write_field(file, field_ptr,
"/VTKHDF/CellData/" + name, cell_slice_base);
211 template<std::ranges::range E, std::ranges::range S>
212 FieldPtr _reshape(
FieldPtr f,
const E& row_major_extents, S&& _slice_end)
const {
213 auto flat = _flatten(f);
214 auto structured = _make_structured(flat, row_major_extents);
215 const auto structured_layout = structured->layout();
217 std::vector<std::size_t> slice_end;
218 std::ranges::copy(_slice_end, std::back_inserter(slice_end));
219 for (std::size_t i = slice_end.size(); i < structured_layout.dimension(); ++i)
220 slice_end.push_back(structured_layout.extent(i));
221 return transform(structured, FieldTransformation::take_slice({
222 .from = std::vector<std::size_t>(slice_end.size(), 0),
228 const auto layout = f->layout();
229 if (layout.dimension() <= 2)
232 auto nl = MDLayout{{layout.extent(0), layout.number_of_entries(1)}};
233 return transform(f, FieldTransformation::reshape_to(std::move(nl)));
236 template<std::ranges::range E>
238 const auto layout = f->layout();
239 std::vector<std::size_t> target_layout(Ranges::size(row_major_extents));
240 std::ranges::copy(row_major_extents, target_layout.begin());
241 if (layout.dimension() > 1)
242 target_layout.push_back(layout.extent(1));
243 return transform(f, FieldTransformation::reshape_to(MDLayout{std::move(target_layout)}));
246 auto _get_image_specs(
const ImageSpecs& piece_specs)
const {
247 using OffsetType = std::ranges::range_value_t<
decltype(piece_specs.extents)>;
248 if (Parallel::size(_comm) > 1) {
250 const auto all_origins = Parallel::gather(_comm, piece_specs.origin, root_rank);
251 const auto all_extents = Parallel::gather(_comm, piece_specs.extents, root_rank);
252 const auto is_negative_axis = VTK::CommonDetail::structured_grid_axis_orientation(piece_specs.spacing);
253 const auto [exts_begin, exts_end, whole_extent, origin] = helper.compute_extents_and_origin(
260 const auto my_whole_extent = Parallel::broadcast(_comm, whole_extent, root_rank);
261 const auto my_whole_origin = Parallel::broadcast(_comm, origin, root_rank);
262 const auto my_extent_offset = Parallel::scatter(_comm, Ranges::flat(exts_begin), root_rank);
263 return std::make_tuple(
264 ImageSpecs{my_whole_origin, piece_specs.spacing, my_whole_extent},
268 return std::make_tuple(piece_specs, std::vector<OffsetType>(dim, 0));
272 template<
bool increment = false,
typename TotalExtents,
typename Extents,
typename Offsets>
273 HDF5::Slice _make_slice(
const TotalExtents& total_extents,
274 const Extents& extents,
275 const Offsets& offsets)
const {
277 result.total_size.emplace();
280 std::vector<bool> is_nonzero;
282 total_extents | std::views::transform([] (
const auto& v) {
return v != 0; }),
283 std::back_inserter(is_nonzero)
285 std::ranges::for_each(total_extents, [&, i=0] (
const auto& value)
mutable {
287 result.total_size.value().push_back(value + (increment ? 1 : 0));
289 std::ranges::for_each(extents, [&, i=0] (
const auto& value)
mutable {
291 result.count.push_back(value);
293 std::ranges::for_each(offsets, [&, i=0] (
const auto& value)
mutable {
295 result.offset.push_back(value);
299 std::ranges::reverse(result.total_size.value());
300 std::ranges::reverse(result.count);
301 std::ranges::reverse(result.offset);
306 auto _get_direction()
const {
307 using T = MDRangeScalar<
decltype(basis(this->grid()))>;
308 std::array<T, vtk_space_dim*vtk_space_dim> coefficients;
309 std::ranges::fill(coefficients, T{0});
310 std::ranges::for_each(
311 Matrix{basis(this->grid())}.transposed(),
312 [it = coefficients.begin()] (
const std::ranges::range
auto& row)
mutable {
313 std::ranges::copy(row, it);
314 std::advance(it, vtk_space_dim);
320 void _write_field(HDF5File& file,
322 const std::string& path,
323 const HDF5::Slice& slice)
const {
324 const std::size_t dimension_offset = is_transient ? 1 : 0;
325 std::vector<std::size_t> size(slice.total_size.value().size() + dimension_offset);
326 std::vector<std::size_t> count(slice.count.size() + dimension_offset);
327 std::vector<std::size_t> offset(slice.offset.size() + dimension_offset);
328 std::ranges::copy(slice.total_size.value(), std::ranges::begin(size | std::views::drop(dimension_offset)));
329 std::ranges::copy(slice.count, std::ranges::begin(count | std::views::drop(dimension_offset)));
330 std::ranges::copy(slice.offset, std::ranges::begin(offset | std::views::drop(dimension_offset)));
332 const auto layout = field->layout();
333 const bool is_vector_field = layout.dimension() > size.size() - dimension_offset;
335 std::ranges::for_each(
336 std::views::iota(size.size() - dimension_offset, layout.dimension()),
337 [&] (
const std::size_t codim) {
338 size.push_back(layout.extent(codim));
339 count.push_back(layout.extent(codim));
344 if constexpr (is_transient) {
348 FieldPtr sub_field = transform(field, FieldTransformation::as_sub_field);
349 file.write(*sub_field, path, HDF5::Slice{
350 .offset = std::move(offset),
351 .count = std::move(count),
352 .total_size = std::move(size)
355 file.write(*field, path, HDF5::Slice{
356 .offset = std::move(offset),
357 .count = std::move(count),
358 .total_size = std::move(size)
363 HDF5::Slice _slice_from(
FieldPtr field)
const {
364 const auto layout = field->layout();
365 std::vector<std::size_t> dims(layout.dimension());
366 std::vector<std::size_t> offset(layout.dimension(), 0);
367 layout.export_to(dims);
369 .offset = std::move(offset),
376 std::string _timeseries_filename =
"";
384template<Concepts::ImageGr
id G, Concepts::Communicator C = NullCommunicator>
388 using ParentType::ParentType;
395template<Concepts::ImageGr
id G, Concepts::Communicator C = NullCommunicator>
399 using ParentType::ParentType;
402template<Concepts::ImageGr
id G>
404template<Concepts::ImageGr
id G, Concepts::Communicator C>
407template<Concepts::ImageGr
id G>
409template<Concepts::ImageGr
id G, Concepts::Communicator C>
410VTKHDFImageGridTimeSeriesWriter(
const G&,
const C&, std::string, VTK::HDFTransientOptions = {}) -> VTKHDFImageGridTimeSeriesWriter<G, C>;
415template<
typename... Args>
418template<
typename... Args>
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 function for writing parallel VTK files.
Helper class to store processor offsets in parallel writes.
Definition: hdf_common.hpp:37
Common functionality for VTK writers.