dune-vtk  0.2
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 
19 #include <dune/vtk/utility/enum.hh>
23 
24 namespace Dune {
25 
26 template <class GV, class DC>
27 std::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 
80 template <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 
107 template <class GV, class 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 
130 template <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 
154 template <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 
178 namespace 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 
189 template <class GV, class DC>
190  template <class T>
192  ::writeValuesAscii (std::ofstream& out, std::vector<T> const& values) const
193 {
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';
200 }
201 
202 template <class GV, class DC>
204  ::writeHeader (std::ofstream& out, std::string const& type) const
205 {
206  out << "<VTKFile"
207  << " type=\"" << type << "\""
208  << " version=\"1.0\""
209  << " header_type=\"" << to_string(headertype_) << "\"";
210  if (format_ != Vtk::FormatTypes::ASCII)
211  out << " byte_order=\"" << getEndian() << "\"";
212  if (compressor_ != Vtk::CompressorTypes::NONE)
213  out << " compressor=\"" << to_string(compressor_) << "\"";
214  out << ">\n";
215 }
216 
217 
218 namespace Impl {
219 
220 template <class T>
221 std::uint64_t writeValuesToBuffer (std::size_t max_num_values, unsigned char* buffer,
222  std::vector<T> const& vec, std::size_t shift)
223 {
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 
230 inline 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 
253 inline 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 
267 inline 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 
283 template <class GV, class DC>
284  template <class HeaderType, class FloatType>
285 std::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 
351 template <class GV, class DC>
352 std::string VtkWriterInterface<GV,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: datacollectorinterface.hh:9
std::string to_string(Vtk::FormatTypes type)
Definition: types.cc:12
Path relative(Path const &a, Path const &b)
Find the path of a relative to directory of b
Definition: filesystem.cc:173
@ APPENDED
Definition: types.hh:25
@ ASCII
Definition: types.hh:22
constexpr bool is_a(E a, Integer b)
Definition: enum.hh:12
@ NONE
Definition: types.hh:140
@ ZLIB
Definition: types.hh:141
@ LZ4
Definition: types.hh:142
@ LZMA
Definition: types.hh:143
@ FLOAT32
Definition: types.hh:56
Wrapper class for functions allowing local evaluations.
Definition: function.hh:27
Vtk::DataTypes dataType() const
Return the VTK Datatype associated with the functions range type.
Definition: function.hh:208
std::string const & name() const
Return a name associated with the function.
Definition: function.hh:175
int numComponents() const
Return the number of components of the Range as it is written to the file.
Definition: function.hh:187
Definition: filesystem.hh:15
std::string string() const
Return the path as string.
Definition: filesystem.cc:27
void writeAppended(std::ofstream &out, std::vector< pos_type > const &offsets) const
Definition: vtkwriterinterface.impl.hh:156
void writeValuesAscii(std::ofstream &out, std::vector< T > const &values) const
Definition: vtkwriterinterface.impl.hh:192
std::string getNames(std::vector< VtkFunction > const &data) const
Return PointData/CellData attributes for the name of the first scalar/vector/tensor DataArray.
Definition: vtkwriterinterface.impl.hh:353
typename std::ostream::pos_type pos_type
Definition: vtkwriterinterface.hh:35
void writeHeader(std::ofstream &out, std::string const &type) const
Definition: vtkwriterinterface.impl.hh:204
void writeDataAppended(std::ofstream &out, std::vector< std::uint64_t > &blocks) const
Definition: vtkwriterinterface.impl.hh:132
void writeData(std::ofstream &out, std::vector< pos_type > &offsets, VtkFunction const &fct, PositionTypes type, std::optional< std::size_t > timestep={}) const
Definition: vtkwriterinterface.impl.hh:82
void writePoints(std::ofstream &out, std::vector< pos_type > &offsets, std::optional< std::size_t > timestep={}) const
Definition: vtkwriterinterface.impl.hh:109
virtual std::string write(std::string const &fn, std::optional< std::string > dir={}) const override
Write the attached data to the file.
Definition: vtkwriterinterface.impl.hh:28
std::uint64_t writeValuesAppended(std::ofstream &out, std::vector< FloatType > const &values) const
Definition: vtkwriterinterface.impl.hh:286
PositionTypes
Definition: vtkwriterinterface.hh:37