added datacollector and quadratic vtk cell types

dune/vtk/datacollector.hh

+#pragma once
+
+#include <algorithm>
+#include <cstdint>
+#include <vector>
+
+#include "vtktypes.hh"
+
+namespace Dune { namespace experimental {
+
+struct Cells
+{
+  std::vector<std::uint8_t> types;
+  std::vector<std::int64_t> offsets;
+  std::vector<std::int64_t> connectivity;
+};
+
+class DataCollectorInterface
+{
+private:
+  DataCollectorInterface () = default;
+
+public:
+  /// \brief Return the number of points in the grid
+  std::uint64_t numPoints () const;
+
+  /// \brief Return the number of cells in the grid
+  std::uint64_t numCells () const;
+
+  /// Update the DataCollector on the current GridView
+  void update ();
+
+  /// \brief Return a flat vector of point coordinates
+  /**
+   * All coordinates are extended to 3 components and concatenated.
+   * [p0_x, p0_y, p0_z, p1_x, p1_y, p1_z, ...]
+   * If the GridView::dimensionworld < 3, the remaining components are
+   * set to 0
+   **/
+  template <class T>
+  std::vector<T> points () const;
+
+  /// \brief Return cell types, offsets, and connectivity. \see Cells
+  Cells cells () const;
+
+  /// \brief Return a flat vector of function values evaluated at the grid points.
+  /**
+   * In case of a vector valued function, flat the vector entries:
+   * [fct(p0)_0, fct(p0)_1, fct(p0)_2, fct(p1)_0, ...]
+   * where the vector dimension must be 3 (possible extended by 0s)
+   * In case of tensor valued function, flat the tensor row-wise:
+   * [fct(p0)_00, fct(p0)_01, fct(p0)_02, fct(p0)_10, fct(p0)_11, fct(p0)_12, fct(p0)_20...]
+   * where the tensor dimension must be 3x3 (possible extended by 0s)
+   **/
+  template <class T, class GlobalFunction>
+  std::vector<T> pointData (GlobalFunction const& fct) const;
+
+  /// \brief Return a flat vector of function values evaluated at the grid cells. \see pointData.
+  template <class T, class GlobalFunction>
+  std::vector<T> cellData (GlobalFunction const& fct) const;
+};
+
+
+/// Implementation of \ref DataCollector for linear cells, with continuous data.
+template <class GridView>
+class DefaultDataCollector
+{
+  enum { dim = GridView::dimension };
+
+public:
+  DefaultDataCollector (GridView const& gridView)
+    : gridView_(gridView)
+  {}
+
+  void update () {}
+
+  std::uint64_t numPoints () const
+  {
+    return gridView_.size(dim);
+  }
+
+  std::uint64_t numCells () const
+  {
+    return gridView_.size(0);
+  }
+
+  template <class T>
+  std::vector<T> points () const
+  {
+    std::vector<T> data(numPoints() * 3);
+    auto const& indexSet = gridView_.indexSet();
+    for (auto const& vertex : vertices(gridView_)) {
+      std::size_t idx = 3 * indexSet.index(vertex);
+      auto v = vertex.geometry().center();
+      for (std::size_t j = 0; j < v.size(); ++j)
+        data[idx + j] = T(v[j]);
+      for (std::size_t j = v.size(); j < 3u; ++j)
+        data[idx + j] = T(0);
+    }
+    return data;
+  }
+
+  Cells cells () const
+  {
+    auto const& indexSet = gridView_.indexSet();
+    auto types = indexSet.types(0);
+    int maxVertices = std::accumulate(types.begin(), types.end(), 1, [](int m, GeometryType t) {
+      auto refElem = referenceElement<double,dim>(t);
+      return std::max(m, refElem.size(dim));
+    });
+
+    Cells cells;
+    cells.connectivity.reserve(numCells() * maxVertices);
+    cells.offsets.reserve(numCells());
+    cells.types.reserve(numCells());
+
+    std::int64_t old_o = 0;
+    for (auto const& c : elements(gridView_)) {
+      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)) );
+      cells.offsets.push_back(old_o += c.subEntities(dim));
+      cells.types.push_back(cellType.type());
+    }
+
+    return cells;
+  }
+
+  template <class T, class GlobalFunction>
+  std::vector<T> pointData (GlobalFunction const& fct) const
+  {
+    std::vector<T> data(numPoints() * fct.ncomps());
+    auto const& indexSet = gridView_.indexSet();
+    auto localFct = localFunction(fct);
+    for (auto const& e : elements(gridView_)) {
+      localFct.bind(e);
+      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.permutation(j),dim);
+        for (int comp = 0; comp < fct.ncomps(); ++comp)
+          data[idx + comp] = T(localFct.evaluate(comp, refElem.position(cellType.permutation(j),dim)));
+      }
+      localFct.unbind();
+    }
+    return data;
+  }
+
+  template <class T, class GlobalFunction>
+  std::vector<T> cellData (GlobalFunction const& fct) const
+  {
+    std::vector<T> data(numCells() * fct.ncomps());
+    auto const& indexSet = gridView_.indexSet();
+    auto localFct = localFunction(fct);
+    for (auto const& e : elements(gridView_)) {
+      localFct.bind(e);
+      auto refElem = referenceElement(e.geometry());
+      std::size_t idx = fct.ncomps() * indexSet.index(e);
+      for (int comp = 0; comp < fct.ncomps(); ++comp)
+        data[idx + comp] = T(localFct.evaluate(comp, refElem.position(0,0)));
+      localFct.unbind();
+    }
+    return data;
+  }
+
+private:
+  GridView gridView_;
+};
+
+
+/// Implementation of \ref DataCollector for linear cells, with discontinuous data.
+template <class GridView>
+class DiscontinuousDataCollector
+{
+  enum { dim = GridView::dimension };
+
+public:
+  DiscontinuousDataCollector (GridView const& gridView)
+    : gridView_(gridView)
+  {}
+
+  void update ()
+  {
+    numPoints_ = 0;
+    indexMap_.resize(gridView_.size(dim));
+    std::int64_t vertex_idx = 0;
+    auto const& indexSet = gridView_.indexSet();
+    for (auto const& c : elements(gridView_)) {
+      numPoints_ += c.subEntities(dim);
+      for (int i = 0; i < c.subEntities(dim); ++i)
+        indexMap_[indexSet.subIndex(c, i, dim)] = vertex_idx++;
+    }
+  }
+
+  std::uint64_t numPoints () const
+  {
+    return numPoints_;
+  }
+
+  std::uint64_t numCells () const
+  {
+    return gridView_.size(0);
+  }
+
+  template <class T>
+  std::vector<T> points () const
+  {
+    std::vector<T> data(numPoints() * 3);
+    auto const& indexSet = gridView_.indexSet();
+    for (auto const& element : elements(gridView_)) {
+      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);
+        for (std::size_t j = 0; j < v.size(); ++j)
+          data[idx + j] = T(v[j]);
+        for (std::size_t j = v.size(); j < 3u; ++j)
+          data[idx + j] = T(0);
+      }
+    }
+    return data;
+  }
+
+  Cells cells () const
+  {
+    Cells cells;
+    cells.connectivity.reserve(numPoints());
+    cells.offsets.reserve(numCells());
+    cells.types.reserve(numCells());
+
+    std::int64_t old_o = 0;
+    auto const& indexSet = gridView_.indexSet();
+    for (auto const& c : elements(gridView_)) {
+      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)];
+        cells.connectivity.push_back(vertex_idx);
+      }
+      cells.offsets.push_back(old_o += c.subEntities(dim));
+      cells.types.push_back(cellType.type());
+    }
+
+    return cells;
+  }
+
+  template <class T, class GlobalFunction>
+  std::vector<T> pointData (GlobalFunction const& fct) const
+  {
+    std::vector<T> data(numPoints() * fct.ncomps());
+    auto const& indexSet = gridView_.indexSet();
+    auto localFct = localFunction(fct);
+    for (auto const& e : elements(gridView_)) {
+      localFct.bind(e);
+      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() * indexMap_[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.permutation(j),dim)));
+      }
+      localFct.unbind();
+    }
+    return data;
+  }
+
+  template <class T, class GlobalFunction>
+  std::vector<T> cellData (GlobalFunction const& fct) const
+  {
+    std::vector<T> data(numCells() * fct.ncomps());
+    auto const& indexSet = gridView_.indexSet();
+    auto localFct = localFunction(fct);
+    for (auto const& e : elements(gridView_)) {
+      localFct.bind(e);
+      auto refElem = referenceElement(e.geometry());
+      std::size_t idx = fct.ncomps() * indexSet.index(e);
+      for (int comp = 0; comp < fct.ncomps(); ++comp)
+        data[idx + comp] = T(localFct.evaluate(comp, refElem.position(0,0)));
+      localFct.unbind();
+    }
+    return data;
+  }
+
+private:
+  GridView gridView_;
+  std::uint64_t numPoints_ = 0;
+  std::vector<std::int64_t> indexMap_;
+};
+
+
+/// Implementation of \ref DataCollector for quadratic cells, with continuous data.
+template <class GridView>
+class QuadraticDataCollector
+{
+  enum { dim = GridView::dimension };
+
+public:
+  QuadraticDataCollector (GridView const& gridView)
+    : gridView_(gridView)
+  {}
+
+  void update () {}
+
+  std::uint64_t numPoints () const
+  {
+    return gridView_.size(dim) + gridView_.size(dim-1);
+  }
+
+  std::uint64_t numCells () const
+  {
+    return gridView_.size(0);
+  }
+
+  template <class T>
+  std::vector<T> points () const
+  {
+    std::vector<T> data(numPoints() * 3);
+    auto const& indexSet = gridView_.indexSet();
+    for (auto const& element : elements(gridView_)) {
+      auto geometry = element.geometry();
+      auto refElem = referenceElement<T,dim>(element.type());
+
+      // vertices
+      for (int i = 0; i < element.subEntities(dim); ++i) {
+        std::size_t idx = 3 * indexSet.subIndex(element, i, dim);
+        auto v = geometry.global(refElem.position(i,dim));
+        for (std::size_t j = 0; j < v.size(); ++j)
+          data[idx + j] = T(v[j]);
+        for (std::size_t j = v.size(); j < 3u; ++j)
+          data[idx + j] = T(0);
+      }
+      // edge centers
+      for (int i = 0; i < element.subEntities(dim-1); ++i) {
+        std::size_t idx = 3 * (indexSet.subIndex(element, i, dim-1) + gridView_.size(dim));
+        auto v = geometry.global(refElem.position(i,dim-1));
+        for (std::size_t j = 0; j < v.size(); ++j)
+          data[idx + j] = T(v[j]);
+        for (std::size_t j = v.size(); j < 3u; ++j)
+          data[idx + j] = T(0);
+      }
+    }
+    return data;
+  }
+
+  Cells cells () const
+  {
+    Cells cells;
+    cells.connectivity.reserve(numPoints());
+    cells.offsets.reserve(numCells());
+    cells.types.reserve(numCells());
+
+    std::int64_t old_o = 0;
+    auto const& indexSet = gridView_.indexSet();
+    for (auto const& c : elements(gridView_)) {
+      Vtk::CellType cellType(c.type(), Vtk::QUADRATIC);
+      for (int j = 0; j < c.subEntities(dim); ++j) {
+        int k = cellType.permutation(j);
+        std::int64_t point_idx = indexSet.subIndex(c,k,dim);
+        cells.connectivity.push_back(point_idx);
+      }
+      for (int j = 0; j < c.subEntities(dim-1); ++j) {
+        int k = cellType.permutation(c.subEntities(dim) + j);
+        std::int64_t point_idx = (indexSet.subIndex(c,k,dim-1) + gridView_.size(dim));
+        cells.connectivity.push_back(point_idx);
+      }
+      cells.offsets.push_back(old_o += c.subEntities(dim)+c.subEntities(dim-1));
+      cells.types.push_back(cellType.type());
+    }
+
+    return cells;
+  }
+
+  template <class T, class GlobalFunction>
+  std::vector<T> pointData (GlobalFunction const& fct) const
+  {
+    std::vector<T> data(numPoints() * fct.ncomps());
+    auto const& indexSet = gridView_.indexSet();
+    auto localFct = localFunction(fct);
+    for (auto const& e : elements(gridView_)) {
+      localFct.bind(e);
+      Vtk::CellType cellType{e.type(), Vtk::QUADRATIC};
+      auto refElem = referenceElement(e.geometry());
+      for (int j = 0; j < e.subEntities(dim); ++j) {
+        int k = cellType.permutation(j);
+        std::size_t idx = fct.ncomps() * indexSet.subIndex(e, k, dim);
+        for (int comp = 0; comp < fct.ncomps(); ++comp)
+          data[idx + comp] = T(localFct.evaluate(comp, refElem.position(k, dim)));
+      }
+      for (int j = 0; j < e.subEntities(dim-1); ++j) {
+        int k = cellType.permutation(e.subEntities(dim) + j);
+        std::size_t idx = fct.ncomps() * (indexSet.subIndex(e, k, dim-1) + gridView_.size(dim));
+        for (int comp = 0; comp < fct.ncomps(); ++comp)
+          data[idx + comp] = T(localFct.evaluate(comp, refElem.position(k, dim-1)));
+      }
+      localFct.unbind();
+    }
+    return data;
+  }
+
+  template <class T, class GlobalFunction>
+  std::vector<T> cellData (GlobalFunction const& fct) const
+  {
+    std::vector<T> data(numCells() * fct.ncomps());
+    auto const& indexSet = gridView_.indexSet();
+    auto localFct = localFunction(fct);
+    for (auto const& e : elements(gridView_)) {
+      localFct.bind(e);
+      auto refElem = referenceElement(e.geometry());
+      std::size_t idx = fct.ncomps() * indexSet.index(e);
+      for (int comp = 0; comp < fct.ncomps(); ++comp)
+        data[idx + comp] = T(localFct.evaluate(comp, refElem.position(0,0)));
+      localFct.unbind();
+    }
+    return data;
+  }
+
+private:
+  GridView gridView_;
+};
+
+}} // end namespace Dune::experimental

dune/vtk/filereader.hh

@@ -4,7 +4,7 @@
 #include <string>
 #include <utility>
 
-namespace Dune
+namespace Dune { namespace experimental
 {
   template <class Grid, class FilerReaderImp>
   class FileReader
@@ -68,5 +68,4 @@ namespace Dune
     }
   };
 
-
-} // end namespace Dune
+}} // end namespace Dune::experimental

dune/vtk/filewriter.hh

@@ -2,7 +2,7 @@
 
 #include <string>
 
-namespace Dune
+namespace Dune { namespace experimental
 {
   class FileWriter
   {
@@ -14,4 +14,4 @@ namespace Dune
     virtual void write (std::string const& filename) = 0;
   };
 
-} // end namespace Dune
+}} // end namespace Dune::experimental

dune/vtk/gridcreator.hh

@@ -10,7 +10,7 @@
 
 #include "vtktypes.hh"
 
-namespace Dune
+namespace Dune { namespace experimental
 {
   // Create a grid where the input points and connectivity is already
   // connected correctly.
@@ -28,9 +28,9 @@ namespace Dune
 
       std::size_t idx = 0;
       for (std::size_t i = 0; i < types.size(); ++i) {
-        if (Vtk::Map::type.count(types[i]) == 0)
+        if (Vtk::Map::from_type.count(types[i]) == 0)
           DUNE_THROW(Exception, "Unknown ElementType: " << types[i]);
-        auto type = Vtk::Map::type[types[i]];
+        auto type = Vtk::Map::from_type[types[i]];
         Vtk::CellType cellType{type};
         auto refElem = referenceElement<double,Grid::dimension>(type);
 
@@ -92,9 +92,9 @@ namespace Dune
 
       idx = 0;
       for (std::size_t i = 0; i < types.size(); ++i) {
-        if (Vtk::Map::type.count(types[i]) == 0)
+        if (Vtk::Map::from_type.count(types[i]) == 0)
           DUNE_THROW(Exception, "Unknown ElementType: " << types[i]);
-        auto type = Vtk::Map::type[types[i]];
+        auto type = Vtk::Map::from_type[types[i]];
         Vtk::CellType cellType{type};
 
         std::size_t nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]);
@@ -120,4 +120,4 @@ namespace Dune
     }
   };
 
-} // end namespace Dune
+}} // end namespace Dune::experimental

dune/vtk/vtkfunction.hh

@@ -5,7 +5,7 @@
 #include <dune/common/std/type_traits.hh>
 #include <dune/functions/common/typeerasure.hh>
 
-namespace Dune
+namespace Dune { namespace experimental
 {
   /// An abstract base class for LocalFunctions
   template <class GridView>
@@ -69,7 +69,7 @@ namespace Dune
       double evaluateImpl (int comp, LocalCoordinate const& xi, std::true_type) const
       {
         auto y = this->get()(xi);
-        return y[comp];
+        return comp < y.size() ? y[comp] : 0.0;
       }
 
       // Return the scalar values
@@ -128,7 +128,7 @@ namespace Dune
   {
   public:
     /// Create a local function
-    virtual VTKLocalFunction<GridView> makelocalFunction () const = 0;
+    virtual VTKLocalFunction<GridView> makeLocalFunction () const = 0;
 
     /// Virtual destructor
     virtual ~VTKFunctionInterface () = default;
@@ -143,12 +143,9 @@ namespace Dune
     {
     public:
       using Wrapper::Wrapper;
-      using Function = typename Wrapper::Wrapped;
-      using Interface = VTKFunctionInterface<GridView>;
-
-      virtual VTKLocalFunction<GridView> makelocalFunction () const override
+      virtual VTKLocalFunction<GridView> makeLocalFunction () const override
       {
-        return {localFunction(this->get())};
+        return VTKLocalFunction<GridView>{localFunction(this->get())};
       }
     };
   };
@@ -163,16 +160,17 @@ namespace Dune
                                              VTKFunctionImpl<GridView>::template Model>;
 
     template <class F>
-    using IsGridFunction = decltype(localFunction(std::declval<F>()));
+    using HasLocalFunction = decltype(localFunction(std::declval<F>()));
 
   public:
-    template <class F, disableCopyMove<VTKFunction, F> = 0>
+    template <class F>
     VTKFunction (F&& f, std::string name, int ncomps = 1)
       : Super(std::forward<F>(f))
       , name_(std::move(name))
-      , ncomps_(ncomps)
+      , ncomps_(ncomps > 1 ? 3 : 1)
     {
-      static_assert(Std::is_detected<IsGridFunction,F>::value,
+      assert(1 <= ncomps && ncomps <= 3);
+      static_assert(Std::is_detected<HasLocalFunction,F>::value,
         "Requires A GridFunction to be passed to the VTKFunction.");
     }
 
@@ -181,13 +179,7 @@ namespace Dune
     /// Create a LocalFunction
     friend VTKLocalFunction<GridView> localFunction (VTKFunction const& self)
     {
-      return self.asInterface().makelocalFunction();
-    }
-
-    /// Return the number of components of the Range
-    int ncomps () const
-    {
-      return ncomps_; // TODO: detect automatically. Is this always possible?
+      return self.asInterface().makeLocalFunction();
     }
 
     /// Return a name associated with the function
@@ -196,10 +188,16 @@ namespace Dune
       return name_;
     }
 
+    /// Return the number of components of the Range
+    int ncomps () const
+    {
+      return ncomps_;
+    }
+
   private:
     std::string name_;
 
     int ncomps_ = 1;
   };
 
-} // end namespace Dune
+}} // end namespace Dune::experimental

dune/vtk/vtkreader.hh

@@ -7,7 +7,7 @@
 #include "gridcreator.hh"
 #include "vtktypes.hh"
 
-namespace Dune
+namespace Dune { namespace experimental
 {
   /// File-Reader for Vtk .vtu files
   /**
@@ -129,6 +129,6 @@ namespace Dune
     std::uint64_t offset0_;
   };
 
-} // end namespace Dune
+}} // end namespace Dune::experimental
 
 #include "vtkreader.impl.hh"

dune/vtk/vtkreader.impl.hh

@@ -12,7 +12,7 @@
 #include "utility/filesystem.hh"
 #include "utility/string.hh"