parallel vtk reader partially implemented

.gitlab-ci.yml

@@ -3,18 +3,26 @@ before_script:
   - duneci-install-module https://gitlab.dune-project.org/martin.nolte/dune-polygongrid.git
   - duneci-install-module https://gitlab.dune-project.org/extensions/dune-spgrid.git
 
-debian:10 gcc-8-17:
-  image: registry.dune-project.org/docker/ci/dune:2.6-debian-10-gcc-8-17
+dune-2.6 gcc-6-14:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-9-gcc-6-14
   script: duneci-standard-test
 
-debian:10 clang-6-libcpp-17:
-  image: registry.dune-project.org/docker/ci/dune:2.6-debian-10-clang-6-libcpp-17
+dune-2.6 gcc-8-17:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-10-gcc-8-17
   script: duneci-standard-test
 
-debian:9 gcc-6-14:
-  image: registry.dune-project.org/docker/ci/dune:2.6-debian-9-gcc-6-14
+dune-2.6 clang-6-17:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:2.6-debian-10-clang-6-libcpp-17
   script: duneci-standard-test
 
-ubuntu:18.04 clang-6-17:
-  image: registry.dune-project.org/docker/ci/dune:2.6-ubuntu-18.04-clang-6-17
+dune-git gcc-7-14:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-gcc-7-14
   script: duneci-standard-test
+
+dune-git gcc-8-17:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-gcc-8-noassert-17
+  script: duneci-standard-test
+
+dune-git clang-6-17:
+  image: registry.dune-project.org/docker/ci/dune-pdelab-deps:git-debian-10-clang-6-libcpp-17
+  script: duneci-standard-test
\ No newline at end of file

dune.module

@@ -4,7 +4,7 @@
 
 #Name of the module
 Module: dune-vtk
-Version: 0.1
+Version: 0.2
 Maintainer: simon.praetorius@tu-dresden.de
 #depending on
 Depends: dune-common (>= 2.6) dune-geometry (>= 2.6) dune-grid

dune/vtk/datacollectorinterface.hh

@@ -19,7 +19,7 @@ public:
     asDerived().updateImpl();
   }
 
-  /// Return the number of overlapping elements
+  /// Return the number of ghost elements
   int ghostLevel () const
   {
     return asDerived().ghostLevelImpl();

dune/vtk/datacollectorinterface.impl.hh

@@ -2,6 +2,7 @@
 
 #include <dune/geometry/referenceelements.hh>
 #include <dune/grid/common/gridenums.hh>
+#include <dune/grid/common/mcmgmapper.hh>
 
 namespace Dune {
 
@@ -11,12 +12,12 @@ std::vector<T> DataCollectorInterface<GridView, Derived>
   ::cellDataImpl (VtkFunction const& fct) const
 {
   std::vector<T> data(gridView_.size(0) * fct.ncomps());
-  auto const& indexSet = gridView_.indexSet();
+  MultipleCodimMultipleGeomTypeMapper<GridView> mapper(gridView_, mcmgElementLayout());
   auto localFct = localFunction(fct);
   for (auto const& e : elements(gridView_, Partitions::all)) {
     localFct.bind(e);
     auto refElem = referenceElement<T,GridView::dimension>(e.type());
-    std::size_t idx = fct.ncomps() * indexSet.index(e);
+    std::size_t idx = fct.ncomps() * mapper.index(e);
     for (int comp = 0; comp < fct.ncomps(); ++comp)
       data[idx + comp] = T(localFct.evaluate(comp, refElem.position(0,0)));
     localFct.unbind();

dune/vtk/datacollectors/continuousdatacollector.hh

@@ -3,6 +3,8 @@
 #include <numeric>
 #include "unstructureddatacollector.hh"
 
+#include <dune/grid/utility/globalindexset.hh>
+
 namespace Dune
 {
 
@@ -44,6 +46,18 @@ public:
     return data;
   }
 
+  /// Return a vector of global unique ids of the points
+  std::vector<std::uint64_t> pointIdsImpl () const
+  {
+    std::vector<std::uint64_t> data(gridView_.size(dim));
+    GlobalIndexSet<GridView> globalIndexSet(gridView_, dim);
+    auto const& indexSet = gridView_.indexSet();
+    for (auto const& vertex : vertices(gridView_, Partitions::all)) {
+      data[indexSet.index(vertex)] = std::uint64_t(globalIndexSet.index(vertex));
+    }
+    return data;
+  }
+
   /// Return number of grid cells
   std::uint64_t numCellsImpl () const
   {

dune/vtk/datacollectors/structureddatacollector.hh

@@ -101,15 +101,13 @@ public: // default implementation:
     subDataCollector_.update();
 
 #if HAVE_MPI
-    int rank = -1;
-    int num_ranks = -1;
-    MPI_Comm_rank(gridView_.comm(), &rank);
-    MPI_Comm_size(gridView_.comm(), &num_ranks);
+    int rank = gridView_.comm().rank();
+    int numRanks = gridView_.comm().size();
 
     if (rank == 0) {
-      extents_.resize(num_ranks);
-      requests_.resize(num_ranks, MPI_REQUEST_NULL);
-      for (int i = 1; i < num_ranks; ++i)
+      extents_.resize(numRanks);
+      requests_.resize(numRanks, MPI_REQUEST_NULL);
+      for (int i = 1; i < numRanks; ++i)
         MPI_Irecv(extents_[i].data(), extents_[i].size(), MPI_INT, i, /*tag=*/6, gridView_.comm(), &requests_[i]);
     }
 
@@ -149,8 +147,7 @@ public: // default implementation:
     MPI_Test(&sendRequest_, &sendFlag, &sendStatus);
 
     if (sendFlag) {
-      int rank = -1;
-      MPI_Comm_rank(gridView_.comm(), &rank);
+      int rank = gridView_.comm().rank();
       if (rank != 0) {
         MPI_Isend(extent.data(), extent.size(), MPI_INT, 0, /*tag=*/6, gridView_.comm(), &sendRequest_);
       } else {
@@ -169,12 +166,11 @@ public: // default implementation:
 #if HAVE_MPI
     writer(0, extents_[0], true);
 
-    int num_ranks = -1;
-    MPI_Comm_size(gridView_.comm(), &num_ranks);
-    for (int p = 1; p < num_ranks; ++p) {
+    int numRanks = gridView_.comm().size();
+    for (int p = 1; p < numRanks; ++p) {
       int idx = -1;
       MPI_Status status;
-      MPI_Waitany(num_ranks, requests_.data(), &idx, &status);
+      MPI_Waitany(numRanks, requests_.data(), &idx, &status);
       if (idx != MPI_UNDEFINED) {
         assert(idx == status.MPI_SOURCE && status.MPI_TAG == 6);
         writer(idx, extents_[idx], true);

dune/vtk/datacollectors/unstructureddatacollector.hh

@@ -37,6 +37,18 @@ public:
     return this->asDerived().cellsImpl();
   }
 
+  std::vector<std::uint64_t> pointIds () const
+  {
+    return this->asDerived().pointIdsImpl();
+  }
+
+protected:
+  // default implementation
+  std::vector<std::uint64_t> pointIdsImpl () const
+  {
+    return {};
+  }
+
 protected:
   using Super::gridView_;
 };

dune/vtk/gridcreator.hh

@@ -6,12 +6,29 @@
 #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
@@ -19,12 +36,23 @@ namespace Dune
     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)
     {
-      for (auto const& p : points)
-        factory.insertVertex(p);
+      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) {
@@ -74,20 +102,38 @@ namespace Dune
     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;
-      for (auto const& p : points) {
-        auto b = unique_points.emplace(std::make_pair(p,idx));
-        if (b.second) {
-          factory.insertVertex(p);
-          ++idx;
+
+      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]);

dune/vtk/vtkreader.hh

@@ -9,7 +9,7 @@
 
 namespace Dune
 {
-  /// File-Reader for Vtk .vtu files
+  /// File-Reader for Vtk unstructured .vtu files
   /**
    * Reads .vtu files and constructs a grid from the cells stored in the file
    * Additionally, stored data can be read.
@@ -37,13 +37,17 @@ namespace Dune
     using GlobalCoordinate = typename Entity::Geometry::GlobalCoordinate;
 
   public:
-    /// Constructor. Stores a pointer to the GridFactory.
-    VtkReader (GridFactory<Grid>& factory)
-      : factory_(&factory)
-    {}
+    /// Constructor. Stores a pointer to the GridFactory and rank and size of the communicator.
+    VtkReader (GridFactory<Grid>& factory);
 
-    /// Read the grid from file with `filename` into the GridFactory `factory`
-    void readFromFile (std::string const& filename);
+    /// Read the grid from file with `filename` into the GridFactory \ref factory_
+    void readFromFile (std::string const& filename, bool create = true);
+
+    /// Read the grid from and input stream into the GridFactory \ref factory_
+    void readSerialFileFromStream (std::ifstream& input, bool create = true);
+
+    /// Read the grid from and input stream into the GridFactory \ref factory_
+    void readParallelFileFromStream (std::ifstream& input, int rank, int size, bool create = true);
 
     /// Implementation of \ref FileReader interface
     static void readFactoryImpl (GridFactory<Grid>& factory, std::string const& filename)
@@ -52,6 +56,12 @@ namespace Dune
       reader.readFromFile(filename);
     }
 
+    /// Return the filenames of parallel pieces
+    std::vector<std::string> const& pieces () const
+    {
+      return pieces_;
+    }
+
   private:
     // Read values stored on the cells with name `name`
     template <class T>
@@ -69,7 +79,7 @@ namespace Dune
     }
 
     // Read vertex coordinates from `input` stream and store in into `factory`
-    Sections readPoints (std::ifstream& input);
+    Sections readPoints (std::ifstream& input, std::string name);
 
     template <class T>
     void readPointsAppended (std::ifstream& input);
@@ -107,13 +117,18 @@ namespace Dune
     void createGrid() const;
 
   private:
-    GridFactory<Grid>* factory_;
+    GridFactory<Grid>* factory_ = nullptr;
+
+    // the rank and size of the factory collective communication
+    int rank_ = 0;
+    int numRanks_ = 1;
 
     /// 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::int64_t> vec_connectivity; //< vertex indices of cell
@@ -125,6 +140,9 @@ namespace Dune
     // map Name -> {DataType,NumberOfComponents,Offset}
     std::map<std::string, DataArrayAttributes> dataArray_;
 
+    // vector of filenames of parallel pieces
+    std::vector<std::string> pieces_;
+
     /// Offset of beginning of appended data
     std::uint64_t offset0_;
   };

dune/vtk/vtkreader.impl.hh

@@ -3,11 +3,16 @@
 #include <iterator>
 #include <string>
 
+<<<<<<< 0d5e1c63fbeb492af8b00cdfbddbcd516e2805af
 #ifdef HAVE_VTK_ZLIB
+=======
+#if HAVE_VTK_ZLIB
+>>>>>>> parallel vtk reader partially implemented
 #include <zlib.h>
 #endif
 
 #include <dune/common/classname.hh>
+#include <dune/common/version.hh>
 
 #include "utility/filesystem.hh"
 #include "utility/string.hh"
@@ -15,7 +20,23 @@
 namespace Dune {
 
 template <class Grid, class Creator>
-void VtkReader<Grid,Creator>::readFromFile (std::string const& filename)
+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)
 {
   // check whether file exists!
   if (!filesystem::exists(filename))
@@ -24,6 +45,20 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename)
   std::ifstream input(filename, std::ios_base::in | std::ios_base::binary);
   assert(input.is_open());
 
+  std::string ext = filesystem::path(filename).extension().string();
+  if (ext == ".vtu") {
+    readSerialFileFromStream(input, create);
+  } else if (ext == ".pvtu") {
+    readParallelFileFromStream(input, rank_, numRanks_, create);
+  } else {
+    DUNE_THROW(IOError, "File has unknown file-extension '" << ext << "'. Allowed are only '.vtu' and '.pvtu'.");
+  }
+}
+
+
+template <class Grid, class Creator>
+void VtkReader<Grid,Creator>::readSerialFileFromStream (std::ifstream& input, bool create)
+{
   std::string compressor = "";
   std::string data_name = "", data_format = "";
   Vtk::DataTypes data_type = Vtk::UNKNOWN;
@@ -50,9 +85,6 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename)
         compressor = attr["compressor"];
         assert(compressor == "vtkZLibDataCompressor"); // only ZLib compression supported
       }
-
-      // std::cout << "<VTKFile type='" << attr["type"] << "' version='" << attr["version"] << "' header_type='" << attr["header_type"] << "' byte_order='" << attr["byte_order"] << "' compressor='" << attr["compressor"] << "'>\n";
-
       section = VTK_FILE;
     }
     else if (isSection(line, "/VTKFile", section, VTK_FILE))
@@ -121,8 +153,6 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename)
       // Store attributes of DataArray
       dataArray_[data_name] = {data_type, data_components, data_offset};
 
-      // std::cout << "<DataArray type='" << attr["type"] << "' Name='" << data_name << "' NumberOfComponents='" << attr["NumberOfComponents"] << "' format='" << data_format << "' offset='" << attr["offset"] << "' " << (closed ? "/" : "") << ">\n";
-
       // Skip section in appended mode
       if (data_format == "appended") {
         if (!closed) {
@@ -187,7 +217,7 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename)
         }
         break;
       case POINTS_DATA_ARRAY:
-        section = readPoints(input);
+        section = readPoints(input, data_name);
         break;
       case CD_DATA_ARRAY:
         if (data_type == Vtk::FLOAT32) {
@@ -213,7 +243,62 @@ void VtkReader<Grid,Creator>::readFromFile (std::string const& filename)
   if (section != NO_SECTION)
     DUNE_THROW(IOError, "VTK-File is incomplete. It must end with </VTKFile>!");
 
-  createGrid();
+  if (create)
+    createGrid();
+}
+
+
+template <class Grid, class Creator>
+void VtkReader<Grid,Creator>::readParallelFileFromStream (std::ifstream& input, int commRank, int commSize, bool create)
+{
+  pieces_.clear();
+
+  Sections section = NO_SECTION;
+  for (std::string line; std::getline(input, line); ) {
+    ltrim(line);
+
+    if (isSection(line, "VTKFile", section)) {
+      bool closed = false;
+      auto attr = parseXml(line, closed);
+
+      if (!attr["type"].empty())
+        assert(attr["type"] == "PUnstructuredGrid");
+      if (!attr["version"].empty())
+        assert(std::stod(attr["version"]) == 1.0);
+      if (!attr["byte_order"].empty())
+        assert(attr["byte_order"] == "LittleEndian");
+      if (!attr["header_type"].empty())
+        assert(attr["header_type"] == "UInt64");
+      if (!attr["compressor"].empty())
+        assert(attr["compressor"] == "vtkZLibDataCompressor"); // only ZLib compression supported
+      section = VTK_FILE;
+    }
+    else if (isSection(line, "/VTKFile", section, VTK_FILE))
+      section = NO_SECTION;
+    else if (isSection(line, "PUnstructuredGrid", section, VTK_FILE))
+      section = UNSTRUCTURED_GRID;
+    else if (isSection(line, "/PUnstructuredGrid", section, UNSTRUCTURED_GRID))
+      section = VTK_FILE;
+    else if (isSection(line, "Piece", section, UNSTRUCTURED_GRID)) {
+      bool closed = false;
+      auto attr = parseXml(line, closed);
+
+      assert(attr.count("Source") > 0);
+      pieces_.push_back(attr["Source"]);
+    }
+
+    if (section == NO_SECTION)
+      break;
+  }
+
+  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();
 }
 
 
@@ -250,14 +335,16 @@ Sections readDataArray (IStream& input, std::vector<T>& values, std::size_t max_
 
 template <class Grid, class Creator>
 typename VtkReader<Grid,Creator>::Sections
-VtkReader<Grid,Creator>::readPoints (std::ifstream& input)
+VtkReader<Grid,Creator>::readPoints (std::ifstream& input, std::string name)
 {
   using T = typename GlobalCoordinate::value_type;
   assert(numberOfPoints_ > 0);
   assert(dataArray_["points"].components == 3u);
 
+  Sections sec;
+
   std::vector<T> point_values;
-  auto sec = readDataArray(input, point_values, 3*numberOfPoints_, POINTS_DATA_ARRAY, POINTS);
+  sec = readDataArray(input, point_values, 3*numberOfPoints_, POINTS_DATA_ARRAY, POINTS);
   assert(sec == POINTS);
   assert(point_values.size() == 3*numberOfPoints_);
 
@@ -282,7 +369,6 @@ void VtkReader<Grid,Creator>::readPointsAppended (std::ifstream& input)
 {
   assert(numberOfPoints_ > 0);
   assert(dataArray_["points"].components == 3u);
-
   std::vector<T> point_values;
   readAppended(input, point_values, dataArray_["points"].offset);
   assert(point_values.size() == 3*numberOfPoints_);
@@ -321,6 +407,9 @@ VtkReader<Grid,Creator>::readCells (std::ifstream& input, std::string name)
     else
       max_size = numberOfCells_ * max_vertices;
     sec = readDataArray(input, vec_connectivity, max_size, CELLS_DATA_ARRAY, CELLS);
+  } else if (name == "global_point_ids") {
+    sec = readDataArray(input, vec_point_ids, numberOfPoints_, CELLS_DATA_ARRAY, CELLS);
+    assert(vec_point_ids.size() == numberOfPoints_);
   }
 
   return sec;
@@ -346,6 +435,13 @@ void VtkReader<Grid,Creator>::readCellsAppended (std::ifstream& input)
   assert(connectivity_data.type == Vtk::INT64);
   readAppended(input, vec_connectivity, connectivity_data.offset);
   assert(vec_connectivity.size() == std::size_t(vec_offsets.back()));
+
+  if (dataArray_.count("global_point_ids") > 0) {
+    auto point_id_data = dataArray_["global_point_ids"];
+    assert(point_id_data.type == Vtk::UINT64);
+    readAppended(input, vec_point_ids, point_id_data.offset);
+    assert(vec_point_ids.size() == numberOfPoints_);
+  }
 }
 
 
@@ -361,7 +457,7 @@ template <class T, class IStream>
 void read_compressed (T* buffer, unsigned char* buffer_in,
                       std::uint64_t bs, std::uint64_t cbs, IStream& input)
 {
-#ifdef HAVE_VTK_ZLIB
+#if HAVE_VTK_ZLIB
   uLongf uncompressed_space = uLongf(bs);
   uLongf compressed_space = uLongf(cbs);
 
@@ -442,7 +538,7 @@ void VtkReader<Grid,Creator>::createGrid () const
   assert(vec_types.size() == numberOfCells_);
   assert(vec_offsets.size() == numberOfCells_);
 
-  Creator::create(*factory_, vec_points, vec_types, vec_offsets, vec_connectivity);
+  Creator::create(*factory_, vec_points, vec_point_ids, vec_types, vec_offsets, vec_connectivity);
 }
 
 // Assume input already read the line <AppendedData ...>

dune/vtk/writers/vtkunstructuredgridwriter.hh

@@ -75,10 +75,14 @@ namespace Dune
     // Write the element connectivity to the output stream `out`. In case
     // of binary format, stores the streampos of XML attributes "offset" in the
     // vector `offsets`.
-    void writeCells (std::ofstream& oust,
+    void writeCells (std::ofstream& out,
                      std::vector<pos_type>& offsets,
                      Std::optional<std::size_t> timestep = {}) const;
 
+    void writePointIds (std::ofstream& out,
+                        std::vector<pos_type>& offsets,
+                        Std::optional<std::size_t> timestep = {}) const;
+
   private:
     using Super::dataCollector_;
     using Super::format_;

dune/vtk/writers/vtkunstructuredgridwriter.impl.hh

@@ -37,6 +37,7 @@ void VtkUnstructuredGridWriter<GV,DC>
   // Write element connectivity, types and offsets
   out << "<Cells>\n";
   writeCells(out, offsets);
+  writePointIds(out, offsets);
   out << "</Cells>\n";
 
   // Write data associated with grid points
@@ -134,14 +135,17 @@ void VtkUnstructuredGridWriter<GV,DC>
 
   // Write point coordinates
   out << "<Points>\n";
-  for (std::size_t i = 0; i < timesteps.size(); ++i)
+  for (std::size_t i = 0; i < timesteps.size(); ++i) {
     this->writePoints(out, offsets[i], i);
+  }
   out << "</Points>\n";