Merge branch 'feature/structuredgridwriter' into 'master'

Feature/structuredgridwriter

See merge request spraetor/dune-vtk!1

dune.module

@@ -8,3 +8,4 @@ Version: 0.1
 Maintainer: simon.praetorius@tu-dresden.de
 #depending on
 Depends: dune-grid dune-functions
+Suggests: dune-spgrid

dune/vtk/datacollector.hh

@@ -113,7 +113,7 @@ protected: // default implementations
     std::vector<T> data(numCells() * fct.ncomps());
     auto const& indexSet = gridView_.indexSet();
     auto localFct = localFunction(fct);
-    for (auto const& e : elements(gridView_)) {
+    for (auto const& e : elements(gridView_, Partitions::all)) {
       localFct.bind(e);
       auto geometry = e.geometry();
       std::size_t idx = fct.ncomps() * indexSet.index(e);

dune/vtk/datacollectors/continuousdatacollector.hh

@@ -33,7 +33,7 @@ public:
   {
     std::vector<T> data(this->numPoints() * 3);
     auto const& indexSet = gridView_.indexSet();
-    for (auto const& vertex : vertices(gridView_)) {
+    for (auto const& vertex : vertices(gridView_, Partitions::all)) {
       std::size_t idx = 3 * indexSet.index(vertex);
       auto v = vertex.geometry().center();
       for (std::size_t j = 0; j < v.size(); ++j)
@@ -61,7 +61,7 @@ public:
     cells.types.reserve(this->numCells());
 
     std::int64_t old_o = 0;
-    for (auto const& c : elements(gridView_)) {
+    for (auto const& c : elements(gridView_, Partitions::all)) {
       Vtk::CellType cellType(c.type());
       for (int j = 0; j < c.subEntities(dim); ++j)
         cells.connectivity.push_back( std::int64_t(indexSet.subIndex(c,cellType.permutation(j),dim)) );
@@ -79,7 +79,7 @@ public:
     std::vector<T> data(this->numPoints() * fct.ncomps());
     auto const& indexSet = gridView_.indexSet();
     auto localFct = localFunction(fct);
-    for (auto const& e : elements(gridView_)) {
+    for (auto const& e : elements(gridView_, Partitions::all)) {
       localFct.bind(e);
       Vtk::CellType cellType{e.type()};
       auto refElem = referenceElement(e.geometry());

dune/vtk/datacollectors/discontinuousdatacollector.hh

@@ -30,7 +30,7 @@ public:
     indexMap_.resize(gridView_.size(dim));
     std::int64_t vertex_idx = 0;
     auto const& indexSet = gridView_.indexSet();
-    for (auto const& c : elements(gridView_)) {
+    for (auto const& c : elements(gridView_, Partitions::interior)) {
       numPoints_ += c.subEntities(dim);
       for (int i = 0; i < c.subEntities(dim); ++i)
         indexMap_[indexSet.subIndex(c, i, dim)] = vertex_idx++;
@@ -49,7 +49,7 @@ public:
   {
     std::vector<T> data(this->numPoints() * 3);
     auto const& indexSet = gridView_.indexSet();
-    for (auto const& element : elements(gridView_)) {
+    for (auto const& element : elements(gridView_, Partitions::interior)) {
       for (int i = 0; i < element.subEntities(dim); ++i) {
         std::size_t idx = 3 * indexMap_[indexSet.subIndex(element, i, dim)];
         auto v = element.geometry().corner(i);
@@ -72,7 +72,7 @@ public:
 
     std::int64_t old_o = 0;
     auto const& indexSet = gridView_.indexSet();
-    for (auto const& c : elements(gridView_)) {
+    for (auto const& c : elements(gridView_, Partitions::interior)) {
       Vtk::CellType cellType(c.type());
       for (int j = 0; j < c.subEntities(dim); ++j) {
         std::int64_t vertex_idx = indexMap_[indexSet.subIndex(c,cellType.permutation(j),dim)];
@@ -92,7 +92,7 @@ public:
     std::vector<T> data(this->numPoints() * fct.ncomps());
     auto const& indexSet = gridView_.indexSet();
     auto localFct = localFunction(fct);
-    for (auto const& e : elements(gridView_)) {
+    for (auto const& e : elements(gridView_, Partitions::interior)) {
       localFct.bind(e);
       Vtk::CellType cellType{e.type()};
       auto refElem = referenceElement(e.geometry());

dune/vtk/datacollectors/quadraticdatacollector.hh

@@ -37,7 +37,7 @@ public:
   {
     std::vector<T> data(this->numPoints() * 3);
     auto const& indexSet = gridView_.indexSet();
-    for (auto const& element : elements(gridView_)) {
+    for (auto const& element : elements(gridView_, Partitions::interior)) {
       auto geometry = element.geometry();
       auto refElem = referenceElement<T,dim>(element.type());
 
@@ -77,7 +77,7 @@ public:
 
     std::int64_t old_o = 0;
     auto const& indexSet = gridView_.indexSet();
-    for (auto const& c : elements(gridView_)) {
+    for (auto const& c : elements(gridView_, Partitions::interior)) {
       Vtk::CellType cellType(c.type(), Vtk::QUADRATIC);
       for (int j = 0; j < c.subEntities(dim); ++j) {
         int k = cellType.permutation(j);
@@ -103,7 +103,7 @@ public:
     std::vector<T> data(this->numPoints() * fct.ncomps());
     auto const& indexSet = gridView_.indexSet();
     auto localFct = localFunction(fct);
-    for (auto const& e : elements(gridView_)) {
+    for (auto const& e : elements(gridView_, Partitions::interior)) {
       localFct.bind(e);
       Vtk::CellType cellType{e.type(), Vtk::QUADRATIC};
       auto refElem = referenceElement(e.geometry());

dune/vtk/datacollectors/spdatacollector.hh

+#pragma once
+
+#if HAVE_DUNE_SPGRID
+#include <dune/grid/spgrid.hh>
+#endif
+#include <dune/vtk/datacollectors/structureddatacollector.hh>
+
+namespace Dune { namespace experimental
+{
+#if HAVE_DUNE_SPGRID
+
+// Specialization for SPGrid
+template <class GridView>
+class SPDataCollector
+    : public StructuredDataCollectorInterface<GridView, SPDataCollector<GridView>>
+{
+  enum { dim = GridView::dimension };
+
+  using Self = SPDataCollector;
+  using Super = StructuredDataCollectorInterface<GridView, Self>;
+  using Super::gridView_;
+  using ctype = typename GridView::ctype;
+
+public:
+  SPDataCollector (GridView const& gridView)
+    : Super(gridView)
+    , wholeExtent_(filledArray<6,int>(0))
+    , origin_(0.0)
+    , spacing_(0.0)
+  {}
+
+  std::array<int, 6> const& wholeExtentImpl () const
+  {
+    return wholeExtent_;
+  }
+
+  auto const& originImpl () const
+  {
+    return origin_;
+  }
+
+  auto const& spacingImpl () const
+  {
+    return spacing_;
+  }
+
+  void updateImpl ()
+  {
+    Super::updateImpl();
+
+    auto const& begin = gridView_.impl().gridLevel().globalMesh().begin();
+    auto const& end = gridView_.impl().gridLevel().globalMesh().end();
+    auto const& cube = gridView_.grid().domain().cube();
+
+    for (int i = 0; i < dim; ++i) {
+      wholeExtent_[2*i] = begin[i];
+      wholeExtent_[2*i+1] = end[i];
+      spacing_[i] = cube.width()[i] / (end[i] - begin[i]);
+      origin_[i] = cube.origin()[i];
+    }
+  }
+
+  template <class Writer>
+  void writeLocalPieceImpl (Writer const& writer) const
+  {
+    auto extent = filledArray<6,int>(0);
+    auto const& localMesh = gridView_.impl().gridLevel().localMesh();
+    for (int i = 0; i < dim; ++i) {
+      extent[2*i] = localMesh.begin()[i];
+      extent[2*i+1] = localMesh.end()[i];
+    }
+    writer(extent);
+  }
+
+  template <class Writer>
+  void writePiecesImpl (Writer const& writer) const
+  {
+    auto extent = filledArray<6,int>(0);
+    auto const& partitionList = gridView_.impl().gridLevel().decomposition_;
+    // for (auto const& part : partitionList)
+    for (std::size_t p = 0; p < partitionList.size(); ++p)
+    {
+      auto const& part = partitionList[p];
+      for (int i = 0; i < dim; ++i) {
+        extent[2*i] = part.begin()[i];
+        extent[2*i+1] = part.begin()[i] + part.width(i);
+      }
+      writer(p, extent, true);
+    }
+  }
+
+private:
+  std::array<int, 6> wholeExtent_;
+  std::array<int, 6> extent_;
+  FieldVector<ctype,3> spacing_;
+  FieldVector<ctype,3> origin_;
+  std::vector<std::size_t> indexMap_;
+};
+
+namespace Impl
+{
+  template<class GridView, class ct, int dim, template< int > class Ref, class Comm>
+  struct StructuredDataCollector<GridView, SPGrid<ct,dim,Ref,Comm>>
+  {
+    using type = SPDataCollector<GridView>;
+  };
+}
+
+#endif // HAVE_DUNE_SPGRID
+}} // end namespace Dune::experimental

dune/vtk/datacollectors/structureddatacollector.hh

+#pragma once
+
+#include <dune/vtk/datacollectors/continuousdatacollector.hh>
+
+namespace Dune { namespace experimental
+{
+
+namespace Impl
+{
+  // Should be specialized for concrete structured grid
+  template <class GridView, class Grid>
+  struct StructuredDataCollector;
+}
+
+template <class GridView>
+using StructuredDataCollector = typename Impl::StructuredDataCollector<GridView, typename GridView::Grid>::type;
+
+
+/// The Interface for structured data-collectors
+template <class GridView, class Derived>
+class StructuredDataCollectorInterface
+    : public DataCollectorInterface<GridView, Derived>
+{
+protected:
+  using Super = DataCollectorInterface<GridView, Derived>;
+  using Super::gridView_;
+  using ctype = typename GridView::ctype;
+
+public:
+  StructuredDataCollectorInterface (GridView const& gridView)
+    : Super(gridView)
+    , defaultDataCollector_(gridView)
+    , ghostLevel_(gridView.overlapSize(0))
+  {}
+
+  /// Return number of grid vertices
+  std::uint64_t numPointsImpl () const
+  {
+    return gridView_.size(GridView::dimension);
+  }
+
+  void updateImpl ()
+  {
+    defaultDataCollector_.update();
+  }
+
+  std::array<int, 6> const& wholeExtent () const
+  {
+    return this->asDerived().wholeExtentImpl();
+  }
+
+  FieldVector<ctype, 3> const& origin () const
+  {
+    return this->asDerived().originImpl();
+  }
+
+  FieldVector<ctype, 3> const& spacing () const
+  {
+    return this->asDerived().spacingImpl();
+  }
+
+  template <class Writer>
+  void writeLocalPiece (Writer const& writer) const
+  {
+    this->asDerived().writeLocalPieceImpl(writer);
+  }
+
+  template <class Writer>
+  void writePieces (Writer const& writer) const
+  {
+    this->asDerived().writePiecesImpl(writer);
+  }
+
+  int ghostLevel () const
+  {
+    return this->asDerived().ghostLevelImpl();
+  }
+
+  /// Return the coordinates along the ordinates x, y, and z
+  template <class T>
+  std::array<std::vector<T>, 3> coordinates () const
+  {
+    return this->asDerived().template coordinatesImpl<T>();
+  }
+
+
+public:
+
+  /// Return the coordinates of all grid vertices in the order given by the indexSet
+  template <class T>
+  std::vector<T> pointsImpl () const
+  {
+    return defaultDataCollector_.template points<T>();
+  }
+
+  /// Evaluate the `fct` at the corners of the elements
+  template <class T, class GlobalFunction>
+  std::vector<T> pointDataImpl (GlobalFunction const& fct) const
+  {
+    return defaultDataCollector_.template pointData<T>(fct);
+  }
+
+
+private:
+
+  int ghostLevelImpl () const
+  {
+    return ghostLevel_;
+  }
+
+  template <class T>
+  std::array<std::vector<T>, 3> coordinatesImpl () const
+  {
+    auto origin = this->origin();
+    auto spacing = this->spacing();
+
+    std::array<std::vector<T>, 3> ordinates{};
+    writeLocalPiece([&ordinates,&origin,&spacing](auto const& extent) {
+      for (std::size_t d = 0; d < GridView::dimension; ++d) {
+        auto s = extent[2*d+1] - extent[2*d] + 1;
+        ordinates[d].resize(s);
+        for (std::size_t i = 0; i < s; ++i)
+          ordinates[d][i] = origin[d] + (extent[2*d] + i)*spacing[d];
+      }
+    });
+
+    for (std::size_t d = GridView::dimension; d < 3; ++d)
+      ordinates[d].resize(1, T(0));
+
+    return ordinates;
+  }
+
+private:
+  DefaultDataCollector<GridView> defaultDataCollector_;
+  int ghostLevel_;
+};
+
+}} // end namespace Dune::experimental

dune/vtk/datacollectors/yaspdatacollector.hh

+#pragma once
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/vtk/datacollectors/structureddatacollector.hh>
+
+namespace Dune { namespace experimental
+{
+// Specialization for YaspGrid
+template <class GridView>
+class YaspDataCollector
+    : public StructuredDataCollectorInterface<GridView, YaspDataCollector<GridView>>
+{
+  enum { dim = GridView::dimension };
+
+  using Self = YaspDataCollector;
+  using Super = StructuredDataCollectorInterface<GridView, Self>;
+  using Super::gridView_;
+  using ctype = typename GridView::ctype;
+
+public:
+  YaspDataCollector (GridView const& gridView)
+    : Super(gridView)
+    , wholeExtent_(filledArray<6,int>(0))
+    , origin_(0.0)
+    , spacing_(0.0)
+    , level_(gridView.template begin<0,All_Partition>()->level())
+  {}
+
+  std::array<int, 6> const& wholeExtentImpl () const
+  {
+    return wholeExtent_;
+  }
+
+  auto const& originImpl () const
+  {
+    return origin_;
+  }
+
+  auto const& spacingImpl () const
+  {
+    return spacing_;
+  }
+
+  void updateImpl ()
+  {
+    Super::updateImpl();
+    localPieceCalled_ = false;
+
+    for (int i = 0; i < dim; ++i) {
+      wholeExtent_[2*i] = 0;
+      wholeExtent_[2*i+1] = gridView_.grid().levelSize(level_,i);
+    }
+
+    auto it = gridView_.grid().begin(level_);
+    initGeometry(it->coords);
+
+#if HAVE_MPI
+    int rank = -1;
+    int num_ranks = -1;
+    MPI_Comm_rank(gridView_.comm(), &rank);
+    MPI_Comm_size(gridView_.comm(), &num_ranks);
+
+    if (rank == 0) {
+      extents_.resize(num_ranks);
+      requests_.resize(num_ranks);
+      for (int i = 1; i < num_ranks; ++i)
+        MPI_Irecv(extents_[i].data(), extents_[i].size(), MPI_INT, i, /*tag=*/6, gridView_.comm(), &requests_[i]);
+    }
+#endif
+  }
+
+  template <class Coords>
+  void initGeometry (Coords const& coords) { DUNE_THROW(NotImplemented, "Coordinate-Type not implemented!"); }
+
+  void initGeometry (EquidistantCoordinates<ctype,dim> const& coords)
+  {
+    for (int i = 0; i < dim; ++i) {
+      spacing_[i] = coords.meshsize(i,0);
+      origin_[i] = 0;
+    }
+  }
+
+  void initGeometry (EquidistantOffsetCoordinates<ctype,dim> const& coords)
+  {
+    for (int i = 0; i < dim; ++i) {
+      spacing_[i] = coords.meshsize(i,0);
+      origin_[i] = coords.origin(i);
+    }
+  }
+
+  void initGeometry (TensorProductCoordinates<ctype,dim> const& coords)
+  {
+    for (int i = 0; i < dim; ++i) {
+      spacing_[i] = coords.meshsize(i,0); // is not constant, but also not used.
+      origin_[i] = coords.coordinate(i,0);
+    }
+  }
+
+  template <class Writer>
+  void writeLocalPieceImpl (Writer const& writer) const
+  {
+    auto const& gl = *gridView_.grid().begin(level_);
+    auto const& g = gl.interior[0];
+    auto extent = filledArray<6,int>(0);
+    auto const& gc = *g.dataBegin();
+    for (int i = 0; i < dim; ++i) {
+      extent[2*i] = gc.min(i);
+      extent[2*i+1] = gc.max(i)+1;
+    }
+
+#if HAVE_MPI
+    if (!localPieceCalled_) {
+      int rank = -1;
+      MPI_Comm_rank(gridView_.comm(), &rank);
+      if (rank != 0) {
+        MPI_Isend(extent.data(), extent.size(), MPI_INT, 0, /*tag=*/6, gridView_.comm(), &sendRequest_);
+      } else {
+        extents_[0] = extent;
+      }
+      localPieceCalled_ = true;
+    }
+#endif
+
+    writer(extent);
+  }
+
+  template <class Writer>
+  void writePiecesImpl (Writer const& writer) const
+  {
+    assert(localPieceCalled_);
+#if HAVE_MPI
+    int num_ranks = -1;
+    MPI_Comm_size(gridView_.comm(), &num_ranks);
+
+    std::vector<MPI_Status> status(num_ranks-1);
+    MPI_Waitall(num_ranks-1, requests_.data()+1, status.data());
+
+    for (int p = 0; p < num_ranks; ++p) {
+      writer(p, extents_[p], true);
+    }
+#endif
+  }
+
+  template <class T>
+  std::array<std::vector<T>, 3> coordinatesImpl () const
+  {
+    auto it = gridView_.grid().begin(level_);
+    auto const& coords = it->coords;
+
+    std::array<std::vector<T>, 3> ordinates{};
+    writeLocalPieceImpl([&ordinates,&coords](auto const& extent)
+    {
+      for (std::size_t d = 0; d < dim; ++d) {
+        auto s = extent[2*d+1] - extent[2*d] + 1;
+        ordinates[d].resize(s);
+        for (std::size_t i = 0; i < s; ++i)
+          ordinates[d][i] = coords.coordinate(d, extent[2*d] + i);
+      }
+    });
+
+    for (std::size_t d = dim; d < 3; ++d)
+      ordinates[d].resize(1, T(0));
+
+    return ordinates;
+  }
+
+private:
+  std::array<int, 6> wholeExtent_;
+  FieldVector<ctype,3> spacing_;
+  FieldVector<ctype,3> origin_;
+  std::vector<std::size_t> indexMap_;
+  int level_;
+
+  mutable std::vector<std::array<int,6>> extents_;
+  mutable std::vector<MPI_Request> requests_;
+  mutable MPI_Request sendRequest_;
+  mutable bool localPieceCalled_ = false;
+};
+
+namespace Impl
+{
+  template<class GridView, int dim, class Coordinates>
+  struct StructuredDataCollector<GridView, YaspGrid<dim,Coordinates>>
+  {
+    using type = YaspDataCollector<GridView>;
+  };
+}
+
+}} // end namespace Dune::experimental

dune/vtk/gridcreator.hh

@@ -28,9 +28,7 @@ namespace Dune { namespace experimental
 
       std::size_t idx = 0;
       for (std::size_t i = 0; i < types.size(); ++i) {
-        if (Vtk::Map::from_type.count(types[i]) == 0)
-          DUNE_THROW(Exception, "Unknown ElementType: " << types[i]);
-        auto type = Vtk::Map::from_type[types[i]];
+        auto type = Vtk::to_geometry(types[i]);
         Vtk::CellType cellType{type};
         auto refElem = referenceElement<double,Grid::dimension>(type);
 
@@ -92,9 +90,7 @@ namespace Dune { namespace experimental
 
       idx = 0;
       for (std::size_t i = 0; i < types.size(); ++i) {
-        if (Vtk::Map::from_type.count(types[i]) == 0)
-          DUNE_THROW(Exception, "Unknown ElementType: " << types[i]);
-        auto type = Vtk::Map::from_type[types[i]];
+        auto type = Vtk::to_geometry(types[i]);
         Vtk::CellType cellType{type};
 
         std::size_t nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]);

dune/vtk/utility/string.hh

@@ -3,6 +3,7 @@
 #include <algorithm>
 #include <cctype>
 #include <locale>
+#include <sstream>
 #include <string>
 
 namespace Dune
@@ -74,8 +75,8 @@ namespace Dune
     }
   }
 
-  template <class InputIter, class SeparaterIter, class Func>
-  void split(InputIter first, InputIter end, SeparaterIter s_first, SeparaterIter s_end, Func f)
+  template <class InputIter, class SeparatorIter, class Func>
+  void split(InputIter first, InputIter end, SeparatorIter s_first, SeparatorIter s_end, Func f)
   {
     if (first == end)
       return;
@@ -102,4 +103,18 @@ namespace Dune