Merge branch 'issue/dune_grid_files' into 'master'

Implement reader for dune-grid vtu files

See merge request extensions/dune-vtk!8

dune/vtk/types.hh

@@ -33,6 +33,13 @@ namespace Dune
     };
     std::string to_string (DataTypes);
 
+    enum CompressorTypes {
+      NONE = 0,
+      ZLIB,
+      LZ4,
+      LZMA
+    };
+
     enum CellParametrization {
       LINEAR,
       QUADRATIC,

dune/vtk/vtkreader.hh

@@ -6,6 +6,7 @@
 #include <vector>
 
 #include <dune/common/shared_ptr.hh>
+#include <dune/common/typelist.hh>
 #include <dune/common/typeutilities.hh>
 
 #include <dune/vtk/filereader.hh>
@@ -165,35 +166,35 @@ namespace Dune
 #endif
 
   private:
-    // Read values stored on the cells with name `name`
-    Sections readCellData (std::ifstream& /*input*/, std::string /*name*/);
+    // Read values stored on the cells with ID `id`
+    Sections readCellData (std::ifstream& input, std::string id);
 
-    template <class T>
-    void readCellDataAppended (std::ifstream& /*input*/, std::string /*name*/);
+    template <class F, class H>
+    void readCellDataAppended (MetaType<F>, MetaType<H>, std::ifstream& input, std::string id);
 
-    // Read values stored on the points with name `name`
-    Sections readPointData (std::ifstream& /*input*/, std::string /*name*/);
+    // Read values stored on the points with ID `id`
+    Sections readPointData (std::ifstream& input, std::string id);
 
-    template <class T>
-    void readPointDataAppended (std::ifstream& /*input*/, std::string /*name*/);
+    template <class F, class H>
+    void readPointDataAppended (MetaType<F>, MetaType<H>, std::ifstream& input, std::string id);
 
 
     // Read vertex coordinates from `input` stream and store in into `factory`
-    Sections readPoints (std::ifstream& input, std::string name);
+    Sections readPoints (std::ifstream& input, std::string id);
 
-    template <class T>
-    void readPointsAppended (std::ifstream& input);
+    template <class F, class H>
+    void readPointsAppended (MetaType<F>, MetaType<H>, std::ifstream& input);
 
 
     // Read cell type, cell offsets and connectivity from `input` stream
-    Sections readCells (std::ifstream& input, std::string name);
-
-    void readCellsAppended (std::ifstream& input);
+    Sections readCells (std::ifstream& input, std::string id);
 
+    template <class H>
+    void readCellsAppended (MetaType<H>, std::ifstream& input);
 
     // Read data from appended section in vtk file, starting from `offset`
-    template <class T>
-    void readAppended (std::ifstream& input, std::vector<T>& values, std::uint64_t offset);
+    template <class FloatType, class HeaderType>
+    void readAppended (std::ifstream& input, std::vector<FloatType>& values, HeaderType offset);
 
     // Test whether line belongs to section
     bool isSection (std::string line,
@@ -230,6 +231,9 @@ namespace Dune
     /// Data format, i.e. ASCII, BINARY or COMPRESSED. Read from xml attributes.
     Vtk::FormatTypes format_;
 
+    /// Type of compression algorithm used for binary data
+    Vtk::CompressorTypes compressor_;
+
     // Temporary data to construct the grid elements
     std::vector<GlobalCoordinate> vec_points;
     std::vector<std::uint64_t> vec_point_ids;   //< Global unique vertex ID

dune/vtk/vtkreader.impl.hh

@@ -42,8 +42,10 @@ template <class Grid, class Creator, class Field>
 void VtkReader<Grid,Creator,Field>::readSerialFileFromStream (std::ifstream& input, bool fillCreator)
 {
   clear();
-  std::string compressor = "";
-  std::string data_id = "", data_format = "";
+  compressor_ = Vtk::NONE;
+  Vtk::DataTypes header_type = Vtk::UINT32;
+  std::string data_id = "";
+  std::string data_format = "";
   Vtk::DataTypes data_type = Vtk::UNKNOWN;
   unsigned int data_components = 0;
   std::uint64_t data_offset = 0;
@@ -58,20 +60,27 @@ void VtkReader<Grid,Creator,Field>::readSerialFileFromStream (std::ifstream& inp
 
       if (!attr["type"].empty())
         VTK_ASSERT_MSG(attr["type"] == "UnstructuredGrid", "VtkReader supports UnstructuredGrid types");
-      if (!attr["version"].empty())
-        VTK_ASSERT_MSG(std::stod(attr["version"]) == 1.0, "File format must be 1.0");
       if (!attr["byte_order"].empty())
         VTK_ASSERT_MSG(attr["byte_order"] == "LittleEndian", "LittleEndian byte order supported");
-      if (!attr["header_type"].empty())
-        VTK_ASSERT_MSG(attr["header_type"] == "UInt64", "Header integer type must be UInt64");
-      if (!attr["compressor"].empty()) {
-        compressor = attr["compressor"];
-        VTK_ASSERT_MSG(compressor == "vtkZLibDataCompressor", "Only ZLib compression supported");
-      }
+
+      if (attr["header_type"] == "UInt32")
+        header_type = Vtk::UINT32;
+      else if (attr["header_type"] == "UInt64")
+        header_type = Vtk::UINT64;
+
+      if (attr["compressor"] == "vtkZLibDataCompressor")
+        compressor_ = Vtk::ZLIB;
+      else if (attr["compressor"] == "vtkLZ4DataCompressor")
+        compressor_ = Vtk::LZ4;
+      else if (attr["compressor"] == "vtkLZMADataCompressor")
+        compressor_ = Vtk::LZMA;
+
       section = VTK_FILE;
     }
-    else if (isSection(line, "/VTKFile", section, VTK_FILE))
+    else if (isSection(line, "/VTKFile", section, VTK_FILE)) {
       section = NO_SECTION;
+      break;
+    }
     else if (isSection(line, "UnstructuredGrid", section, VTK_FILE))
       section = UNSTRUCTURED_GRID;
     else if (isSection(line, "/UnstructuredGrid", section, UNSTRUCTURED_GRID))
@@ -124,7 +133,7 @@ void VtkReader<Grid,Creator,Field>::readSerialFileFromStream (std::ifstream& inp
       // determine FormatType
       data_format = Vtk::to_lower(attr["format"]);
       if (data_format == "appended") {
-        format_ = !compressor.empty() ? Vtk::COMPRESSED : Vtk::BINARY;
+        format_ = compressor_ != Vtk::NONE ? Vtk::COMPRESSED : Vtk::BINARY;
       } else {
         format_ = Vtk::ASCII;
       }
@@ -181,30 +190,30 @@ void VtkReader<Grid,Creator,Field>::readSerialFileFromStream (std::ifstream& inp
         VTK_ASSERT_MSG(attr["encoding"] == "raw", "base64 encoding not supported");
 
       offset0_ = findAppendedDataPosition(input);
-      if (dataArray_["Points"].type == Vtk::FLOAT32)
-        readPointsAppended<float>(input);
-      else
-        readPointsAppended<double>(input);
-
-      readCellsAppended(input);
-
-      // read point an cell data
-      for (auto const& [id,data] : dataArray_) {
-        if (data.section == POINT_DATA) {
-          if (data.type == Vtk::FLOAT32)
-            readPointDataAppended<float>(input, id);
-          else
-            readPointDataAppended<double>(input, id);
+      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);
+          });
         }
-        else if (data.section == CELL_DATA) {
-          if (data.type == Vtk::FLOAT32)
-            readCellDataAppended<float>(input, id);
-          else
-            readCellDataAppended<double>(input, id);
+        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);
+          });
         }
       }
 
       section = NO_SECTION; // finish reading after appended section
+      break;
     }
     else if (isSection(line, "/AppendedData", section, APPENDED_DATA))
       section = VTK_FILE;
@@ -223,12 +232,8 @@ void VtkReader<Grid,Creator,Field>::readSerialFileFromStream (std::ifstream& inp
         section = readCells(input, data_id);
         break;
       default:
-        // do nothing
         break;
     }
-
-    if (section == NO_SECTION)
-      break;
   }
 
   if (section != NO_SECTION)
@@ -240,10 +245,13 @@ void VtkReader<Grid,Creator,Field>::readSerialFileFromStream (std::ifstream& inp
 
 
 template <class Grid, class Creator, class Field>
-void VtkReader<Grid,Creator,Field>::readParallelFileFromStream (std::ifstream& input, int commRank, int commSize, bool fillCreator)
+void VtkReader<Grid,Creator,Field>::readParallelFileFromStream (std::ifstream& input, int /* commRank */, int /* commSize */, bool fillCreator)
 {
   clear();
 
+  [[maybe_unused]] Vtk::DataTypes header_type = Vtk::UINT32;
+  compressor_ = Vtk::NONE;
+
   Sections section = NO_SECTION;
   for (std::string line; std::getline(input, line); ) {
     Vtk::ltrim(line);
@@ -258,15 +266,25 @@ void VtkReader<Grid,Creator,Field>::readParallelFileFromStream (std::ifstream& i
         VTK_ASSERT_MSG(std::stod(attr["version"]) == 1.0, "File format must be 1.0");
       if (!attr["byte_order"].empty())
         VTK_ASSERT_MSG(attr["byte_order"] == "LittleEndian", "LittleEndian byte order supported");
-      if (!attr["header_type"].empty())
-        VTK_ASSERT_MSG(attr["header_type"] == "UInt64", "Header integer type must be UInt64");
-      if (!attr["compressor"].empty())
-        VTK_ASSERT_MSG(attr["compressor"] == "vtkZLibDataCompressor", "Only ZLib compression supported");
+
+      if (attr["header_type"] == "UInt32")
+        header_type = Vtk::UINT32;
+      else if (attr["header_type"] == "UInt64")
+        header_type = Vtk::UINT64;
+
+      if (attr["compressor"] == "vtkZLibDataCompressor")
+        compressor_ = Vtk::ZLIB;
+      else if (attr["compressor"] == "vtkLZ4DataCompressor")
+        compressor_ = Vtk::LZ4;
+      else if (attr["compressor"] == "vtkLZMADataCompressor")
+        compressor_ = Vtk::LZMA;
+
       section = VTK_FILE;
     }
-    else if (isSection(line, "/VTKFile", section, VTK_FILE))
+    else if (isSection(line, "/VTKFile", section, VTK_FILE)) {
       section = NO_SECTION;
-    else if (isSection(line, "PUnstructuredGrid", section, VTK_FILE))
+      break;
+    } else if (isSection(line, "PUnstructuredGrid", section, VTK_FILE))
       section = UNSTRUCTURED_GRID;
     else if (isSection(line, "/PUnstructuredGrid", section, UNSTRUCTURED_GRID))
       section = VTK_FILE;
@@ -277,9 +295,6 @@ void VtkReader<Grid,Creator,Field>::readParallelFileFromStream (std::ifstream& i
       VTK_ASSERT_MSG(attr.count("Source") > 0, "No source files for partitions provided");
       pieces_.push_back(attr["Source"]);
     }
-
-    if (section == NO_SECTION)
-      break;
   }
 
   VTK_ASSERT_MSG(section == NO_SECTION, "VTK-File is incomplete. It must end with </VTKFile>!");
@@ -357,14 +372,14 @@ VtkReader<Grid,Creator,Field>::readCellData (std::ifstream& input, std::string i
 
 
 template <class Grid, class Creator, class Field>
-  template <class T>
-void VtkReader<Grid,Creator,Field>::readCellDataAppended (std::ifstream& input, std::string id)
+  template <class FloatType, class HeaderType>
+void VtkReader<Grid,Creator,Field>::readCellDataAppended (MetaType<FloatType>, MetaType<HeaderType>, std::ifstream& input, std::string id)
 {
   VTK_ASSERT(numberOfCells_ > 0);
   unsigned int components = dataArray_[id].components;
 
-  std::vector<T> values;
-  readAppended(input, values, dataArray_[id].offset);
+  std::vector<FloatType> values;
+  readAppended(input, values, HeaderType(dataArray_[id].offset));
   VTK_ASSERT(values.size() == components*numberOfCells_);
 
   cellData_[id].resize(values.size());
@@ -392,14 +407,14 @@ VtkReader<Grid,Creator,Field>::readPointData (std::ifstream& input, std::string
 
 
 template <class Grid, class Creator, class Field>
-  template <class T>
-void VtkReader<Grid,Creator,Field>::readPointDataAppended (std::ifstream& input, std::string id)
+  template <class FloatType, class HeaderType>
+void VtkReader<Grid,Creator,Field>::readPointDataAppended (MetaType<FloatType>, MetaType<HeaderType>, std::ifstream& input, std::string id)
 {
   VTK_ASSERT(numberOfPoints_ > 0);
   unsigned int components = dataArray_[id].components;
 
-  std::vector<T> values;
-  readAppended(input, values, dataArray_[id].offset);
+  std::vector<FloatType> values;
+  readAppended(input, values, HeaderType(dataArray_[id].offset));
   VTK_ASSERT(values.size() == components*numberOfPoints_);
 
   pointData_[id].resize(values.size());
@@ -413,6 +428,7 @@ VtkReader<Grid,Creator,Field>::readPoints (std::ifstream& input, std::string id)
 {
   using T = typename GlobalCoordinate::value_type;
   VTK_ASSERT(numberOfPoints_ > 0);
+  VTK_ASSERT(id == "Points");
   VTK_ASSERT(dataArray_["Points"].components == 3u);
 
   Sections sec;
@@ -440,13 +456,13 @@ VtkReader<Grid,Creator,Field>::readPoints (std::ifstream& input, std::string id)
 
 
 template <class Grid, class Creator, class Field>
-  template <class T>
-void VtkReader<Grid,Creator,Field>::readPointsAppended (std::ifstream& input)
+  template <class FloatType, class HeaderType>
+void VtkReader<Grid,Creator,Field>::readPointsAppended (MetaType<FloatType>, MetaType<HeaderType>, std::ifstream& input)
 {
   VTK_ASSERT(numberOfPoints_ > 0);
   VTK_ASSERT(dataArray_["Points"].components == 3u);
-  std::vector<T> point_values;
-  readAppended(input, point_values, dataArray_["Points"].offset);
+  std::vector<FloatType> point_values;
+  readAppended(input, point_values, HeaderType(dataArray_["Points"].offset));
   VTK_ASSERT(point_values.size() == 3*numberOfPoints_);
 
   // extract points from continuous values
@@ -455,7 +471,7 @@ void VtkReader<Grid,Creator,Field>::readPointsAppended (std::ifstream& input)
   std::size_t idx = 0;
   for (std::size_t i = 0; i < numberOfPoints_; ++i) {
     for (std::size_t j = 0; j < p.size(); ++j)
-      p[j] = T(point_values[idx++]);
+      p[j] = FloatType(point_values[idx++]);
     idx += (3u - p.size());
     vec_points.push_back(p);
   }
@@ -487,7 +503,8 @@ VtkReader<Grid,Creator,Field>::readCells (std::ifstream& input, std::string id)
 
 
 template <class Grid, class Creator, class Field>
-void VtkReader<Grid,Creator,Field>::readCellsAppended (std::ifstream& input)
+  template <class HeaderType>
+void VtkReader<Grid,Creator,Field>::readCellsAppended (MetaType<HeaderType>, std::ifstream& input)
 {
   VTK_ASSERT(numberOfCells_ > 0);
   auto types_data = dataArray_["Cells.types"];
@@ -495,14 +512,14 @@ void VtkReader<Grid,Creator,Field>::readCellsAppended (std::ifstream& input)
   auto connectivity_data = dataArray_["Cells.connectivity"];
 
   VTK_ASSERT(types_data.type == Vtk::UINT8);
-  readAppended(input, vec_types, types_data.offset);
+  readAppended(input, vec_types, HeaderType(types_data.offset));
   VTK_ASSERT(vec_types.size() == numberOfCells_);
 
   if (offsets_data.type == Vtk::INT64)
-    readAppended(input, vec_offsets, offsets_data.offset);
+    readAppended(input, vec_offsets, HeaderType(offsets_data.offset));
   else if (offsets_data.type == Vtk::INT32) {
     std::vector<std::int32_t> offsets;
-    readAppended(input, offsets, offsets_data.offset);
+    readAppended(input, offsets, HeaderType(offsets_data.offset));
     vec_offsets.resize(offsets.size());
     std::copy(offsets.begin(), offsets.end(), vec_offsets.begin());
   }
@@ -510,10 +527,10 @@ void VtkReader<Grid,Creator,Field>::readCellsAppended (std::ifstream& input)
   VTK_ASSERT(vec_offsets.size() == numberOfCells_);
 
   if (connectivity_data.type == Vtk::INT64)
-    readAppended(input, vec_connectivity, connectivity_data.offset);
+    readAppended(input, vec_connectivity, HeaderType(connectivity_data.offset));
   else if (connectivity_data.type == Vtk::INT32) {
     std::vector<std::int32_t> connectivity;
-    readAppended(input, connectivity, connectivity_data.offset);
+    readAppended(input, connectivity, HeaderType(connectivity_data.offset));
     vec_connectivity.resize(connectivity.size());
     std::copy(connectivity.begin(), connectivity.end(), vec_connectivity.begin());
   }
@@ -523,13 +540,15 @@ void VtkReader<Grid,Creator,Field>::readCellsAppended (std::ifstream& input)
   if (dataArray_.count("Cells.global_point_ids") > 0) {
     auto point_id_data = dataArray_["Cells.global_point_ids"];
     VTK_ASSERT(point_id_data.type == Vtk::UINT64);
-    readAppended(input, vec_point_ids, point_id_data.offset);
+    readAppended(input, vec_point_ids, HeaderType(point_id_data.offset));
     VTK_ASSERT(vec_point_ids.size() == numberOfPoints_);
   }
 }
 
 
 // @{ implementation detail
+namespace {
+
 /**
  * Read compressed data into `buffer_in`, uncompress it and store the result in
  * the concrete-data-type `buffer`
@@ -538,8 +557,8 @@ void VtkReader<Grid,Creator,Field>::readCellsAppended (std::ifstream& input)
  * \param input  Stream to read from.
  **/
 template <class T, class IStream>
-void read_compressed (T* buffer, unsigned char* buffer_in,
-                      std::uint64_t bs, std::uint64_t cbs, IStream& input)
+void read_compressed_zlib (T* buffer, unsigned char* buffer_in,
+                           std::uint64_t bs, std::uint64_t cbs, IStream& input)
 {
 #if HAVE_VTK_ZLIB
   uLongf uncompressed_space = uLongf(bs);
@@ -557,56 +576,98 @@ void read_compressed (T* buffer, unsigned char* buffer_in,
   }
   VTK_ASSERT(uLongf(bs) == uncompressed_space);
 #else
-  std::cerr << "Can not call read_compressed without compression enabled!\n";
+  std::cerr << "ZLib Compression not supported. Provide the ZLIB package to CMake." << std::endl;
   std::abort();
 #endif
+}
+
+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
+}
+
 }
 // @}
 
 
 template <class Grid, class Creator, class Field>
-  template <class T>
-void VtkReader<Grid,Creator,Field>::readAppended (std::ifstream& input, std::vector<T>& values, std::uint64_t offset)
+  template <class FloatType, class HeaderType>
+void VtkReader<Grid,Creator,Field>::readAppended (std::ifstream& input, std::vector<FloatType>& values, HeaderType offset)
 {
   input.seekg(offset0_ + offset);
 
-  std::uint64_t size = 0;
+  HeaderType 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
+  HeaderType num_blocks = 0;
+  HeaderType block_size = 0;
+  HeaderType last_block_size = 0;
+  std::vector<HeaderType> 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));
+  if (compressor_ != Vtk::NONE) {
+    input.read((char*)&num_blocks, sizeof(HeaderType));
+    input.read((char*)&block_size, sizeof(HeaderType));
+    input.read((char*)&last_block_size, sizeof(HeaderType));
 
-    VTK_ASSERT(block_size % sizeof(T) == 0);
+    VTK_ASSERT(block_size % sizeof(FloatType) == 0);
 
     // 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));
+    input.read((char*)cbs.data(), num_blocks*sizeof(HeaderType));
   } else {
-    input.read((char*)&size, sizeof(std::uint64_t));
+    input.read((char*)&size, sizeof(HeaderType));
   }
-  VTK_ASSERT(size > 0 && (size % sizeof(T)) == 0);
-  values.resize(size / sizeof(T));
+  VTK_ASSERT(size > 0 && (size % sizeof(FloatType)) == 0);
+  values.resize(size / sizeof(FloatType));
 
-  if (format_ == Vtk::COMPRESSED) {
+  if (compressor_ != Vtk::NONE) {
     // upper bound for compressed block-size
-    std::uint64_t compressed_block_size = block_size + (block_size + 999)/1000 + 12;
+    HeaderType 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::size_t num_values = block_size / sizeof(FloatType);
 
     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);
+      HeaderType bs = i < std::size_t(num_blocks-1) ? block_size : last_block_size;
+
+      switch (compressor_) {
+        case Vtk::ZLIB:
+          read_compressed_zlib(values.data() + i*num_values, buffer_in.data(), bs, cbs[i], input);
+          break;
+        case Vtk::LZ4:
+          read_compressed_lz4(values.data() + i*num_values, buffer_in.data(), bs, cbs[i], input);
+          break;
+        case Vtk::LZMA:
+          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;
+      }
     }
   } else {
     input.read((char*)(values.data()), size);

src/CMakeLists.txt

@@ -23,6 +23,10 @@ dune_add_test(SOURCES datareader.cc
   LINK_LIBRARIES dunevtk
   CMAKE_GUARD dune-functions_FOUND)
 
+dune_add_test(SOURCES dunegridvtkreader.cc
+  LINK_LIBRARIES dunevtk
+  CMAKE_GUARD "dune-uggrid_FOUND OR dune-alugrid_FOUND")
+
 dune_add_test(SOURCES lagrangepoints.cc
   LINK_LIBRARIES dunevtk)
 

src/dunegridvtkreader.cc

+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/grid/io/file/vtk.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/vtk/vtkreader.hh>
+#include <dune/vtk/gridcreators/discontinuousgridcreator.hh>
+#include <dune/vtk/utility/errors.hh>
+
+#if HAVE_DUNE_UGGRID
+#include <dune/grid/uggrid.hh>
+#elif HAVE_DUNE_ALUGRID
+#include <dune/alugrid/grid.hh>
+#endif
+
+using namespace Dune;
+
+template <class GridType>
+void test(std::string prefix, GridType const& grid)
+{
+  { // write conforming binary file
+    Dune::VTKWriter vtkWriter{grid.leafGridView(), Dune::VTK::DataMode::conforming, Dune::VTK::Precision::float32};
+    auto fn = vtkWriter.write(prefix + "_conforming_binary", Dune::VTK::OutputType::appendedraw);
+
+    // Read conforming binary file
+    std::cout << "read '" << fn << "'..." << std::endl;
+    auto gridPtr = VtkReader<GridType>::createGridFromFile(fn);
+
+    VTK_ASSERT(gridPtr->size(0) == grid.size(0));
+    VTK_ASSERT(gridPtr->size(GridType::dimension) == grid.size(GridType::dimension));
+    std::remove(fn.c_str());
+  }
+
+  { // write non-conforming ascii file
+    Dune::VTKWriter vtkWriter{grid.leafGridView(), Dune::VTK::DataMode::nonconforming, Dune::VTK::Precision::float32};
+    auto fn = vtkWriter.write(prefix + "_nonconforming_ascii", Dune::VTK::OutputType::ascii);
+
+    // Read non-conforming binary file
+    std::cout << "read '" << fn << "'..." << std::endl;
+    using DGC = Dune::Vtk::DiscontinuousGridCreator<GridType>;
+    auto gridPtr = VtkReader<GridType, DGC>::createGridFromFile(fn);
+
+    VTK_ASSERT(gridPtr->size(0) == grid.size(0));
+    VTK_ASSERT(gridPtr->size(GridType::dimension) == grid.size(GridType::dimension));
+    std::remove(fn.c_str());
+  }
+}
+
+#if HAVE_DUNE_UGGRID
+  template <int dim> using CubeGrid = Dune::UGGrid<dim>;
+  template <int dim> using SimplexGrid = Dune::UGGrid<dim>;
+#elif HAVE_DUNE_ALUGRID
+  template <int dim> using CubeGrid = Dune::ALUGrid<dim,dim,Dune::cube,Dune::nonconforming>;
+  template <int dim> using SimplexGrid = Dune::ALUGrid<dim,dim,Dune::simplex,Dune::conforming>;
+#endif
+
+int main(int argc, char** argv)
+{
+  Dune::MPIHelper::instance(argc, argv);