GridFormat 0.2.1
I/O-Library for grid-like data structures
Loading...
Searching...
No Matches
hdf_image_grid_writer.hpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2022-2023 Dennis Gläser <dennis.glaeser@iws.uni-stuttgart.de>
2// SPDX-License-Identifier: MIT
8#ifndef GRIDFORMAT_VTK_HDF_IMAGE_GRID_WRITER_HPP_
9#define GRIDFORMAT_VTK_HDF_IMAGE_GRID_WRITER_HPP_
10#if GRIDFORMAT_HAVE_HIGH_FIVE
11
12#include <ranges>
13#include <iterator>
14#include <algorithm>
15#include <utility>
16#include <tuple>
17
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>
26
27#include <gridformat/parallel/communication.hpp>
28#include <gridformat/parallel/concepts.hpp>
29
32
36
37namespace GridFormat {
38
39template<bool is_transient, Concepts::ImageGrid Grid, Concepts::Communicator Communicator = NullCommunicator>
40class VTKHDFImageGridWriterImpl : public GridDetail::WriterBase<is_transient, Grid>::type {
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};
45
46 using CT = CoordinateType<Grid>;
48 using HDF5File = HDF5::File<Communicator>;
49
50 struct ImageSpecs {
51 const std::array<CT, dim> origin;
52 const std::array<CT, dim> spacing;
53 const std::array<std::size_t, dim> extents;
54
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))}
60 {}
61 };
62
63 static constexpr WriterOptions writer_opts{
65 .append_null_terminator_to_strings = true
66 };
67
68 public:
69 explicit VTKHDFImageGridWriterImpl(LValueReferenceOf<const Grid> grid)
70 requires(std::is_same_v<Communicator, NullCommunicator>)
71 : GridWriter<Grid>(grid.get(), ".hdf", writer_opts)
72 {}
73
74 explicit VTKHDFImageGridWriterImpl(LValueReferenceOf<const Grid> grid, const Communicator& comm)
75 requires(std::is_copy_constructible_v<Communicator>)
76 : GridWriter<Grid>(grid.get(), ".hdf", writer_opts)
77 , _comm{comm}
78 {}
79
80 explicit VTKHDFImageGridWriterImpl(LValueReferenceOf<const Grid> grid,
81 std::string filename_without_extension,
83 .static_grid = true,
84 .static_meta_data = false
85 })
86 requires(is_transient && std::is_same_v<Communicator, NullCommunicator>)
87 : VTKHDFImageGridWriterImpl(grid.get(), NullCommunicator{}, filename_without_extension, std::move(opts))
88 {}
89
90 explicit VTKHDFImageGridWriterImpl(LValueReferenceOf<const Grid> grid,
91 const Communicator& comm,
92 std::string filename_without_extension,
94 .static_grid = true,
95 .static_meta_data = false
96 })
97 requires(is_transient && std::is_copy_constructible_v<Communicator>)
98 : TimeSeriesGridWriter<Grid>(grid.get(), writer_opts)
99 , _comm{comm}
100 , _timeseries_filename{std::move(filename_without_extension) + ".hdf"}
101 , _transient_opts{std::move(opts)} {
102 if (!_transient_opts.static_grid)
103 throw ValueError("Transient VTK-HDF ImageData files do not support evolving grids");
104 }
105
106 const Communicator& communicator() const {
107 return _comm;
108 }
109
110 private:
111 void _write(std::ostream&) const {
112 throw InvalidState("VTKHDFImageGridWriter does not support export into stream");
113 }
114
115 std::string _write([[maybe_unused]] double t) {
116 if constexpr (!is_transient)
117 throw InvalidState("This overload only works for transient output");
118
119 if (this->_step_count == 0)
120 HDF5File::clear(_timeseries_filename, _comm);
121
122 HDF5File file{_timeseries_filename, _comm, HDF5File::Mode::append};
123 _write_to(file);
124 file.write_attribute(this->_step_count+1, "/VTKHDF/Steps/NSteps");
125 file.write(std::array{t}, "/VTKHDF/Steps/Values");
126 return _timeseries_filename;
127 }
128
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};
133 _write_to(file);
134 }
135
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,
141 my_specs.extents,
142 my_offset
143 );
144
145 // in order to avoid overlapping hyperslabs for point data in parallel I/O,
146 // we only write the last entries of the slab (per direction) when our portion
147 // of the image is the last portion of the overall image in that direction.
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;
154 }));
155 const auto point_slice_base = _make_slice<true>(
156 overall_specs.extents,
157 my_point_extents,
158 my_offset
159 );
160
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");
167
168 std::ranges::for_each(this->_meta_data_field_names(), [&] (const std::string& name) {
169 if constexpr (is_transient) {
170 if (this->_step_count > 0 && _transient_opts.static_meta_data) {
171 file.write(std::array{0}, "/VTKHDF/Steps/FieldDataOffsets/" + name);
172 return;
173 } else {
174 file.write(std::array{this->_step_count}, "/VTKHDF/Steps/FieldDataOffsets/" + name);
175 }
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));
179 } else {
180 auto field_ptr = this->_get_meta_data_field_ptr(name);
181 _write_field(file, field_ptr, "/VTKHDF/FieldData/" + name, _slice_from(field_ptr));
182 }
183
184 });
185
186 std::vector<std::size_t> non_zero_extents;
187 std::ranges::copy(
188 my_specs.extents | std::views::filter([] (auto e) { return e != 0; }),
189 std::back_inserter(non_zero_extents)
190 );
191
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
197 );
198 _write_field(file, field_ptr, "/VTKHDF/PointData/" + name, point_slice_base);
199 });
200
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
206 );
207 _write_field(file, field_ptr, "/VTKHDF/CellData/" + name, cell_slice_base);
208 });
209 }
210
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();
216
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),
223 .to = slice_end
224 }));
225 }
226
227 FieldPtr _flatten(FieldPtr f) const {
228 const auto layout = f->layout();
229 if (layout.dimension() <= 2)
230 return f;
231 // vtk requires tensors to be made flat
232 auto nl = MDLayout{{layout.extent(0), layout.number_of_entries(1)}};
233 return transform(f, FieldTransformation::reshape_to(std::move(nl)));
234 }
235
236 template<std::ranges::range E>
237 FieldPtr _make_structured(FieldPtr f, const E& row_major_extents) const {
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)}));
244 }
245
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(
254 all_origins,
255 all_extents,
256 is_negative_axis,
257 basis(this->grid())
258 );
259
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},
265 my_extent_offset
266 );
267 } else {
268 return std::make_tuple(piece_specs, std::vector<OffsetType>(dim, 0));
269 }
270 }
271
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 {
276 HDF5::Slice result;
277 result.total_size.emplace();
278
279 // use only those dimensions that are not zero
280 std::vector<bool> is_nonzero;
281 std::ranges::copy(
282 total_extents | std::views::transform([] (const auto& v) { return v != 0; }),
283 std::back_inserter(is_nonzero)
284 );
285 std::ranges::for_each(total_extents, [&, i=0] (const auto& value) mutable {
286 if (is_nonzero[i++])
287 result.total_size.value().push_back(value + (increment ? 1 : 0));
288 });
289 std::ranges::for_each(extents, [&, i=0] (const auto& value) mutable {
290 if (is_nonzero[i++])
291 result.count.push_back(value);
292 });
293 std::ranges::for_each(offsets, [&, i=0] (const auto& value) mutable {
294 if (is_nonzero[i++])
295 result.offset.push_back(value);
296 });
297
298 // slices in VTK are accessed with the last coordinate first (i.e. values[z][y][x])
299 std::ranges::reverse(result.total_size.value());
300 std::ranges::reverse(result.count);
301 std::ranges::reverse(result.offset);
302
303 return result;
304 }
305
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);
315 }
316 );
317 return coefficients;
318 }
319
320 void _write_field(HDF5File& file,
321 FieldPtr field,
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)));
331
332 const auto layout = field->layout();
333 const bool is_vector_field = layout.dimension() > size.size() - dimension_offset;
334 if (is_vector_field)
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));
340 offset.push_back(0);
341 }
342 );
343
344 if constexpr (is_transient) {
345 size.at(0) = 1;
346 count.at(0) = 1;
347 offset.at(0) = 0;
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)
353 });
354 } else {
355 file.write(*field, path, HDF5::Slice{
356 .offset = std::move(offset),
357 .count = std::move(count),
358 .total_size = std::move(size)
359 });
360 }
361 }
362
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);
368 return {
369 .offset = std::move(offset),
370 .count = dims,
371 .total_size = dims
372 };
373 }
374
375 Communicator _comm;
376 std::string _timeseries_filename = "";
377 VTK::HDFTransientOptions _transient_opts;
378};
379
384template<Concepts::ImageGrid G, Concepts::Communicator C = NullCommunicator>
387 public:
388 using ParentType::ParentType;
389};
390
395template<Concepts::ImageGrid G, Concepts::Communicator C = NullCommunicator>
398 public:
399 using ParentType::ParentType;
400};
401
402template<Concepts::ImageGrid G>
404template<Concepts::ImageGrid G, Concepts::Communicator C>
406
407template<Concepts::ImageGrid G>
408VTKHDFImageGridTimeSeriesWriter(const G&, std::string, VTK::HDFTransientOptions = {}) -> VTKHDFImageGridTimeSeriesWriter<G, NullCommunicator>;
409template<Concepts::ImageGrid G, Concepts::Communicator C>
410VTKHDFImageGridTimeSeriesWriter(const G&, const C&, std::string, VTK::HDFTransientOptions = {}) -> VTKHDFImageGridTimeSeriesWriter<G, C>;
411
412
413namespace Traits {
414
415template<typename... Args>
416struct WritesConnectivity<VTKHDFImageGridWriter<Args...>> : public std::false_type {};
417
418template<typename... Args>
419struct WritesConnectivity<VTKHDFImageGridTimeSeriesWriter<Args...>> : public std::false_type {};
420
421} // namespace Traits
422} // namespace GridFormat
423
424#endif // GRIDFORMAT_HAVE_HIGH_FIVE
425#endif // GRIDFORMAT_VTK_HDF_IMAGE_GRID_WRITER_HPP_
Abstract base class for grid file writers.
Definition: writer.hpp:303
Abstract base class for time series file writers.
Definition: writer.hpp:344
Writer for the transient VTK HDF file format for image grids.
Definition: hdf_image_grid_writer.hpp:396
Definition: hdf_image_grid_writer.hpp:40
Writer for the VTK HDF file format for image grids.
Definition: hdf_image_grid_writer.hpp:385
Grid concepts.
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.
Can be specialized by writers in case the file format does not contain connectivity information.
Definition: writer.hpp:36
Helper class to store processor offsets in parallel writes.
Definition: hdf_common.hpp:37
Options for transient vtk-hdf file formats.
Definition: hdf_common.hpp:26
bool static_grid
Set to true the grid is the same for all time steps (will only be written once)
Definition: hdf_common.hpp:27
bool static_meta_data
Set to true if the metadata is same for all time steps (will only be written once)
Definition: hdf_common.hpp:28
Options that writer implementations can pass to the base class.
Definition: writer.hpp:59
bool use_structured_grid_ordering
Use row-major structured grid ordering.
Definition: writer.hpp:60
Common functionality for VTK writers.