10#include <dune/common/classname.hh>
11#include <dune/common/version.hh>
19template <
class Gr
id,
class Creator,
class Field>
24 DUNE_THROW(IOError,
"File " << filename <<
" does not exist!");
26 std::ifstream input(filename, std::ios_base::in | std::ios_base::binary);
31 readSerialFileFromStream(input, fillCreator);
32 pieces_.push_back(filename);
33 }
else if (ext ==
".pvtu") {
34 readParallelFileFromStream(input, comm().rank(), comm().size(), fillCreator);
36 DUNE_THROW(
Dune::VtkError,
"File has unknown file-extension '" << ext <<
"'. Allowed are only '.vtu' and '.pvtu'.");
42template <
class Gr
id,
class Creator,
class Field>
48 std::string data_id =
"";
49 std::string data_format =
"";
51 unsigned int data_components = 0;
52 std::uint64_t data_offset = 0;
54 Sections section = NO_SECTION;
55 for (std::string line; std::getline(input, line); ) {
58 if (isSection(line,
"VTKFile", section)) {
60 auto attr = parseXml(line, closed);
62 if (!attr[
"type"].empty())
63 VTK_ASSERT_MSG(attr[
"type"] ==
"UnstructuredGrid",
"VtkReader supports UnstructuredGrid types");
64 if (!attr[
"byte_order"].empty())
65 VTK_ASSERT_MSG(attr[
"byte_order"] ==
"LittleEndian",
"LittleEndian byte order supported");
67 if (attr[
"header_type"] ==
"UInt32")
69 else if (attr[
"header_type"] ==
"UInt64")
72 if (attr[
"compressor"] ==
"vtkZLibDataCompressor")
74 else if (attr[
"compressor"] ==
"vtkLZ4DataCompressor")
76 else if (attr[
"compressor"] ==
"vtkLZMADataCompressor")
81 else if (isSection(line,
"/VTKFile", section, VTK_FILE)) {
85 else if (isSection(line,
"UnstructuredGrid", section, VTK_FILE))
86 section = UNSTRUCTURED_GRID;
87 else if (isSection(line,
"/UnstructuredGrid", section, UNSTRUCTURED_GRID))
89 else if (isSection(line,
"Piece", section, UNSTRUCTURED_GRID)) {
91 auto attr = parseXml(line, closed);
93 VTK_ASSERT_MSG(attr.count(
"NumberOfPoints") > 0 && attr.count(
"NumberOfCells") > 0,
94 "Number of points or cells in file must be > 0");
95 numberOfPoints_ = std::stoul(attr[
"NumberOfPoints"]);
96 numberOfCells_ = std::stoul(attr[
"NumberOfCells"]);
99 else if (isSection(line,
"/Piece", section, PIECE))
100 section = UNSTRUCTURED_GRID;
101 else if (isSection(line,
"PointData", section, PIECE))
102 section = POINT_DATA;
103 else if (isSection(line,
"/PointData", section, POINT_DATA))
105 else if (isSection(line,
"CellData", section, PIECE))
107 else if (isSection(line,
"/CellData", section, CELL_DATA))
109 else if (isSection(line,
"Points", section, PIECE))
111 else if (isSection(line,
"/Points", section, POINTS))
113 else if (isSection(line,
"Cells", section, PIECE))
115 else if (isSection(line,
"/Cells", section, CELLS))
117 else if (line.substr(1,9) ==
"DataArray") {
119 auto attr = parseXml(line, closed);
124 data_id = toString(section) +
"." + attr[
"Name"];
126 if (section == POINTS)
131 if (!attr[
"NumberOfComponents"].empty())
132 data_components = std::stoul(attr[
"NumberOfComponents"]);
136 if (data_format ==
"appended") {
144 if (!attr[
"offset"].empty()) {
145 data_offset = std::stoul(attr[
"offset"]);
146 VTK_ASSERT_MSG(data_format ==
"appended",
"Attribute 'offset' only supported by appended mode");
150 dataArray_[data_id] = {attr[
"Name"], data_type, data_components, data_offset, section};
153 if (data_format ==
"appended") {
155 while (std::getline(input, line)) {
157 if (line.substr(1,10) ==
"/DataArray")
164 if (section == POINT_DATA)
165 section = PD_DATA_ARRAY;
166 else if (section == POINTS)
167 section = POINTS_DATA_ARRAY;
168 else if (section == CELL_DATA)
169 section = CD_DATA_ARRAY;
170 else if (section == CELLS)
171 section = CELLS_DATA_ARRAY;
175 else if (line.substr(1,10) ==
"/DataArray") {
176 if (section == PD_DATA_ARRAY)
177 section = POINT_DATA;
178 else if (section == POINTS_DATA_ARRAY)
180 else if (section == CD_DATA_ARRAY)
182 else if (section == CELLS_DATA_ARRAY)
187 else if (isSection(line,
"AppendedData", section, VTK_FILE)) {
189 auto attr = parseXml(line, closed);
190 if (!attr[
"encoding"].empty())
191 VTK_ASSERT_MSG(attr[
"encoding"] ==
"raw",
"base64 encoding not supported");
193 offset0_ = findAppendedDataPosition(input);
194 Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(dataArray_[
"Points"].type, header_type,
195 [&](
auto f,
auto h) {
196 this->readPointsAppended(f,h,input);
197 this->readCellsAppended(h,input);
201 for (
auto const& d : dataArray_) {
202 if (d.second.section == POINT_DATA) {
203 Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(d.second.type, header_type,
204 [&](
auto f,
auto h) {
205 this->readPointDataAppended(f,h,input,d.first);
208 else if (d.second.section == CELL_DATA) {
209 Vtk::mapDataTypes<std::is_floating_point, std::is_integral>(d.second.type, header_type,
210 [&](
auto f,
auto h) {
211 this->readCellDataAppended(f,h,input,d.first);
216 section = NO_SECTION;
219 else if (isSection(line,
"/AppendedData", section, APPENDED_DATA))
224 section = readPointData(input, data_id);
226 case POINTS_DATA_ARRAY:
227 section = readPoints(input, data_id);
230 section = readCellData(input, data_id);
232 case CELLS_DATA_ARRAY:
233 section = readCells(input, data_id);
240 if (section != NO_SECTION)
241 DUNE_THROW(
Dune::VtkError,
"VTK-File is incomplete. It must end with </VTKFile>!");
248template <
class Gr
id,
class Creator,
class Field>
256 Sections section = NO_SECTION;
257 for (std::string line; std::getline(input, line); ) {
260 if (isSection(line,
"VTKFile", section)) {
262 auto attr = parseXml(line, closed);
264 if (!attr[
"type"].empty())
265 VTK_ASSERT_MSG(attr[
"type"] ==
"PUnstructuredGrid",
"VtkReader supports PUnstructuredGrid types");
266 if (!attr[
"version"].empty())
267 VTK_ASSERT_MSG(std::stod(attr[
"version"]) == 1.0,
"File format must be 1.0");
268 if (!attr[
"byte_order"].empty())
269 VTK_ASSERT_MSG(attr[
"byte_order"] ==
"LittleEndian",
"LittleEndian byte order supported");
271 if (attr[
"header_type"] ==
"UInt32")
273 else if (attr[
"header_type"] ==
"UInt64")
276 if (attr[
"compressor"] ==
"vtkZLibDataCompressor")
278 else if (attr[
"compressor"] ==
"vtkLZ4DataCompressor")
280 else if (attr[
"compressor"] ==
"vtkLZMADataCompressor")
285 else if (isSection(line,
"/VTKFile", section, VTK_FILE)) {
286 section = NO_SECTION;
288 }
else if (isSection(line,
"PUnstructuredGrid", section, VTK_FILE))
289 section = UNSTRUCTURED_GRID;
290 else if (isSection(line,
"/PUnstructuredGrid", section, UNSTRUCTURED_GRID))
292 else if (isSection(line,
"Piece", section, UNSTRUCTURED_GRID)) {
294 auto attr = parseXml(line, closed);
296 VTK_ASSERT_MSG(attr.count(
"Source") > 0,
"No source files for partitions provided");
297 pieces_.push_back(attr[
"Source"]);
301 VTK_ASSERT_MSG(section == NO_SECTION,
"VTK-File is incomplete. It must end with </VTKFile>!");
315template <
class IStream,
class T,
class Sections>
316Sections
readDataArray (IStream& input, std::vector<T>& values, std::size_t max_size,
317 Sections section, Sections parent_section)
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>;
323 for (std::string line; std::getline(input, line);) {
325 if (line.substr(1,10) ==
"/DataArray")
326 return parent_section;
330 std::istringstream stream(line);
332 for (; stream >> value; idx++)
333 values.push_back(T(value));
341template <
class IStream,
class Sections>
344 for (std::string line; std::getline(input, line);) {
346 if (line.substr(1,10) ==
"/DataArray")
347 return parent_section;
356template <
class Gr
id,
class Creator,
class Field>
357typename VtkReader<Grid,Creator,Field>::Sections
358VtkReader<Grid,Creator,Field>::readCellData (std::ifstream& input, std::string
id)
361 unsigned int components = dataArray_[id].components;
364 std::vector<Field>& values = cellData_[id];
365 sec =
readDataArray(input, values, components*numberOfCells_, CD_DATA_ARRAY, CELL_DATA);
366 if (sec != CELL_DATA)
369 VTK_ASSERT(values.size() == components*numberOfCells_);
375template <
class Gr
id,
class Creator,
class Field>
376 template <
class FloatType,
class HeaderType>
380 unsigned int components = dataArray_[id].components;
382 std::vector<FloatType> values;
383 readAppended(input, values, HeaderType(dataArray_[
id].offset));
384 VTK_ASSERT(values.size() == components*numberOfCells_);
386 cellData_[id].resize(values.size());
387 std::copy(values.begin(), values.end(), cellData_[
id].begin());
391template <
class Gr
id,
class Creator,
class Field>
392typename VtkReader<Grid,Creator,Field>::Sections
396 unsigned int components = dataArray_[id].components;
399 std::vector<Field>& values = pointData_[id];
400 sec =
readDataArray(input, values, components*numberOfPoints_, PD_DATA_ARRAY, POINT_DATA);
401 if (sec != POINT_DATA)
404 VTK_ASSERT(values.size() == components*numberOfPoints_);
410template <
class Gr
id,
class Creator,
class Field>
411 template <
class FloatType,
class HeaderType>
415 unsigned int components = dataArray_[id].components;
417 std::vector<FloatType> values;
418 readAppended(input, values, HeaderType(dataArray_[
id].offset));
419 VTK_ASSERT(values.size() == components*numberOfPoints_);
421 pointData_[id].resize(values.size());
422 std::copy(values.begin(), values.end(), pointData_[
id].begin());
426template <
class Gr
id,
class Creator,
class Field>
427typename VtkReader<Grid,Creator,Field>::Sections
430 using T =
typename GlobalCoordinate::value_type;
433 VTK_ASSERT(dataArray_[
"Points"].components == 3u);
437 std::vector<T> point_values;
438 sec =
readDataArray(input, point_values, 3*numberOfPoints_, POINTS_DATA_ARRAY, POINTS);
442 VTK_ASSERT(point_values.size() == 3*numberOfPoints_);
446 vec_points.reserve(numberOfPoints_);
448 for (std::size_t i = 0; i < numberOfPoints_; ++i) {
449 for (std::size_t j = 0; j < p.size(); ++j)
450 p[j] = point_values[idx++];
451 idx += (3u - p.size());
452 vec_points.push_back(p);
459template <
class Gr
id,
class Creator,
class Field>
460 template <
class FloatType,
class HeaderType>
464 VTK_ASSERT(dataArray_[
"Points"].components == 3u);
465 std::vector<FloatType> point_values;
466 readAppended(input, point_values, HeaderType(dataArray_[
"Points"].offset));
467 VTK_ASSERT(point_values.size() == 3*numberOfPoints_);
471 vec_points.reserve(numberOfPoints_);
473 for (std::size_t i = 0; i < numberOfPoints_; ++i) {
474 for (std::size_t j = 0; j < p.size(); ++j)
475 p[j] = FloatType(point_values[idx++]);
476 idx += (3u - p.size());
477 vec_points.push_back(p);
482template <
class Gr
id,
class Creator,
class Field>
483typename VtkReader<Grid,Creator,Field>::Sections
486 Sections sec = CELLS_DATA_ARRAY;
489 if (
id ==
"Cells.types") {
490 sec =
readDataArray(input, vec_types, numberOfCells_, CELLS_DATA_ARRAY, CELLS);
491 VTK_ASSERT(vec_types.size() == numberOfCells_);
492 }
else if (
id ==
"Cells.offsets") {
493 sec =
readDataArray(input, vec_offsets, numberOfCells_, CELLS_DATA_ARRAY, CELLS);
494 VTK_ASSERT(vec_offsets.size() == numberOfCells_);
495 }
else if (
id ==
"Cells.connectivity") {
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);
499 VTK_ASSERT(vec_point_ids.size() == numberOfPoints_);
506template <
class Gr
id,
class Creator,
class Field>
507 template <
class HeaderType>
511 auto types_data = dataArray_[
"Cells.types"];
512 auto offsets_data = dataArray_[
"Cells.offsets"];
513 auto connectivity_data = dataArray_[
"Cells.connectivity"];
516 readAppended(input, vec_types, HeaderType(types_data.offset));
517 VTK_ASSERT(vec_types.size() == numberOfCells_);
520 readAppended(input, vec_offsets, HeaderType(offsets_data.offset));
522 std::vector<std::int32_t> offsets;
523 readAppended(input, offsets, HeaderType(offsets_data.offset));
524 vec_offsets.resize(offsets.size());
525 std::copy(offsets.begin(), offsets.end(), vec_offsets.begin());
527 else { DUNE_THROW(Dune::NotImplemented,
"Unsupported DataType in Cell offsets."); }
528 VTK_ASSERT(vec_offsets.size() == numberOfCells_);
531 readAppended(input, vec_connectivity, HeaderType(connectivity_data.offset));
533 std::vector<std::int32_t> connectivity;
534 readAppended(input, connectivity, HeaderType(connectivity_data.offset));
535 vec_connectivity.resize(connectivity.size());
536 std::copy(connectivity.begin(), connectivity.end(), vec_connectivity.begin());
538 else { DUNE_THROW(Dune::NotImplemented,
"Unsupported DataType in Cell connectivity."); }
539 VTK_ASSERT(vec_connectivity.size() == std::size_t(vec_offsets.back()));
541 if (dataArray_.count(
"Cells.global_point_ids") > 0) {
542 auto point_id_data = dataArray_[
"Cells.global_point_ids"];
544 readAppended(input, vec_point_ids, HeaderType(point_id_data.offset));
545 VTK_ASSERT(vec_point_ids.size() == numberOfPoints_);
560template <
class T,
class IStream>
561void read_compressed_zlib (T* buffer,
unsigned char* buffer_in,
562 std::uint64_t bs, std::uint64_t cbs, IStream& input)
565 uLongf uncompressed_space = uLongf(bs);
566 uLongf compressed_space = uLongf(cbs);
568 Bytef* compressed_buffer =
reinterpret_cast<Bytef*
>(buffer_in);
569 Bytef* uncompressed_buffer =
reinterpret_cast<Bytef*
>(buffer);
571 input.read((
char*)(compressed_buffer), compressed_space);
572 VTK_ASSERT(uLongf(input.gcount()) == compressed_space);
574 if (uncompress(uncompressed_buffer, &uncompressed_space, compressed_buffer, compressed_space) != Z_OK) {
575 std::cerr <<
"Zlib error while uncompressing data.\n";
580 std::cerr <<
"ZLib Compression not supported. Provide the ZLIB package to CMake." << std::endl;
585template <
class T,
class IStream>
586void read_compressed_lz4 (T* ,
unsigned char* ,
587 std::uint64_t , std::uint64_t , IStream& )
590 std::cerr <<
"LZ4 Compression not yet implemented" << std::endl;
593 std::cerr <<
"LZ4 Compression not supported. Provide the LZ4 package to CMake." << std::endl;
598template <
class T,
class IStream>
599void read_compressed_lzma (T* ,
unsigned char* ,
600 std::uint64_t , std::uint64_t , IStream& )
603 std::cerr <<
"LZMA Compression not yet implemented" << std::endl;
606 std::cerr <<
"LZMA Compression not supported. Provide the LZMA package to CMake." << std::endl;
615template <
class Gr
id,
class Creator,
class Field>
616 template <
class FloatType,
class HeaderType>
617void VtkReader<Grid,Creator,Field>::readAppended (std::ifstream& input, std::vector<FloatType>& values, HeaderType offset)
619 input.seekg(offset0_ + offset);
623 HeaderType num_blocks = 0;
624 HeaderType block_size = 0;
625 HeaderType last_block_size = 0;
626 std::vector<HeaderType> cbs;
630 input.read((
char*)&num_blocks,
sizeof(HeaderType));
631 input.read((
char*)&block_size,
sizeof(HeaderType));
632 input.read((
char*)&last_block_size,
sizeof(HeaderType));
634 VTK_ASSERT(block_size %
sizeof(FloatType) == 0);
637 size = block_size * (num_blocks-1) + last_block_size;
640 cbs.resize(num_blocks);
641 input.read((
char*)cbs.data(), num_blocks*
sizeof(HeaderType));
643 input.read((
char*)&size,
sizeof(HeaderType));
645 VTK_ASSERT(size > 0 && (size %
sizeof(FloatType)) == 0);
646 values.resize(size /
sizeof(FloatType));
650 HeaderType compressed_block_size = block_size + (block_size + 999)/1000 + 12;
652 std::size_t num_values = block_size /
sizeof(FloatType);
654 std::vector<unsigned char> buffer_in(compressed_block_size);
655 for (std::size_t i = 0; i < std::size_t(num_blocks); ++i) {
656 HeaderType bs = i < std::size_t(num_blocks-1) ? block_size : last_block_size;
658 switch (compressor_) {
660 read_compressed_zlib(values.data() + i*num_values, buffer_in.data(), bs, cbs[i], input);
663 read_compressed_lz4(values.data() + i*num_values, buffer_in.data(), bs, cbs[i], input);
666 read_compressed_lzma(values.data() + i*num_values, buffer_in.data(), bs, cbs[i], input);
674 input.read((
char*)(values.data()), size);
675 VTK_ASSERT(input.gcount() == std::streamsize(size));
680template <
class Gr
id,
class Creator,
class Field>
683 VTK_ASSERT(vec_points.size() == numberOfPoints_);
684 VTK_ASSERT(vec_types.size() == numberOfCells_);
685 VTK_ASSERT(vec_offsets.size() == numberOfCells_);
687 if (!vec_points.empty())
688 creator_->insertVertices(vec_points, vec_point_ids);
689 if (!vec_types.empty())
690 creator_->insertElements(vec_types, vec_offsets, vec_connectivity);
692 creator_->insertPieces(pieces_);
697template <
class Gr
id,
class Creator,
class Field>
703 case UNSTRUCTURED_GRID:
704 return "UnstructuredGrid";
716 return "AppendedData";
719 case POINTS_DATA_ARRAY:
720 case CELLS_DATA_ARRAY:
729template <
class Gr
id,
class Creator,
class Field>
730std::uint64_t VtkReader<Grid,Creator,Field>::findAppendedDataPosition (std::ifstream& input)
const
733 while (input.get(c) && std::isblank(c)) { }
735 std::uint64_t offset = input.tellg();
743template <
class Gr
id,
class Creator,
class Field>
744std::map<std::string, std::string> VtkReader<Grid,Creator,Field>::parseXml (std::string
const& line,
bool& closed)
747 std::map<std::string, std::string> attr;
749 Sections sec = NO_SECTION;
752 std::string name =
"";
753 std::string value =
"";
754 for (
auto c : line) {
757 if (std::isalpha(c) || c ==
'_') {
761 }
else if (c ==
'/') {
766 if (std::isalpha(c) || c ==
'_')
769 sec = (c ==
'=' ? XML_NAME_ASSIGN : NO_SECTION);
771 case XML_NAME_ASSIGN:
778 if (c ==
'"' && !escape) {
781 }
else if (c ==
'\\' && !escape) {
797template <
class Gr
id,
class Creator,
class Field>
798void VtkReader<Grid,Creator,Field>::clear ()
801 vec_point_ids.clear();
804 vec_connectivity.clear();
Macro for wrapping error checks and throwing exceptions.
#define VTK_ASSERT_MSG(cond, text)
check if condition cond holds; otherwise, throw a VtkError with a message.
Definition: errors.hh:19
#define VTK_ASSERT(cond)
check if condition cond holds; otherwise, throw a VtkError.
Definition: errors.hh:29
Sections skipRestOfDataArray(IStream &input, Sections section, Sections parent_section)
Definition: vtkreader.impl.hh:342
Sections readDataArray(IStream &input, std::vector< T > &values, std::size_t max_size, Sections section, Sections parent_section)
Definition: vtkreader.impl.hh:316
std::string & ltrim(std::string &str)
trim a string from the left
Definition: string.hh:30
std::string & trim(std::string &str)
trim a string from both sides
Definition: string.hh:52
Vtk::DataTypes dataTypeOf()
Definition: types.hh:79
std::string to_lower(std::string input)
convert all characters in a string to upper case
Definition: string.hh:22
DataTypes
Definition: types.hh:52
bool exists(Path const &p)
Test whether the path is a valid (existing and accessible) file / directory.
Definition: filesystem.cc:134
Definition: filesystem.hh:15
Path extension() const
Returns the file extension path component.
Definition: filesystem.cc:77
std::string string() const
Return the path as string.
Definition: filesystem.cc:27
File-Reader for Vtk unstructured .vtu files.
Definition: vtkreader.hh:37
void read(std::string const &filename, bool fillCreator=true)
Read the grid from file with filename into the GridCreator.
Definition: vtkreader.impl.hh:20
void readParallelFileFromStream(std::ifstream &input, int rank, int size, bool create=true)
Read the grid from and input stream, referring to a .pvtu file, into the GridFactory factory_.
Definition: vtkreader.impl.hh:249
void fillGridCreator(bool insertPieces=true)
Definition: vtkreader.impl.hh:681
void readSerialFileFromStream(std::ifstream &input, bool create=true)
Read the grid from an input stream, referring to a .vtu file, into the GridFactory factory_.
Definition: vtkreader.impl.hh:43