From 938c9f1e1f9e3250f003f948ab12fd4da0bd3d69 Mon Sep 17 00:00:00 2001
From: Simon Praetorius <simon.praetorius@tu-dresden.de>
Date: Mon, 3 Dec 2018 20:37:12 -0500
Subject: [PATCH] restructured gridcreators

---
 dune/vtk/CMakeLists.txt                       |   2 +-
 .../datacollectors/structureddatacollector.hh |  12 --
 dune/vtk/filereader.hh                        |   2 +
 dune/vtk/filewriter.hh                        |   2 +
 dune/vtk/forward.hh                           | 104 ++++++++++
 dune/vtk/gridcreator.hh                       | 165 ----------------
 dune/vtk/gridcreatorinterface.hh              | 114 +++++++++++
 dune/vtk/gridcreators/CMakeLists.txt          |   9 +
 dune/vtk/gridcreators/common.hh               |  22 +++
 .../vtk/gridcreators/continuousgridcreator.hh |  69 +++++++
 dune/vtk/gridcreators/derivedgridcreator.hh   |  51 +++++
 .../gridcreators/discontinuousgridcreator.hh  |  97 +++++++++
 dune/vtk/gridcreators/parallelgridcreator.hh  |  50 +++++
 dune/vtk/gridcreators/serialgridcreator.hh    |  71 +++++++
 dune/vtk/pvdwriter.hh                         |   1 +
 dune/vtk/vtkreader.hh                         |  52 +++--
 dune/vtk/vtkreader.impl.hh                    |  49 ++---
 dune/vtk/vtktimeserieswriter.hh               |   1 +
 dune/vtk/vtkwriter.hh                         |   1 +
 dune/vtk/vtkwriterinterface.hh                |   8 +-
 dune/vtk/writers/vtkimagedatawriter.hh        |   3 +-
 dune/vtk/writers/vtkrectilineargridwriter.hh  |   3 +-
 dune/vtk/writers/vtkstructuredgridwriter.hh   |   3 +-
 dune/vtk/writers/vtkunstructuredgridwriter.hh |   3 +-
 src/test/CMakeLists.txt                       |   3 +
 src/test/parallel_reader_writer_test.cc       | 185 ++++++++++++++++++
 src/test/reader_writer_test.cc                |  39 ++--
 src/vtkreader.cc                              |   4 +-
 28 files changed, 880 insertions(+), 245 deletions(-)
 create mode 100644 dune/vtk/forward.hh
 delete mode 100644 dune/vtk/gridcreator.hh
 create mode 100644 dune/vtk/gridcreatorinterface.hh
 create mode 100644 dune/vtk/gridcreators/CMakeLists.txt
 create mode 100644 dune/vtk/gridcreators/common.hh
 create mode 100644 dune/vtk/gridcreators/continuousgridcreator.hh
 create mode 100644 dune/vtk/gridcreators/derivedgridcreator.hh
 create mode 100644 dune/vtk/gridcreators/discontinuousgridcreator.hh
 create mode 100644 dune/vtk/gridcreators/parallelgridcreator.hh
 create mode 100644 dune/vtk/gridcreators/serialgridcreator.hh
 create mode 100644 src/test/parallel_reader_writer_test.cc

diff --git a/dune/vtk/CMakeLists.txt b/dune/vtk/CMakeLists.txt
index 7c6d8a5..64ee710 100644
--- a/dune/vtk/CMakeLists.txt
+++ b/dune/vtk/CMakeLists.txt
@@ -8,7 +8,6 @@ install(FILES
   defaultvtkfunction.hh
   filereader.hh
   filewriter.hh
-  gridcreator.hh
   legacyvtkfunction.hh
   pvdwriter.hh
   pvdwriter.impl.hh
@@ -26,5 +25,6 @@ install(FILES
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vtkwriter)
 
 add_subdirectory(datacollectors)
+add_subdirectory(gridcreators)
 add_subdirectory(utility)
 add_subdirectory(writers)
diff --git a/dune/vtk/datacollectors/structureddatacollector.hh b/dune/vtk/datacollectors/structureddatacollector.hh
index 01c6760..c7c1167 100644
--- a/dune/vtk/datacollectors/structureddatacollector.hh
+++ b/dune/vtk/datacollectors/structureddatacollector.hh
@@ -9,18 +9,6 @@
 
 namespace Dune
 {
-
-namespace Impl
-{
-  // Should be specialized for concrete structured grid
-  template <class GridView, class Grid>
-  struct StructuredDataCollectorImpl;
-}
-
-template <class GridView>
-using StructuredDataCollector = typename Impl::StructuredDataCollectorImpl<GridView, typename GridView::Grid>::type;
-
-
 /// The Interface for structured data-collectors
 template <class GridView, class Derived>
 class StructuredDataCollectorInterface
diff --git a/dune/vtk/filereader.hh b/dune/vtk/filereader.hh
index abc9d0c..ff7e36c 100644
--- a/dune/vtk/filereader.hh
+++ b/dune/vtk/filereader.hh
@@ -4,6 +4,8 @@
 #include <string>
 #include <utility>
 
+#include <dune/vtk/forward.hh>
+
 namespace Dune
 {
   template <class Grid, class FilerReaderImp>
diff --git a/dune/vtk/filewriter.hh b/dune/vtk/filewriter.hh
index 8873e75..3760185 100644
--- a/dune/vtk/filewriter.hh
+++ b/dune/vtk/filewriter.hh
@@ -3,6 +3,8 @@
 #include <string>
 #include <dune/common/std/optional.hh>
 
+#include <dune/vtk/forward.hh>
+
 namespace Dune
 {
   class FileWriter
diff --git a/dune/vtk/forward.hh b/dune/vtk/forward.hh
new file mode 100644
index 0000000..efc8538
--- /dev/null
+++ b/dune/vtk/forward.hh
@@ -0,0 +1,104 @@
+#pragma once
+
+namespace Dune
+{
+  // forward declaration of all classes in dune-vtk
+
+  template <class GridView, class Derived>
+  class DataCollectorInterface;
+
+  // @{ datacollectors
+  template <class GridView, class Derived>
+  class UnstructuredDataCollectorInterface;
+
+  // @{ unstructured-datacollectors
+  template <class GridView>
+  class ContinuousDataCollector;
+
+  template <class GridView>
+  class DiscontinuousDataCollector;
+
+  template <class GridView>
+  class QuadraticDataCollector;
+  // @} unstructured-datacollectors
+
+  template <class GridView, class Derived>
+  class StructuredDataCollectorInterface;
+
+  namespace Impl
+  {
+    // Should be specialized for concrete structured grid
+    template <class GridView, class Grid>
+    struct StructuredDataCollectorImpl;
+  }
+
+  template <class GridView>
+  using StructuredDataCollector = typename Impl::StructuredDataCollectorImpl<GridView, typename GridView::Grid>::type;
+
+  // @{ structured-datacollectors
+  template <class GridView>
+  class SPDataCollector;
+
+  template <class GridView>
+  class YaspDataCollector;
+  // @} structured-datacollectors
+
+  // @} datacollectors
+
+  template <class Grid, class Derived>
+  class GridCreatorInterface;
+
+  template <class GridCreator, class Derived>
+  struct DerivedGridCreator;
+
+  // @{ gridcreators
+  template <class Grid>
+  struct ContinuousGridCreator;
+
+  template <class Grid>
+  struct DiscontinuousGridCreator;
+
+  template <class Grid>
+  struct ParallelGridCreator;
+
+  template <class Grid>
+  struct SerialGridCreator;
+  // @} gridcreators
+
+
+  template <class Grid, class FilerReaderImp>
+  class FileReader;
+
+  template <class Grid, class GridCreator = ContinuousGridCreator<Grid>>
+  class VtkReader;
+
+
+  class FileWriter;
+
+  // @{ filewriters
+  template <class VtkWriter>
+  class PvdWriter;
+
+  template <class VtkWriter>
+  class VtkTimeseriesWriter;
+
+  template <class GridView, class DataCollector>
+  class VtkWriterInterface;
+
+  // @{ vtkwriters
+  template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+  class VtkImageDataWriter;
+
+  template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+  class VtkRectilinearGridWriter;
+
+  template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+  class VtkStructuredGridWriter;
+
+  template <class GridView, class DataCollector = ContinuousDataCollector<GridView>>
+  class VtkUnstructuredGridWriter;
+  // @} vtkwriters
+
+  // @} filewriters
+
+} // end namespace Dune
diff --git a/dune/vtk/gridcreator.hh b/dune/vtk/gridcreator.hh
deleted file mode 100644
index f7f4e43..0000000
--- a/dune/vtk/gridcreator.hh
+++ /dev/null
@@ -1,165 +0,0 @@
-#pragma once
-
-#include <cassert>
-#include <cstdint>
-#include <limits>
-#include <vector>
-
-#include <dune/common/exceptions.hh>
-#include <dune/common/hybridutilities.hh>
-#include <dune/grid/common/gridfactory.hh>
-
-#include "vtktypes.hh"
-
-namespace Dune
-{
-  template <class Factory, class... Args>
-  using HasInsertVertex = decltype( std::declval<Factory>().insertVertex(std::declval<Args>()...) );
-
-  namespace Impl
-  {
-    template <class GF, class = void>
-    struct VertexIdType { using type = unsigned int; };
-
-    template <class GF>
-    struct VertexIdType<GF, typename GF::VertexId> { using type = typename GF::VertexId; };
-  }
-
-  template <class GF>
-  using VertexId_t = typename Impl::VertexIdType<GF>::type;
-
-
-  // Create a grid where the input points and connectivity is already
-  // connected correctly.
-  struct DefaultGridCreator
-  {
-    template <class Grid, class Coord>
-    static void create (GridFactory<Grid>& factory,
-                        std::vector<Coord> const& points,
-                        std::vector<std::uint64_t> const& point_ids,
-                        std::vector<std::uint8_t> const& types,
-                        std::vector<std::int64_t> const& offsets,
-                        std::vector<std::int64_t> const& connectivity)
-    {
-      const auto hasInsertVertex = Std::is_detected<HasInsertVertex, GridFactory<Grid>, Coord, unsigned int>{};
-      if (point_ids.empty() || !hasInsertVertex.value) {
-        for (auto const& p : points)
-          factory.insertVertex(p);
-      } else {
-        Hybrid::ifElse(hasInsertVertex,
-        [&](auto id) {
-          using VertexId = VertexId_t<GridFactory<Grid>>;
-          for (std::size_t i = 0; i < points.size(); ++i)
-            id(factory).insertVertex(points[i], VertexId(point_ids[i]));
-        });
-      }
-
-      std::size_t idx = 0;
-      for (std::size_t i = 0; i < types.size(); ++i) {
-        auto type = Vtk::to_geometry(types[i]);
-        Vtk::CellType cellType{type};
-        auto refElem = referenceElement<double,Grid::dimension>(type);
-
-        int nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]);
-        assert(nNodes == refElem.size(Grid::dimension));
-        std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes);
-        for (int j = 0; j < nNodes; ++j)
-          vtk_cell.push_back( connectivity[idx++] );
-
-        if (cellType.noPermutation())
-          factory.insertElement(type,vtk_cell);
-        else {
-          // apply index permutation
-          std::vector<unsigned int> cell(nNodes);
-          for (int j = 0; j < nNodes; ++j)
-            cell[j] = vtk_cell[cellType.permutation(j)];
-
-          factory.insertElement(type,cell);
-        }
-      }
-    }
-  };
-
-
-  // Create a grid where the input points are not connected and the connectivity
-  // describes separated elements.
-  struct ConnectedGridCreator
-  {
-    struct CoordLess
-    {
-      template <class T, int N>
-      bool operator() (FieldVector<T,N> const& lhs, FieldVector<T,N> const& rhs) const
-      {
-        for (int i = 0; i < N; ++i) {
-          if (std::abs(lhs[i] - rhs[i]) < std::numeric_limits<T>::epsilon())
-            continue;
-          return lhs[i] < rhs[i];
-        }
-        return false;
-      }
-    };
-
-    template <class Grid, class Coord>
-    static void create (GridFactory<Grid>& factory,
-                        std::vector<Coord> const& points,
-                        std::vector<std::uint64_t> const& point_ids,
-                        std::vector<std::uint8_t> const& types,
-                        std::vector<std::int64_t> const& offsets,
-                        std::vector<std::int64_t> const& connectivity)
-    {
-      std::size_t idx = 0;
-      std::map<Coord, std::size_t, CoordLess> unique_points;
-
-      const auto hasInsertVertex = Std::is_detected<HasInsertVertex, GridFactory<Grid>, Coord, unsigned int>{};
-      if (point_ids.empty() || !hasInsertVertex.value) {
-        for (auto const& p : points) {
-          auto b = unique_points.emplace(std::make_pair(p,idx));
-          if (b.second) {
-            factory.insertVertex(p);
-            ++idx;
-          }
-        }
-      } else {
-        Hybrid::ifElse(hasInsertVertex,
-        [&](auto id) {
-          using VertexId = VertexId_t<GridFactory<Grid>>;
-          for (std::size_t i = 0; i < points.size(); ++i) {
-            auto b = unique_points.emplace(std::make_pair(points[i],idx));
-            if (b.second) {
-              id(factory).insertVertex(points[i], VertexId(point_ids[i]));
-              ++idx;
-            }
-          }
-        });
-      }
-
-
-      idx = 0;
-      for (std::size_t i = 0; i < types.size(); ++i) {
-        auto type = Vtk::to_geometry(types[i]);
-        Vtk::CellType cellType{type};
-
-        int nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]);
-        assert(nNodes > 0);
-        std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes);
-        for (int j = 0; j < nNodes; ++j) {
-          std::size_t v_j = connectivity[idx++];
-          std::size_t new_idx = unique_points[points[v_j]];
-          vtk_cell.push_back(new_idx);
-        }
-
-        if (cellType.noPermutation())
-          factory.insertElement(type,vtk_cell);
-        else {
-          // apply index permutation
-          std::vector<unsigned int> cell(nNodes);
-          for (int j = 0; j < nNodes; ++j)
-            cell[j] = vtk_cell[cellType.permutation(j)];
-
-          factory.insertElement(type,cell);
-        }
-      }
-    }
-  };
-
-} // end namespace Dune
diff --git a/dune/vtk/gridcreatorinterface.hh b/dune/vtk/gridcreatorinterface.hh
new file mode 100644
index 0000000..634a248
--- /dev/null
+++ b/dune/vtk/gridcreatorinterface.hh
@@ -0,0 +1,114 @@
+#pragma once
+
+#include <cstdint>
+#include <string>
+#include <vector>
+
+#include <dune/common/version.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/vtk/forward.hh>
+
+namespace Dune {
+
+template <class G, class Derived>
+class GridCreatorInterface
+{
+public:
+  using Grid = G;
+  using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
+
+public:
+  GridCreatorInterface (GridFactory<Grid>& factory)
+    : factory_(&factory)
+  {
+#if DUNE_VERSION_LT(DUNE_GRID,2,7)
+  // old GridFactory implementation does not provide access to its collective communication
+  #if HAVE_MPI
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank_);
+    MPI_Comm_size(MPI_COMM_WORLD, &numRanks_);
+  #endif
+#else
+    rank_ = factory.comm().rank();
+    numRanks_ = factory.comm().size();
+#endif
+  }
+
+  /// Insert all points as vertices into the factory
+  void insertVertices (std::vector<GlobalCoordinate> const& points,
+                       std::vector<std::uint64_t> const& point_ids)
+  {
+    asDerived().insertVerticesImpl(points, point_ids);
+  }
+
+  /// Create elements based on type and connectivity description
+  void insertElements (std::vector<std::uint8_t> const& types,
+                       std::vector<std::int64_t> const& offsets,
+                       std::vector<std::int64_t> const& connectivity)
+  {
+    asDerived().insertElementsImpl(types, offsets, connectivity);
+  }
+
+  /// Insert part of a grid stored in file into factory
+  void insertPieces (std::vector<std::string> const& pieces)
+  {
+    asDerived().insertPiecesImpl(pieces);
+  }
+
+  /// Return the associated GridFactory
+  GridFactory<Grid>& factory ()
+  {
+    return *factory_;
+  }
+
+  /// Return the MPI_Comm_rank of the factory, (or of MPI_COMM_WORLD)
+  int rank () const
+  {
+    return rank_;
+  }
+
+  /// Return the MPI_Comm_size of the factory, (or of MPI_COMM_WORLD)
+  int size () const
+  {
+    return numRanks_;
+  }
+
+protected: // cast to derived type
+
+  Derived& asDerived ()
+  {
+    return static_cast<Derived&>(*this);
+  }
+
+  const Derived& asDerived () const
+  {
+    return static_cast<const Derived&>(*this);
+  }
+
+public: // default implementations
+
+  void insertVerticesImpl (std::vector<GlobalCoordinate> const&,
+                           std::vector<std::uint64_t> const&)
+  {
+    /* do nothing */
+  }
+
+  void insertElementsImpl (std::vector<std::uint8_t> const&,
+                           std::vector<std::int64_t> const&,
+                           std::vector<std::int64_t> const&)
+  {
+    /* do nothing */
+  }
+
+  void insertPiecesImpl (std::vector<std::string> const&)
+  {
+    /* do nothing */;
+  }
+
+protected:
+  GridFactory<Grid>* factory_;
+  int rank_ = 0;
+  int numRanks_ = 1;
+};
+
+} // end namespace Dune
diff --git a/dune/vtk/gridcreators/CMakeLists.txt b/dune/vtk/gridcreators/CMakeLists.txt
new file mode 100644
index 0000000..299b4d3
--- /dev/null
+++ b/dune/vtk/gridcreators/CMakeLists.txt
@@ -0,0 +1,9 @@
+#install headers
+install(FILES
+  common.hh
+  continuousgridcreator.hh
+  derivedgridcreator.hh
+  discontinuousgridcreator.hh
+  parallelgridcreator.hh
+  serialgridcreator.hh
+  DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/vtkwriter/gridcreators)
diff --git a/dune/vtk/gridcreators/common.hh b/dune/vtk/gridcreators/common.hh
new file mode 100644
index 0000000..e47eb4d
--- /dev/null
+++ b/dune/vtk/gridcreators/common.hh
@@ -0,0 +1,22 @@
+#pragma once
+
+#include <type_traits>
+
+namespace Dune
+{
+  template <class Factory, class... Args>
+  using HasInsertVertex = decltype( std::declval<Factory>().insertVertex(std::declval<Args>()...) );
+
+  namespace Impl
+  {
+    template <class GF, class = void>
+    struct VertexIdType { using type = unsigned int; };
+
+    template <class GF>
+    struct VertexIdType<GF, typename GF::VertexId> { using type = typename GF::VertexId; };
+  }
+
+  template <class GF>
+  using VertexId_t = typename Impl::VertexIdType<GF>::type;
+
+} // end namespace Dune
diff --git a/dune/vtk/gridcreators/continuousgridcreator.hh b/dune/vtk/gridcreators/continuousgridcreator.hh
new file mode 100644
index 0000000..03270a2
--- /dev/null
+++ b/dune/vtk/gridcreators/continuousgridcreator.hh
@@ -0,0 +1,69 @@
+#pragma once
+
+#include <cassert>
+#include <cstdint>
+#include <limits>
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/vtk/vtktypes.hh>
+#include <dune/vtk/gridcreatorinterface.hh>
+
+namespace Dune
+{
+  // Create a grid where the input points and connectivity is already
+  // connected correctly.
+  template <class Grid>
+  struct ContinuousGridCreator
+      : public GridCreatorInterface<Grid, ContinuousGridCreator<Grid>>
+  {
+    using Super = GridCreatorInterface<Grid, ContinuousGridCreator<Grid>>;
+    using GlobalCoordinate = typename Super::GlobalCoordinate;
+
+    ContinuousGridCreator (GridFactory<Grid>& factory)
+      : Super(factory)
+    {}
+
+    using Super::factory;
+
+    void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
+                             std::vector<std::uint64_t> const& /*point_ids*/)
+    {
+      for (auto const& p : points)
+        factory().insertVertex(p);
+    }
+
+    void insertElementsImpl (std::vector<std::uint8_t> const& types,
+                             std::vector<std::int64_t> const& offsets,
+                             std::vector<std::int64_t> const& connectivity)
+    {
+      std::size_t idx = 0;
+      for (std::size_t i = 0; i < types.size(); ++i) {
+        auto type = Vtk::to_geometry(types[i]);
+        Vtk::CellType cellType{type};
+        auto refElem = referenceElement<double,Grid::dimension>(type);
+
+        int nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]);
+        assert(nNodes == refElem.size(Grid::dimension));
+        std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes);
+        for (int j = 0; j < nNodes; ++j)
+          vtk_cell.push_back( connectivity[idx++] );
+
+        if (cellType.noPermutation())
+          factory().insertElement(type,vtk_cell);
+        else {
+          // apply index permutation
+          std::vector<unsigned int> cell(nNodes);
+          for (int j = 0; j < nNodes; ++j)
+            cell[j] = vtk_cell[cellType.permutation(j)];
+
+          factory().insertElement(type,cell);
+        }
+      }
+    }
+  };
+
+} // end namespace Dune
diff --git a/dune/vtk/gridcreators/derivedgridcreator.hh b/dune/vtk/gridcreators/derivedgridcreator.hh
new file mode 100644
index 0000000..91fa846
--- /dev/null
+++ b/dune/vtk/gridcreators/derivedgridcreator.hh
@@ -0,0 +1,51 @@
+#pragma once
+
+#include <cstdint>
+#include <string>
+#include <vector>
+
+#include <dune/grid/common/gridfactory.hh>
+#include <dune/vtk/forward.hh>
+#include <dune/vtk/gridcreatorinterface.hh>
+#include <dune/vtk/gridcreators/common.hh>
+#include <dune/vtk/gridcreators/continuousgridcreator.hh>
+
+namespace Dune
+{
+  template <class GridCreator, class Derived>
+  struct DerivedGridCreator
+      : public GridCreatorInterface<typename GridCreator::Grid, Derived>
+  {
+    using Self = DerivedGridCreator;
+    using Super = GridCreatorInterface<typename GridCreator::Grid, Derived>;
+    using Grid = typename GridCreator::Grid;
+    using GlobalCoordinate = typename Super::GlobalCoordinate;
+
+    DerivedGridCreator (GridFactory<Grid>& factory)
+      : Super(factory)
+      , gridCreator_(factory)
+    {}
+
+    void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
+                             std::vector<std::uint64_t> const& point_ids)
+    {
+      gridCreator_.insertVertices(points, point_ids);
+    }
+
+    void insertElementsImpl (std::vector<std::uint8_t> const& types,
+                             std::vector<std::int64_t> const& offsets,
+                             std::vector<std::int64_t> const& connectivity)
+    {
+      gridCreator_.insertElements(types, offsets, connectivity);
+    }
+
+    void insertPiecesImpl (std::vector<std::string> const& pieces)
+    {
+      gridCreator_.insertPieces(pieces);
+    }
+
+  private:
+    GridCreator gridCreator_;
+  };
+
+} // end namespace Dune
diff --git a/dune/vtk/gridcreators/discontinuousgridcreator.hh b/dune/vtk/gridcreators/discontinuousgridcreator.hh
new file mode 100644
index 0000000..c392d4f
--- /dev/null
+++ b/dune/vtk/gridcreators/discontinuousgridcreator.hh
@@ -0,0 +1,97 @@
+#pragma once
+
+#include <cassert>
+#include <cstdint>
+#include <limits>
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/vtk/vtktypes.hh>
+#include <dune/vtk/gridcreatorinterface.hh>
+namespace Dune
+{
+  // Create a grid where the input points are not connected and the connectivity
+  // describes separated elements.
+  template <class Grid>
+  struct DiscontinuousGridCreator
+      : public GridCreatorInterface<Grid, DiscontinuousGridCreator<Grid>>
+  {
+    using Super = GridCreatorInterface<Grid, DiscontinuousGridCreator<Grid>>;
+    using GlobalCoordinate = typename Super::GlobalCoordinate;
+
+    struct CoordLess
+    {
+      template <class T, int N>
+      bool operator() (FieldVector<T,N> const& lhs, FieldVector<T,N> const& rhs) const
+      {
+        for (int i = 0; i < N; ++i) {
+          if (std::abs(lhs[i] - rhs[i]) < std::numeric_limits<T>::epsilon())
+            continue;
+          return lhs[i] < rhs[i];
+        }
+        return false;
+      }
+    };
+
+    DiscontinuousGridCreator (GridFactory<Grid>& factory)
+      : Super(factory)
+    {}
+
+    using Super::factory;
+    void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
+                             std::vector<std::uint64_t> const& /*point_ids*/)
+    {
+      points_ = &points;
+      uniquePoints_.clear();
+      std::size_t idx = 0;
+
+      for (auto const& p : points) {
+        auto b = uniquePoints_.emplace(std::make_pair(p,idx));
+        if (b.second) {
+          factory().insertVertex(p);
+          ++idx;
+        }
+      }
+    }
+
+    void insertElementsImpl (std::vector<std::uint8_t> const& types,
+                             std::vector<std::int64_t> const& offsets,
+                             std::vector<std::int64_t> const& connectivity)
+    {
+      assert(points_ != nullptr);
+      std::size_t idx = 0;
+      for (std::size_t i = 0; i < types.size(); ++i) {
+        auto type = Vtk::to_geometry(types[i]);
+        Vtk::CellType cellType{type};
+
+        int nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]);
+        assert(nNodes > 0);
+        std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes);
+        for (int j = 0; j < nNodes; ++j) {
+          std::size_t v_j = connectivity[idx++];
+          std::size_t new_idx = uniquePoints_[(*points_)[v_j]];
+          vtk_cell.push_back(new_idx);
+        }
+
+        if (cellType.noPermutation())
+          factory().insertElement(type,vtk_cell);
+        else {
+          // apply index permutation
+          std::vector<unsigned int> cell(nNodes);
+          for (int j = 0; j < nNodes; ++j)
+            cell[j] = vtk_cell[cellType.permutation(j)];
+
+          factory().insertElement(type,cell);
+        }
+      }
+    }
+
+  private:
+    std::vector<GlobalCoordinate> const* points_ = nullptr;
+    std::map<GlobalCoordinate, std::size_t, CoordLess> uniquePoints_;
+  };
+
+} // end namespace Dune
diff --git a/dune/vtk/gridcreators/parallelgridcreator.hh b/dune/vtk/gridcreators/parallelgridcreator.hh
new file mode 100644
index 0000000..4b69554
--- /dev/null
+++ b/dune/vtk/gridcreators/parallelgridcreator.hh
@@ -0,0 +1,50 @@
+#pragma once
+
+#include <cstdint>
+#include <string>
+#include <vector>
+
+#include <dune/grid/common/gridfactory.hh>
+#include <dune/vtk/forward.hh>
+#include <dune/vtk/gridcreatorinterface.hh>
+#include <dune/vtk/gridcreators/common.hh>
+#include <dune/vtk/gridcreators/derivedgridcreator.hh>
+#include <dune/vtk/gridcreators/continuousgridcreator.hh>
+
+namespace Dune
+{
+  template <class Grid>
+  struct ParallelGridCreator
+      : public DerivedGridCreator<ContinuousGridCreator<Grid>, ParallelGridCreator<Grid>>
+  {
+    using Self = ParallelGridCreator;
+    using Super = DerivedGridCreator<ContinuousGridCreator<Grid>, Self>;
+    using GlobalCoordinate = typename Super::GlobalCoordinate;
+    using VertexId = VertexId_t<GridFactory<Grid>>;
+
+    // The GridFactory must support insertion of global vertex IDs
+    static_assert(Std::is_detected<HasInsertVertex, GridFactory<Grid>, GlobalCoordinate, VertexId>{}, "");
+
+    ParallelGridCreator (GridFactory<Grid>& factory)
+      : Super(factory)
+    {}
+
+    using Super::factory;
+    void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
+                             std::vector<std::uint64_t> const& point_ids)
+    {
+      assert(point_ids.size() == points.size());
+      for (std::size_t i = 0; i < points.size(); ++i)
+        factory().insertVertex(points[i], VertexId(point_ids[i]));
+    }
+
+    void insertPiecesImpl (std::vector<std::string> const& pieces)
+    {
+      if (int(pieces.size()) == this->size()) {
+        VtkReader<Grid, Self> pieceReader(factory());
+        pieceReader.readFromFile(pieces[this->rank()], true);
+      }
+    }
+  };
+
+} // end namespace Dune
diff --git a/dune/vtk/gridcreators/serialgridcreator.hh b/dune/vtk/gridcreators/serialgridcreator.hh
new file mode 100644
index 0000000..c19d33f
--- /dev/null
+++ b/dune/vtk/gridcreators/serialgridcreator.hh
@@ -0,0 +1,71 @@
+#pragma once
+
+#include <cstdint>
+#include <string>
+#include <vector>
+
+#include <dune/grid/common/gridfactory.hh>
+#include <dune/vtk/forward.hh>
+#include <dune/vtk/gridcreatorinterface.hh>
+#include <dune/vtk/gridcreators/discontinuousgridcreator.hh>
+
+namespace Dune
+{
+  template <class Grid>
+  struct SerialGridCreator
+      : public GridCreatorInterface<Grid, SerialGridCreator<Grid>>
+  {
+    using Self = SerialGridCreator;
+    using Super = GridCreatorInterface<Grid, Self>;
+    using GlobalCoordinate = typename Super::GlobalCoordinate;
+
+    SerialGridCreator (GridFactory<Grid>& factory)
+      : Super(factory)
+    {}
+
+    void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
+                             std::vector<std::uint64_t> const& /*point_ids*/)
+    {
+      shift_.push_back(points_.size());
+      points_.reserve(points_.size() + points.size());
+      points_.insert(points_.end(), points.begin(), points.end());
+    }
+
+    void insertElementsImpl (std::vector<std::uint8_t> const& types,
+                             std::vector<std::int64_t> const& offsets,
+                             std::vector<std::int64_t> const& connectivity)
+    {
+      types_.reserve(types_.size() + types.size());
+      types_.insert(types_.end(), types.begin(), types.end());
+
+      offsets_.reserve(offsets_.size() + offsets.size());
+      std::transform(offsets.begin(), offsets.end(), std::back_inserter(offsets_),
+        [shift=offsets_.empty() ? 0 : offsets_.back()](std::int64_t o) { return o + shift; });
+
+      connectivity_.reserve(connectivity_.size() + connectivity.size());
+      std::transform(connectivity.begin(), connectivity.end(), std::back_inserter(connectivity_),
+        [shift=shift_.back()](std::int64_t idx) { return idx + shift; });
+    }
+
+    void insertPiecesImpl (std::vector<std::string> const& pieces)
+    {
+      if (this->rank() == 0) {
+        VtkReader<Grid, Self> pieceReader(*this);
+        for (std::string const& piece : pieces)
+          pieceReader.readFromFile(piece);
+
+        DiscontinuousGridCreator<Grid> creator(this->factory());
+        creator.insertVertices(points_, {});
+        creator.insertElements(types_, offsets_, connectivity_);
+      }
+    }
+
+  private:
+    std::vector<GlobalCoordinate> points_;
+    std::vector<std::uint8_t> types_;
+    std::vector<std::int64_t> offsets_;
+    std::vector<std::int64_t> connectivity_;
+    std::vector<std::int64_t> shift_;
+  };
+
+} // end namespace Dune
diff --git a/dune/vtk/pvdwriter.hh b/dune/vtk/pvdwriter.hh
index 46343bd..8fe864e 100644
--- a/dune/vtk/pvdwriter.hh
+++ b/dune/vtk/pvdwriter.hh
@@ -7,6 +7,7 @@
 
 #include <dune/vtk/vtktypes.hh>
 #include <dune/vtk/filewriter.hh>
+#include <dune/vtk/forward.hh>
 
 namespace Dune
 {
diff --git a/dune/vtk/vtkreader.hh b/dune/vtk/vtkreader.hh
index fb514bd..82a2b53 100644
--- a/dune/vtk/vtkreader.hh
+++ b/dune/vtk/vtkreader.hh
@@ -3,9 +3,12 @@
 #include <iosfwd>
 #include <map>
 
-#include "filereader.hh"
-#include "gridcreator.hh"
-#include "vtktypes.hh"
+#include <dune/vtk/filereader.hh>
+#include <dune/vtk/forward.hh>
+#include <dune/vtk/vtktypes.hh>
+
+// default GridCreator
+#include <dune/vtk/gridcreators/continuousgridcreator.hh>
 
 namespace Dune
 {
@@ -16,7 +19,7 @@ namespace Dune
    *
    * Assumption on the file structure: Each XML tag must be on a separate line.
    **/
-  template <class Grid, class GridCreator = DefaultGridCreator>
+  template <class Grid, class GridCreator>
   class VtkReader
       : public FileReader<Grid, VtkReader<Grid, GridCreator>>
   {
@@ -37,8 +40,17 @@ namespace Dune
     using GlobalCoordinate = typename Entity::Geometry::GlobalCoordinate;
 
   public:
-    /// Constructor. Stores a pointer to the GridFactory and rank and size of the communicator.
-    VtkReader (GridFactory<Grid>& factory);
+    /// Constructor. Creates a new GridCreator with the passed factory
+    template <class... Args>
+    explicit VtkReader (Args&&... args)
+      : creatorStorage_(std::make_unique<GridCreator>(std::forward<Args>(args)...))
+      , creator_(*creatorStorage_)
+    {}
+
+    /// Constructor. Stores a references to the passed creator
+    VtkReader (GridCreator& creator)
+      : creator_(creator)
+    {}
 
     /// Read the grid from file with `filename` into the GridFactory \ref factory_
     void readFromFile (std::string const& filename, bool create = true);
@@ -56,6 +68,9 @@ namespace Dune
       reader.readFromFile(filename);
     }
 
+    /// Construct a grid using the GridCreator
+    void createGrid();
+
     /// Return the filenames of parallel pieces
     std::vector<std::string> const& pieces () const
     {
@@ -111,30 +126,25 @@ namespace Dune
     // Read attributes from current xml tag
     std::map<std::string, std::string> parseXml(std::string const& line, bool& closed);
 
-    // Construct a grid using the GridFactory `factory` and the read vectors
-    // \ref vec_points, \ref vec_types, \ref vec_offsets, and \ref vec_connectivity,
-    // by passing to \ref GridCreator.
-    void createGrid() const;
+    // clear all vectors
+    void clear ();
 
   private:
-    GridFactory<Grid>* factory_ = nullptr;
-
-    // the rank and size of the factory collective communication
-    int rank_ = 0;
-    int numRanks_ = 1;
+    std::unique_ptr<GridCreator> creatorStorage_ = nullptr;
+    GridCreator& creator_;
 
     /// Data format, i.e. ASCII, BINARY or COMPRESSED. Read from xml attributes.
     Vtk::FormatTypes format_;
 
     // Temporary data to construct the grid elements
     std::vector<GlobalCoordinate> vec_points;
-    std::vector<std::uint64_t> vec_point_ids; //< Global unique vertex ID
-    std::vector<std::uint8_t> vec_types; //< VTK cell type ID
-    std::vector<std::int64_t> vec_offsets; //< offset of vertices of cell
+    std::vector<std::uint64_t> vec_point_ids;   //< Global unique vertex ID
+    std::vector<std::uint8_t> vec_types;        //< VTK cell type ID
+    std::vector<std::int64_t> vec_offsets;      //< offset of vertices of cell
     std::vector<std::int64_t> vec_connectivity; //< vertex indices of cell
 
-    std::size_t numberOfCells_; //< Number of cells in the grid
-    std::size_t numberOfPoints_; // Number of vertices in the grid
+    std::size_t numberOfCells_ = 0;   //< Number of cells in the grid
+    std::size_t numberOfPoints_ = 0;  //< Number of vertices in the grid
 
     // offset information for appended data
     // map Name -> {DataType,NumberOfComponents,Offset}
@@ -144,7 +154,7 @@ namespace Dune
     std::vector<std::string> pieces_;
 
     /// Offset of beginning of appended data
-    std::uint64_t offset0_;
+    std::uint64_t offset0_ = 0;
   };
 
 } // end namespace Dune
diff --git a/dune/vtk/vtkreader.impl.hh b/dune/vtk/vtkreader.impl.hh
index e94905f..0f8099a 100644
--- a/dune/vtk/vtkreader.impl.hh
+++ b/dune/vtk/vtkreader.impl.hh
@@ -15,22 +15,6 @@
 
 namespace Dune {
 
-template <class Grid, class Creator>
-VtkReader<Grid,Creator>::VtkReader (GridFactory<Grid>& factory)
-  : factory_(&factory)
-{
-#if DUNE_VERSION_LT(DUNE_GRID,2,7)
-  #if HAVE_MPI
-  MPI_Comm_rank(MPI_COMM_WORLD, &rank_);
-  MPI_Comm_size(MPI_COMM_WORLD, &numRanks_);
-  #endif
-#else
-  rank_ = factory.comm().rank();
-  numRanks_ = factory.comm().size();
-#endif
-}
-
-
 template <class Grid, class Creator>
 void VtkReader<Grid,Creator>::readFromFile (std::string const& filename, bool create)
 {
@@ -45,7 +29,7 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename, bool cr
   if (ext == ".vtu") {
     readSerialFileFromStream(input, create);
   } else if (ext == ".pvtu") {
-    readParallelFileFromStream(input, rank_, numRanks_, create);
+    readParallelFileFromStream(input, creator_.rank(), creator_.size(), create);
   } else {
     DUNE_THROW(IOError, "File has unknown file-extension '" << ext << "'. Allowed are only '.vtu' and '.pvtu'.");
   }
@@ -55,6 +39,7 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename, bool cr
 template <class Grid, class Creator>
 void VtkReader<Grid,Creator>::readSerialFileFromStream (std::ifstream& input, bool create)
 {
+  clear();
   std::string compressor = "";
   std::string data_name = "", data_format = "";
   Vtk::DataTypes data_type = Vtk::UNKNOWN;
@@ -247,7 +232,7 @@ void VtkReader<Grid,Creator>::readSerialFileFromStream (std::ifstream& input, bo
 template <class Grid, class Creator>
 void VtkReader<Grid,Creator>::readParallelFileFromStream (std::ifstream& input, int commRank, int commSize, bool create)
 {
-  pieces_.clear();
+  clear();
 
   Sections section = NO_SECTION;
   for (std::string line; std::getline(input, line); ) {
@@ -290,9 +275,6 @@ void VtkReader<Grid,Creator>::readParallelFileFromStream (std::ifstream& input,
   if (section != NO_SECTION)
     DUNE_THROW(IOError, "VTK-File is incomplete. It must end with </VTKFile>!");
 
-  assert(pieces_.size() == commSize);
-  readFromFile(pieces_[commRank], false);
-
   if (create)
     createGrid();
 }
@@ -528,13 +510,17 @@ void VtkReader<Grid,Creator>::readAppended (std::ifstream& input, std::vector<T>
 
 
 template <class Grid, class Creator>
-void VtkReader<Grid,Creator>::createGrid () const
+void VtkReader<Grid,Creator>::createGrid ()
 {
   assert(vec_points.size() == numberOfPoints_);
   assert(vec_types.size() == numberOfCells_);
   assert(vec_offsets.size() == numberOfCells_);
 
-  Creator::create(*factory_, vec_points, vec_point_ids, vec_types, vec_offsets, vec_connectivity);
+  if (!vec_points.empty())
+    creator_.insertVertices(vec_points, vec_point_ids);
+  if (!vec_types.empty())
+    creator_.insertElements(vec_types, vec_offsets, vec_connectivity);
+  creator_.insertPieces(pieces_);
 }
 
 // Assume input already read the line <AppendedData ...>
@@ -605,4 +591,21 @@ std::map<std::string, std::string> VtkReader<Grid,Creator>::parseXml (std::strin
   return attr;
 }
 
+
+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;
+}
+
 } // end namespace Dune
diff --git a/dune/vtk/vtktimeserieswriter.hh b/dune/vtk/vtktimeserieswriter.hh
index 5874677..a2b022e 100644
--- a/dune/vtk/vtktimeserieswriter.hh
+++ b/dune/vtk/vtktimeserieswriter.hh
@@ -5,6 +5,7 @@
 #include <vector>
 
 #include <dune/vtk/filewriter.hh>
+#include <dune/vtk/forward.hh>
 #include <dune/vtk/vtktypes.hh>
 #include <dune/vtk/utility/filesystem.hh>
 #include <dune/vtk/utility/uid.hh>
diff --git a/dune/vtk/vtkwriter.hh b/dune/vtk/vtkwriter.hh
index 91befcf..ff44145 100644
--- a/dune/vtk/vtkwriter.hh
+++ b/dune/vtk/vtkwriter.hh
@@ -12,6 +12,7 @@
 
 #include <dune/grid/geometrygrid.hh>
 #include <dune/grid/yaspgrid.hh>
+#include <dune/vtk/forward.hh>
 #include <dune/vtk/datacollectors/yaspdatacollector.hh>
 
 namespace Dune
diff --git a/dune/vtk/vtkwriterinterface.hh b/dune/vtk/vtkwriterinterface.hh
index 674bde1..7f6ea25 100644
--- a/dune/vtk/vtkwriterinterface.hh
+++ b/dune/vtk/vtkwriterinterface.hh
@@ -7,18 +7,12 @@
 
 #include <dune/common/std/optional.hh>
 #include <dune/vtk/filewriter.hh>
+#include <dune/vtk/forward.hh>
 #include <dune/vtk/vtkfunction.hh>
 #include <dune/vtk/vtktypes.hh>
 
 namespace Dune
 {
-  // forward declaration
-  template <class VtkWriter>
-  class VtkTimeseriesWriter;
-
-  template <class VtkWriter>
-  class PvdWriter;
-
   /// Interface for file writers for the Vtk XML file formats
   template <class GridView, class DataCollector>
   class VtkWriterInterface
diff --git a/dune/vtk/writers/vtkimagedatawriter.hh b/dune/vtk/writers/vtkimagedatawriter.hh
index 50db5b7..e5260e3 100644
--- a/dune/vtk/writers/vtkimagedatawriter.hh
+++ b/dune/vtk/writers/vtkimagedatawriter.hh
@@ -5,6 +5,7 @@
 #include <map>
 
 #include <dune/vtk/filewriter.hh>
+#include <dune/vtk/forward.hh>
 #include <dune/vtk/vtkfunction.hh>
 #include <dune/vtk/vtktypes.hh>
 #include <dune/vtk/datacollectors/structureddatacollector.hh>
@@ -18,7 +19,7 @@ namespace Dune
    * Requirement:
    * - DataCollector must be a model of \ref StructuredDataCollector
    **/
-  template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+  template <class GridView, class DataCollector>
   class VtkImageDataWriter
       : public VtkWriterInterface<GridView, DataCollector>
   {
diff --git a/dune/vtk/writers/vtkrectilineargridwriter.hh b/dune/vtk/writers/vtkrectilineargridwriter.hh
index f7c0317..5bb2da2 100644
--- a/dune/vtk/writers/vtkrectilineargridwriter.hh
+++ b/dune/vtk/writers/vtkrectilineargridwriter.hh
@@ -5,6 +5,7 @@
 #include <map>
 
 #include <dune/vtk/filewriter.hh>
+#include <dune/vtk/forward.hh>
 #include <dune/vtk/vtkfunction.hh>
 #include <dune/vtk/vtktypes.hh>
 #include <dune/vtk/datacollectors/structureddatacollector.hh>
@@ -18,7 +19,7 @@ namespace Dune
    * Requirement:
    * - DataCollector must be a model of \ref StructuredDataCollector
    **/
-  template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+  template <class GridView, class DataCollector>
   class VtkRectilinearGridWriter
       : public VtkWriterInterface<GridView, DataCollector>
   {
diff --git a/dune/vtk/writers/vtkstructuredgridwriter.hh b/dune/vtk/writers/vtkstructuredgridwriter.hh
index 522d472..d5c240f 100644
--- a/dune/vtk/writers/vtkstructuredgridwriter.hh
+++ b/dune/vtk/writers/vtkstructuredgridwriter.hh
@@ -5,6 +5,7 @@
 #include <map>
 
 #include <dune/vtk/filewriter.hh>
+#include <dune/vtk/forward.hh>
 #include <dune/vtk/vtkfunction.hh>
 #include <dune/vtk/vtktypes.hh>
 #include <dune/vtk/datacollectors/structureddatacollector.hh>
@@ -18,7 +19,7 @@ namespace Dune
    * Requirement:
    * - DataCollector must be a model of \ref StructuredDataCollector
    **/
-  template <class GridView, class DataCollector = StructuredDataCollector<GridView>>
+  template <class GridView, class DataCollector>
   class VtkStructuredGridWriter
       : public VtkWriterInterface<GridView, DataCollector>
   {
diff --git a/dune/vtk/writers/vtkunstructuredgridwriter.hh b/dune/vtk/writers/vtkunstructuredgridwriter.hh
index b37399a..ee88bf7 100644
--- a/dune/vtk/writers/vtkunstructuredgridwriter.hh
+++ b/dune/vtk/writers/vtkunstructuredgridwriter.hh
@@ -5,6 +5,7 @@
 #include <map>
 
 #include <dune/vtk/filewriter.hh>
+#include <dune/vtk/forward.hh>
 #include <dune/vtk/vtkfunction.hh>
 #include <dune/vtk/vtktypes.hh>
 #include <dune/vtk/datacollectors/continuousdatacollector.hh>
@@ -18,7 +19,7 @@ namespace Dune
    * Requirement:
    * - DataCollector must be a model of \ref DataCollector
    **/
-  template <class GridView, class DataCollector = ContinuousDataCollector<GridView>>
+  template <class GridView, class DataCollector>
   class VtkUnstructuredGridWriter
       : public VtkWriterInterface<GridView, DataCollector>
   {
diff --git a/src/test/CMakeLists.txt b/src/test/CMakeLists.txt
index 77409fe..86a7f5f 100644
--- a/src/test/CMakeLists.txt
+++ b/src/test/CMakeLists.txt
@@ -1,6 +1,9 @@
 dune_add_test(SOURCES reader_writer_test.cc
               LINK_LIBRARIES dunevtk)
 
+dune_add_test(SOURCES parallel_reader_writer_test.cc
+              LINK_LIBRARIES dunevtk)
+
 dune_add_test(SOURCES mixed_element_test.cc
               LINK_LIBRARIES dunevtk
               CMAKE_GUARD HAVE_UG)
\ No newline at end of file
diff --git a/src/test/parallel_reader_writer_test.cc b/src/test/parallel_reader_writer_test.cc
new file mode 100644
index 0000000..cb126ba
--- /dev/null
+++ b/src/test/parallel_reader_writer_test.cc
@@ -0,0 +1,185 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifdef HAVE_CONFIG_H
+# include "config.h"
+#endif
+
+#include <cstring>
+#include <iostream>
+#include <vector>
+
+#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
+#include <dune/common/filledarray.hh>
+#include <dune/common/test/testsuite.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#if HAVE_DUNE_ALUGRID
+#include <dune/alugrid/grid.hh>
+#endif
+
+#include <dune/vtk/vtkreader.hh>
+#include <dune/vtk/writers/vtkunstructuredgridwriter.hh>
+#include <dune/vtk/gridcreators/parallelgridcreator.hh>
+#include <dune/vtk/gridcreators/serialgridcreator.hh>
+
+using namespace Dune;
+
+// see https://stackoverflow.com/questions/6163611/compare-two-files
+bool compare_files (std::string const& fn1, std::string const& fn2)
+{
+  std::ifstream in1(fn1, std::ios::binary);
+  std::ifstream in2(fn2, std::ios::binary);
+  if (!in1 || !in2) {
+    std::cout << "can not find file " << fn1 << " or file " << fn2 << "\n";
+    return false;
+  }
+
+  std::ifstream::pos_type size1 = in1.seekg(0, std::ifstream::end).tellg();
+  in1.seekg(0, std::ifstream::beg);
+
+  std::ifstream::pos_type size2 = in2.seekg(0, std::ifstream::end).tellg();
+  in2.seekg(0, std::ifstream::beg);
+
+  if (size1 != size2)
+    return false;
+
+  static const std::size_t BLOCKSIZE = 4096;
+  std::size_t remaining = size1;
+
+  while (remaining) {
+    char buffer1[BLOCKSIZE], buffer2[BLOCKSIZE];
+    std::size_t size = std::min(BLOCKSIZE, remaining);
+
+    in1.read(buffer1, size);
+    in2.read(buffer2, size);
+
+    if (0 != std::memcmp(buffer1, buffer2, size))
+      return false;
+
+    remaining -= size;
+  }
+
+  return true;
+}
+
+using TestCase = std::tuple<std::string,Vtk::FormatTypes,Vtk::DataTypes>;
+using TestCases = std::set<TestCase>;
+static TestCases test_cases = {
+  {"ascii32", Vtk::ASCII, Vtk::FLOAT32},
+  {"bin32", Vtk::BINARY, Vtk::FLOAT32},
+  {"zlib32", Vtk::COMPRESSED, Vtk::FLOAT32},
+  {"ascii64", Vtk::ASCII, Vtk::FLOAT64},
+  {"bin64", Vtk::BINARY, Vtk::FLOAT64},
+  {"zlib64", Vtk::COMPRESSED, Vtk::FLOAT64},
+};
+
+template <class Test>
+void compare (Test& test, filesystem::path const& dir, filesystem::path const& name)
+{
+  test.check(compare_files(dir.string() + '/' + name.string() + ".vtu",
+                           dir.string() + '/' + name.string() + "_2.vtu"));
+}
+
+template <class GridView>
+void writer_test (GridView const& gridView)
+{
+  for (auto const& test_case : test_cases) {
+    VtkUnstructuredGridWriter<GridView> vtkWriter(gridView, std::get<1>(test_case), std::get<2>(test_case));
+    vtkWriter.write("parallel_reader_writer_test_" + std::get<0>(test_case) + ".vtu");
+  }
+}
+
+template <class G> struct IsALUGrid : std::false_type {};
+#if DUNE_VERSION_GT(DUNE_GRID,2,6) && HAVE_DUNE_ALUGRID
+template<int dim, int dimworld, Dune::ALUGridElementType elType, Dune::ALUGridRefinementType refineType, class Comm>
+struct IsALUGrid<Dune::ALUGrid<dim,dimworld,elType,refineType,Comm>> : std::true_type {};
+#endif
+
+template <class Grid, class Creator>
+void reader_test (MPIHelper& mpi, TestSuite& test)
+{
+  std::string ext = ".vtu";
+  if (mpi.size() > 1)
+    ext = ".pvtu";
+
+  TestCase test_case = {"ascii32", Vtk::ASCII, Vtk::FLOAT32};
+
+  GridFactory<Grid> factory;
+  VtkReader<Grid, Creator> reader{factory};
+  reader.readFromFile("parallel_reader_writer_test_" + std::get<0>(test_case) + ext);
+
+  std::unique_ptr<Grid> grid = Hybrid::ifElse(IsALUGrid<Grid>{},
+    [&](auto id) { return id(factory).createGrid(std::true_type{}); },
+    [&](auto id) { return id(factory).createGrid(); });
+  std::vector<std::string> pieces1 = grid->comm().size() > 1 ?
+    reader.pieces() :
+    std::vector<std::string>{"parallel_reader_writer_test_" + std::get<0>(test_case) + ".vtu"};
+
+  VtkUnstructuredGridWriter<typename Grid::LeafGridView> vtkWriter(grid->leafGridView(),
+    std::get<1>(test_case), std::get<2>(test_case));
+  vtkWriter.write("parallel_reader_writer_test_" + std::get<0>(test_case) + "_2.vtu");
+
+  GridFactory<Grid> factory2;
+  VtkReader<Grid, Creator> reader2{factory2};
+  reader2.readFromFile("parallel_reader_writer_test_" + std::get<0>(test_case) + "_2" + ext, false);
+  std::vector<std::string> pieces2 = grid->comm().size() > 1 ?
+    reader.pieces() :
+    std::vector<std::string>{"parallel_reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"};
+
+  test.check(pieces1.size() == pieces2.size(), "pieces1.size == pieces2.size");
+  for (std::size_t i = 0; i < pieces1.size(); ++i)
+    test.check(compare_files(pieces1[i], pieces2[i]));
+}
+
+
+template <class Grid, class Creator>
+void reader_writer_test(MPIHelper& mpi, TestSuite& test, std::string const& testName)
+{
+  std::cout << "== " << testName << "\n";
+  const int dim = Grid::dimension;
+  FieldVector<double,dim> lowerLeft; lowerLeft = 0.0;
+  FieldVector<double,dim> upperRight; upperRight = 1.0;
+  auto numElements = filledArray<dim,unsigned int>(4);
+  auto gridPtr = StructuredGridFactory<Grid>::createSimplexGrid(lowerLeft, upperRight, numElements);
+  gridPtr->loadBalance();
+
+  writer_test(gridPtr->leafGridView());
+  MPI_Barrier(MPI_COMM_WORLD);
+  reader_test<Grid, Creator>(mpi,test);
+}
+
+#if HAVE_DUNE_ALUGRID
+  template <int dim>
+  using ALUGridType = Dune::ALUGrid<dim, dim, Dune::simplex, Dune::conforming>;
+#endif
+
+template <int I>
+using int_ = std::integral_constant<int,I>;
+
+int main (int argc, char** argv)
+{
+  auto& mpi = Dune::MPIHelper::instance(argc, argv);
+
+  TestSuite test{};
+
+#if HAVE_UG
+  reader_writer_test<UGGrid<2>, SerialGridCreator<UGGrid<2>>>(mpi, test, "UGGrid<2>");
+  reader_writer_test<UGGrid<3>, SerialGridCreator<UGGrid<3>>>(mpi, test, "UGGrid<3>");
+#endif
+
+#if HAVE_DUNE_ALUGRID
+  reader_writer_test<ALUGridType<2>, SerialGridCreator<ALUGridType<2>>>(mpi, test, "ALUGridType<2>");
+  reader_writer_test<ALUGridType<3>, SerialGridCreator<ALUGridType<3>>>(mpi, test, "ALUGridType<3>");
+
+// #if DUNE_VERSION_LT(DUNE_GRID,2,7)
+  reader_writer_test<ALUGridType<2>, ParallelGridCreator<ALUGridType<2>>>(mpi, test, "ALUGridType<2, Parallel>");
+  reader_writer_test<ALUGridType<3>, ParallelGridCreator<ALUGridType<3>>>(mpi, test, "ALUGridType<3, Parallel>");
+// #endif
+#endif
+
+  return test.exit();
+}
diff --git a/src/test/reader_writer_test.cc b/src/test/reader_writer_test.cc
index 840d9c3..860ff56 100644
--- a/src/test/reader_writer_test.cc
+++ b/src/test/reader_writer_test.cc
@@ -23,6 +23,7 @@
 
 #include <dune/vtk/vtkreader.hh>
 #include <dune/vtk/writers/vtkunstructuredgridwriter.hh>
+#include <dune/vtk/gridcreators/continuousgridcreator.hh>
 
 using namespace Dune;
 
@@ -90,18 +91,30 @@ void writer_test (GridView const& gridView)
   }
 }
 
+template <class G> struct IsALUGrid : std::false_type {};
+#if DUNE_VERSION_GT(DUNE_GRID,2,6) && HAVE_DUNE_ALUGRID
+template<int dim, int dimworld, Dune::ALUGridElementType elType, Dune::ALUGridRefinementType refineType, class Comm>
+struct IsALUGrid<Dune::ALUGrid<dim,dimworld,elType,refineType,Comm>> : std::true_type {};
+#endif
+
 template <class Grid, class Test>
-void reader_test (Test& test)
+void reader_test (MPIHelper& mpi, Test& test)
 {
   std::string ext = ".vtu";
+  if (mpi.size() > 1)
+    ext = ".pvtu";
 
   for (auto const& test_case : test_cases) {
     GridFactory<Grid> factory;
     VtkReader<Grid> reader{factory};
     reader.readFromFile("reader_writer_test_" + std::get<0>(test_case) + ext);
-    std::unique_ptr<Grid> grid(factory.createGrid());
-    // auto pieces1 = reader.pieces();
-    std::vector<std::string> pieces1 = {"reader_writer_test_" + std::get<0>(test_case) + ".vtu"};
+
+    std::unique_ptr<Grid> grid = Hybrid::ifElse(IsALUGrid<Grid>{},
+      [&](auto id) { return id(factory).createGrid(std::true_type{}); },
+      [&](auto id) { return id(factory).createGrid(); });
+    std::vector<std::string> pieces1 = grid->comm().size() > 1 ?
+      reader.pieces() :
+      std::vector<std::string>{"reader_writer_test_" + std::get<0>(test_case) + ".vtu"};
 
     VtkUnstructuredGridWriter<typename Grid::LeafGridView> vtkWriter(grid->leafGridView(),
       std::get<1>(test_case), std::get<2>(test_case));
@@ -110,8 +123,9 @@ void reader_test (Test& test)
     GridFactory<Grid> factory2;
     VtkReader<Grid> reader2{factory2};
     reader2.readFromFile("reader_writer_test_" + std::get<0>(test_case) + "_2" + ext, false);
-    // auto pieces2 = reader.pieces();
-    std::vector<std::string> pieces2 = {"reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"};
+    std::vector<std::string> pieces2 = grid->comm().size() > 1 ?
+      reader.pieces() :
+      std::vector<std::string>{"reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"};
 
     test.check(pieces1.size() == pieces2.size(), "pieces1.size == pieces2.size");
     for (std::size_t i = 0; i < pieces1.size(); ++i)
@@ -135,7 +149,8 @@ int main (int argc, char** argv)
 
 #if HAVE_UG
   // Test VtkWriter for UGGrid
-  Hybrid::forEach(std::make_tuple(int_<2>{}, int_<3>{}), [&test](auto dim)
+  if (mpi.size() == 1) {
+  Hybrid::forEach(std::make_tuple(int_<2>{}, int_<3>{}), [&test,&mpi](auto dim)
   {
     using GridType = UGGrid<dim.value>;
     {
@@ -148,13 +163,15 @@ int main (int argc, char** argv)
       writer_test(gridPtr->leafGridView());
     }
 
-    reader_test<GridType>(test);
+    reader_test<GridType>(mpi,test);
   });
+  }
 #endif
 
-#if DUNE_VERSION_LT(DUNE_GRID,2,7) && HAVE_DUNE_ALUGRID
+// DUNE_VERSION_LT(DUNE_GRID,2,7) &&
+#if  HAVE_DUNE_ALUGRID
   // Test VtkWriter for ALUGrid. Currently the 2.7 branch is not working.
-  Hybrid::forEach(std::make_tuple(int_<2>{}, int_<3>{}), [&test](auto dim)
+  Hybrid::forEach(std::make_tuple(int_<2>{}, int_<3>{}), [&test,&mpi](auto dim)
   {
     using GridType = Dune::ALUGrid<dim.value, dim.value, Dune::simplex, Dune::conforming>;
     {
@@ -167,7 +184,7 @@ int main (int argc, char** argv)
       writer_test(gridPtr->leafGridView());
     }
 
-    reader_test<GridType>(test);
+    reader_test<GridType>(mpi,test);
   });
 #endif
 
diff --git a/src/vtkreader.cc b/src/vtkreader.cc
index d82e96b..277e47a 100644
--- a/src/vtkreader.cc
+++ b/src/vtkreader.cc
@@ -16,6 +16,8 @@
 #include <dune/grid/utility/structuredgridfactory.hh>
 
 #include <dune/vtk/vtkreader.hh>
+#include <dune/vtk/gridcreators/continuousgridcreator.hh>
+#include <dune/vtk/gridcreators/discontinuousgridcreator.hh>
 #include <dune/vtk/writers/vtkunstructuredgridwriter.hh>
 
 using namespace Dune;
@@ -72,7 +74,7 @@ int main(int argc, char** argv)
   }
 
   {
-    auto gridPtr = VtkReader<GridType,ConnectedGridCreator>::read("test_ascii_float32.vtu");
+    auto gridPtr = VtkReader<GridType,DiscontinuousGridCreator<GridType>>::read("test_ascii_float32.vtu");
     auto& grid = *gridPtr;
 
     VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::ASCII);
-- 
GitLab