vtkreader.impl.hh 26.4 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
#include <fstream>
#include <iterator>
3
#include <sstream>
Praetorius, Simon's avatar
Praetorius, Simon committed
4
5
#include <string>

6
#if HAVE_VTK_ZLIB
Praetorius, Simon's avatar
Praetorius, Simon committed
7
8
9
#include <zlib.h>
#endif

10
#include <dune/common/classname.hh>
11
#include <dune/common/version.hh>
12

13
#include "utility/errors.hh"
Praetorius, Simon's avatar
Praetorius, Simon committed
14
15
16
#include "utility/filesystem.hh"
#include "utility/string.hh"

17
namespace Dune {
Praetorius, Simon's avatar
Praetorius, Simon committed
18

19
20
template <class Grid, class Creator, class Field>
void VtkReader<Grid,Creator,Field>::read (std::string const& filename, bool fillCreator)
Praetorius, Simon's avatar
Praetorius, Simon committed
21
22
{
  // check whether file exists!
Stenger, Florian's avatar
Stenger, Florian committed
23
  if (!Vtk::exists(filename))
Praetorius, Simon's avatar
Praetorius, Simon committed
24
25
26
    DUNE_THROW(IOError, "File " << filename << " does not exist!");

  std::ifstream input(filename, std::ios_base::in | std::ios_base::binary);
Praetorius, Simon's avatar
Praetorius, Simon committed
27
  VTK_ASSERT(input.is_open());
Praetorius, Simon's avatar
Praetorius, Simon committed
28

29
  std::string ext = Vtk::Path(filename).extension().string();
30
  if (ext == ".vtu") {
31
    readSerialFileFromStream(input, fillCreator);
32
    pieces_.push_back(filename);
33
  } else if (ext == ".pvtu") {
34
    readParallelFileFromStream(input, comm().rank(), comm().size(), fillCreator);
35
  } else {
36
    DUNE_THROW(Dune::VtkError, "File has unknown file-extension '" << ext << "'. Allowed are only '.vtu' and '.pvtu'.");
37
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
38
  read_ = true;
39
40
41
}


42
43
template <class Grid, class Creator, class Field>
void VtkReader<Grid,Creator,Field>::readSerialFileFromStream (std::ifstream& input, bool fillCreator)
44
{
Praetorius, Simon's avatar
Praetorius, Simon committed
45
  clear();
46
47
  compressor_ = Vtk::CompressorTypes::NONE;
  Vtk::DataTypes header_type = Vtk::DataTypes::UINT32;
48
49
  std::string data_id = "";
  std::string data_format = "";
50
  Vtk::DataTypes data_type = Vtk::DataTypes::UNKNOWN;
51
  unsigned int data_components = 0;
52
  std::uint64_t data_offset = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
53
54
55

  Sections section = NO_SECTION;
  for (std::string line; std::getline(input, line); ) {
Stenger, Florian's avatar
Stenger, Florian committed
56
    Vtk::ltrim(line);
Praetorius, Simon's avatar
Praetorius, Simon committed
57
58

    if (isSection(line, "VTKFile", section)) {
59
60
      bool closed = false;
      auto attr = parseXml(line, closed);
61
62

      if (!attr["type"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
63
        VTK_ASSERT_MSG(attr["type"] == "UnstructuredGrid", "VtkReader supports UnstructuredGrid types");
64
      if (!attr["byte_order"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
65
        VTK_ASSERT_MSG(attr["byte_order"] == "LittleEndian", "LittleEndian byte order supported");
66
67

      if (attr["header_type"] == "UInt32")
68
        header_type = Vtk::DataTypes::UINT32;
69
      else if (attr["header_type"] == "UInt64")
70
        header_type = Vtk::DataTypes::UINT64;
71
72

      if (attr["compressor"] == "vtkZLibDataCompressor")
73
        compressor_ = Vtk::CompressorTypes::ZLIB;
74
      else if (attr["compressor"] == "vtkLZ4DataCompressor")
75
        compressor_ = Vtk::CompressorTypes::LZ4;
76
      else if (attr["compressor"] == "vtkLZMADataCompressor")
77
        compressor_ = Vtk::CompressorTypes::LZMA;
78

Praetorius, Simon's avatar
Praetorius, Simon committed
79
80
      section = VTK_FILE;
    }
81
    else if (isSection(line, "/VTKFile", section, VTK_FILE)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
82
      section = NO_SECTION;
83
84
      break;
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
85
86
87
88
89
    else if (isSection(line, "UnstructuredGrid", section, VTK_FILE))
      section = UNSTRUCTURED_GRID;
    else if (isSection(line, "/UnstructuredGrid", section, UNSTRUCTURED_GRID))
      section = VTK_FILE;
    else if (isSection(line, "Piece", section, UNSTRUCTURED_GRID)) {
90
91
      bool closed = false;
      auto attr = parseXml(line, closed);
92

Praetorius, Simon's avatar
Praetorius, Simon committed
93
94
      VTK_ASSERT_MSG(attr.count("NumberOfPoints") > 0 && attr.count("NumberOfCells") > 0,
        "Number of points or cells in file must be > 0");
95
96
      numberOfPoints_ = std::stoul(attr["NumberOfPoints"]);
      numberOfCells_ = std::stoul(attr["NumberOfCells"]);
Praetorius, Simon's avatar
Praetorius, Simon committed
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
      section = PIECE;
    }
    else if (isSection(line, "/Piece", section, PIECE))
      section = UNSTRUCTURED_GRID;
    else if (isSection(line, "PointData", section, PIECE))
      section = POINT_DATA;
    else if (isSection(line, "/PointData", section, POINT_DATA))
      section = PIECE;
    else if (isSection(line, "CellData", section, PIECE))
      section = CELL_DATA;
    else if (isSection(line, "/CellData", section, CELL_DATA))
      section = PIECE;
    else if (isSection(line, "Points", section, PIECE))
      section = POINTS;
    else if (isSection(line, "/Points", section, POINTS))
      section = PIECE;
    else if (isSection(line, "Cells", section, PIECE))
      section = CELLS;
    else if (isSection(line, "/Cells", section, CELLS))
      section = PIECE;
    else if (line.substr(1,9) == "DataArray") {
118
119
      bool closed = false;
      auto attr = parseXml(line, closed);
Praetorius, Simon's avatar
Praetorius, Simon committed
120

121
      data_type = Vtk::dataTypeOf(attr["type"]);
122

Praetorius, Simon's avatar
Praetorius, Simon committed
123
124
125
      // Use Section.Name as id
      data_id = toString(section) + "." + attr["Name"];

126
127
      if (section == POINTS)
        // In the Points section must only be one DataArray with id=Points
Praetorius, Simon's avatar
Praetorius, Simon committed
128
        data_id = "Points";
129
130

      data_components = 1;
Praetorius, Simon's avatar
Praetorius, Simon committed
131
132
133
134
      if (!attr["NumberOfComponents"].empty())
        data_components = std::stoul(attr["NumberOfComponents"]);

      // determine FormatType
Stenger, Florian's avatar
Stenger, Florian committed
135
      data_format = Vtk::to_lower(attr["format"]);
Praetorius, Simon's avatar
Praetorius, Simon committed
136
      if (data_format == "appended") {
137
        format_ = compressor_ != Vtk::CompressorTypes::NONE ? Vtk::FormatTypes::COMPRESSED : Vtk::FormatTypes::BINARY;
Praetorius, Simon's avatar
Praetorius, Simon committed
138
      } else {
139
        format_ = Vtk::FormatTypes::ASCII;
Praetorius, Simon's avatar
Praetorius, Simon committed
140
141
      }

142
143
      // Offset makes sense in appended mode only
      data_offset = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
144
145
      if (!attr["offset"].empty()) {
        data_offset = std::stoul(attr["offset"]);
Praetorius, Simon's avatar
Praetorius, Simon committed
146
        VTK_ASSERT_MSG(data_format == "appended", "Attribute 'offset' only supported by appended mode");
Praetorius, Simon's avatar
Praetorius, Simon committed
147
148
      }

149
      // Store attributes of DataArray
Praetorius, Simon's avatar
Praetorius, Simon committed
150
      dataArray_[data_id] = {attr["Name"], data_type, data_components, data_offset, section};
151
152
153
154
155

      // Skip section in appended mode
      if (data_format == "appended") {
        if (!closed) {
          while (std::getline(input, line)) {
Stenger, Florian's avatar
Stenger, Florian committed
156
            Vtk::ltrim(line);
157
            if (line.substr(1,10) == "/DataArray")
158
159
160
161
162
163
164
              break;
          }
        }
        continue;
      }

      if (section == POINT_DATA)
Praetorius, Simon's avatar
Praetorius, Simon committed
165
166
167
168
169
170
171
172
        section = PD_DATA_ARRAY;
      else if (section == POINTS)
        section = POINTS_DATA_ARRAY;
      else if (section == CELL_DATA)
        section = CD_DATA_ARRAY;
      else if (section == CELLS)
        section = CELLS_DATA_ARRAY;
      else
173
        DUNE_THROW(Dune::VtkError, "Wrong section for <DataArray>");
Praetorius, Simon's avatar
Praetorius, Simon committed
174
175
176
177
178
179
180
181
182
183
184
    }
    else if (line.substr(1,10) == "/DataArray") {
      if (section == PD_DATA_ARRAY)
        section = POINT_DATA;
      else if (section == POINTS_DATA_ARRAY)
        section = POINTS;
      else if (section == CD_DATA_ARRAY)
        section = CELL_DATA;
      else if (section == CELLS_DATA_ARRAY)
        section = CELLS;
      else
185
        DUNE_THROW(Dune::VtkError, "Wrong section for </DataArray>");
Praetorius, Simon's avatar
Praetorius, Simon committed
186
187
    }
    else if (isSection(line, "AppendedData", section, VTK_FILE)) {
188
189
      bool closed = false;
      auto attr = parseXml(line, closed);
190
      if (!attr["encoding"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
191
        VTK_ASSERT_MSG(attr["encoding"] == "raw", "base64 encoding not supported");
192

193
      offset0_ = findAppendedDataPosition(input);
194
195
196
197
198
199
200
201
202
203
204
205
206
      Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(dataArray_["Points"].type, header_type,
      [&](auto f, auto h) {
        this->readPointsAppended(f,h,input);
        this->readCellsAppended(h,input);
      });

      // read point and cell data
      for (auto const& d : dataArray_) {
        if (d.second.section == POINT_DATA) {
          Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(d.second.type, header_type,
          [&](auto f, auto h) {
            this->readPointDataAppended(f,h,input,d.first);
          });
207
        }
208
209
210
211
212
        else if (d.second.section == CELL_DATA) {
          Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(d.second.type, header_type,
          [&](auto f, auto h) {
            this->readCellDataAppended(f,h,input,d.first);
          });
213
214
215
        }
      }

216
      section = NO_SECTION; // finish reading after appended section
217
      break;
Praetorius, Simon's avatar
Praetorius, Simon committed
218
219
220
221
222
223
    }
    else if (isSection(line, "/AppendedData", section, APPENDED_DATA))
      section = VTK_FILE;

    switch (section) {
      case PD_DATA_ARRAY:
Praetorius, Simon's avatar
Praetorius, Simon committed
224
        section = readPointData(input, data_id);
Praetorius, Simon's avatar
Praetorius, Simon committed
225
226
        break;
      case POINTS_DATA_ARRAY:
Praetorius, Simon's avatar
Praetorius, Simon committed
227
        section = readPoints(input, data_id);
Praetorius, Simon's avatar
Praetorius, Simon committed
228
229
        break;
      case CD_DATA_ARRAY:
Praetorius, Simon's avatar
Praetorius, Simon committed
230
        section = readCellData(input, data_id);
Praetorius, Simon's avatar
Praetorius, Simon committed
231
232
        break;
      case CELLS_DATA_ARRAY:
Praetorius, Simon's avatar
Praetorius, Simon committed
233
        section = readCells(input, data_id);
Praetorius, Simon's avatar
Praetorius, Simon committed
234
235
236
237
238
239
240
        break;
      default:
        break;
    }
  }

  if (section != NO_SECTION)
241
    DUNE_THROW(Dune::VtkError, "VTK-File is incomplete. It must end with </VTKFile>!");
Praetorius, Simon's avatar
Praetorius, Simon committed
242

243
244
  if (fillCreator)
    fillGridCreator();
245
246
247
}


248
template <class Grid, class Creator, class Field>
249
void VtkReader<Grid,Creator,Field>::readParallelFileFromStream (std::ifstream& input, int /* commRank */, int /* commSize */, bool fillCreator)
250
{
Praetorius, Simon's avatar
Praetorius, Simon committed
251
  clear();
252

253
254
  [[maybe_unused]] Vtk::DataTypes header_type = Vtk::DataTypes::UINT32;
  compressor_ = Vtk::CompressorTypes::NONE;
255

256
257
  Sections section = NO_SECTION;
  for (std::string line; std::getline(input, line); ) {
Stenger, Florian's avatar
Stenger, Florian committed
258
    Vtk::ltrim(line);
259
260
261
262
263
264

    if (isSection(line, "VTKFile", section)) {
      bool closed = false;
      auto attr = parseXml(line, closed);

      if (!attr["type"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
265
        VTK_ASSERT_MSG(attr["type"] == "PUnstructuredGrid", "VtkReader supports PUnstructuredGrid types");
266
      if (!attr["version"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
267
        VTK_ASSERT_MSG(std::stod(attr["version"]) == 1.0, "File format must be 1.0");
268
      if (!attr["byte_order"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
269
        VTK_ASSERT_MSG(attr["byte_order"] == "LittleEndian", "LittleEndian byte order supported");
270
271

      if (attr["header_type"] == "UInt32")
272
        header_type = Vtk::DataTypes::UINT32;
273
      else if (attr["header_type"] == "UInt64")
274
        header_type = Vtk::DataTypes::UINT64;
275
276

      if (attr["compressor"] == "vtkZLibDataCompressor")
277
        compressor_ = Vtk::CompressorTypes::ZLIB;
278
      else if (attr["compressor"] == "vtkLZ4DataCompressor")
279
        compressor_ = Vtk::CompressorTypes::LZ4;
280
      else if (attr["compressor"] == "vtkLZMADataCompressor")
281
        compressor_ = Vtk::CompressorTypes::LZMA;
282

283
284
      section = VTK_FILE;
    }
285
    else if (isSection(line, "/VTKFile", section, VTK_FILE)) {
286
      section = NO_SECTION;
287
288
      break;
    } else if (isSection(line, "PUnstructuredGrid", section, VTK_FILE))
289
290
291
292
293
294
295
      section = UNSTRUCTURED_GRID;
    else if (isSection(line, "/PUnstructuredGrid", section, UNSTRUCTURED_GRID))
      section = VTK_FILE;
    else if (isSection(line, "Piece", section, UNSTRUCTURED_GRID)) {
      bool closed = false;
      auto attr = parseXml(line, closed);

Praetorius, Simon's avatar
Praetorius, Simon committed
296
      VTK_ASSERT_MSG(attr.count("Source") > 0, "No source files for partitions provided");
297
298
299
300
      pieces_.push_back(attr["Source"]);
    }
  }

301
  VTK_ASSERT_MSG(section == NO_SECTION, "VTK-File is incomplete. It must end with </VTKFile>!");
302

303
304
  if (fillCreator)
    fillGridCreator();
Praetorius, Simon's avatar
Praetorius, Simon committed
305
306
307
}


308
309
310
311
312
313
314
// @{ implementation detail
/**
 * Read ASCII data from `input` stream into vector `values`
 * \param max_size  Upper bound for the number of values
 * \param section   Current XML section you are reading in
 * \param parent_section   XML Section to return when current `section` is finished.
 **/
Praetorius, Simon's avatar
Praetorius, Simon committed
315
template <class IStream, class T, class Sections>
316
317
Sections readDataArray (IStream& input, std::vector<T>& values, std::size_t max_size,
                        Sections section, Sections parent_section)
Praetorius, Simon's avatar
Praetorius, Simon committed
318
{
Praetorius, Simon's avatar
Praetorius, Simon committed
319
  values.reserve(max_size < std::size_t(-1) ? max_size : 0);
320
  using S = std::conditional_t<(sizeof(T) <= 1), std::uint16_t, T>; // problem when reading chars as ints
Praetorius, Simon's avatar
Praetorius, Simon committed
321
322

  std::size_t idx = 0;
323
  for (std::string line; std::getline(input, line);) {
Stenger, Florian's avatar
Stenger, Florian committed
324
    Vtk::trim(line);
325
    if (line.substr(1,10) == "/DataArray")
Praetorius, Simon's avatar
Praetorius, Simon committed
326
      return parent_section;
Praetorius, Simon's avatar
Praetorius, Simon committed
327
328
    if (line[0] == '<')
      break;
Praetorius, Simon's avatar
Praetorius, Simon committed
329
330
331
332
333

    std::istringstream stream(line);
    S value;
    for (; stream >> value; idx++)
      values.push_back(T(value));
Praetorius, Simon's avatar
Praetorius, Simon committed
334
335
336
337
338
339
340
341
342
343
344
    if (idx >= max_size)
      break;
  }

  return section;
}

template <class IStream, class Sections>
Sections skipRestOfDataArray (IStream& input, Sections section, Sections parent_section)
{
  for (std::string line; std::getline(input, line);) {
Stenger, Florian's avatar
Stenger, Florian committed
345
    Vtk::ltrim(line);
Praetorius, Simon's avatar
Praetorius, Simon committed
346
347
    if (line.substr(1,10) == "/DataArray")
      return parent_section;
Praetorius, Simon's avatar
Praetorius, Simon committed
348
349
350
351
  }

  return section;
}
352
// @}
Praetorius, Simon's avatar
Praetorius, Simon committed
353
354


355
356
357
// Read values stored on the cells with name `name`
template <class Grid, class Creator, class Field>
typename VtkReader<Grid,Creator,Field>::Sections
358
VtkReader<Grid,Creator,Field>::readCellData (std::ifstream& input, std::string id)
359
360
{
  VTK_ASSERT(numberOfCells_ > 0);
361
  unsigned int components = dataArray_[id].components;
362
363

  Sections sec;
364
  std::vector<Field>& values = cellData_[id];
365
366
367
368
369
370
371
372
373
374
375
  sec = readDataArray(input, values, components*numberOfCells_, CD_DATA_ARRAY, CELL_DATA);
  if (sec != CELL_DATA)
    sec = skipRestOfDataArray(input, CD_DATA_ARRAY, CELL_DATA);
  VTK_ASSERT(sec == CELL_DATA);
  VTK_ASSERT(values.size() == components*numberOfCells_);

  return sec;
}


template <class Grid, class Creator, class Field>
376
377
  template <class FloatType, class HeaderType>
void VtkReader<Grid,Creator,Field>::readCellDataAppended (MetaType<FloatType>, MetaType<HeaderType>, std::ifstream& input, std::string id)
378
379
{
  VTK_ASSERT(numberOfCells_ > 0);
380
  unsigned int components = dataArray_[id].components;
381

382
383
  std::vector<FloatType> values;
  readAppended(input, values, HeaderType(dataArray_[id].offset));
384
385
  VTK_ASSERT(values.size() == components*numberOfCells_);

386
387
  cellData_[id].resize(values.size());
  std::copy(values.begin(), values.end(), cellData_[id].begin());
388
389
390
391
392
}


template <class Grid, class Creator, class Field>
typename VtkReader<Grid,Creator,Field>::Sections
393
VtkReader<Grid,Creator,Field>::readPointData (std::ifstream& input, std::string id)
394
395
{
  VTK_ASSERT(numberOfPoints_ > 0);
396
  unsigned int components = dataArray_[id].components;
397
398

  Sections sec;
399
  std::vector<Field>& values = pointData_[id];
400
401
402
403
404
405
406
407
408
409
410
  sec = readDataArray(input, values, components*numberOfPoints_, PD_DATA_ARRAY, POINT_DATA);
  if (sec != POINT_DATA)
    sec = skipRestOfDataArray(input, PD_DATA_ARRAY, POINT_DATA);
  VTK_ASSERT(sec == POINT_DATA);
  VTK_ASSERT(values.size() == components*numberOfPoints_);

  return sec;
}


template <class Grid, class Creator, class Field>
411
412
  template <class FloatType, class HeaderType>
void VtkReader<Grid,Creator,Field>::readPointDataAppended (MetaType<FloatType>, MetaType<HeaderType>, std::ifstream& input, std::string id)
413
414
{
  VTK_ASSERT(numberOfPoints_ > 0);
415
  unsigned int components = dataArray_[id].components;
416

417
418
  std::vector<FloatType> values;
  readAppended(input, values, HeaderType(dataArray_[id].offset));
419
420
  VTK_ASSERT(values.size() == components*numberOfPoints_);

421
422
  pointData_[id].resize(values.size());
  std::copy(values.begin(), values.end(), pointData_[id].begin());
423
424
425
426
427
}


template <class Grid, class Creator, class Field>
typename VtkReader<Grid,Creator,Field>::Sections
428
VtkReader<Grid,Creator,Field>::readPoints (std::ifstream& input, std::string id)
Praetorius, Simon's avatar
Praetorius, Simon committed
429
430
{
  using T = typename GlobalCoordinate::value_type;
Praetorius, Simon's avatar
Praetorius, Simon committed
431
  VTK_ASSERT(numberOfPoints_ > 0);
432
  VTK_ASSERT(id == "Points");
433
  VTK_ASSERT(dataArray_["Points"].components == 3u);
434

435
436
  Sections sec;

Praetorius, Simon's avatar
Praetorius, Simon committed
437
  std::vector<T> point_values;
438
  sec = readDataArray(input, point_values, 3*numberOfPoints_, POINTS_DATA_ARRAY, POINTS);
Praetorius, Simon's avatar
Praetorius, Simon committed
439
440
  if (sec != POINTS)
    sec = skipRestOfDataArray(input, POINTS_DATA_ARRAY, POINTS);
Praetorius, Simon's avatar
Praetorius, Simon committed
441
442
  VTK_ASSERT(sec == POINTS);
  VTK_ASSERT(point_values.size() == 3*numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
443
444
445

  // extract points from continuous values
  GlobalCoordinate p;
446
  vec_points.reserve(numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
447
  std::size_t idx = 0;
448
  for (std::size_t i = 0; i < numberOfPoints_; ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
449
450
451
    for (std::size_t j = 0; j < p.size(); ++j)
      p[j] = point_values[idx++];
    idx += (3u - p.size());
452
    vec_points.push_back(p);
Praetorius, Simon's avatar
Praetorius, Simon committed
453
454
455
456
457
458
  }

  return sec;
}


459
template <class Grid, class Creator, class Field>
460
461
  template <class FloatType, class HeaderType>
void VtkReader<Grid,Creator,Field>::readPointsAppended (MetaType<FloatType>, MetaType<HeaderType>, std::ifstream& input)
Praetorius, Simon's avatar
Praetorius, Simon committed
462
{
Praetorius, Simon's avatar
Praetorius, Simon committed
463
  VTK_ASSERT(numberOfPoints_ > 0);
464
  VTK_ASSERT(dataArray_["Points"].components == 3u);
465
466
  std::vector<FloatType> point_values;
  readAppended(input, point_values, HeaderType(dataArray_["Points"].offset));
Praetorius, Simon's avatar
Praetorius, Simon committed
467
  VTK_ASSERT(point_values.size() == 3*numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
468
469
470

  // extract points from continuous values
  GlobalCoordinate p;
471
  vec_points.reserve(numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
472
  std::size_t idx = 0;
473
  for (std::size_t i = 0; i < numberOfPoints_; ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
474
    for (std::size_t j = 0; j < p.size(); ++j)
475
      p[j] = FloatType(point_values[idx++]);
Praetorius, Simon's avatar
Praetorius, Simon committed
476
    idx += (3u - p.size());
477
    vec_points.push_back(p);
Praetorius, Simon's avatar
Praetorius, Simon committed
478
479
480
481
  }
}


482
483
template <class Grid, class Creator, class Field>
typename VtkReader<Grid,Creator,Field>::Sections
484
VtkReader<Grid,Creator,Field>::readCells (std::ifstream& input, std::string id)
Praetorius, Simon's avatar
Praetorius, Simon committed
485
486
487
{
  Sections sec = CELLS_DATA_ARRAY;

Praetorius, Simon's avatar
Praetorius, Simon committed
488
  VTK_ASSERT(numberOfCells_ > 0);
489
  if (id == "Cells.types") {
490
    sec = readDataArray(input, vec_types, numberOfCells_, CELLS_DATA_ARRAY, CELLS);
Praetorius, Simon's avatar
Praetorius, Simon committed
491
    VTK_ASSERT(vec_types.size() == numberOfCells_);
492
  } else if (id == "Cells.offsets") {
493
    sec = readDataArray(input, vec_offsets, numberOfCells_, CELLS_DATA_ARRAY, CELLS);
Praetorius, Simon's avatar
Praetorius, Simon committed
494
    VTK_ASSERT(vec_offsets.size() == numberOfCells_);
495
  } else if (id == "Cells.connectivity") {
Praetorius, Simon's avatar
Praetorius, Simon committed
496
    sec = readDataArray(input, vec_connectivity, std::size_t(-1), CELLS_DATA_ARRAY, CELLS);
497
  } else if (id == "Cells.global_point_ids") {
498
    sec = readDataArray(input, vec_point_ids, numberOfPoints_, CELLS_DATA_ARRAY, CELLS);
Praetorius, Simon's avatar
Praetorius, Simon committed
499
    VTK_ASSERT(vec_point_ids.size() == numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
500
501
502
503
504
505
  }

  return sec;
}


506
template <class Grid, class Creator, class Field>
507
508
  template <class HeaderType>
void VtkReader<Grid,Creator,Field>::readCellsAppended (MetaType<HeaderType>, std::ifstream& input)
Praetorius, Simon's avatar
Praetorius, Simon committed
509
{
Praetorius, Simon's avatar
Praetorius, Simon committed
510
  VTK_ASSERT(numberOfCells_ > 0);
511
512
513
  auto types_data = dataArray_["Cells.types"];
  auto offsets_data = dataArray_["Cells.offsets"];
  auto connectivity_data = dataArray_["Cells.connectivity"];
Praetorius, Simon's avatar
Praetorius, Simon committed
514

515
  VTK_ASSERT(types_data.type == Vtk::DataTypes::UINT8);
516
  readAppended(input, vec_types, HeaderType(types_data.offset));
Praetorius, Simon's avatar
Praetorius, Simon committed
517
  VTK_ASSERT(vec_types.size() == numberOfCells_);
Praetorius, Simon's avatar
Praetorius, Simon committed
518

Praetorius, Simon's avatar
Praetorius, Simon committed
519
  if (offsets_data.type == Vtk::DataTypes::INT64)
520
    readAppended(input, vec_offsets, HeaderType(offsets_data.offset));
Praetorius, Simon's avatar
Praetorius, Simon committed
521
  else if (offsets_data.type == Vtk::DataTypes::INT32) {
522
    std::vector<std::int32_t> offsets;
523
    readAppended(input, offsets, HeaderType(offsets_data.offset));
524
525
526
527
    vec_offsets.resize(offsets.size());
    std::copy(offsets.begin(), offsets.end(), vec_offsets.begin());
  }
  else { DUNE_THROW(Dune::NotImplemented, "Unsupported DataType in Cell offsets."); }
Praetorius, Simon's avatar
Praetorius, Simon committed
528
  VTK_ASSERT(vec_offsets.size() == numberOfCells_);
Praetorius, Simon's avatar
Praetorius, Simon committed
529

Praetorius, Simon's avatar
Praetorius, Simon committed
530
  if (connectivity_data.type == Vtk::DataTypes::INT64)
531
    readAppended(input, vec_connectivity, HeaderType(connectivity_data.offset));
Praetorius, Simon's avatar
Praetorius, Simon committed
532
  else if (connectivity_data.type == Vtk::DataTypes::INT32) {
533
    std::vector<std::int32_t> connectivity;
534
    readAppended(input, connectivity, HeaderType(connectivity_data.offset));
535
536
537
538
    vec_connectivity.resize(connectivity.size());
    std::copy(connectivity.begin(), connectivity.end(), vec_connectivity.begin());
  }
  else { DUNE_THROW(Dune::NotImplemented, "Unsupported DataType in Cell connectivity."); }
Praetorius, Simon's avatar
Praetorius, Simon committed
539
  VTK_ASSERT(vec_connectivity.size() == std::size_t(vec_offsets.back()));
540

541
542
  if (dataArray_.count("Cells.global_point_ids") > 0) {
    auto point_id_data = dataArray_["Cells.global_point_ids"];
543
    VTK_ASSERT(point_id_data.type == Vtk::DataTypes::UINT64);
544
    readAppended(input, vec_point_ids, HeaderType(point_id_data.offset));
Praetorius, Simon's avatar
Praetorius, Simon committed
545
    VTK_ASSERT(vec_point_ids.size() == numberOfPoints_);
546
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
547
548
549
}


550
// @{ implementation detail
551
552
namespace {

553
554
555
556
557
558
559
/**
 * Read compressed data into `buffer_in`, uncompress it and store the result in
 * the concrete-data-type `buffer`
 * \param bs     Size of the uncompressed data
 * \param cbs    Size of the compressed data
 * \param input  Stream to read from.
 **/
Praetorius, Simon's avatar
Praetorius, Simon committed
560
template <class T, class IStream>
561
562
void read_compressed_zlib (T* buffer, unsigned char* buffer_in,
                           std::uint64_t bs, std::uint64_t cbs, IStream& input)
Praetorius, Simon's avatar
Praetorius, Simon committed
563
{
564
#if HAVE_VTK_ZLIB
Praetorius, Simon's avatar
Praetorius, Simon committed
565
566
567
568
569
570
  uLongf uncompressed_space = uLongf(bs);
  uLongf compressed_space = uLongf(cbs);

  Bytef* compressed_buffer = reinterpret_cast<Bytef*>(buffer_in);
  Bytef* uncompressed_buffer = reinterpret_cast<Bytef*>(buffer);

571
  input.read((char*)(compressed_buffer), compressed_space);
Praetorius, Simon's avatar
Praetorius, Simon committed
572
  VTK_ASSERT(uLongf(input.gcount()) == compressed_space);
Praetorius, Simon's avatar
Praetorius, Simon committed
573
574
575
576
577

  if (uncompress(uncompressed_buffer, &uncompressed_space, compressed_buffer, compressed_space) != Z_OK) {
    std::cerr << "Zlib error while uncompressing data.\n";
    std::abort();
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
578
  VTK_ASSERT(uLongf(bs) == uncompressed_space);
Praetorius, Simon's avatar
Praetorius, Simon committed
579
#else
580
  std::cerr << "ZLib Compression not supported. Provide the ZLIB package to CMake." << std::endl;
Praetorius, Simon's avatar
Praetorius, Simon committed
581
582
  std::abort();
#endif
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
}

template <class T, class IStream>
void read_compressed_lz4 (T* /* buffer */, unsigned char* /* buffer_in */,
                          std::uint64_t /* bs */, std::uint64_t /* cbs */, IStream& /* input */)
{
#if HAVE_VTK_LZ4
  std::cerr << "LZ4 Compression not yet implemented" << std::endl;
  std::abort();
#else
  std::cerr << "LZ4 Compression not supported. Provide the LZ4 package to CMake." << std::endl;
  std::abort();
#endif
}

template <class T, class IStream>
void read_compressed_lzma (T* /* buffer */, unsigned char* /* buffer_in */,
                           std::uint64_t /* bs */, std::uint64_t /* cbs */, IStream& /* input */)
{
#if HAVE_VTK_LZMA
  std::cerr << "LZMA Compression not yet implemented" << std::endl;
  std::abort();
#else
  std::cerr << "LZMA Compression not supported. Provide the LZMA package to CMake." << std::endl;
  std::abort();
#endif
}

Praetorius, Simon's avatar
Praetorius, Simon committed
611
}
612
// @}
Praetorius, Simon's avatar
Praetorius, Simon committed
613
614


615
template <class Grid, class Creator, class Field>
616
617
  template <class FloatType, class HeaderType>
void VtkReader<Grid,Creator,Field>::readAppended (std::ifstream& input, std::vector<FloatType>& values, HeaderType offset)
Praetorius, Simon's avatar
Praetorius, Simon committed
618
619
620
{
  input.seekg(offset0_ + offset);

621
  HeaderType size = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
622

623
624
625
626
  HeaderType num_blocks = 0;
  HeaderType block_size = 0;
  HeaderType last_block_size = 0;
  std::vector<HeaderType> cbs; // compressed block sizes
Praetorius, Simon's avatar
Praetorius, Simon committed
627
628

  // read total size / block-size(s)
629
  if (compressor_ != Vtk::CompressorTypes::NONE) {
630
631
632
    input.read((char*)&num_blocks, sizeof(HeaderType));
    input.read((char*)&block_size, sizeof(HeaderType));
    input.read((char*)&last_block_size, sizeof(HeaderType));
Praetorius, Simon's avatar
Praetorius, Simon committed
633

634
    VTK_ASSERT(block_size % sizeof(FloatType) == 0);
Praetorius, Simon's avatar
Praetorius, Simon committed
635
636
637
638
639
640

    // total size of the uncompressed data
    size = block_size * (num_blocks-1) + last_block_size;

    // size of the compressed blocks
    cbs.resize(num_blocks);
641
    input.read((char*)cbs.data(), num_blocks*sizeof(HeaderType));
Praetorius, Simon's avatar
Praetorius, Simon committed
642
  } else {
643
    input.read((char*)&size, sizeof(HeaderType));
Praetorius, Simon's avatar
Praetorius, Simon committed
644
  }
645
646
  VTK_ASSERT(size > 0 && (size % sizeof(FloatType)) == 0);
  values.resize(size / sizeof(FloatType));
Praetorius, Simon's avatar
Praetorius, Simon committed
647

648
  if (compressor_ != Vtk::CompressorTypes::NONE) {
Praetorius, Simon's avatar
Praetorius, Simon committed
649
    // upper bound for compressed block-size
650
    HeaderType compressed_block_size = block_size + (block_size + 999)/1000 + 12;
Praetorius, Simon's avatar
Praetorius, Simon committed
651
    // number of values in the full blocks
652
    std::size_t num_values = block_size / sizeof(FloatType);
Praetorius, Simon's avatar
Praetorius, Simon committed
653
654
655

    std::vector<unsigned char> buffer_in(compressed_block_size);
    for (std::size_t i = 0; i < std::size_t(num_blocks); ++i) {
656
657
658
      HeaderType bs = i < std::size_t(num_blocks-1) ? block_size : last_block_size;

      switch (compressor_) {
659
        case Vtk::CompressorTypes::ZLIB:
660
661
          read_compressed_zlib(values.data() + i*num_values, buffer_in.data(), bs, cbs[i], input);
          break;
662
        case Vtk::CompressorTypes::LZ4:
663
664
          read_compressed_lz4(values.data() + i*num_values, buffer_in.data(), bs, cbs[i], input);
          break;
665
        case Vtk::CompressorTypes::LZMA:
666
667
668
669
670
671
          read_compressed_lzma(values.data() + i*num_values, buffer_in.data(), bs, cbs[i], input);
          break;
        default:
          VTK_ASSERT_MSG(false, "Unsupported Compressor type.");
          break;
      }
Praetorius, Simon's avatar
Praetorius, Simon committed
672
673
    }
  } else {
674
    input.read((char*)(values.data()), size);
Praetorius, Simon's avatar
Praetorius, Simon committed
675
    VTK_ASSERT(input.gcount() == std::streamsize(size));
Praetorius, Simon's avatar
Praetorius, Simon committed
676
677
678
679
  }
}


680
681
template <class Grid, class Creator, class Field>
void VtkReader<Grid,Creator,Field>::fillGridCreator (bool insertPieces)
Praetorius, Simon's avatar
Praetorius, Simon committed
682
{
Praetorius, Simon's avatar
Praetorius, Simon committed
683
684
685
  VTK_ASSERT(vec_points.size() == numberOfPoints_);
  VTK_ASSERT(vec_types.size() == numberOfCells_);
  VTK_ASSERT(vec_offsets.size() == numberOfCells_);
Praetorius, Simon's avatar
Praetorius, Simon committed
686

Praetorius, Simon's avatar
Praetorius, Simon committed
687
  if (!vec_points.empty())
688
    creator_->insertVertices(vec_points, vec_point_ids);
Praetorius, Simon's avatar
Praetorius, Simon committed
689
  if (!vec_types.empty())
690
    creator_->insertElements(vec_types, vec_offsets, vec_connectivity);
691
  if (insertPieces)
692
    creator_->insertPieces(pieces_);
693
}
Praetorius, Simon's avatar
Praetorius, Simon committed
694

695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727

// Convert section into string
template <class Grid, class Creator, class Field>
std::string VtkReader<Grid,Creator,Field>::toString(Sections s) const
{
  switch (s) {
    case VTK_FILE:
      return "VTKFile";
    case UNSTRUCTURED_GRID:
      return "UnstructuredGrid";
    case PIECE:
      return "Piece";
    case POINT_DATA:
      return "PointData";
    case CELL_DATA:
      return "CellData";
    case POINTS:
      return "Points";
    case CELLS:
      return "Cells";
    case APPENDED_DATA:
      return "AppendedData";
    case PD_DATA_ARRAY:
    case CD_DATA_ARRAY:
    case POINTS_DATA_ARRAY:
    case CELLS_DATA_ARRAY:
      return "DataArray";
    default:
      return "Unknown";
  }
}


728
// Assume input already read the line <AppendedData ...>
729
730
template <class Grid, class Creator, class Field>
std::uint64_t VtkReader<Grid,Creator,Field>::findAppendedDataPosition (std::ifstream& input) const
731
732
733
734
735
736
737
738
739
740
741
742
{
  char c;
  while (input.get(c) && std::isblank(c)) { /*do nothing*/ }

  std::uint64_t offset = input.tellg();
  if (c != '_')
    --offset; // if char is not '_', assume it is part of the data.

  return offset;
}


743
744
template <class Grid, class Creator, class Field>
std::map<std::string, std::string> VtkReader<Grid,Creator,Field>::parseXml (std::string const& line, bool& closed)
Praetorius, Simon's avatar
Praetorius, Simon committed
745
{
746
  closed = false;
Praetorius, Simon's avatar
Praetorius, Simon committed
747
748
749
750
751
752
753
754
755
756
757
758
759
760
  std::map<std::string, std::string> attr;

  Sections sec = NO_SECTION;
  bool escape = false;

  std::string name = "";
  std::string value = "";
  for (auto c : line) {
    switch (sec) {
    case NO_SECTION:
      if (std::isalpha(c) || c == '_') {
        name.clear();
        sec = XML_NAME;
        name.push_back(c);
761
762
      } else if (c == '/') {
        closed = true;
Praetorius, Simon's avatar
Praetorius, Simon committed
763
764
765
766
767
768
769
770
771
772
773
      }
      break;
    case XML_NAME:
      if (std::isalpha(c) || c == '_')
        name.push_back(c);
      else
        sec = (c == '=' ? XML_NAME_ASSIGN : NO_SECTION);
      break;
    case XML_NAME_ASSIGN:
      value.clear();
      escape = false;
774
      VTK_ASSERT_MSG( c == '"', "Format error!" );
Praetorius, Simon's avatar
Praetorius, Simon committed
775
776
777
778
779
780
781
782
783
784
785
786
787
788
      sec = XML_VALUE;
      break;
    case XML_VALUE:
      if (c == '"' && !escape) {
        attr[name] = value;
        sec = NO_SECTION;
      } else if (c == '\\' && !escape) {
        escape = true;
      }  else {
        value.push_back(c);
        escape = false;
      }
      break;
    default:
789
      VTK_ASSERT_MSG(false, "Format error!");
Praetorius, Simon's avatar
Praetorius, Simon committed
790
791
792
793
794
795
    }
  }

  return attr;
}

Praetorius, Simon's avatar
Praetorius, Simon committed
796

797
798
template <class Grid, class Creator, class Field>
void VtkReader<Grid,Creator,Field>::clear ()
Praetorius, Simon's avatar
Praetorius, Simon committed
799
800
801
802
803
804
805
806
807
808
809
810
{
  vec_points.clear();
  vec_point_ids.clear();
  vec_types.clear();
  vec_offsets.clear();
  vec_connectivity.clear();
  dataArray_.clear();
  pieces_.clear();

  numberOfCells_ = 0;
  numberOfPoints_ = 0;
  offset0_ = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
811
  read_ = false;
Praetorius, Simon's avatar
Praetorius, Simon committed
812
813
}

814
} // end namespace Dune