continuousdatacollector.hh 3.32 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
6
7
8
9

namespace Dune { namespace experimental
{

/// Implementation of \ref DataCollector for linear cells, with continuous data.
template <class GridView>
10
11
class ContinuousDataCollector
    : public UnstructuredDataCollectorInterface<GridView, ContinuousDataCollector<GridView>>
Praetorius, Simon's avatar
Praetorius, Simon committed
12
13
14
{
  enum { dim = GridView::dimension };

15
16
  using Self = ContinuousDataCollector;
  using Super = UnstructuredDataCollectorInterface<GridView, Self>;
Praetorius, Simon's avatar
Praetorius, Simon committed
17
18

public:
19
  ContinuousDataCollector (GridView const& gridView)
Praetorius, Simon's avatar
Praetorius, Simon committed
20
21
22
23
24
25
26
27
28
29
30
31
32
    : Super(gridView)
  {}

  /// Return number of grid vertices
  std::uint64_t numPointsImpl () const
  {
    return gridView_.size(dim);
  }

  /// Return the coordinates of all grid vertices in the order given by the indexSet
  template <class T>
  std::vector<T> pointsImpl () const
  {
33
    std::vector<T> data(gridView_.size(dim) * 3);
Praetorius, Simon's avatar
Praetorius, Simon committed
34
    auto const& indexSet = gridView_.indexSet();
35
    for (auto const& vertex : vertices(gridView_, Partitions::all)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
36
37
38
39
40
41
42
43
44
45
      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;
  }

46
47
48
49
50
51
  /// Return number of grid cells
  std::uint64_t numCellsImpl () const
  {
    return gridView_.size(0);
  }

Praetorius, Simon's avatar
Praetorius, Simon committed
52
53
54
55
56
57
58
59
60
61
62
63
  /// 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;
64
65
66
    cells.connectivity.reserve(gridView_.size(0) * maxVertices);
    cells.offsets.reserve(gridView_.size(0));
    cells.types.reserve(gridView_.size(0));
Praetorius, Simon's avatar
Praetorius, Simon committed
67
68

    std::int64_t old_o = 0;
69
    for (auto const& c : elements(gridView_, Partitions::all)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
70
      Vtk::CellType cellType(c.type());
71
      for (unsigned int j = 0; j < c.subEntities(dim); ++j)
Praetorius, Simon's avatar
Praetorius, Simon committed
72
73
74
75
76
77
78
79
80
81
82
        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;
  }

  /// Evaluate the `fct` at the corners of the elements
  template <class T, class GlobalFunction>
  std::vector<T> pointDataImpl (GlobalFunction const& fct) const
  {
83
    std::vector<T> data(gridView_.size(dim) * fct.ncomps());
Praetorius, Simon's avatar
Praetorius, Simon committed
84
85
    auto const& indexSet = gridView_.indexSet();
    auto localFct = localFunction(fct);
86
    for (auto const& e : elements(gridView_, Partitions::all)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
87
88
89
      localFct.bind(e);
      Vtk::CellType cellType{e.type()};
      auto refElem = referenceElement(e.geometry());
90
      for (unsigned int j = 0; j < e.subEntities(dim); ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
91
92
93
94
95
96
97
98
        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;
  }
99
100
101

protected:
  using Super::gridView_;
Praetorius, Simon's avatar
Praetorius, Simon committed
102
103
104
};

}} // end namespace Dune::experimental