vtkreader.impl.hh 23.8 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
38
39
40
  }
}


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

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

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

      if (!attr["type"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
60
        VTK_ASSERT_MSG(attr["type"] == "UnstructuredGrid", "VtkReader supports UnstructuredGrid types");
Praetorius, Simon's avatar
Praetorius, Simon committed
61
      if (!attr["version"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
62
        VTK_ASSERT_MSG(std::stod(attr["version"]) == 1.0, "File format must be 1.0");
63
      if (!attr["byte_order"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
64
        VTK_ASSERT_MSG(attr["byte_order"] == "LittleEndian", "LittleEndian byte order supported");
65
      if (!attr["header_type"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
66
        VTK_ASSERT_MSG(attr["header_type"] == "UInt64", "Header integer type must be UInt64");
67
      if (!attr["compressor"].empty()) {
Praetorius, Simon's avatar
Praetorius, Simon committed
68
        compressor = attr["compressor"];
Praetorius, Simon's avatar
Praetorius, Simon committed
69
        VTK_ASSERT_MSG(compressor == "vtkZLibDataCompressor", "Only ZLib compression supported");
70
      }
Praetorius, Simon's avatar
Praetorius, Simon committed
71
72
73
74
75
76
77
78
79
      section = VTK_FILE;
    }
    else if (isSection(line, "/VTKFile", section, VTK_FILE))
      section = NO_SECTION;
    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)) {
80
81
      bool closed = false;
      auto attr = parseXml(line, closed);
82

Praetorius, Simon's avatar
Praetorius, Simon committed
83
84
      VTK_ASSERT_MSG(attr.count("NumberOfPoints") > 0 && attr.count("NumberOfCells") > 0,
        "Number of points or cells in file must be > 0");
85
86
      numberOfPoints_ = std::stoul(attr["NumberOfPoints"]);
      numberOfCells_ = std::stoul(attr["NumberOfCells"]);
Praetorius, Simon's avatar
Praetorius, Simon committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
      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") {
108
109
      bool closed = false;
      auto attr = parseXml(line, closed);
Praetorius, Simon's avatar
Praetorius, Simon committed
110

111
      data_type = Vtk::Map::to_datatype[attr["type"]];
112

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

116
117
      if (section == POINTS)
        // In the Points section must only be one DataArray with id=Points
Praetorius, Simon's avatar
Praetorius, Simon committed
118
        data_id = "Points";
119
120

      data_components = 1;
Praetorius, Simon's avatar
Praetorius, Simon committed
121
122
123
124
      if (!attr["NumberOfComponents"].empty())
        data_components = std::stoul(attr["NumberOfComponents"]);

      // determine FormatType
Stenger, Florian's avatar
Stenger, Florian committed
125
      data_format = Vtk::to_lower(attr["format"]);
Praetorius, Simon's avatar
Praetorius, Simon committed
126
      if (data_format == "appended") {
127
        format_ = !compressor.empty() ? Vtk::COMPRESSED : Vtk::BINARY;
Praetorius, Simon's avatar
Praetorius, Simon committed
128
129
130
131
      } else {
        format_ = Vtk::ASCII;
      }

132
133
      // Offset makes sense in appended mode only
      data_offset = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
134
135
      if (!attr["offset"].empty()) {
        data_offset = std::stoul(attr["offset"]);
Praetorius, Simon's avatar
Praetorius, Simon committed
136
        VTK_ASSERT_MSG(data_format == "appended", "Attribute 'offset' only supported by appended mode");
Praetorius, Simon's avatar
Praetorius, Simon committed
137
138
      }

139
      // Store attributes of DataArray
Praetorius, Simon's avatar
Praetorius, Simon committed
140
      dataArray_[data_id] = {attr["Name"], data_type, data_components, data_offset, section};
141
142
143
144
145

      // Skip section in appended mode
      if (data_format == "appended") {
        if (!closed) {
          while (std::getline(input, line)) {
Stenger, Florian's avatar
Stenger, Florian committed
146
            Vtk::ltrim(line);
147
            if (line.substr(1,10) == "/DataArray")
148
149
150
151
152
153
154
              break;
          }
        }
        continue;
      }

      if (section == POINT_DATA)
Praetorius, Simon's avatar
Praetorius, Simon committed
155
156
157
158
159
160
161
162
        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
163
        DUNE_THROW(Dune::VtkError, "Wrong section for <DataArray>");
Praetorius, Simon's avatar
Praetorius, Simon committed
164
165
166
167
168
169
170
171
172
173
174
    }
    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
175
        DUNE_THROW(Dune::VtkError, "Wrong section for </DataArray>");
Praetorius, Simon's avatar
Praetorius, Simon committed
176
177
    }
    else if (isSection(line, "AppendedData", section, VTK_FILE)) {
178
179
      bool closed = false;
      auto attr = parseXml(line, closed);
180
      if (!attr["encoding"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
181
        VTK_ASSERT_MSG(attr["encoding"] == "raw", "base64 encoding not supported");
182

183
      offset0_ = findAppendedDataPosition(input);
184
      if (dataArray_["Points"].type == Vtk::FLOAT32)
185
186
187
188
189
        readPointsAppended<float>(input);
      else
        readPointsAppended<double>(input);

      readCellsAppended(input);
190
191

      // read point an cell data
192
      for (auto const& [id,data] : dataArray_) {
193
194
        if (data.section == POINT_DATA) {
          if (data.type == Vtk::FLOAT32)
195
            readPointDataAppended<float>(input, id);
196
          else
197
            readPointDataAppended<double>(input, id);
198
199
200
        }
        else if (data.section == CELL_DATA) {
          if (data.type == Vtk::FLOAT32)
201
            readCellDataAppended<float>(input, id);
202
          else
203
            readCellDataAppended<double>(input, id);
204
205
206
        }
      }

207
      section = NO_SECTION; // finish reading after appended section
Praetorius, Simon's avatar
Praetorius, Simon committed
208
209
210
211
212
213
    }
    else if (isSection(line, "/AppendedData", section, APPENDED_DATA))
      section = VTK_FILE;

    switch (section) {
      case PD_DATA_ARRAY:
Praetorius, Simon's avatar
Praetorius, Simon committed
214
        section = readPointData(input, data_id);
Praetorius, Simon's avatar
Praetorius, Simon committed
215
216
        break;
      case POINTS_DATA_ARRAY:
Praetorius, Simon's avatar
Praetorius, Simon committed
217
        section = readPoints(input, data_id);
Praetorius, Simon's avatar
Praetorius, Simon committed
218
219
        break;
      case CD_DATA_ARRAY:
Praetorius, Simon's avatar
Praetorius, Simon committed
220
        section = readCellData(input, data_id);
Praetorius, Simon's avatar
Praetorius, Simon committed
221
222
        break;
      case CELLS_DATA_ARRAY:
Praetorius, Simon's avatar
Praetorius, Simon committed
223
        section = readCells(input, data_id);
Praetorius, Simon's avatar
Praetorius, Simon committed
224
225
226
227
228
229
230
231
232
233
234
        break;
      default:
        // do nothing
        break;
    }

    if (section == NO_SECTION)
      break;
  }

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

237
238
  if (fillCreator)
    fillGridCreator();
239
240
241
}


242
243
template <class Grid, class Creator, class Field>
void VtkReader<Grid,Creator,Field>::readParallelFileFromStream (std::ifstream& input, int commRank, int commSize, bool fillCreator)
244
{
Praetorius, Simon's avatar
Praetorius, Simon committed
245
  clear();
246
247
248

  Sections section = NO_SECTION;
  for (std::string line; std::getline(input, line); ) {
Stenger, Florian's avatar
Stenger, Florian committed
249
    Vtk::ltrim(line);
250
251
252
253
254
255

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

      if (!attr["type"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
256
        VTK_ASSERT_MSG(attr["type"] == "PUnstructuredGrid", "VtkReader supports PUnstructuredGrid types");
257
      if (!attr["version"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
258
        VTK_ASSERT_MSG(std::stod(attr["version"]) == 1.0, "File format must be 1.0");
259
      if (!attr["byte_order"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
260
        VTK_ASSERT_MSG(attr["byte_order"] == "LittleEndian", "LittleEndian byte order supported");
261
      if (!attr["header_type"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
262
        VTK_ASSERT_MSG(attr["header_type"] == "UInt64", "Header integer type must be UInt64");
263
      if (!attr["compressor"].empty())
Praetorius, Simon's avatar
Praetorius, Simon committed
264
        VTK_ASSERT_MSG(attr["compressor"] == "vtkZLibDataCompressor", "Only ZLib compression supported");
265
266
267
268
269
270
271
272
273
274
275
276
      section = VTK_FILE;
    }
    else if (isSection(line, "/VTKFile", section, VTK_FILE))
      section = NO_SECTION;
    else if (isSection(line, "PUnstructuredGrid", section, VTK_FILE))
      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
277
      VTK_ASSERT_MSG(attr.count("Source") > 0, "No source files for partitions provided");
278
279
280
281
282
283
284
      pieces_.push_back(attr["Source"]);
    }

    if (section == NO_SECTION)
      break;
  }

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

287
288
  if (fillCreator)
    fillGridCreator();
Praetorius, Simon's avatar
Praetorius, Simon committed
289
290
291
}


292
293
294
295
296
297
298
// @{ 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
299
template <class IStream, class T, class Sections>
300
301
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
302
{
Praetorius, Simon's avatar
Praetorius, Simon committed
303
  values.reserve(max_size < std::size_t(-1) ? max_size : 0);
304
  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
305
306

  std::size_t idx = 0;
307
  for (std::string line; std::getline(input, line);) {
Stenger, Florian's avatar
Stenger, Florian committed
308
    Vtk::trim(line);
309
    if (line.substr(1,10) == "/DataArray")
Praetorius, Simon's avatar
Praetorius, Simon committed
310
      return parent_section;
Praetorius, Simon's avatar
Praetorius, Simon committed
311
312
    if (line[0] == '<')
      break;
Praetorius, Simon's avatar
Praetorius, Simon committed
313
314
315
316
317

    std::istringstream stream(line);
    S value;
    for (; stream >> value; idx++)
      values.push_back(T(value));
Praetorius, Simon's avatar
Praetorius, Simon committed
318
319
320
321
322
323
324
325
326
327
328
    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
329
    Vtk::ltrim(line);
Praetorius, Simon's avatar
Praetorius, Simon committed
330
331
    if (line.substr(1,10) == "/DataArray")
      return parent_section;
Praetorius, Simon's avatar
Praetorius, Simon committed
332
333
334
335
  }

  return section;
}
336
// @}
Praetorius, Simon's avatar
Praetorius, Simon committed
337
338


339
340
341
// Read values stored on the cells with name `name`
template <class Grid, class Creator, class Field>
typename VtkReader<Grid,Creator,Field>::Sections
342
VtkReader<Grid,Creator,Field>::readCellData (std::ifstream& input, std::string id)
343
344
{
  VTK_ASSERT(numberOfCells_ > 0);
345
  unsigned int components = dataArray_[id].components;
346
347

  Sections sec;
348
  std::vector<Field>& values = cellData_[id];
349
350
351
352
353
354
355
356
357
358
359
360
  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>
  template <class T>
361
void VtkReader<Grid,Creator,Field>::readCellDataAppended (std::ifstream& input, std::string id)
362
363
{
  VTK_ASSERT(numberOfCells_ > 0);
364
  unsigned int components = dataArray_[id].components;
365
366

  std::vector<T> values;
367
  readAppended(input, values, dataArray_[id].offset);
368
369
  VTK_ASSERT(values.size() == components*numberOfCells_);

370
371
  cellData_[id].resize(values.size());
  std::copy(values.begin(), values.end(), cellData_[id].begin());
372
373
374
375
376
}


template <class Grid, class Creator, class Field>
typename VtkReader<Grid,Creator,Field>::Sections
377
VtkReader<Grid,Creator,Field>::readPointData (std::ifstream& input, std::string id)
378
379
{
  VTK_ASSERT(numberOfPoints_ > 0);
380
  unsigned int components = dataArray_[id].components;
381
382

  Sections sec;
383
  std::vector<Field>& values = pointData_[id];
384
385
386
387
388
389
390
391
392
393
394
395
  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>
  template <class T>
396
void VtkReader<Grid,Creator,Field>::readPointDataAppended (std::ifstream& input, std::string id)
397
398
{
  VTK_ASSERT(numberOfPoints_ > 0);
399
  unsigned int components = dataArray_[id].components;
400
401

  std::vector<T> values;
402
  readAppended(input, values, dataArray_[id].offset);
403
404
  VTK_ASSERT(values.size() == components*numberOfPoints_);

405
406
  pointData_[id].resize(values.size());
  std::copy(values.begin(), values.end(), pointData_[id].begin());
407
408
409
410
411
}


template <class Grid, class Creator, class Field>
typename VtkReader<Grid,Creator,Field>::Sections
412
VtkReader<Grid,Creator,Field>::readPoints (std::ifstream& input, std::string id)
Praetorius, Simon's avatar
Praetorius, Simon committed
413
414
{
  using T = typename GlobalCoordinate::value_type;
Praetorius, Simon's avatar
Praetorius, Simon committed
415
  VTK_ASSERT(numberOfPoints_ > 0);
416
  VTK_ASSERT(dataArray_["Points"].components == 3u);
417

418
419
  Sections sec;

Praetorius, Simon's avatar
Praetorius, Simon committed
420
  std::vector<T> point_values;
421
  sec = readDataArray(input, point_values, 3*numberOfPoints_, POINTS_DATA_ARRAY, POINTS);
Praetorius, Simon's avatar
Praetorius, Simon committed
422
423
  if (sec != POINTS)
    sec = skipRestOfDataArray(input, POINTS_DATA_ARRAY, POINTS);
Praetorius, Simon's avatar
Praetorius, Simon committed
424
425
  VTK_ASSERT(sec == POINTS);
  VTK_ASSERT(point_values.size() == 3*numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
426
427
428

  // extract points from continuous values
  GlobalCoordinate p;
429
  vec_points.reserve(numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
430
  std::size_t idx = 0;
431
  for (std::size_t i = 0; i < numberOfPoints_; ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
432
433
434
    for (std::size_t j = 0; j < p.size(); ++j)
      p[j] = point_values[idx++];
    idx += (3u - p.size());
435
    vec_points.push_back(p);
Praetorius, Simon's avatar
Praetorius, Simon committed
436
437
438
439
440
441
  }

  return sec;
}


442
template <class Grid, class Creator, class Field>
Praetorius, Simon's avatar
Praetorius, Simon committed
443
  template <class T>
444
void VtkReader<Grid,Creator,Field>::readPointsAppended (std::ifstream& input)
Praetorius, Simon's avatar
Praetorius, Simon committed
445
{
Praetorius, Simon's avatar
Praetorius, Simon committed
446
  VTK_ASSERT(numberOfPoints_ > 0);
447
  VTK_ASSERT(dataArray_["Points"].components == 3u);
Praetorius, Simon's avatar
Praetorius, Simon committed
448
  std::vector<T> point_values;
449
  readAppended(input, point_values, dataArray_["Points"].offset);
Praetorius, Simon's avatar
Praetorius, Simon committed
450
  VTK_ASSERT(point_values.size() == 3*numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
451
452
453

  // extract points from continuous values
  GlobalCoordinate p;
454
  vec_points.reserve(numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
455
  std::size_t idx = 0;
456
  for (std::size_t i = 0; i < numberOfPoints_; ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
457
458
459
    for (std::size_t j = 0; j < p.size(); ++j)
      p[j] = T(point_values[idx++]);
    idx += (3u - p.size());
460
    vec_points.push_back(p);
Praetorius, Simon's avatar
Praetorius, Simon committed
461
462
463
464
  }
}


465
466
template <class Grid, class Creator, class Field>
typename VtkReader<Grid,Creator,Field>::Sections
467
VtkReader<Grid,Creator,Field>::readCells (std::ifstream& input, std::string id)
Praetorius, Simon's avatar
Praetorius, Simon committed
468
469
470
{
  Sections sec = CELLS_DATA_ARRAY;

Praetorius, Simon's avatar
Praetorius, Simon committed
471
  VTK_ASSERT(numberOfCells_ > 0);
472
  if (id == "Cells.types") {
473
    sec = readDataArray(input, vec_types, numberOfCells_, CELLS_DATA_ARRAY, CELLS);
Praetorius, Simon's avatar
Praetorius, Simon committed
474
    VTK_ASSERT(vec_types.size() == numberOfCells_);
475
  } else if (id == "Cells.offsets") {
476
    sec = readDataArray(input, vec_offsets, numberOfCells_, CELLS_DATA_ARRAY, CELLS);
Praetorius, Simon's avatar
Praetorius, Simon committed
477
    VTK_ASSERT(vec_offsets.size() == numberOfCells_);
478
  } else if (id == "Cells.connectivity") {
Praetorius, Simon's avatar
Praetorius, Simon committed
479
    sec = readDataArray(input, vec_connectivity, std::size_t(-1), CELLS_DATA_ARRAY, CELLS);
480
  } else if (id == "Cells.global_point_ids") {
481
    sec = readDataArray(input, vec_point_ids, numberOfPoints_, CELLS_DATA_ARRAY, CELLS);
Praetorius, Simon's avatar
Praetorius, Simon committed
482
    VTK_ASSERT(vec_point_ids.size() == numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
483
484
485
486
487
488
  }

  return sec;
}


489
490
template <class Grid, class Creator, class Field>
void VtkReader<Grid,Creator,Field>::readCellsAppended (std::ifstream& input)
Praetorius, Simon's avatar
Praetorius, Simon committed
491
{
Praetorius, Simon's avatar
Praetorius, Simon committed
492
  VTK_ASSERT(numberOfCells_ > 0);
493
494
495
  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
496

Praetorius, Simon's avatar
Praetorius, Simon committed
497
  VTK_ASSERT(types_data.type == Vtk::UINT8);
498
  readAppended(input, vec_types, types_data.offset);
Praetorius, Simon's avatar
Praetorius, Simon committed
499
  VTK_ASSERT(vec_types.size() == numberOfCells_);
Praetorius, Simon's avatar
Praetorius, Simon committed
500

501
502
503
504
505
506
507
508
509
  if (offsets_data.type == Vtk::INT64)
    readAppended(input, vec_offsets, offsets_data.offset);
  else if (offsets_data.type == Vtk::INT32) {
    std::vector<std::int32_t> offsets;
    readAppended(input, offsets, offsets_data.offset);
    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
510
  VTK_ASSERT(vec_offsets.size() == numberOfCells_);
Praetorius, Simon's avatar
Praetorius, Simon committed
511

512
513
514
515
516
517
518
519
520
  if (connectivity_data.type == Vtk::INT64)
    readAppended(input, vec_connectivity, connectivity_data.offset);
  else if (connectivity_data.type == Vtk::INT32) {
    std::vector<std::int32_t> connectivity;
    readAppended(input, connectivity, connectivity_data.offset);
    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
521
  VTK_ASSERT(vec_connectivity.size() == std::size_t(vec_offsets.back()));
522

523
524
  if (dataArray_.count("Cells.global_point_ids") > 0) {
    auto point_id_data = dataArray_["Cells.global_point_ids"];
Praetorius, Simon's avatar
Praetorius, Simon committed
525
    VTK_ASSERT(point_id_data.type == Vtk::UINT64);
526
    readAppended(input, vec_point_ids, point_id_data.offset);
Praetorius, Simon's avatar
Praetorius, Simon committed
527
    VTK_ASSERT(vec_point_ids.size() == numberOfPoints_);
528
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
529
530
531
}


532
533
534
535
536
537
538
539
// @{ implementation detail
/**
 * 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
540
541
542
543
template <class T, class IStream>
void read_compressed (T* buffer, unsigned char* buffer_in,
                      std::uint64_t bs, std::uint64_t cbs, IStream& input)
{
544
#if HAVE_VTK_ZLIB
Praetorius, Simon's avatar
Praetorius, Simon committed
545
546
547
548
549
550
  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);

551
  input.read((char*)(compressed_buffer), compressed_space);
Praetorius, Simon's avatar
Praetorius, Simon committed
552
  VTK_ASSERT(uLongf(input.gcount()) == compressed_space);
Praetorius, Simon's avatar
Praetorius, Simon committed
553
554
555
556
557

  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
558
  VTK_ASSERT(uLongf(bs) == uncompressed_space);
Praetorius, Simon's avatar
Praetorius, Simon committed
559
560
561
562
563
#else
  std::cerr << "Can not call read_compressed without compression enabled!\n";
  std::abort();
#endif
}
564
// @}
Praetorius, Simon's avatar
Praetorius, Simon committed
565
566


567
template <class Grid, class Creator, class Field>
Praetorius, Simon's avatar
Praetorius, Simon committed
568
  template <class T>
569
void VtkReader<Grid,Creator,Field>::readAppended (std::ifstream& input, std::vector<T>& values, std::uint64_t offset)
Praetorius, Simon's avatar
Praetorius, Simon committed
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
{
  input.seekg(offset0_ + offset);

  std::uint64_t size = 0;

  std::uint64_t num_blocks = 0;
  std::uint64_t block_size = 0;
  std::uint64_t last_block_size = 0;
  std::vector<std::uint64_t> cbs; // compressed block sizes

  // read total size / block-size(s)
  if (format_ == Vtk::COMPRESSED) {
    input.read((char*)&num_blocks, sizeof(std::uint64_t));
    input.read((char*)&block_size, sizeof(std::uint64_t));
    input.read((char*)&last_block_size, sizeof(std::uint64_t));

Praetorius, Simon's avatar
Praetorius, Simon committed
586
    VTK_ASSERT(block_size % sizeof(T) == 0);
Praetorius, Simon's avatar
Praetorius, Simon committed
587
588
589
590
591
592
593
594
595
596

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

    // size of the compressed blocks
    cbs.resize(num_blocks);
    input.read((char*)cbs.data(), num_blocks*sizeof(std::uint64_t));
  } else {
    input.read((char*)&size, sizeof(std::uint64_t));
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
597
  VTK_ASSERT(size > 0 && (size % sizeof(T)) == 0);
Praetorius, Simon's avatar
Praetorius, Simon committed
598
599
600
601
602
603
604
605
606
607
608
609
610
611
  values.resize(size / sizeof(T));

  if (format_ == Vtk::COMPRESSED) {
    // upper bound for compressed block-size
    std::uint64_t compressed_block_size = block_size + (block_size + 999)/1000 + 12;
    // number of values in the full blocks
    std::size_t num_values = block_size / sizeof(T);

    std::vector<unsigned char> buffer_in(compressed_block_size);
    for (std::size_t i = 0; i < std::size_t(num_blocks); ++i) {
      std::uint64_t bs = i < std::size_t(num_blocks-1) ? block_size : last_block_size;
      read_compressed(values.data() + i*num_values, buffer_in.data(), bs, cbs[i], input);
    }
  } else {
612
    input.read((char*)(values.data()), size);
Praetorius, Simon's avatar
Praetorius, Simon committed
613
    VTK_ASSERT(input.gcount() == std::streamsize(size));
Praetorius, Simon's avatar
Praetorius, Simon committed
614
615
616
617
  }
}


618
619
template <class Grid, class Creator, class Field>
void VtkReader<Grid,Creator,Field>::fillGridCreator (bool insertPieces)
Praetorius, Simon's avatar
Praetorius, Simon committed
620
{
Praetorius, Simon's avatar
Praetorius, Simon committed
621
622
623
  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
624

Praetorius, Simon's avatar
Praetorius, Simon committed
625
  if (!vec_points.empty())
626
    creator_->insertVertices(vec_points, vec_point_ids);
Praetorius, Simon's avatar
Praetorius, Simon committed
627
  if (!vec_types.empty())
628
    creator_->insertElements(vec_types, vec_offsets, vec_connectivity);
629
  if (insertPieces)
630
    creator_->insertPieces(pieces_);
631
}
Praetorius, Simon's avatar
Praetorius, Simon committed
632

633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665

// 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";
  }
}


666
// Assume input already read the line <AppendedData ...>
667
668
template <class Grid, class Creator, class Field>
std::uint64_t VtkReader<Grid,Creator,Field>::findAppendedDataPosition (std::ifstream& input) const
669
670
671
672
673
674
675
676
677
678
679
680
{
  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;
}


681
682
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
683
{
684
  closed = false;
Praetorius, Simon's avatar
Praetorius, Simon committed
685
686
687
688
689
690
691
692
693
694
695
696
697
698
  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);
699
700
      } else if (c == '/') {
        closed = true;
Praetorius, Simon's avatar
Praetorius, Simon committed
701
702
703
704
705
706
707
708
709
710
711
      }
      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;
712
      VTK_ASSERT_MSG( c == '"', "Format error!" );
Praetorius, Simon's avatar
Praetorius, Simon committed
713
714
715
716
717
718
719
720
721
722
723
724
725
726
      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:
727
      VTK_ASSERT_MSG(false, "Format error!");
Praetorius, Simon's avatar
Praetorius, Simon committed
728
729
730
731
732
733
    }
  }

  return attr;
}

Praetorius, Simon's avatar
Praetorius, Simon committed
734

735
736
template <class Grid, class Creator, class Field>
void VtkReader<Grid,Creator,Field>::clear ()
Praetorius, Simon's avatar
Praetorius, Simon committed
737
738
739
740
741
742
743
744
745
746
747
748
749
750
{
  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;
}

751
} // end namespace Dune