GridFormat 0.4.0
I/O-Library for grid-like data structures
Loading...
Searching...
No Matches
common.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_COMMON_HPP_
9#define GRIDFORMAT_VTK_COMMON_HPP_
10
11#include <ranges>
12#include <cassert>
13#include <utility>
14#include <type_traits>
15#include <algorithm>
16#include <array>
17#include <cmath>
18
20#include <gridformat/common/concepts.hpp>
21#include <gridformat/common/exceptions.hpp>
22#include <gridformat/common/precision.hpp>
23#include <gridformat/common/serialization.hpp>
24#include <gridformat/common/md_layout.hpp>
25#include <gridformat/common/ranges.hpp>
26#include <gridformat/common/matrix.hpp>
27#include <gridformat/common/type_traits.hpp>
28#include <gridformat/common/string_conversion.hpp>
29#include <gridformat/common/flat_index_mapper.hpp>
30#include <gridformat/common/field_transformations.hpp>
32
33#include <gridformat/grid/entity_fields.hpp>
34#include <gridformat/grid/cell_type.hpp>
36#include <gridformat/grid/_detail.hpp>
37#include <gridformat/grid/grid.hpp>
38
39#ifndef DOXYGEN
40namespace GridFormat::Encoding { struct Ascii; struct Base64; struct RawBinary; }
41#endif // DOXYGEN
42
43namespace GridFormat::VTK {
44
47
48namespace DataFormat {
49
50struct Inlined {};
51struct Appended {};
52
53inline constexpr Inlined inlined;
54inline constexpr Appended appended;
55
56} // namespace DataFormat
57
58#ifndef DOXYGEN
59namespace Traits {
60
61template<typename T> struct ProducesValidXML;
62template<> struct ProducesValidXML<Encoding::Ascii> : public std::true_type {};
63template<> struct ProducesValidXML<Encoding::Base64> : public std::true_type {};
64template<> struct ProducesValidXML<Encoding::RawBinary> : public std::false_type {};
65
66} // namespace Traits
67
68template<typename Encoder>
69inline constexpr bool produces_valid_xml(const Encoder&) {
70 static_assert(
71 is_complete<Traits::ProducesValidXML<Encoder>>,
72 "Traits::ProducesValidXML was not specialized for the given encoder"
73 );
74 return Traits::ProducesValidXML<Encoder>::value;
75}
76
77inline constexpr std::uint8_t cell_type_number(CellType t) {
78 switch (t) {
79 case (CellType::vertex): return 1;
80 case (CellType::segment): return 3;
81 case (CellType::triangle): return 5;
82 case (CellType::pixel): return 8;
83 case (CellType::quadrilateral): return 9;
84 case (CellType::polygon): return 7;
85 case (CellType::tetrahedron): return 10;
86 case (CellType::hexahedron): return 12;
87 case (CellType::voxel): return 11;
88 case (CellType::lagrange_segment): return 68;
89 case (CellType::lagrange_triangle): return 69;
90 case (CellType::lagrange_quadrilateral): return 70;
91 case (CellType::lagrange_tetrahedron): return 71;
92 case (CellType::lagrange_hexahedron): return 72;
93 }
94
95 throw NotImplemented("VTK cell type number for the given cell type");
96}
97
98inline constexpr CellType cell_type(std::uint8_t vtk_id) {
99 switch (vtk_id) {
100 case 1: return CellType::vertex;
101 case 3: return CellType::segment;
102 case 5: return CellType::triangle;
103 case 8: return CellType::pixel;
104 case 9: return CellType::quadrilateral;
105 case 7: return CellType::polygon;
106 case 10: return CellType::tetrahedron;
107 case 12: return CellType::hexahedron;
108 case 11: return CellType::voxel;
109 case 68: return CellType::lagrange_segment;
110 case 69: return CellType::lagrange_triangle;
111 case 70: return CellType::lagrange_quadrilateral;
112 case 71: return CellType::lagrange_tetrahedron;
113 case 72: return CellType::lagrange_hexahedron;
114 }
115
116 throw NotImplemented("Cell type for the given VTK cell type number: " + std::to_string(vtk_id));
117}
118
119FieldPtr make_vtk_field(FieldPtr field) {
120 const auto layout = field->layout();
121 if (layout.dimension() < 2)
122 return field;
123 // (maybe) make vector/tensor fields 3d
124 if (std::ranges::all_of(
125 std::views::iota(std::size_t{1}, layout.dimension()),
126 [&] (const std::size_t codim) { return layout.extent(codim) < 3; }
127 ))
128 return transform(field, FieldTransformation::extend_all_to(3));
129 return field;
130}
131
132template<std::derived_from<Field> F>
133 requires(!std::is_lvalue_reference_v<F>)
134FieldPtr make_vtk_field(F&& field) {
135 return make_vtk_field(make_field_ptr(std::forward<F>(field)));
136}
137
138template<typename ctype, GridDetail::ExposesPointRange Grid>
139auto make_coordinates_field(const Grid& grid, bool structured_grid_ordering) {
140 return make_vtk_field(PointField{
141 grid,
142 [&] (const auto& point) { return coordinates(grid, point); },
143 structured_grid_ordering,
144 Precision<ctype>{}
145 });
146}
147
148template<typename HeaderType = std::size_t,
149 Concepts::UnstructuredGrid Grid,
150 std::ranges::forward_range Cells,
151 typename PointMap>
152 requires(std::is_lvalue_reference_v<PointMap>)
153auto make_connectivity_field(const Grid& grid,
154 Cells&& cells,
155 PointMap&& map) {
156 class ConnectivityField : public Field {
157 public:
158 explicit ConnectivityField(const Grid& g,
159 Cells&& cells,
160 PointMap&& map)
161 : _grid(g)
162 , _cells{std::forward<Cells>(cells)}
163 , _point_map{std::forward<PointMap>(map)} {
164 _num_values = 0;
165 std::ranges::for_each(_cells, [&] (const auto& cell) {
166 _num_values += number_of_points(_grid, cell);
167 });
168 }
169
170 private:
171 MDLayout _layout() const override { return MDLayout{{_num_values}}; }
172 DynamicPrecision _precision() const override { return Precision<HeaderType>{}; }
173 Serialization _serialized() const override {
174 Serialization serialization(sizeof(HeaderType)*_num_values);
175 HeaderType* data = serialization.as_span_of<HeaderType>().data();
176
177 std::size_t i = 0;
178 std::ranges::for_each(_cells, [&] (const auto& cell) {
179 std::ranges::for_each(points(_grid, cell), [&] (const auto& point) {
180 data[i] = _point_map.at(id(_grid, point));
181 i++;
182 });
183 });
184 return serialization;
185 }
186
187 const Grid& _grid;
188 LVReferenceOrValue<Cells> _cells;
189 LVReferenceOrValue<PointMap> _point_map;
190 HeaderType _num_values;
191 } _field{grid, std::forward<Cells>(cells), std::forward<PointMap>(map)};
192
193 return make_vtk_field(std::move(_field));
194}
195
196template<typename HeaderType = std::size_t,
197 Concepts::UnstructuredGrid Grid,
198 typename PointMap>
199auto make_connectivity_field(const Grid& grid, PointMap&& map) {
200 return make_connectivity_field<HeaderType>(grid, cells(grid), std::forward<PointMap>(map));
201}
202
203template<typename HeaderType = std::size_t,
204 Concepts::UnstructuredGrid Grid,
205 std::ranges::range Cells>
206auto make_offsets_field(const Grid& grid, Cells&& cells) {
207 class OffsetField : public Field {
208 public:
209 explicit OffsetField(const Grid& g, Cells&& cells)
210 : _grid(g)
211 , _cells{std::forward<Cells>(cells)}
212 , _num_cells{static_cast<HeaderType>(Ranges::size(_cells))}
213 {}
214
215 private:
216 MDLayout _layout() const override { return MDLayout{{_num_cells}}; }
217 DynamicPrecision _precision() const override { return Precision<HeaderType>{}; }
218 Serialization _serialized() const override {
219 Serialization serialization(sizeof(HeaderType)*_num_cells);
220 HeaderType* data = serialization.as_span_of<HeaderType>().data();
221
222 std::size_t i = 0;
223 HeaderType offset = 0;
224 std::ranges::for_each(_cells, [&] (const auto& cell) {
225 offset += number_of_points(_grid, cell);
226 data[i] = offset;
227 i++;
228 });
229 return serialization;
230 }
231
232 const Grid& _grid;
233 LVReferenceOrValue<Cells> _cells;
234 HeaderType _num_cells;
235 } _field{grid, std::forward<Cells>(cells)};
236
237 return make_vtk_field(std::move(_field));
238}
239
240template<typename HeaderType = std::size_t, Concepts::UnstructuredGrid Grid>
241auto make_offsets_field(const Grid& grid) {
242 return make_offsets_field<HeaderType>(grid, cells(grid));
243}
244
245template<Concepts::UnstructuredGrid Grid>
246auto make_cell_types_field(const Grid& grid) {
247 return make_vtk_field(CellField{
248 grid,
249 [&] (const auto& cell) {
250 return VTK::cell_type_number(type(grid, cell));
251 },
252 false
253 });
254}
255
256inline auto active_array_attribute_for_rank(unsigned int rank) {
257 if (rank > 2)
258 throw ValueError("Rank must be <= 2");
259 static constexpr std::array attributes{"Scalars", "Vectors", "Tensors"};
260 return attributes[rank];
261}
262
263namespace CommonDetail {
264
265 template<Concepts::StaticallySizedRange R>
266 std::string number_string_3d(const R& r) {
267 static_assert(static_size<R> >= 1 || static_size<R> <= 3);
268 if constexpr (static_size<R> == 3)
269 return as_string(r);
270 if constexpr (static_size<R> == 2)
271 return as_string(r) + " 0";
272 if constexpr (static_size<R> == 1)
273 return as_string(r) + " 0 0";
274 }
275
276 template<Concepts::StaticallySizedMDRange<2> R>
277 std::string direction_string(const R& basis) {
278 Matrix m{basis};
279 m.transpose();
280
281 std::string result = "";
282 std::ranges::for_each(m, [&] (const auto& row) {
283 if (result != "")
284 result += " ";
285 result += number_string_3d(row);
286 });
287 for (int i = static_size<R>; i < 3; ++i)
288 result += " 0 0 0";
289 return result;
290 }
291
292 template<Concepts::StaticallySizedRange R1,
293 Concepts::StaticallySizedRange R2>
294 requires(std::integral<std::ranges::range_value_t<R1>> and
295 std::integral<std::ranges::range_value_t<R2>>)
296 std::array<std::size_t, 6> get_extents(const R1& from, const R2& to) {
297 static_assert(static_size<R1> == static_size<R2>);
298 static_assert(static_size<R1> <= 3);
299
300 int i = 0;
301 auto result = Ranges::filled_array<6>(std::size_t{0});
302 auto it1 = std::ranges::begin(from);
303 auto it2 = std::ranges::begin(to);
304 for (; it1 != std::ranges::end(from); ++it1, ++it2, ++i) {
305 result[i*2 + 0] = *it1;
306 result[i*2 + 1] = *it2;
307 }
308 return result;
309 }
310
311 template<Concepts::StaticallySizedRange R>
312 std::array<std::size_t, 6> get_extents(const R& to) {
313 using T = std::ranges::range_value_t<R>;
314 return get_extents(Ranges::filled_array<static_size<R>>(T{0}), to);
315 }
316
317 template<Concepts::StaticallySizedRange R1,
318 Concepts::StaticallySizedRange R2>
319 std::string extents_string(const R1& from, const R2& to) {
320 return as_string(get_extents(from, to));
321 }
322
323 template<Concepts::StaticallySizedRange R>
324 std::string extents_string(const R& r) {
325 using T = std::ranges::range_value_t<R>;
326 return extents_string(Ranges::filled_array<static_size<R>>(T{0}), r);
327 }
328
329 template<Concepts::StructuredEntitySet Grid>
330 requires(!Concepts::StaticallySizedRange<Grid>)
331 std::string extents_string(const Grid& grid) {
332 return extents_string(extents(grid));
333 }
334
335 template<Concepts::StaticallySizedRange Spacing>
336 auto structured_grid_axis_orientation(const Spacing& spacing) {
337 std::array<bool, static_size<Spacing>> result;
338 std::ranges::copy(
339 spacing | std::views::transform([&] <typename CT> (const CT dx) { return dx <= CT{0}; }),
340 result.begin()
341 );
342 return result;
343 }
344
345 std::size_t number_of_entities(const std::array<std::size_t, 6>& extents) {
346 return std::max(extents[1] - extents[0], std::size_t{1})
347 *std::max(extents[3] - extents[2], std::size_t{1})
348 *std::max(extents[5] - extents[4], std::size_t{1});
349 }
350
351 unsigned int structured_grid_dimension(const std::array<std::size_t, 3>& cells_per_direction) {
352 return std::ranges::count_if(cells_per_direction, [] (const std::size_t e) { return e > 0; });
353 }
354
355 template<typename T>
356 std::array<T, 3> compute_location(const std::array<T, 3>& origin,
357 const std::array<T, 3>& coordinate,
358 const std::array<T, 9>& direction) {
359 const auto& [x, y, z] = coordinate;
360 return {
361 origin[0] + x*direction[0] + y*direction[1] + z*direction[2],
362 origin[1] + x*direction[3] + y*direction[4] + z*direction[5],
363 origin[2] + x*direction[6] + y*direction[7] + z*direction[8]
364 };
365 }
366
367 template<typename T>
368 std::array<T, 3> compute_piece_origin(const std::array<T, 3>& global_origin,
369 const std::array<T, 3>& spacing,
370 const std::array<std::size_t, 3>& extents_begin,
371 const std::array<T, 9>& direction) {
372 return compute_location(
373 global_origin,
374 {
375 spacing[0]*static_cast<T>(extents_begin[0]),
376 spacing[1]*static_cast<T>(extents_begin[1]),
377 spacing[2]*static_cast<T>(extents_begin[2])
378 },
379 direction
380 );
381 }
382
383 template<typename T>
384 Serialization serialize_structured_points(const std::array<std::size_t, 6>& extents,
385 const std::array<T, 3>& origin,
386 const std::array<T, 3>& spacing,
387 const std::array<T, 9>& direction) {
388 const MDLayout layout{{
389 extents[1] - extents[0],
390 extents[3] - extents[2],
391 extents[5] - extents[4]
392 }};
393 const FlatIndexMapper mapper{layout};
394 const auto piece_origin = compute_piece_origin(
395 origin, spacing, {extents[0], extents[2], extents[4]}, direction
396 );
397
398 static constexpr unsigned int vtk_space_dim = 3;
399 Serialization result(layout.number_of_entries()*sizeof(T)*vtk_space_dim);
400 auto span_out = result.as_span_of(Precision<T>{});
401 for (const auto& md_index : MDIndexRange{layout}) {
402 const auto offset = mapper.map(md_index)*vtk_space_dim;
403 assert(offset + 2 < span_out.size());
404 std::ranges::copy(
405 compute_location(piece_origin, {
406 static_cast<T>(md_index.get(0))*spacing[0],
407 static_cast<T>(md_index.get(1))*spacing[1],
408 static_cast<T>(md_index.get(2))*spacing[2]
409 }, direction),
410 span_out.data() + offset
411 );
412 }
413 return result;
414 }
415
416 template<typename Visitor>
417 void visit_structured_cells(const Visitor& visitor,
418 const std::array<std::size_t, 6>& extents,
419 const bool is_axis_aligned = true) {
420 std::array<CellType, 4> grid_dim_to_cell_type{
421 CellType::vertex,
422 CellType::segment,
423 (is_axis_aligned ? CellType::pixel : CellType::quadrilateral),
424 (is_axis_aligned ? CellType::voxel : CellType::hexahedron)
425 };
426
427 std::array<std::size_t, 3> counts{
428 extents[1] - extents[0],
429 extents[3] - extents[2],
430 extents[5] - extents[4]
431 };
432
433 const std::size_t grid_dim = structured_grid_dimension(counts);
434 if (grid_dim == 0)
435 throw ValueError("Grid must be at least 1d");
436
437 const MDLayout point_layout{Ranges::incremented(counts, 1)};
438 const FlatIndexMapper point_mapper{point_layout};
439 const auto x_offset = grid_dim > 1 ? point_layout.extent(0) : std::size_t{0};
440 const auto y_offset = grid_dim > 2 ? point_layout.extent(0)*point_layout.extent(1) : std::size_t{0};
441
442 // avoid zero counts s.t. the index range does not degenerate
443 std::ranges::for_each(counts, [] (std::size_t& count) { count = std::max(count, std::size_t{1}); });
444 std::vector<std::size_t> corners(std::pow(2, grid_dim), 0);
445
446 const MDIndexRange index_range{MDLayout{counts}};
447 if (grid_dim == 1) {
448 std::ranges::for_each(index_range, [&] (const auto& md_index) {
449 const auto p0 = point_mapper.map(md_index);
450 corners = {p0, p0 + 1};
451 visitor(grid_dim_to_cell_type[grid_dim], corners);
452 });
453 } else if (grid_dim == 2) {
454 std::ranges::for_each(index_range, [&] (const auto& md_index) {
455 const auto p0 = point_mapper.map(md_index);
456 corners = {p0, p0 + 1, p0 + x_offset, p0 + 1 + x_offset};
457 visitor(grid_dim_to_cell_type[grid_dim], corners);
458 });
459 } else {
460 std::ranges::for_each(index_range, [&] (const auto& md_index) {
461 const auto p0 = point_mapper.map(md_index);
462 corners = {
463 p0, p0 + 1, p0 + x_offset, p0 + 1 + x_offset,
464 p0 + y_offset, p0 + y_offset + 1, p0 + y_offset + x_offset, p0 + 1 + y_offset + x_offset
465 };
466 visitor(grid_dim_to_cell_type[grid_dim], corners);
467 });
468 }
469 }
470
471} // namespace CommonDetail
472#endif // DOXYGEN
473
475} // namespace GridFormat::VTK
476
477
478#endif // GRIDFORMAT_VTK_COMMON_HPP_
Grid concepts.
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
Appended data format (all data is appended at the end of the xml file)
Definition: common.hpp:51
Inline data format (inside xml elements)
Definition: common.hpp:50
constexpr Inlined inlined
Instance of the inline data format.
Definition: common.hpp:53
constexpr Appended appended
Instance of the appended data format.
Definition: common.hpp:54