discontinuousdatacollector.hh 3.78 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
#pragma once

3
#include "unstructureddatacollector.hh"
Praetorius, Simon's avatar
Praetorius, Simon committed
4

5
namespace Dune
Praetorius, Simon's avatar
Praetorius, Simon committed
6
7
8
{

/// Implementation of \ref DataCollector for linear cells, with discontinuous data.
9
template <class GridView, class Partition>
Praetorius, Simon's avatar
Praetorius, Simon committed
10
class DiscontinuousDataCollector
11
    : public UnstructuredDataCollectorInterface<GridView, DiscontinuousDataCollector<GridView,Partition>, Partition>
Praetorius, Simon's avatar
Praetorius, Simon committed
12
13
{
  using Self = DiscontinuousDataCollector;
14
15
16
17
18
  using Super = UnstructuredDataCollectorInterface<GridView, Self, Partition>;

public:
  using Super::dim;
  using Super::partition;
Praetorius, Simon's avatar
Praetorius, Simon committed
19
20
21
22

public:
  DiscontinuousDataCollector (GridView const& gridView)
    : Super(gridView)
23
  {}
Praetorius, Simon's avatar
Praetorius, Simon committed
24
25
26
27
28

  /// Create an index map the uniquely assignes an index to each pair (element,corner)
  void updateImpl ()
  {
    numPoints_ = 0;
29
    numCells_ = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
30
31
32
    indexMap_.resize(gridView_.size(dim));
    std::int64_t vertex_idx = 0;
    auto const& indexSet = gridView_.indexSet();
33
34
    for (auto const& c : elements(gridView_, partition)) {
      numCells_++;
Praetorius, Simon's avatar
Praetorius, Simon committed
35
      numPoints_ += c.subEntities(dim);
36
      for (unsigned int i = 0; i < c.subEntities(dim); ++i)
Praetorius, Simon's avatar
Praetorius, Simon committed
37
38
39
40
41
42
43
44
45
46
47
48
49
50
        indexMap_[indexSet.subIndex(c, i, dim)] = vertex_idx++;
    }
  }

  /// The number of pointsi approx. #cell * #corners-per-cell
  std::uint64_t numPointsImpl () const
  {
    return numPoints_;
  }

  /// Return the coordinates of the corners of all cells
  template <class T>
  std::vector<T> pointsImpl () const
  {
51
    std::vector<T> data(numPoints_ * 3);
Praetorius, Simon's avatar
Praetorius, Simon committed
52
    auto const& indexSet = gridView_.indexSet();
53
    for (auto const& element : elements(gridView_, partition)) {
54
      for (unsigned int i = 0; i < element.subEntities(dim); ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
55
56
57
58
59
60
61
62
63
64
65
        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;
  }

66
67
68
  /// Return number of grid cells
  std::uint64_t numCellsImpl () const
  {
69
    return numCells_;
70
71
  }

Praetorius, Simon's avatar
Praetorius, Simon committed
72
73
74
75
  /// Connect the corners of each cell. The leads to a global discontinuous grid
  Cells cellsImpl () const
  {
    Cells cells;
76
    cells.connectivity.reserve(numPoints_);
77
78
    cells.offsets.reserve(numCells_);
    cells.types.reserve(numCells_);
Praetorius, Simon's avatar
Praetorius, Simon committed
79
80
81

    std::int64_t old_o = 0;
    auto const& indexSet = gridView_.indexSet();
82
    for (auto const& c : elements(gridView_, partition)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
83
      Vtk::CellType cellType(c.type());
84
      for (unsigned int j = 0; j < c.subEntities(dim); ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
85
86
87
88
89
90
91
92
93
94
95
96
97
98
        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;
  }

  /// Evaluate the `fct` in the corners of each cell
  template <class T, class GlobalFunction>
  std::vector<T> pointDataImpl (GlobalFunction const& fct) const
  {
99
    std::vector<T> data(numPoints_ * fct.ncomps());
Praetorius, Simon's avatar
Praetorius, Simon committed
100
101
    auto const& indexSet = gridView_.indexSet();
    auto localFct = localFunction(fct);
102
    for (auto const& e : elements(gridView_, partition)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
103
104
105
      localFct.bind(e);
      Vtk::CellType cellType{e.type()};
      auto refElem = referenceElement(e.geometry());
106
      for (unsigned int j = 0; j < e.subEntities(dim); ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
107
108
109
110
111
112
113
114
115
        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;
  }

116
117
protected:
  using Super::gridView_;
118
  std::uint64_t numCells_ = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
119
120
121
122
  std::uint64_t numPoints_ = 0;
  std::vector<std::int64_t> indexMap_;
};

123
} // end namespace Dune