quadraticdatacollector.hh 4.81 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
9
10
{

/// Implementation of \ref DataCollector for quadratic cells, with continuous data.
template <class GridView>
class QuadraticDataCollector
11
    : public UnstructuredDataCollectorInterface<GridView, QuadraticDataCollector<GridView>, Partitions::All>
Praetorius, Simon's avatar
Praetorius, Simon committed
12
13
{
  using Self = QuadraticDataCollector;
14
15
16
17
18
  using Super = UnstructuredDataCollectorInterface<GridView, Self, Partitions::All>;

public:
  using Super::dim;
  using Super::partition; // NOTE: quadratic data-collector currently implemented for the All partition only
Praetorius, Simon's avatar
Praetorius, Simon committed
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

public:
  QuadraticDataCollector (GridView const& gridView)
    : Super(gridView)
  {}

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

  /// Return a vector of point coordinates.
  /**
   * The vector of point coordinates is composed of vertex coordinates first and second
   * edge center coordinates.
   **/
  template <class T>
  std::vector<T> pointsImpl () const
  {
    std::vector<T> data(this->numPoints() * 3);
    auto const& indexSet = gridView_.indexSet();
41
    for (auto const& element : elements(gridView_, partition)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
42
43
44
45
      auto geometry = element.geometry();
      auto refElem = referenceElement<T,dim>(element.type());

      // vertices
46
      for (unsigned int i = 0; i < element.subEntities(dim); ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
47
48
49
50
51
52
53
54
        std::size_t idx = 3 * indexSet.subIndex(element, i, dim);
        auto v = geometry.global(refElem.position(i,dim));
        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);
      }
      // edge centers
55
      for (unsigned int i = 0; i < element.subEntities(dim-1); ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
56
57
58
59
60
61
62
63
64
65
66
        std::size_t idx = 3 * (indexSet.subIndex(element, i, dim-1) + gridView_.size(dim));
        auto v = geometry.global(refElem.position(i,dim-1));
        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;
  }

67
68
69
70
71
72
  /// Return number of grid cells
  std::uint64_t numCellsImpl () const
  {
    return gridView_.size(0);
  }

Praetorius, Simon's avatar
Praetorius, Simon committed
73
74
75
76
77
78
79
80
81
82
83
84
85
86
  /// \brief Return cell types, offsets, and connectivity. \see Cells
  /**
   * The cell connectivity is composed of cell vertices first and second cell edges,
   * where the indices are grouped [vertex-indices..., (#vertices)+edge-indices...]
   **/
  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;
    auto const& indexSet = gridView_.indexSet();
87
    for (auto const& c : elements(gridView_, partition)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
88
      Vtk::CellType cellType(c.type(), Vtk::QUADRATIC);
89
      for (unsigned int j = 0; j < c.subEntities(dim); ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
90
91
92
93
        int k = cellType.permutation(j);
        std::int64_t point_idx = indexSet.subIndex(c,k,dim);
        cells.connectivity.push_back(point_idx);
      }
94
      for (unsigned int j = 0; j < c.subEntities(dim-1); ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
        int k = cellType.permutation(c.subEntities(dim) + j);
        std::int64_t point_idx = (indexSet.subIndex(c,k,dim-1) + gridView_.size(dim));
        cells.connectivity.push_back(point_idx);
      }
      cells.offsets.push_back(old_o += c.subEntities(dim)+c.subEntities(dim-1));
      cells.types.push_back(cellType.type());
    }
    return cells;
  }

  /// Evaluate the `fct` at element vertices and edge centers in the same order as the point coords.
  template <class T, class GlobalFunction>
  std::vector<T> pointDataImpl (GlobalFunction const& fct) const
  {
    std::vector<T> data(this->numPoints() * fct.ncomps());
    auto const& indexSet = gridView_.indexSet();
    auto localFct = localFunction(fct);
112
    for (auto const& e : elements(gridView_, partition)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
113
114
115
      localFct.bind(e);
      Vtk::CellType cellType{e.type(), Vtk::QUADRATIC};
      auto refElem = referenceElement(e.geometry());
116
      for (unsigned int j = 0; j < e.subEntities(dim); ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
117
118
119
120
121
        int k = cellType.permutation(j);
        std::size_t idx = fct.ncomps() * indexSet.subIndex(e, k, dim);
        for (int comp = 0; comp < fct.ncomps(); ++comp)
          data[idx + comp] = T(localFct.evaluate(comp, refElem.position(k, dim)));
      }
122
      for (unsigned int j = 0; j < e.subEntities(dim-1); ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
123
124
125
126
127
128
129
130
131
        int k = cellType.permutation(e.subEntities(dim) + j);
        std::size_t idx = fct.ncomps() * (indexSet.subIndex(e, k, dim-1) + gridView_.size(dim));
        for (int comp = 0; comp < fct.ncomps(); ++comp)
          data[idx + comp] = T(localFct.evaluate(comp, refElem.position(k, dim-1)));
      }
      localFct.unbind();
    }
    return data;
  }
132
133
134

protected:
  using Super::gridView_;
Praetorius, Simon's avatar
Praetorius, Simon committed
135
136
};

137
} // end namespace Dune