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

3
#include <numeric>
4
#include <vector>
Praetorius, Simon's avatar
Praetorius, Simon committed
5

6
#include <dune/geometry/referenceelements.hh>
7
8
#include <dune/grid/utility/globalindexset.hh>

9
10
11
12
13
#include <dune/vtk/forward.hh>
#include <dune/vtk/vtktypes.hh>

#include "unstructureddatacollector.hh"

14
namespace Dune
Praetorius, Simon's avatar
Praetorius, Simon committed
15
16
17
{

/// Implementation of \ref DataCollector for linear cells, with continuous data.
18
template <class GridView, class Partition>
19
class ContinuousDataCollector
20
    : public UnstructuredDataCollectorInterface<GridView, ContinuousDataCollector<GridView,Partition>, Partition>
Praetorius, Simon's avatar
Praetorius, Simon committed
21
{
22
  using Self = ContinuousDataCollector;
23
24
25
26
27
  using Super = UnstructuredDataCollectorInterface<GridView, Self, Partition>;

public:
  using Super::dim;
  using Super::partition;
Praetorius, Simon's avatar
Praetorius, Simon committed
28
29

public:
30
  ContinuousDataCollector (GridView const& gridView)
Praetorius, Simon's avatar
Praetorius, Simon committed
31
32
33
    : Super(gridView)
  {}

34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
  /// Collect the vertex indices
  void updateImpl ()
  {
    numPoints_ = 0;
    indexMap_.resize(gridView_.size(dim));
    auto const& indexSet = gridView_.indexSet();
    for (auto const& vertex : vertices(gridView_, partition))
      indexMap_[indexSet.index(vertex)] = std::int64_t(numPoints_++);

    if (gridView_.comm().size() > 1) {
      auto&& e = elements(gridView_, partition);
      numCells_ = std::distance(std::begin(e), std::end(e));
    } else {
      numCells_ = gridView_.size(0);
    }
  }

Praetorius, Simon's avatar
Praetorius, Simon committed
51
52
53
  /// Return number of grid vertices
  std::uint64_t numPointsImpl () const
  {
54
    return numPoints_;
Praetorius, Simon's avatar
Praetorius, Simon committed
55
56
57
58
59
60
  }

  /// Return the coordinates of all grid vertices in the order given by the indexSet
  template <class T>
  std::vector<T> pointsImpl () const
  {
61
62
63
    std::vector<T> data;
    data.reserve(numPoints_ * 3);
    for (auto const& vertex : vertices(gridView_, partition)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
64
65
      auto v = vertex.geometry().center();
      for (std::size_t j = 0; j < v.size(); ++j)
66
        data.emplace_back(v[j]);
Praetorius, Simon's avatar
Praetorius, Simon committed
67
      for (std::size_t j = v.size(); j < 3u; ++j)
68
        data.emplace_back(0);
Praetorius, Simon's avatar
Praetorius, Simon committed
69
70
71
72
    }
    return data;
  }

73
74
75
  /// Return a vector of global unique ids of the points
  std::vector<std::uint64_t> pointIdsImpl () const
  {
76
77
    std::vector<std::uint64_t> data;
    data.reserve(numPoints_);
78
    GlobalIndexSet<GridView> globalIndexSet(gridView_, dim);
79
80
    for (auto const& vertex : vertices(gridView_, partition)) {
      data.emplace_back(globalIndexSet.index(vertex));
81
82
83
84
    }
    return data;
  }

85
86
87
  /// Return number of grid cells
  std::uint64_t numCellsImpl () const
  {
88
    return numCells_;
89
90
  }

Praetorius, Simon's avatar
Praetorius, Simon committed
91
92
93
94
95
96
97
98
99
100
101
102
  /// Return the types, offsets and connectivity of the cells, using the same connectivity as
  /// given by the grid.
  Cells cellsImpl () 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;
103
104
105
    cells.connectivity.reserve(numCells_ * maxVertices);
    cells.offsets.reserve(numCells_);
    cells.types.reserve(numCells_);
Praetorius, Simon's avatar
Praetorius, Simon committed
106
107

    std::int64_t old_o = 0;
108
    for (auto const& c : elements(gridView_, partition)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
109
      Vtk::CellType cellType(c.type());
110
      for (unsigned int j = 0; j < c.subEntities(dim); ++j)
111
        cells.connectivity.emplace_back(indexMap_[indexSet.subIndex(c,cellType.permutation(j),dim)]);
Praetorius, Simon's avatar
Praetorius, Simon committed
112
113
114
115
116
117
118
119
120
121
      cells.offsets.push_back(old_o += c.subEntities(dim));
      cells.types.push_back(cellType.type());
    }
    return cells;
  }

  /// Evaluate the `fct` at the corners of the elements
  template <class T, class GlobalFunction>
  std::vector<T> pointDataImpl (GlobalFunction const& fct) const
  {
122
    std::vector<T> data(numPoints_ * fct.ncomps());
Praetorius, Simon's avatar
Praetorius, Simon committed
123
124
    auto const& indexSet = gridView_.indexSet();
    auto localFct = localFunction(fct);
125
    for (auto const& e : elements(gridView_, partition)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
126
127
128
      localFct.bind(e);
      Vtk::CellType cellType{e.type()};
      auto refElem = referenceElement(e.geometry());
129
      for (unsigned int j = 0; j < e.subEntities(dim); ++j) {
130
        std::size_t idx = fct.ncomps() * indexMap_[indexSet.subIndex(e,cellType.permutation(j),dim)];
Praetorius, Simon's avatar
Praetorius, Simon committed
131
132
133
134
135
136
137
        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;
  }
138

139
protected:
140
  using Super::gridView_;
141
142
143
  std::uint64_t numPoints_ = 0;
  std::uint64_t numCells_ = 0;
  std::vector<std::int64_t> indexMap_;
Praetorius, Simon's avatar
Praetorius, Simon committed
144
145
};

146
} // end namespace Dune