8#ifndef GRIDFORMAT_VTK_PVTR_WRITER_HPP_
9#define GRIDFORMAT_VTK_PVTR_WRITER_HPP_
20#include <gridformat/common/ranges.hpp>
21#include <gridformat/common/exceptions.hpp>
22#include <gridformat/common/lvalue_reference.hpp>
24#include <gridformat/parallel/communication.hpp>
25#include <gridformat/parallel/helpers.hpp>
27#include <gridformat/grid/grid.hpp>
28#include <gridformat/xml/element.hpp>
38template<Concepts::RectilinearGrid Grid,
39 Concepts::Communicator Communicator>
42 using CT = CoordinateType<Grid>;
44 static constexpr std::size_t space_dim = 3;
45 static constexpr std::size_t dim = dimension<Grid>;
46 static constexpr int root_rank = 0;
49 explicit PVTRWriter(LValueReferenceOf<const Grid> grid,
52 :
ParentType(grid.get(),
".pvtr",
true, xml_opts)
56 const Communicator& communicator()
const {
64 return PVTRWriter{this->grid(), _comm, std::move(xml_opts)};
67 void _write(std::ostream&)
const override {
69 "PVTRWriter does not support direct export into stream. "
70 "Use overload with filename instead!"
74 virtual void _write(
const std::string& filename_with_ext)
const override {
75 const auto& local_extents = extents(this->grid());
76 const auto [origin, is_negative_axis] = _get_origin_and_orientations();
79 const auto all_origins = Parallel::gather(_comm, origin, root_rank);
80 const auto all_extents = Parallel::gather(_comm, local_extents, root_rank);
81 const auto [exts_begin, exts_end, whole_extent, _] = helper.compute_extents_and_origin(
87 const auto my_whole_extent = Parallel::broadcast(_comm, whole_extent, root_rank);
88 const auto my_extent_offset = Parallel::scatter(_comm, Ranges::flat(exts_begin), root_rank);
90 _write_piece(filename_with_ext, Ranges::to_array<dim>(my_extent_offset), {my_whole_extent});
91 Parallel::barrier(_comm);
92 if (Parallel::rank(_comm) == 0)
93 _write_pvtr_file(filename_with_ext, my_whole_extent, exts_begin, exts_end);
94 Parallel::barrier(_comm);
97 auto _get_origin_and_orientations()
const {
98 std::array<CT, dim> origin;
99 std::array<bool, dim> is_negative_axis;
100 for (
unsigned dir = 0; dir < dim; ++dir) {
101 std::array<CT, 2> ordinates_01{0, 0};
102 const auto& dir_ordinates = ordinates(this->grid(), dir);
103 std::ranges::copy(std::views::take(dir_ordinates, 2), ordinates_01.begin());
104 origin[dir] = ordinates_01[0];
105 is_negative_axis[dir] = ordinates_01[1] - ordinates_01[0] < CT{0};
107 return std::make_tuple(origin, is_negative_axis);
110 void _write_piece(
const std::string& par_filename,
111 const std::array<std::size_t, dim>& offset,
113 auto writer =
VTRWriter{this->grid(), this->_xml_opts}
114 .as_piece_for(std::move(domain))
115 .with_offset(offset);
116 this->copy_fields(writer);
117 writer.write(PVTK::piece_basefilename(par_filename, Parallel::rank(_comm)));
120 void _write_pvtr_file(
const std::string& filename_with_ext,
121 const std::array<std::size_t, dim>& extents,
122 const std::vector<std::array<std::size_t, dim>>& proc_extents_begin,
123 const std::vector<std::array<std::size_t, dim>>& proc_extents_end)
const {
124 std::ofstream file_stream(filename_with_ext, std::ios::out);
126 XMLElement pvtk_xml(
"VTKFile");
127 pvtk_xml.set_attribute(
"type",
"PRectilinearGrid");
129 XMLElement& grid = pvtk_xml.add_child(
"PRectilinearGrid");
130 grid.set_attribute(
"WholeExtent", VTK::CommonDetail::extents_string(extents));
132 XMLElement& ppoint_data = grid.add_child(
"PPointData");
133 XMLElement& pcell_data = grid.add_child(
"PCellData");
134 XMLElement& pcoords = grid.add_child(
"PCoordinates");
135 std::visit([&] (
const auto& encoder) {
136 std::visit([&] (
const auto& data_format) {
139 std::ranges::for_each(this->_point_field_names(), [&] (
const std::string& name) {
140 pdata_helper.add(name, this->_get_point_field(name));
142 std::ranges::for_each(this->_cell_field_names(), [&] (
const std::string& name) {
143 cdata_helper.add(name, this->_get_cell_field(name));
146 std::visit([&] <
typename T> (
const Precision<T>& prec) {
147 for (
unsigned int i = 0; i < space_dim; ++i) {
148 XMLElement& pdata_array = pcoords.add_child(
"PDataArray");
149 pdata_array.set_attribute(
"NumberOfComponents",
"1");
150 pdata_array.set_attribute(
"Name",
"X_" + std::to_string(i));
151 pdata_array.set_attribute(
"format", VTK::data_format_name(encoder, data_format));
152 pdata_array.set_attribute(
"type", VTK::attribute_name(prec));
154 }, this->_xml_settings.coordinate_precision);
155 }, this->_xml_settings.data_format);
156 }, this->_xml_settings.encoder);
158 std::ranges::for_each(Parallel::ranks(_comm), [&] (
int rank) {
159 std::array<std::size_t, dim> extents_begin;
160 std::array<std::size_t, dim> extents_end;
161 for (
unsigned dir = 0; dir < dim; ++dir) {
162 extents_begin[dir] = proc_extents_begin[rank][dir];
163 extents_end[dir] = proc_extents_end[rank][dir];
165 auto& piece = grid.add_child(
"Piece");
166 piece.set_attribute(
"Extent", VTK::CommonDetail::extents_string(extents_begin, extents_end));
167 piece.set_attribute(
"Source", std::filesystem::path{
168 PVTK::piece_basefilename(filename_with_ext, rank) +
".vtr"
172 this->_set_default_active_fields(pvtk_xml.get_child(
"PRectilinearGrid"));
173 write_xml_with_version_header(pvtk_xml, file_stream, Indentation{{.width = 2}});
177template<
typename G, Concepts::Communicator C>
182template<
typename... Args>
Helper function for writing parallel VTK files.
Definition: vtr_writer.hpp:43