dune-vtk 2.8
Loading...
Searching...
No Matches
vtkwriterinterface.impl.hh
Go to the documentation of this file.
1#pragma once
2
3#include <algorithm>
4#include <iomanip>
5#include <iostream>
6#include <iterator>
7#include <fstream>
8#include <limits>
9#include <sstream>
10#include <string>
11
12#if HAVE_VTK_ZLIB
13#include <zlib.h>
14#endif
15
16#include <dune/geometry/referenceelements.hh>
17#include <dune/geometry/type.hh>
18
23
24namespace Dune {
25
26template <class GV, class DC>
27std::string VtkWriterInterface<GV,DC>
28 ::write (std::string const& fn, std::optional<std::string> dir) const
29{
30 dataCollector_->update();
31
32 auto p = Vtk::Path(fn);
33 auto name = p.stem();
34 p.removeFilename();
35
36 Vtk::Path fn_dir = p;
37 Vtk::Path data_dir = dir ? Vtk::Path(*dir) : fn_dir;
38 Vtk::Path rel_dir = Vtk::relative(data_dir, fn_dir);
39
40 std::string serial_fn = data_dir.string() + '/' + name.string();
41 std::string parallel_fn = fn_dir.string() + '/' + name.string();
42 std::string rel_fn = rel_dir.string() + '/' + name.string();
43
44 if (comm().size() > 1)
45 serial_fn += "_p" + std::to_string(comm().rank());
46
47 std::string outputFilename;
48
49 { // write serial file
50 outputFilename = serial_fn + "." + fileExtension();
51 std::ofstream serial_out(outputFilename, std::ios_base::ate | std::ios::binary);
52 assert(serial_out.is_open());
53
54 serial_out.imbue(std::locale::classic());
55 serial_out << std::setprecision(datatype_ == Vtk::DataTypes::FLOAT32
56 ? std::numeric_limits<float>::max_digits10
57 : std::numeric_limits<double>::max_digits10);
58
59 writeSerialFile(serial_out);
60 }
61
62 if (comm().size() > 1 && comm().rank() == 0) {
63 // write parallel filee
64 outputFilename = parallel_fn + ".p" + fileExtension();
65 std::ofstream parallel_out(outputFilename, std::ios_base::ate | std::ios::binary);
66 assert(parallel_out.is_open());
67
68 parallel_out.imbue(std::locale::classic());
69 parallel_out << std::setprecision(datatype_ == Vtk::DataTypes::FLOAT32
70 ? std::numeric_limits<float>::max_digits10
71 : std::numeric_limits<double>::max_digits10);
72
73 writeParallelFile(parallel_out, rel_fn, comm().size());
74 }
75
76 return outputFilename;
77}
78
79
80template <class GV, class DC>
82 ::writeData (std::ofstream& out, std::vector<pos_type>& offsets,
83 VtkFunction const& fct, PositionTypes type,
84 std::optional<std::size_t> timestep) const
85{
86 out << "<DataArray Name=\"" << fct.name() << "\" type=\"" << to_string(fct.dataType()) << "\""
87 << " NumberOfComponents=\"" << fct.numComponents() << "\" format=\"" << (format_ == Vtk::FormatTypes::ASCII ? "ascii\"" : "appended\"");
88 if (timestep)
89 out << " TimeStep=\"" << *timestep << "\"";
90
91 if (format_ == Vtk::FormatTypes::ASCII) {
92 out << ">\n";
93 if (type == POINT_DATA)
94 writeValuesAscii(out, dataCollector_->template pointData<double>(fct));
95 else
96 writeValuesAscii(out, dataCollector_->template cellData<double>(fct));
97 out << "</DataArray>\n";
98 } else {
99 out << " offset=";
100 offsets.push_back(out.tellp());
101 out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' ');
102 out << "/>\n";
103 }
104}
105
106
107template <class GV, class DC>
108void VtkWriterInterface<GV,DC>
109 ::writePoints (std::ofstream& out, std::vector<pos_type>& offsets,
110 std::optional<std::size_t> timestep) const
111{
112 out << "<DataArray type=\"" << to_string(datatype_) << "\""
113 << " NumberOfComponents=\"3\" format=\"" << (format_ == Vtk::FormatTypes::ASCII ? "ascii\"" : "appended\"");
114 if (timestep)
115 out << " TimeStep=\"" << *timestep << "\"";
116
117 if (format_ == Vtk::FormatTypes::ASCII) {
118 out << ">\n";
119 writeValuesAscii(out, dataCollector_->template points<double>());
120 out << "</DataArray>\n";
121 } else {
122 out << " offset=";
123 offsets.push_back(out.tellp());
124 out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' ');
125 out << "/>\n";
126 }
127}
128
129
130template <class GV, class DC>
132 ::writeDataAppended (std::ofstream& out, std::vector<std::uint64_t>& blocks) const
133{
134 for (auto const& v : pointData_) {
135 Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(v.dataType(), headertype_,
136 [&](auto f, auto h) {
137 using F = typename decltype(f)::type;
138 using H = typename decltype(h)::type;
139 blocks.push_back(this->template writeValuesAppended<H>(out, dataCollector_->template pointData<F>(v)));
140 });
141 }
142
143 for (auto const& v : cellData_) {
144 Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(v.dataType(), headertype_,
145 [&](auto f, auto h) {
146 using F = typename decltype(f)::type;
147 using H = typename decltype(h)::type;
148 blocks.push_back(this->template writeValuesAppended<H>(out, dataCollector_->template cellData<F>(v)));
149 });
150 }
151}
152
153
154template <class GV, class DC>
156 ::writeAppended (std::ofstream& out, std::vector<pos_type> const& offsets) const
157{
158 if (is_a(format_, Vtk::FormatTypes::APPENDED)) {
159 out << "<AppendedData encoding=\"raw\">\n_";
160 std::vector<std::uint64_t> blocks;
161 writeGridAppended(out, blocks);
162 writeDataAppended(out, blocks);
163 out << "</AppendedData>\n";
164 pos_type appended_pos = out.tellp();
165
166 pos_type offset = 0;
167 for (std::size_t i = 0; i < offsets.size(); ++i) {
168 out.seekp(offsets[i]);
169 out << '"' << offset << '"';
170 offset += pos_type(blocks[i]);
171 }
172
173 out.seekp(appended_pos);
174 }
175}
176
177
178namespace Impl {
179
180 template <class T, std::enable_if_t<(sizeof(T)>1), int> = 0>
181 inline T const& printable (T const& t) { return t; }
182
183 inline std::int16_t printable (std::int8_t c) { return std::int16_t(c); }
184 inline std::uint16_t printable (std::uint8_t c) { return std::uint16_t(c); }
185
186} // end namespace Impl
187
188
189template <class GV, class DC>
190 template <class T>
191void VtkWriterInterface<GV,DC>
192 ::writeValuesAscii (std::ofstream& out, std::vector<T> const& values) const
194 assert(is_a(format_, Vtk::FormatTypes::ASCII) && "Function should by called only in ascii mode!\n");
195 std::size_t i = 0;
196 for (auto const& v : values)
197 out << Impl::printable(v) << (++i % 6 != 0 ? ' ' : '\n');
198 if (i % 6 != 0)
199 out << '\n';
201
202template <class GV, class DC>
204 ::writeHeader (std::ofstream& out, std::string const& type) const
206 out << "<VTKFile"
207 << " type=\"" << type << "\""
208 << " version=\"1.0\""
209 << " header_type=\"" << to_string(headertype_) << "\"";
211 out << " byte_order=\"" << getEndian() << "\"";
212 if (compressor_ != Vtk::CompressorTypes::NONE)
213 out << " compressor=\"" << to_string(compressor_) << "\"";
214 out << ">\n";
216
217
218namespace Impl {
219
220template <class T>
221std::uint64_t writeValuesToBuffer (std::size_t max_num_values, unsigned char* buffer,
222 std::vector<T> const& vec, std::size_t shift)
224 std::size_t num_values = std::min(max_num_values, vec.size()-shift);
225 std::uint64_t bs = num_values*sizeof(T);
226 std::memcpy(buffer, (unsigned char*)(vec.data()+shift), std::size_t(bs));
227 return bs;
228}
229
230inline std::uint64_t compressBuffer_zlib (unsigned char const* buffer, unsigned char* buffer_out,
231 std::uint64_t bs, std::uint64_t cbs, int level)
232{
233#if HAVE_VTK_ZLIB
234 uLongf uncompressed_space = uLongf(bs);
235 uLongf compressed_space = uLongf(cbs);
236
237 Bytef* out = reinterpret_cast<Bytef*>(buffer_out);
238 Bytef const* in = reinterpret_cast<Bytef const*>(buffer);
239
240 if (compress2(out, &compressed_space, in, uncompressed_space, level) != Z_OK) {
241 std::cerr << "Zlib error while compressing data.\n";
242 std::abort();
243 }
244
245 return compressed_space;
246#else
247 std::cerr << "Can not call writeCompressed without compression enabled!\n";
248 std::abort();
249 return 0;
250#endif
251}
252
253inline std::uint64_t compressBuffer_lz4 (unsigned char const* /* buffer */, unsigned char* /* buffer_out */,
254 std::uint64_t /* bs */, std::uint64_t /* cbs */, int /* level */)
255{
256#if HAVE_VTK_LZ4
257 std::cerr << "LZ4 Compression not yet implemented" << std::endl;
258 std::abort();
259 return 0;
260#else
261 std::cerr << "LZ4 Compression not supported. Provide the LZ4 package to CMake." << std::endl;
262 std::abort();
263 return 0;
264#endif
265}
266
267inline std::uint64_t compressBuffer_lzma (unsigned char const* /* buffer */, unsigned char* /* buffer_out */,
268 std::uint64_t /* bs */, std::uint64_t /* cbs */, int /* level */)
269{
270#if HAVE_VTK_LZMA
271 std::cerr << "LZMA Compression not yet implemented" << std::endl;
272 std::abort();
273 return 0;
274#else
275 std::cerr << "LZMA Compression not supported. Provide the LZMA package to CMake." << std::endl;
276 std::abort();
277 return 0;
278#endif
279}
280
281} // end namespace Impl
282
283template <class GV, class DC>
284 template <class HeaderType, class FloatType>
285std::uint64_t VtkWriterInterface<GV,DC>
286 ::writeValuesAppended (std::ofstream& out, std::vector<FloatType> const& values) const
287{
288 assert(is_a(format_, Vtk::FormatTypes::APPENDED) && "Function should by called only in appended mode!\n");
289 pos_type begin_pos = out.tellp();
290
291 HeaderType size = values.size() * sizeof(FloatType);
292
293 HeaderType num_full_blocks = size / block_size;
294 HeaderType last_block_size = size % block_size;
295 HeaderType num_blocks = num_full_blocks + (last_block_size > 0 ? 1 : 0);
296
297 // write block-size(s)
298 HeaderType zero = 0;
299 if (compressor_ != Vtk::CompressorTypes::NONE) {
300 out.write((char*)&num_blocks, sizeof(HeaderType));
301 out.write((char*)&block_size, sizeof(HeaderType));
302 out.write((char*)&last_block_size, sizeof(HeaderType));
303 for (HeaderType i = 0; i < num_blocks; ++i)
304 out.write((char*)&zero, sizeof(HeaderType));
305 } else {
306 out.write((char*)&size, sizeof(HeaderType));
307 }
308
309 HeaderType compressed_block_size = block_size + (block_size + 999)/1000 + 12;
310 std::vector<unsigned char> buffer(block_size);
311 std::vector<unsigned char> buffer_out;
312 if (compressor_ != Vtk::CompressorTypes::NONE)
313 buffer_out.resize(std::size_t(compressed_block_size));
314
315 std::size_t num_values = block_size / sizeof(FloatType);
316
317 std::vector<HeaderType> cbs(std::size_t(num_blocks), 0); // compressed block sizes
318 for (std::size_t i = 0; i < std::size_t(num_blocks); ++i) {
319 HeaderType bs = Impl::writeValuesToBuffer<FloatType>(num_values, buffer.data(), values, i*num_values);
320
321 switch (compressor_) {
323 out.write((char*)buffer.data(), bs);
324 break;
326 cbs[i] = Impl::compressBuffer_zlib(buffer.data(), buffer_out.data(), bs, compressed_block_size, compression_level);
327 out.write((char*)buffer_out.data(), cbs[i]);
328 break;
330 cbs[i] = Impl::compressBuffer_zlib(buffer.data(), buffer_out.data(), bs, compressed_block_size, compression_level);
331 out.write((char*)buffer_out.data(), cbs[i]);
332 break;
334 cbs[i] = Impl::compressBuffer_zlib(buffer.data(), buffer_out.data(), bs, compressed_block_size, compression_level);
335 out.write((char*)buffer_out.data(), cbs[i]);
336 break;
337 }
338 }
339
340 pos_type end_pos = out.tellp();
341 if (compressor_ != Vtk::CompressorTypes::NONE) {
342 out.seekp(begin_pos + std::streamoff(3*sizeof(HeaderType)));
343 out.write((char*)cbs.data(), std::streamsize(num_blocks*sizeof(HeaderType)));
344 out.seekp(end_pos);
345 }
346
347 return std::uint64_t(end_pos - begin_pos);
348}
349
350
351template <class GV, class DC>
353 ::getNames (std::vector<VtkFunction> const& data) const
354{
355 auto scalar = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.rangeType() == Vtk::RangeTypes::SCALAR; });
356 auto vector = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.rangeType() == Vtk::RangeTypes::VECTOR; });
357 auto tensor = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.rangeType() == Vtk::RangeTypes::TENSOR; });
358 return (scalar != data.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
359 + (vector != data.end() ? " Vectors=\"" + vector->name() + "\"" : "")
360 + (tensor != data.end() ? " Tensors=\"" + tensor->name() + "\"" : "");
361}
362
363} // end namespace Dune
Macro for wrapping error checks and throwing exceptions.
Definition: writer.hh:13
Path relative(Path const &a, Path const &b)
Find the path of a relative to directory of b
Definition: filesystem.cc:173
Wrapper class for functions allowing local evaluations.
Definition: function.hh:28
Vtk::DataTypes dataType() const
Return the VTK Datatype associated with the functions range type.
Definition: function.hh:211
std::string const & name() const
Return a name associated with the function.
Definition: function.hh:178
int numComponents() const
Return the number of components of the Range as it is written to the file.
Definition: function.hh:190
Definition: filesystem.hh:15
std::string string() const
Return the path as string.
Definition: filesystem.cc:27
Interface for file writers for the Vtk XML file formats.
Definition: vtkwriterinterface.hh:25
typename std::ostream::pos_type pos_type
Definition: vtkwriterinterface.hh:35
PositionTypes
Definition: vtkwriterinterface.hh:37