vtkreader.impl.hh 20.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
template <class Grid, class Creator>
20
void VtkReader<Grid,Creator>::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
  }
}


template <class Grid, class Creator>
42
void VtkReader<Grid,Creator>::readSerialFileFromStream (std::ifstream& input, bool fillCreator)
43
{
Praetorius, Simon's avatar
Praetorius, Simon committed
44
  clear();
45
46
  std::string compressor = "";
  std::string data_name = "", data_format = "";
Praetorius, Simon's avatar
Praetorius, Simon committed
47
48
  Vtk::DataTypes data_type = Vtk::UNKNOWN;
  std::size_t 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
113

      if (!attr["Name"].empty())
Stenger, Florian's avatar
Stenger, Florian committed
114
        data_name = Vtk::to_lower(attr["Name"]);
115
116
117
118
      else if (section == POINTS)
        data_name = "points";

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

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

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

137
138
139
140
141
142
143
      // Store attributes of DataArray
      dataArray_[data_name] = {data_type, data_components, data_offset};

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

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

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

      readCellsAppended(input);
      section = NO_SECTION; // finish reading after appended section
Praetorius, Simon's avatar
Praetorius, Simon committed
189
190
191
192
193
194
    }
    else if (isSection(line, "/AppendedData", section, APPENDED_DATA))
      section = VTK_FILE;

    switch (section) {
      case PD_DATA_ARRAY:
195
196
        if (data_type == Vtk::FLOAT32) {
          std::vector<float> values;
197
198
          section = readPointData(input, values, data_name);
        } else if (data_type == Vtk::FLOAT64) {
199
          std::vector<double> values;
200
          section = readPointData(input, values, data_name);
201
        }
Praetorius, Simon's avatar
Praetorius, Simon committed
202
203
        break;
      case POINTS_DATA_ARRAY:
204
        section = readPoints(input, data_name);
Praetorius, Simon's avatar
Praetorius, Simon committed
205
206
        break;
      case CD_DATA_ARRAY:
207
208
        if (data_type == Vtk::FLOAT32) {
          std::vector<float> values;
209
210
          section = readCellData(input, values, data_name);
        } else if (data_type == Vtk::FLOAT64) {
211
          std::vector<double> values;
212
          section = readCellData(input, values, data_name);
213
        }
Praetorius, Simon's avatar
Praetorius, Simon committed
214
215
        break;
      case CELLS_DATA_ARRAY:
216
        section = readCells(input, data_name);
Praetorius, Simon's avatar
Praetorius, Simon committed
217
218
219
220
221
222
223
224
225
226
227
        break;
      default:
        // do nothing
        break;
    }

    if (section == NO_SECTION)
      break;
  }

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

230
231
  if (fillCreator)
    fillGridCreator();
232
233
234
235
}


template <class Grid, class Creator>
236
void VtkReader<Grid,Creator>::readParallelFileFromStream (std::ifstream& input, int commRank, int commSize, bool fillCreator)
237
{
Praetorius, Simon's avatar
Praetorius, Simon committed
238
  clear();
239
240
241

  Sections section = NO_SECTION;
  for (std::string line; std::getline(input, line); ) {
Stenger, Florian's avatar
Stenger, Florian committed
242
    Vtk::ltrim(line);
243
244
245
246
247
248

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

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

    if (section == NO_SECTION)
      break;
  }

  if (section != NO_SECTION)
279
    DUNE_THROW(Dune::VtkError, "VTK-File is incomplete. It must end with </VTKFile>!");
280

281
282
  if (fillCreator)
    fillGridCreator();
Praetorius, Simon's avatar
Praetorius, Simon committed
283
284
285
}


286
287
288
289
290
291
292
// @{ 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
293
template <class IStream, class T, class Sections>
294
295
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
296
{
Praetorius, Simon's avatar
Praetorius, Simon committed
297
  values.reserve(max_size < std::size_t(-1) ? max_size : 0);
298
  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
299
300

  std::size_t idx = 0;
301
  for (std::string line; std::getline(input, line);) {
Stenger, Florian's avatar
Stenger, Florian committed
302
    Vtk::trim(line);
303
    if (line.substr(1,10) == "/DataArray")
Praetorius, Simon's avatar
Praetorius, Simon committed
304
      return parent_section;
Praetorius, Simon's avatar
Praetorius, Simon committed
305
306
    if (line[0] == '<')
      break;
Praetorius, Simon's avatar
Praetorius, Simon committed
307
308
309
310
311

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

  return section;
}
330
// @}
Praetorius, Simon's avatar
Praetorius, Simon committed
331
332


333
334
template <class Grid, class Creator>
typename VtkReader<Grid,Creator>::Sections
335
VtkReader<Grid,Creator>::readPoints (std::ifstream& input, std::string name)
Praetorius, Simon's avatar
Praetorius, Simon committed
336
337
{
  using T = typename GlobalCoordinate::value_type;
Praetorius, Simon's avatar
Praetorius, Simon committed
338
339
  VTK_ASSERT(numberOfPoints_ > 0);
  VTK_ASSERT(dataArray_["points"].components == 3u);
340

341
342
  Sections sec;

Praetorius, Simon's avatar
Praetorius, Simon committed
343
  std::vector<T> point_values;
344
  sec = readDataArray(input, point_values, 3*numberOfPoints_, POINTS_DATA_ARRAY, POINTS);
Praetorius, Simon's avatar
Praetorius, Simon committed
345
346
  if (sec != POINTS)
    sec = skipRestOfDataArray(input, POINTS_DATA_ARRAY, POINTS);
Praetorius, Simon's avatar
Praetorius, Simon committed
347
348
  VTK_ASSERT(sec == POINTS);
  VTK_ASSERT(point_values.size() == 3*numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
349
350
351

  // extract points from continuous values
  GlobalCoordinate p;
352
  vec_points.reserve(numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
353
  std::size_t idx = 0;
354
  for (std::size_t i = 0; i < numberOfPoints_; ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
355
356
357
    for (std::size_t j = 0; j < p.size(); ++j)
      p[j] = point_values[idx++];
    idx += (3u - p.size());
358
    vec_points.push_back(p);
Praetorius, Simon's avatar
Praetorius, Simon committed
359
360
361
362
363
364
  }

  return sec;
}


365
template <class Grid, class Creator>
Praetorius, Simon's avatar
Praetorius, Simon committed
366
  template <class T>
367
void VtkReader<Grid,Creator>::readPointsAppended (std::ifstream& input)
Praetorius, Simon's avatar
Praetorius, Simon committed
368
{
Praetorius, Simon's avatar
Praetorius, Simon committed
369
370
  VTK_ASSERT(numberOfPoints_ > 0);
  VTK_ASSERT(dataArray_["points"].components == 3u);
Praetorius, Simon's avatar
Praetorius, Simon committed
371
  std::vector<T> point_values;
372
  readAppended(input, point_values, dataArray_["points"].offset);
Praetorius, Simon's avatar
Praetorius, Simon committed
373
  VTK_ASSERT(point_values.size() == 3*numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
374
375
376

  // extract points from continuous values
  GlobalCoordinate p;
377
  vec_points.reserve(numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
378
  std::size_t idx = 0;
379
  for (std::size_t i = 0; i < numberOfPoints_; ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
380
381
382
    for (std::size_t j = 0; j < p.size(); ++j)
      p[j] = T(point_values[idx++]);
    idx += (3u - p.size());
383
    vec_points.push_back(p);
Praetorius, Simon's avatar
Praetorius, Simon committed
384
385
386
387
  }
}


388
389
390
template <class Grid, class Creator>
typename VtkReader<Grid,Creator>::Sections
VtkReader<Grid,Creator>::readCells (std::ifstream& input, std::string name)
Praetorius, Simon's avatar
Praetorius, Simon committed
391
392
393
{
  Sections sec = CELLS_DATA_ARRAY;

Praetorius, Simon's avatar
Praetorius, Simon committed
394
  VTK_ASSERT(numberOfCells_ > 0);
Praetorius, Simon's avatar
Praetorius, Simon committed
395
  if (name == "types") {
396
    sec = readDataArray(input, vec_types, numberOfCells_, CELLS_DATA_ARRAY, CELLS);
Praetorius, Simon's avatar
Praetorius, Simon committed
397
    VTK_ASSERT(vec_types.size() == numberOfCells_);
Praetorius, Simon's avatar
Praetorius, Simon committed
398
  } else if (name == "offsets") {
399
    sec = readDataArray(input, vec_offsets, numberOfCells_, CELLS_DATA_ARRAY, CELLS);
Praetorius, Simon's avatar
Praetorius, Simon committed
400
    VTK_ASSERT(vec_offsets.size() == numberOfCells_);
Praetorius, Simon's avatar
Praetorius, Simon committed
401
  } else if (name == "connectivity") {
Praetorius, Simon's avatar
Praetorius, Simon committed
402
    sec = readDataArray(input, vec_connectivity, std::size_t(-1), CELLS_DATA_ARRAY, CELLS);
403
404
  } else if (name == "global_point_ids") {
    sec = readDataArray(input, vec_point_ids, numberOfPoints_, CELLS_DATA_ARRAY, CELLS);
Praetorius, Simon's avatar
Praetorius, Simon committed
405
    VTK_ASSERT(vec_point_ids.size() == numberOfPoints_);
Praetorius, Simon's avatar
Praetorius, Simon committed
406
407
408
409
410
411
  }

  return sec;
}


412
413
template <class Grid, class Creator>
void VtkReader<Grid,Creator>::readCellsAppended (std::ifstream& input)
Praetorius, Simon's avatar
Praetorius, Simon committed
414
{
Praetorius, Simon's avatar
Praetorius, Simon committed
415
  VTK_ASSERT(numberOfCells_ > 0);
416
  auto types_data = dataArray_["types"];
417
  auto offsets_data = dataArray_["offsets"];
418
  auto connectivity_data = dataArray_["connectivity"];
Praetorius, Simon's avatar
Praetorius, Simon committed
419

Praetorius, Simon's avatar
Praetorius, Simon committed
420
  VTK_ASSERT(types_data.type == Vtk::UINT8);
421
  readAppended(input, vec_types, types_data.offset);
Praetorius, Simon's avatar
Praetorius, Simon committed
422
  VTK_ASSERT(vec_types.size() == numberOfCells_);
Praetorius, Simon's avatar
Praetorius, Simon committed
423

424
425
426
427
428
429
430
431
432
  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
433
  VTK_ASSERT(vec_offsets.size() == numberOfCells_);
Praetorius, Simon's avatar
Praetorius, Simon committed
434

435
436
437
438
439
440
441
442
443
  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
444
  VTK_ASSERT(vec_connectivity.size() == std::size_t(vec_offsets.back()));
445
446
447

  if (dataArray_.count("global_point_ids") > 0) {
    auto point_id_data = dataArray_["global_point_ids"];
Praetorius, Simon's avatar
Praetorius, Simon committed
448
    VTK_ASSERT(point_id_data.type == Vtk::UINT64);
449
    readAppended(input, vec_point_ids, point_id_data.offset);
Praetorius, Simon's avatar
Praetorius, Simon committed
450
    VTK_ASSERT(vec_point_ids.size() == numberOfPoints_);
451
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
452
453
454
}


455
456
457
458
459
460
461
462
// @{ 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
463
464
465
466
template <class T, class IStream>
void read_compressed (T* buffer, unsigned char* buffer_in,
                      std::uint64_t bs, std::uint64_t cbs, IStream& input)
{
467
#if HAVE_VTK_ZLIB
Praetorius, Simon's avatar
Praetorius, Simon committed
468
469
470
471
472
473
  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);

474
  input.read((char*)(compressed_buffer), compressed_space);
Praetorius, Simon's avatar
Praetorius, Simon committed
475
  VTK_ASSERT(uLongf(input.gcount()) == compressed_space);
Praetorius, Simon's avatar
Praetorius, Simon committed
476
477
478
479
480

  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
481
  VTK_ASSERT(uLongf(bs) == uncompressed_space);
Praetorius, Simon's avatar
Praetorius, Simon committed
482
483
484
485
486
#else
  std::cerr << "Can not call read_compressed without compression enabled!\n";
  std::abort();
#endif
}
487
// @}
Praetorius, Simon's avatar
Praetorius, Simon committed
488
489


490
template <class Grid, class Creator>
Praetorius, Simon's avatar
Praetorius, Simon committed
491
  template <class T>
492
void VtkReader<Grid,Creator>::readAppended (std::ifstream& input, std::vector<T>& values, std::uint64_t offset)
Praetorius, Simon's avatar
Praetorius, Simon committed
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
{
  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
509
    VTK_ASSERT(block_size % sizeof(T) == 0);
Praetorius, Simon's avatar
Praetorius, Simon committed
510
511
512
513
514
515
516
517
518
519

    // 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
520
  VTK_ASSERT(size > 0 && (size % sizeof(T)) == 0);
Praetorius, Simon's avatar
Praetorius, Simon committed
521
522
523
524
525
526
527
528
529
530
531
532
533
534
  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 {
535
    input.read((char*)(values.data()), size);
Praetorius, Simon's avatar
Praetorius, Simon committed
536
    VTK_ASSERT(input.gcount() == std::streamsize(size));
Praetorius, Simon's avatar
Praetorius, Simon committed
537
538
539
540
  }
}


541
template <class Grid, class Creator>
542
void VtkReader<Grid,Creator>::fillGridCreator (bool insertPieces)
Praetorius, Simon's avatar
Praetorius, Simon committed
543
{
Praetorius, Simon's avatar
Praetorius, Simon committed
544
545
546
  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
547

Praetorius, Simon's avatar
Praetorius, Simon committed
548
  if (!vec_points.empty())
549
    creator_->insertVertices(vec_points, vec_point_ids);
Praetorius, Simon's avatar
Praetorius, Simon committed
550
  if (!vec_types.empty())
551
    creator_->insertElements(vec_types, vec_offsets, vec_connectivity);
552
  if (insertPieces)
553
    creator_->insertPieces(pieces_);
554
}
Praetorius, Simon's avatar
Praetorius, Simon committed
555

556
557
558
// Assume input already read the line <AppendedData ...>
template <class Grid, class Creator>
std::uint64_t VtkReader<Grid,Creator>::findAppendedDataPosition (std::ifstream& input) const
559
560
561
562
563
564
565
566
567
568
569
570
{
  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;
}


571
572
template <class Grid, class Creator>
std::map<std::string, std::string> VtkReader<Grid,Creator>::parseXml (std::string const& line, bool& closed)
Praetorius, Simon's avatar
Praetorius, Simon committed
573
{
574
  closed = false;
Praetorius, Simon's avatar
Praetorius, Simon committed
575
576
577
578
579
580
581
582
583
584
585
586
587
588
  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);
589
590
      } else if (c == '/') {
        closed = true;
Praetorius, Simon's avatar
Praetorius, Simon committed
591
592
593
594
595
596
597
598
599
600
601
      }
      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;
Praetorius, Simon's avatar
Praetorius, Simon committed
602
      VTK_ASSERT( c == '"' && "Format error!" );
Praetorius, Simon's avatar
Praetorius, Simon committed
603
604
605
606
607
608
609
610
611
612
613
614
615
616
      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:
Praetorius, Simon's avatar
Praetorius, Simon committed
617
      VTK_ASSERT(false && "Format error!");
Praetorius, Simon's avatar
Praetorius, Simon committed
618
619
620
621
622
623
    }
  }

  return attr;
}

Praetorius, Simon's avatar
Praetorius, Simon committed
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640

template <class Grid, class Creator>
void VtkReader<Grid,Creator>::clear ()
{
  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;
}

641
} // end namespace Dune