grid creators and tests added

dune/vtk/gridcreator.hh

+#pragma once
+
+#include <cassert>
+#include <cstdint>
+#include <limits>
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#include "vtktypes.hh"
+
+namespace Dune
+{
+  // 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::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);
+
+      std::size_t idx = 0;
+      for (std::size_t i = 0; i < types.size(); ++i) {
+        if (Vtk::Map::type.count(types[i]) == 0)
+          DUNE_THROW(Exception, "Unknown ElementType: " << types[i]);
+        auto type = Vtk::Map::type[types[i]];
+        Vtk::CellType cellType{type};
+        auto refElem = referenceElement<double,Grid::dimension>(type);
+
+        std::size_t 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 (std::size_t 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 (std::size_t 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::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;
+        }
+      }
+
+      idx = 0;
+      for (std::size_t i = 0; i < types.size(); ++i) {
+        if (Vtk::Map::type.count(types[i]) == 0)
+          DUNE_THROW(Exception, "Unknown ElementType: " << types[i]);
+        auto type = Vtk::Map::type[types[i]];
+        Vtk::CellType cellType{type};
+
+        std::size_t nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]);
+        assert(nNodes > 0);
+        std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes);
+        for (std::size_t 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 (std::size_t j = 0; j < nNodes; ++j)
+            cell[j] = vtk_cell[cellType.permutation(j)];
+
+          factory.insertElement(type,cell);
+        }
+      }
+    }
+  };
+
+} // end namespace Dune

dune/vtk/vtkreader.hh

@@ -4,6 +4,7 @@
 #include <map>
 
 #include "filereader.hh"
+#include "gridcreator.hh"
 #include "vtktypes.hh"
 
 namespace Dune
@@ -15,9 +16,9 @@ namespace Dune
    *
    * Assumption on the file structure: Each XML tag must be on a separate line.
    **/
-  template <class Grid>
+  template <class Grid, class GridCreator = DefaultGridCreator>
   class VtkReader
-      : public FileReader<Grid, VtkReader<Grid>>
+      : public FileReader<Grid, VtkReader<Grid, GridCreator>>
   {
     // Sections visited during the xml parsing
     enum Sections {
@@ -101,7 +102,8 @@ namespace Dune
     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_types, \ref vec_offsets, and \ref vec_connectivity
+    // \ref vec_points, \ref vec_types, \ref vec_offsets, and \ref vec_connectivity,
+    // by passing to \ref GridCreator.
     void createGrid() const;
 
   private:
@@ -111,6 +113,7 @@ namespace Dune
     Vtk::FormatTypes format_;
 
     // Temporary data to construct the grid elements
+    std::vector<GlobalCoordinate> vec_points;
     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

dune/vtk/vtkreader.impl.hh

@@ -14,8 +14,8 @@
 
 namespace Dune {
 
-template <class Grid>
-void VtkReader<Grid>::readFromFile (std::string const& filename)
+template <class Grid, class Creator>
+void VtkReader<Grid,Creator>::readFromFile (std::string const& filename)
 {
   // check whether file exists!
   if (!filesystem::exists(filename))
@@ -127,10 +127,8 @@ void VtkReader<Grid>::readFromFile (std::string const& filename)
         if (!closed) {
           while (std::getline(input, line)) {
             ltrim(line);
-            if (line.substr(1,10) == "/DataArray") {
-              // std::cout << "</DataArray>\n";
+            if (line.substr(1,10) == "/DataArray")
               break;
-            }
           }
         }
         continue;
@@ -249,9 +247,9 @@ Sections readDataArray (IStream& input, std::vector<T>& values, std::size_t max_
 // @}
 
 
-template <class Grid>
-typename VtkReader<Grid>::Sections
-VtkReader<Grid>::readPoints (std::ifstream& input)
+template <class Grid, class Creator>
+typename VtkReader<Grid,Creator>::Sections
+VtkReader<Grid,Creator>::readPoints (std::ifstream& input)
 {
   using T = typename GlobalCoordinate::value_type;
   assert(numberOfPoints_ > 0);
@@ -264,21 +262,22 @@ VtkReader<Grid>::readPoints (std::ifstream& input)
 
   // extract points from continuous values
   GlobalCoordinate p;
+  vec_points.reserve(numberOfPoints_);
   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] = point_values[idx++];
     idx += (3u - p.size());
-    factory_->insertVertex(p);
+    vec_points.push_back(p);
   }
 
   return sec;
 }
 
 
-template <class Grid>
+template <class Grid, class Creator>
   template <class T>
-void VtkReader<Grid>::readPointsAppended (std::ifstream& input)
+void VtkReader<Grid,Creator>::readPointsAppended (std::ifstream& input)
 {
   assert(numberOfPoints_ > 0);
   assert(dataArray_["points"].components == 3u);
@@ -289,20 +288,20 @@ void VtkReader<Grid>::readPointsAppended (std::ifstream& input)
 
   // extract points from continuous values
   GlobalCoordinate p;
+  vec_points.reserve(numberOfPoints_);
   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++]);
     idx += (3u - p.size());
-
-    factory_->insertVertex(p);
+    vec_points.push_back(p);
   }
 }
 
 
-template <class Grid>
-typename VtkReader<Grid>::Sections
-VtkReader<Grid>::readCells (std::ifstream& input, std::string name)
+template <class Grid, class Creator>
+typename VtkReader<Grid,Creator>::Sections
+VtkReader<Grid,Creator>::readCells (std::ifstream& input, std::string name)
 {
   Sections sec = CELLS_DATA_ARRAY;
 
@@ -327,8 +326,8 @@ VtkReader<Grid>::readCells (std::ifstream& input, std::string name)
 }
 
 
-template <class Grid>
-void VtkReader<Grid>::readCellsAppended (std::ifstream& input)
+template <class Grid, class Creator>
+void VtkReader<Grid,Creator>::readCellsAppended (std::ifstream& input)
 {
   assert(numberOfCells_ > 0);
   auto types_data = dataArray_["types"];
@@ -384,9 +383,9 @@ void read_compressed (T* buffer, unsigned char* buffer_in,
 // @}
 
 
-template <class Grid>
+template <class Grid, class Creator>
   template <class T>
-void VtkReader<Grid>::readAppended (std::ifstream& input, std::vector<T>& values, std::uint64_t offset)
+void VtkReader<Grid,Creator>::readAppended (std::ifstream& input, std::vector<T>& values, std::uint64_t offset)
 {
   input.seekg(offset0_ + offset);
 
@@ -435,35 +434,19 @@ void VtkReader<Grid>::readAppended (std::ifstream& input, std::vector<T>& values
 }
 
 
-template <class Grid>
-void VtkReader<Grid>::createGrid () const
+template <class Grid, class Creator>
+void VtkReader<Grid,Creator>::createGrid () const
 {
-  assert(vec_types.size() == vec_offsets.size());
-  std::size_t idx = 0;
-  for (std::size_t i = 0; i < vec_types.size(); ++i) {
-    if (Vtk::Map::type.count(vec_types[i]) == 0)
-      DUNE_THROW(Exception, "Unknown ElementType: " << vec_types[i]);
-    auto type = Vtk::Map::type[vec_types[i]];
-    Vtk::CellType cellType{type};
-
-    std::size_t nNodes = vec_offsets[i] - (i == 0 ? 0 : vec_offsets[i-1]);
-    assert(nNodes > 0);
-    std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes);
-    for (std::size_t j = 0; j < nNodes; ++j)
-      vtk_cell.push_back( vec_connectivity[idx++] );
-
-    // apply index permutation
-    std::vector<unsigned int> cell(nNodes);
-    for (std::size_t j = 0; j < nNodes; ++j)
-      cell[j] = vtk_cell[cellType.localIndex(j)];
-
-    factory_->insertElement(type,cell);
-  }
-}
+  assert(vec_points.size() == numberOfPoints_);
+  assert(vec_types.size() == numberOfCells_);
+  assert(vec_offsets.size() == numberOfCells_);
 
+  Creator::create(*factory_, vec_points, vec_types, vec_offsets, vec_connectivity);
+}
 
-template <class Grid>
-std::uint64_t VtkReader<Grid>::findAppendedDataPosition (std::ifstream& input) const
+// Assume input already read the line <AppendedData ...>
+template <class Grid, class Creator>
+std::uint64_t VtkReader<Grid,Creator>::findAppendedDataPosition (std::ifstream& input) const
 {
   char c;
   while (input.get(c) && std::isblank(c)) { /*do nothing*/ }
@@ -476,8 +459,8 @@ std::uint64_t VtkReader<Grid>::findAppendedDataPosition (std::ifstream& input) c
 }
 
 
-template <class Grid>
-std::map<std::string, std::string> VtkReader<Grid>::parseXml (std::string const& line, bool& closed)
+template <class Grid, class Creator>
+std::map<std::string, std::string> VtkReader<Grid,Creator>::parseXml (std::string const& line, bool& closed)
 {
   closed = false;
   std::map<std::string, std::string> attr;

dune/vtk/vtktypes.cc

@@ -31,6 +31,7 @@ std::map<std::string, DataTypes> Map::datatype = {
 
 
 CellType::CellType (GeometryType const& t)
+  : noPermutation_(true)
 {
   if (t.isVertex()) {
     type_ = 1;
@@ -55,14 +56,17 @@ CellType::CellType (GeometryType const& t)
   else if (t.isHexahedron()) {
     type_ = 12;
     permutation_ = {0,1,3,2,4,5,7,6};
+    noPermutation_ = false;
   }
   else if (t.isPrism()) {
     type_ = 13;
     permutation_ = {0,2,1,3,5,4};
+    noPermutation_ = false;
   }
   else if (t.isPyramid()) {
     type_ = 14;
     permutation_ = {0,1,3,2,4};
+    noPermutation_ = false;
   }
   else {
     std::cerr << "Geometry Type not supported by VTK!\n";

dune/vtk/vtktypes.hh

@@ -58,14 +58,20 @@ namespace Dune
       }
 
       /// Return a permutation of Dune elemenr vertices to conform to VTK element numbering
-      int localIndex (int idx) const
+      int permutation (int idx) const
       {
         return permutation_[idx];
       }
 
+      bool noPermutation () const
+      {
+        return noPermutation_;
+      }
+
     private:
       std::uint8_t type_;
       std::vector<int> permutation_;
+      bool noPermutation_;
     };
 
   } // end namespace Vtk

dune/vtk/vtkwriter.impl.hh

@@ -165,9 +165,9 @@ std::vector<T> getVertexData (GridView const& gridView, GlobalFunction const& fc
     Vtk::CellType cellType{e.type()};
     auto refElem = referenceElement(e.geometry());
     for (int j = 0; j < e.subEntities(dim); ++j) {
-      std::size_t idx = fct.ncomps() * indexSet.subIndex(e,cellType.localIndex(j),dim);
+      std::size_t idx = fct.ncomps() * indexSet.subIndex(e,cellType.permutation(j),dim);
       for (int comp = 0; comp < fct.ncomps(); ++comp)
-        data[idx + comp] = T(localFct.evaluate(comp, refElem.position(cellType.localIndex(j),dim)));
+        data[idx + comp] = T(localFct.evaluate(comp, refElem.position(cellType.permutation(j),dim)));
     }
     localFct.unbind();
   }
@@ -254,7 +254,7 @@ void VtkWriter<GridView>::writeCells (std::ofstream& out, std::vector<pos_type>&
     for (auto const& c : elements(gridView_)) {
       auto cellType = getType(c.type());
       for (int j = 0; j < c.subEntities(dimension); ++j)
-        out << indexSet.subIndex(c,cellType.localIndex(j),dimension) << (++i % 6 != 0 ? ' ' : '\n');
+        out << indexSet.subIndex(c,cellType.permutation(j),dimension) << (++i % 6 != 0 ? ' ' : '\n');
     }
     out << (i % 6 != 0 ? "\n" : "") << "</DataArray>\n";
 
@@ -438,7 +438,7 @@ std::array<std::uint64_t,3> VtkWriter<GridView>::writeCellsAppended (std::ofstre
   for (auto const& c : elements(gridView_)) {
     auto cellType = getType(c.type());
     for (int j = 0; j < c.subEntities(dimension); ++j)
-      connectivity.push_back( std::int64_t(indexSet.subIndex(c,cellType.localIndex(j),dimension)) );
+      connectivity.push_back( std::int64_t(indexSet.subIndex(c,cellType.permutation(j),dimension)) );
     cell_offsets.push_back(old_o += c.subEntities(dimension));
     cell_types.push_back(cellType.type());
   }

src/CMakeLists.txt

@@ -4,4 +4,6 @@ target_link_libraries("vtkwriter" dunevtk)
 
 add_executable("vtkreader" vtkreader.cc)
 target_link_dune_default_libraries("vtkreader")
-target_link_libraries("vtkreader" dunevtk)
\ No newline at end of file
+target_link_libraries("vtkreader" dunevtk)
+
+add_subdirectory(test)
\ No newline at end of file

src/test/CMakeLists.txt

+dune_add_test(SOURCES 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

src/test/mixed_element_test.cc

+// -*- 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>
+
+#include <dune/vtk/vtkreader.hh>
+#include <dune/vtk/vtkwriter.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);
+
+  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 TestCases = std::set<std::tuple<std::string,Vtk::FormatTypes,Vtk::DataTypes>>;
+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 GridView>
+void writer_test (GridView const& gridView)
+{
+  VtkWriter<GridView> vtkWriter(gridView);
+  for (auto const& test_case : test_cases) {
+    vtkWriter.write("/tmp/reader_writer_test_" + std::get<0>(test_case) + ".vtu",
+      std::get<1>(test_case), std::get<2>(test_case));
+  }
+}
+
+template <class Grid, class Test>
+void reader_test (Test& test)
+{
+  for (auto const& test_case : test_cases) {
+    auto grid = VtkReader<Grid>::read("/tmp/reader_writer_test_" + std::get<0>(test_case) + ".vtu");
+    VtkWriter<typename Grid::LeafGridView> vtkWriter(grid->leafGridView());
+    vtkWriter.write("/tmp/reader_writer_test_" + std::get<0>(test_case) + "_2.vtu",
+      std::get<1>(test_case), std::get<2>(test_case));
+    test.check(compare_files("/tmp/reader_writer_test_" + std::get<0>(test_case) + ".vtu",
+                             "/tmp/reader_writer_test_" + std::get<0>(test_case) + "_2.vtu"));
+  }
+}
+
+
+template <int I>
+using int_ = std::integral_constant<int,I>;
+
+int main (int argc, char** argv)
+{
+  Dune::MPIHelper::instance(argc, argv);
+
+  TestSuite test{};
+
+#ifdef HAVE_UG
+  // Test VtkWriter for UGGrid
+  using GridType = UGGrid<2>;
+  GridFactory<GridType> factory;
+
+  using X = FieldVector<double,2>;
+  using E = std::vector<unsigned int>;
+  factory.insertVertex(X{0.0, 0.0}); // 0
+  factory.insertVertex(X{1.0, 0.0}); // 1
+  factory.insertVertex(X{1.0, 1.0}); // 2
+  factory.insertVertex(X{0.0, 1.0}); // 3
+  factory.insertVertex(X{1.5, 0.5}); // 4
+
+  factory.insertElement(GeometryTypes::quadrilateral, E{0,1,3,2}); // quadrilateral
+  factory.insertElement(GeometryTypes::triangle, E{1,4,2}); // triangle
+
+  {
+    std::unique_ptr<GridType> gridPtr{ factory.createGrid() };
+    writer_test(gridPtr->leafGridView());
+  }
+  reader_test<GridType>(test);
+#endif
+
+  return test.exit();
+}

src/test/reader_writer_test.cc

+// -*- 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>
+
+#include <dune/vtk/vtkreader.hh>
+#include <dune/vtk/vtkwriter.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);
+
+  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;