celltypedatacollector.hh

#pragma once
 
#include <vector>
 
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/partitionset.hh>
#include <dune/vtk/types.hh>
 
#include "unstructureddatacollector.hh"
 
namespace Dune::Vtk
{
  template <class GridView>
  class CellTypeDataCollector
      : public UnstructuredDataCollectorInterface<GridView, CellTypeDataCollector<GridView>, Partitions::All>
  {
    using Self = CellTypeDataCollector;
    using Super = UnstructuredDataCollectorInterface<GridView, Self, Partitions::All>;
 
    static auto nodeMapper (const GridView& gridView, Vtk::CellType::Parametrization param)
    {
      auto layout = [param](Dune::GeometryType gt, int dimgrid) -> unsigned int {
        auto cellType = Dune::Vtk::CellType{gt,param};
        return cellType.layout(gt.dim());
      };
      return Dune::MultipleCodimMultipleGeomTypeMapper(gridView, layout);
    }
 
    Vtk::CellType::Parametrization parametrization_;
    Dune::MultipleCodimMultipleGeomTypeMapper<GridView> mapper_;
 
  public:
    using Super::dim;
    using Super::partition; // NOTE: cell-type data-collector currently implemented for the All partition only
    using Super::gridView;
 
  public:
    CellTypeDataCollector (GridView const& gridView, Vtk::CellType::Parametrization param = Vtk::CellType::LINEAR)
      : Super(gridView)
      , parametrization_(param)
      , mapper_(nodeMapper(Super::gridView(), param))
    {}
 
    void updateImpl ()
    {
      mapper_.update(Super::gridView());
    }
 
    std::uint64_t numPointsImpl () const
    {
      return mapper_.size();
    }
 
 
    template <class T>
    std::vector<T> pointsImpl () const
    {
      std::vector<T> data(this->numPoints() * 3);
      for (auto const& c : elements(gridView(), partition)) {
        auto geometry = c.geometry();
        Vtk::CellType cellType(c.type(), parametrization_);
        auto refElem = referenceElement<T,dim>(c.type());
 
        for (int codim = dim, j = 0; codim >= 0 && j < cellType.size(); --codim) {
          for (int i = 0; i < refElem.size(codim) && j < cellType.size(); ++i,++j) {
            std::int64_t idx = 3 * mapper_.subIndex(c,i,codim);
            auto v = geometry.global(refElem.position(i,codim));
            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 gridView().size(0);
    }
 
 
    Cells cellsImpl () const
    {
      Cells cells;
      cells.connectivity.reserve(this->numPoints());
      cells.offsets.reserve(this->numCells());
      cells.types.reserve(this->numCells());
 
      std::int64_t old_o = 0;
      for (auto const& c : elements(gridView(), partition)) {
        Vtk::CellType cellType(c.type(), parametrization_);
        auto refElem = referenceElement(c);
 
        for (int codim = dim, j = 0; codim >= 0 && j < cellType.size(); --codim) {
          for (int i = 0; i < refElem.size(codim) && j < cellType.size(); ++i,++j) {
            int k = cellType.permutation(j);
            cells.connectivity.push_back(mapper_.subIndex(c,k,codim));
          }
        }
        cells.offsets.push_back(old_o += cellType.size());
        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(this->numPoints() * fct.numComponents());
      auto localFct = localFunction(fct);
      for (auto const& c : elements(gridView(), partition)) {
        localFct.bind(c);
        Vtk::CellType cellType{c.type(), parametrization_};
        auto refElem = referenceElement(c);
 
        for (int codim = dim, j = 0; codim >= 0 && j < cellType.size(); --codim) {
          for (int i = 0; i < refElem.size(codim) && j < cellType.size(); ++i,++j) {
            int k = cellType.permutation(j);
            std::int64_t idx = fct.numComponents() * mapper_.subIndex(c,k,codim);
            auto v = refElem.position(k,codim);
            for (int comp = 0; comp < fct.numComponents(); ++comp)
              data[idx + comp] = T(localFct.evaluate(comp, v));
          }
        }
        localFct.unbind();
      }
      return data;
    }
  };
 
} // end namespace Dune::Vtk