vtkwriterinterface.impl.hh 11.7 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
#pragma once

3
#include <algorithm>
Praetorius, Simon's avatar
Praetorius, Simon committed
4
5
6
7
#include <iomanip>
#include <iostream>
#include <iterator>
#include <fstream>
8
#include <limits>
Praetorius, Simon's avatar
Praetorius, Simon committed
9
10
11
#include <sstream>
#include <string>

12
#if HAVE_VTK_ZLIB
13
14
15
#include <zlib.h>
#endif

Praetorius, Simon's avatar
Praetorius, Simon committed
16
17
18
#include <dune/geometry/referenceelements.hh>
#include <dune/geometry/type.hh>

19
#include <dune/vtk/utility/enum.hh>
20
#include <dune/vtk/utility/errors.hh>
21
22
#include <dune/vtk/utility/filesystem.hh>
#include <dune/vtk/utility/string.hh>
Praetorius, Simon's avatar
Praetorius, Simon committed
23

24
namespace Dune {
Praetorius, Simon's avatar
Praetorius, Simon committed
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
Praetorius, Simon's avatar
Praetorius, Simon committed
29
{
30
  dataCollector_->update();
31

32
  auto p = Vtk::Path(fn);
Praetorius, Simon's avatar
Praetorius, Simon committed
33
  auto name = p.stem();
34
  p.removeFilename();
Praetorius, Simon's avatar
Praetorius, Simon committed
35

36
37
38
  Vtk::Path fn_dir = p;
  Vtk::Path data_dir = dir ? Vtk::Path(*dir) : fn_dir;
  Vtk::Path rel_dir = Vtk::relative(data_dir, fn_dir);
Praetorius, Simon's avatar
Praetorius, Simon committed
39

40
41
  std::string serial_fn = data_dir.string() + '/' + name.string();
  std::string parallel_fn = fn_dir.string() + '/' + name.string();
42
  std::string rel_fn = rel_dir.string() + '/' + name.string();
43

44
45
  if (comm().size() > 1)
    serial_fn += "_p" + std::to_string(comm().rank());
Praetorius, Simon's avatar
Praetorius, Simon committed
46

47
  std::string outputFilename;
Stenger, Florian's avatar
Stenger, Florian committed
48

49
  { // write serial file
50
51
    outputFilename = serial_fn + "." + fileExtension();
    std::ofstream serial_out(outputFilename, std::ios_base::ate | std::ios::binary);
52
53
54
55
    assert(serial_out.is_open());

    serial_out.imbue(std::locale::classic());
    serial_out << std::setprecision(datatype_ == Vtk::FLOAT32
56
57
      ? std::numeric_limits<float>::max_digits10
      : std::numeric_limits<double>::max_digits10);
58
59
60

    writeSerialFile(serial_out);
  }
61

62
  if (comm().size() > 1 && comm().rank() == 0) {
63
64
65
    // write parallel filee
    outputFilename = parallel_fn + ".p" + fileExtension();
    std::ofstream parallel_out(outputFilename, std::ios_base::ate | std::ios::binary);
66
67
68
69
    assert(parallel_out.is_open());

    parallel_out.imbue(std::locale::classic());
    parallel_out << std::setprecision(datatype_ == Vtk::FLOAT32
70
71
      ? std::numeric_limits<float>::max_digits10
      : std::numeric_limits<double>::max_digits10);
72

73
    writeParallelFile(parallel_out, rel_fn, comm().size());
74
  }
75
76

  return outputFilename;
Praetorius, Simon's avatar
Praetorius, Simon committed
77
78
79
}


80
template <class GV, class DC>
81
void VtkWriterInterface<GV,DC>
Praetorius, Simon's avatar
Praetorius, Simon committed
82
  ::writeData (std::ofstream& out, std::vector<pos_type>& offsets,
83
               VtkFunction const& fct, PositionTypes type,
84
               std::optional<std::size_t> timestep) const
Praetorius, Simon's avatar
Praetorius, Simon committed
85
{
Praetorius, Simon's avatar
Praetorius, Simon committed
86
  out << "<DataArray Name=\"" << fct.name() << "\" type=\"" << to_string(fct.type()) << "\""
87
88
89
      << " NumberOfComponents=\"" << fct.ncomps() << "\" format=\"" << (format_ == Vtk::ASCII ? "ascii\"" : "appended\"");
  if (timestep)
    out << " TimeStep=\"" << *timestep << "\"";
Praetorius, Simon's avatar
Praetorius, Simon committed
90
91

  if (format_ == Vtk::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";
Praetorius, Simon's avatar
Praetorius, Simon committed
98
99
100
101
102
103
104
105
106
  } else {
    out << " offset=";
    offsets.push_back(out.tellp());
    out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' ');
    out << "/>\n";
  }
}


107
template <class GV, class DC>
108
void VtkWriterInterface<GV,DC>
109
  ::writePoints (std::ofstream& out, std::vector<pos_type>& offsets,
110
                std::optional<std::size_t> timestep) const
Praetorius, Simon's avatar
Praetorius, Simon committed
111
{
Praetorius, Simon's avatar
Praetorius, Simon committed
112
  out << "<DataArray type=\"" << to_string(datatype_) << "\""
113
114
115
      << " NumberOfComponents=\"3\" format=\"" << (format_ == Vtk::ASCII ? "ascii\"" : "appended\"");
  if (timestep)
    out << " TimeStep=\"" << *timestep << "\"";
Praetorius, Simon's avatar
Praetorius, Simon committed
116
117

  if (format_ == Vtk::ASCII) {
118
    out << ">\n";
119
    writeValuesAscii(out, dataCollector_->template points<double>());
120
    out << "</DataArray>\n";
Praetorius, Simon's avatar
Praetorius, Simon committed
121
122
123
124
125
126
127
128
  } else {
    out << " offset=";
    offsets.push_back(out.tellp());
    out << std::string(std::numeric_limits<std::uint64_t>::digits10 + 2, ' ');
    out << "/>\n";
  }
}

129

Praetorius, Simon's avatar
Praetorius, Simon committed
130
template <class GV, class DC>
131
132
void VtkWriterInterface<GV,DC>
  ::writeDataAppended (std::ofstream& out, std::vector<std::uint64_t>& blocks) const
Praetorius, Simon's avatar
Praetorius, Simon committed
133
{
134
  for (auto const& v : pointData_) {
135
136
137
138
139
140
    Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(v.type(), headertype_,
    [&](auto f, auto h) {
      using F = typename decltype(f)::type;
      using H = typename decltype(h)::type;
      blocks.push_back(this->template writeValuesAppended<H>(out, dataCollector_->template pointData<F>(v)));
    });
141
  }
142

143
  for (auto const& v : cellData_) {
144
145
146
147
148
149
    Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(v.type(), headertype_,
    [&](auto f, auto h) {
      using F = typename decltype(f)::type;
      using H = typename decltype(h)::type;
      blocks.push_back(this->template writeValuesAppended<H>(out, dataCollector_->template cellData<F>(v)));
    });
Praetorius, Simon's avatar
Praetorius, Simon committed
150
151
152
153
154
  }
}


template <class GV, class DC>
155
156
void VtkWriterInterface<GV,DC>
  ::writeAppended (std::ofstream& out, std::vector<pos_type> const& offsets) const
Praetorius, Simon's avatar
Praetorius, Simon committed
157
{
158
  if (is_a(format_, Vtk::FormatTypes::APPENDED)) {
159
160
161
162
163
164
165
166
167
168
169
170
171
    out << "<AppendedData encoding=\"raw\">\n_";
    std::vector<std::uint64_t> blocks;
    writeGridAppended(out, blocks);
    writeDataAppended(out, blocks);
    out << "</AppendedData>\n";
    pos_type appended_pos = out.tellp();

    pos_type offset = 0;
    for (std::size_t i = 0; i < offsets.size(); ++i) {
      out.seekp(offsets[i]);
      out << '"' << offset << '"';
      offset += pos_type(blocks[i]);
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
172

173
174
    out.seekp(appended_pos);
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
175
176
177
}


178
179
180
namespace Impl {

  template <class T, std::enable_if_t<(sizeof(T)>1), int> = 0>
Praetorius, Simon's avatar
Praetorius, Simon committed
181
  inline T const& printable (T const& t) { return t; }
182

Praetorius, Simon's avatar
Praetorius, Simon committed
183
184
  inline std::int16_t printable (std::int8_t c) { return std::int16_t(c); }
  inline std::uint16_t printable (std::uint8_t c) { return std::uint16_t(c); }
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208

} // end namespace Impl


template <class GV, class DC>
  template <class T>
void VtkWriterInterface<GV,DC>
  ::writeValuesAscii (std::ofstream& out, std::vector<T> const& values) const
{
  assert(is_a(format_, Vtk::ASCII) && "Function should by called only in ascii mode!\n");
  std::size_t i = 0;
  for (auto const& v : values)
    out << Impl::printable(v) << (++i % 6 != 0 ? ' ' : '\n');
  if (i % 6 != 0)
    out << '\n';
}

template <class GV, class DC>
void VtkWriterInterface<GV,DC>
  ::writeHeader (std::ofstream& out, std::string const& type) const
{
  out << "<VTKFile"
      << " type=\"" << type << "\""
      << " version=\"1.0\""
209
      << " header_type=\"" << to_string(headertype_) << "\"";
210
211
  if (format_ != Vtk::ASCII)
    out << " byte_order=\"" << getEndian() << "\"";
212
213
  if (compressor_ != Vtk::NONE)
    out << " compressor=\"" << to_string(compressor_) << "\"";
214
215
216
217
  out << ">\n";
}


218
namespace Impl {
Praetorius, Simon's avatar
Praetorius, Simon committed
219
220

template <class T>
221
222
std::uint64_t writeValuesToBuffer (std::size_t max_num_values, unsigned char* buffer,
                                   std::vector<T> const& vec, std::size_t shift)
Praetorius, Simon's avatar
Praetorius, Simon committed
223
224
225
226
227
228
229
{
  std::size_t num_values = std::min(max_num_values, vec.size()-shift);
  std::uint64_t bs = num_values*sizeof(T);
  std::memcpy(buffer, (unsigned char*)(vec.data()+shift), std::size_t(bs));
  return bs;
}

230
231
inline std::uint64_t compressBuffer_zlib (unsigned char const* buffer, unsigned char* buffer_out,
                                          std::uint64_t bs, std::uint64_t cbs, int level)
Praetorius, Simon's avatar
Praetorius, Simon committed
232
{
233
#if HAVE_VTK_ZLIB
Praetorius, Simon's avatar
Praetorius, Simon committed
234
235
236
237
238
239
240
241
242
243
244
245
246
  uLongf uncompressed_space = uLongf(bs);
  uLongf compressed_space = uLongf(cbs);

  Bytef* out = reinterpret_cast<Bytef*>(buffer_out);
  Bytef const* in = reinterpret_cast<Bytef const*>(buffer);

  if (compress2(out, &compressed_space, in, uncompressed_space, level) != Z_OK) {
    std::cerr << "Zlib error while compressing data.\n";
    std::abort();
  }

  return compressed_space;
#else
247
  std::cerr << "Can not call writeCompressed without compression enabled!\n";
Praetorius, Simon's avatar
Praetorius, Simon committed
248
249
250
251
252
  std::abort();
  return 0;
#endif
}

253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
inline std::uint64_t compressBuffer_lz4 (unsigned char const* /* buffer */, unsigned char* /* buffer_out */,
                                         std::uint64_t /* bs */, std::uint64_t /* cbs */, int /* level */)
{
#if HAVE_VTK_LZ4
  std::cerr << "LZ4 Compression not yet implemented" << std::endl;
  std::abort();
  return 0;
#else
  std::cerr << "LZ4 Compression not supported. Provide the LZ4 package to CMake." << std::endl;
  std::abort();
  return 0;
#endif
}

inline std::uint64_t compressBuffer_lzma (unsigned char const* /* buffer */, unsigned char* /* buffer_out */,
                                          std::uint64_t /* bs */, std::uint64_t /* cbs */, int /* level */)
{
#if HAVE_VTK_LZMA
  std::cerr << "LZMA Compression not yet implemented" << std::endl;
  std::abort();
  return 0;
#else
  std::cerr << "LZMA Compression not supported. Provide the LZMA package to CMake." << std::endl;
  std::abort();
  return 0;
#endif
}

281
282
283
} // end namespace Impl

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
Praetorius, Simon's avatar
Praetorius, Simon committed
287
288
289
290
{
  assert(is_a(format_, Vtk::APPENDED) && "Function should by called only in appended mode!\n");
  pos_type begin_pos = out.tellp();

291
  HeaderType size = values.size() * sizeof(FloatType);
Praetorius, Simon's avatar
Praetorius, Simon committed
292

293
294
295
  HeaderType num_full_blocks = size / block_size;
  HeaderType last_block_size = size % block_size;
  HeaderType num_blocks = num_full_blocks + (last_block_size > 0 ? 1 : 0);
Praetorius, Simon's avatar
Praetorius, Simon committed
296
297

  // write block-size(s)
298
299
300
301
302
303
304
  HeaderType zero = 0;
  if (compressor_ != Vtk::NONE) {
    out.write((char*)&num_blocks, sizeof(HeaderType));
    out.write((char*)&block_size, sizeof(HeaderType));
    out.write((char*)&last_block_size, sizeof(HeaderType));
    for (HeaderType i = 0; i < num_blocks; ++i)
      out.write((char*)&zero, sizeof(HeaderType));
Praetorius, Simon's avatar
Praetorius, Simon committed
305
  } else {
306
    out.write((char*)&size, sizeof(HeaderType));
Praetorius, Simon's avatar
Praetorius, Simon committed
307
308
  }

309
  HeaderType compressed_block_size = block_size + (block_size + 999)/1000 + 12;
Praetorius, Simon's avatar
Praetorius, Simon committed
310
311
  std::vector<unsigned char> buffer(block_size);
  std::vector<unsigned char> buffer_out;
312
313
314
315
  if (compressor_ != Vtk::NONE)
    buffer_out.resize(std::size_t(compressed_block_size));

  std::size_t num_values = block_size / sizeof(FloatType);
Praetorius, Simon's avatar
Praetorius, Simon committed
316

317
  std::vector<HeaderType> cbs(std::size_t(num_blocks), 0); // compressed block sizes
Praetorius, Simon's avatar
Praetorius, Simon committed
318
  for (std::size_t i = 0; i < std::size_t(num_blocks); ++i) {
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
    HeaderType bs = Impl::writeValuesToBuffer<FloatType>(num_values, buffer.data(), values, i*num_values);

    switch (compressor_) {
      case Vtk::NONE:
        out.write((char*)buffer.data(), bs);
        break;
      case Vtk::ZLIB:
        cbs[i] = Impl::compressBuffer_zlib(buffer.data(), buffer_out.data(), bs, compressed_block_size, compression_level);
        out.write((char*)buffer_out.data(), cbs[i]);
        break;
      case Vtk::LZ4:
        cbs[i] = Impl::compressBuffer_zlib(buffer.data(), buffer_out.data(), bs, compressed_block_size, compression_level);
        out.write((char*)buffer_out.data(), cbs[i]);
        break;
      case Vtk::LZMA:
        cbs[i] = Impl::compressBuffer_zlib(buffer.data(), buffer_out.data(), bs, compressed_block_size, compression_level);
        out.write((char*)buffer_out.data(), cbs[i]);
        break;
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
338
339
340
  }

  pos_type end_pos = out.tellp();
341
342
343
  if (compressor_ != Vtk::NONE) {
    out.seekp(begin_pos + std::streamoff(3*sizeof(HeaderType)));
    out.write((char*)cbs.data(), std::streamsize(num_blocks*sizeof(HeaderType)));
Praetorius, Simon's avatar
Praetorius, Simon committed
344
345
346
347
348
349
    out.seekp(end_pos);
  }

  return std::uint64_t(end_pos - begin_pos);
}

350
351
352
353
354
355
356
357
358

template <class GV, class DC>
std::string VtkWriterInterface<GV,DC>
  ::getNames (std::vector<VtkFunction> const& data) const
{
  auto scalar = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.ncomps() == 1; });
  auto vector = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.ncomps() == 3; });
  auto tensor = std::find_if(data.begin(), data.end(), [](auto const& v) { return v.ncomps() == 9; });
  return (scalar != data.end() ? " Scalars=\"" + scalar->name() + "\"" : "")
359
360
       + (vector != data.end() ? " Vectors=\"" + vector->name() + "\"" : "")
       + (tensor != data.end() ? " Tensors=\"" + tensor->name() + "\"" : "");
361
362
}

363
} // end namespace Dune