discontinuousdatacollector.hh

#pragma once
 
#include <vector>
 
#include <dune/grid/common/partitionset.hh>
#include <dune/vtk/types.hh>
 
#include "unstructureddatacollector.hh"
 
namespace Dune::Vtk
{
  template <class GridView, class Partition = Partitions::InteriorBorder>
  class DiscontinuousDataCollector
      : public UnstructuredDataCollectorInterface<GridView, DiscontinuousDataCollector<GridView,Partition>, Partition>
  {
    using Self = DiscontinuousDataCollector;
    using Super = UnstructuredDataCollectorInterface<GridView, Self, Partition>;
 
  public:
    using Super::dim;
    using Super::partition;
    using Super::gridView;
 
  public:
    DiscontinuousDataCollector (GridView const& gridView)
      : Super(gridView)
    {}
 
    void updateImpl ()
    {
      numPoints_ = 0;
      numCells_ = 0;
      for (auto const& c : elements(gridView(), partition)) {
        numCells_++;
        numPoints_ += c.subEntities(dim);
      }
    }
 
    std::uint64_t numPointsImpl () const
    {
      return numPoints_;
    }
 
    template <class T>
    std::vector<T> pointsImpl () const
    {
      std::vector<T> data(numPoints_ * 3);
      std::size_t pointIdx = 0;
      for (auto const& element : elements(gridView(), partition)) {
        Vtk::CellType cellType(element.type());
        for (unsigned int i = 0; i < element.subEntities(dim); ++i,++pointIdx) {
          std::size_t idx = 3 * pointIdx;
          auto v = element.geometry().corner(cellType.permutation(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;
    }
 
    std::uint64_t numCellsImpl () const
    {
      return numCells_;
    }
 
    Cells cellsImpl () const
    {
      Cells cells;
      cells.connectivity.reserve(numPoints_);
      cells.offsets.reserve(numCells_);
      cells.types.reserve(numCells_);
 
      std::int64_t old_o = 0;
      std::size_t pointIdx = 0;
      for (auto const& c : elements(gridView(), partition)) {
        Vtk::CellType cellType(c.type());
        for (unsigned int j = 0; j < c.subEntities(dim); ++j,++pointIdx) {
          cells.connectivity.push_back(pointIdx);
        }
        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> pointDataImpl (GlobalFunction const& fct) const
    {
      std::vector<T> data(numPoints_ * fct.numComponents());
      auto localFct = localFunction(fct);
      std::size_t pointIdx = 0;
      for (auto const& e : elements(gridView(), partition)) {
        localFct.bind(e);
        Vtk::CellType cellType{e.type()};
        auto refElem = referenceElement(e.geometry());
        for (unsigned int j = 0; j < e.subEntities(dim); ++j,++pointIdx) {
          std::size_t idx = fct.numComponents() * pointIdx;
          for (int comp = 0; comp < fct.numComponents(); ++comp)
            data[idx + comp] = T(localFct.evaluate(comp, refElem.position(cellType.permutation(j),dim)));
        }
        localFct.unbind();
      }
      return data;
    }
 
  private:
    std::uint64_t numCells_ = 0;
    std::uint64_t numPoints_ = 0;
    std::vector<std::int64_t> indexMap_;
  };
 
} // end namespace Dune::Vtk