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

3
4
5
6
#include <vector>

#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/partitionset.hh>
Stenger, Florian's avatar
Stenger, Florian committed
7
#include <dune/vtk/types.hh>
8

9
#include "unstructureddatacollector.hh"
Praetorius, Simon's avatar
Praetorius, Simon committed
10

11
namespace Dune
Praetorius, Simon's avatar
Praetorius, Simon committed
12
{
Stenger, Florian's avatar
Stenger, Florian committed
13
  namespace Vtk
Praetorius, Simon's avatar
Praetorius, Simon committed
14
  {
Stenger, Florian's avatar
Stenger, Florian committed
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
    /// Implementation of \ref DataCollector for quadratic cells, with continuous data.
    template <class GridView>
    class QuadraticDataCollector
        : public UnstructuredDataCollectorInterface<GridView, QuadraticDataCollector<GridView>, Partitions::All>
    {
      using Self = QuadraticDataCollector;
      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

    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);
Praetorius, Simon's avatar
Praetorius, Simon committed
36
37
      }

Stenger, Florian's avatar
Stenger, Florian committed
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
      /// 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();
        for (auto const& element : elements(gridView_, partition)) {
          auto geometry = element.geometry();
          auto refElem = referenceElement<T,dim>(element.type());

          // vertices
          for (unsigned int i = 0; i < element.subEntities(dim); ++i) {
            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
          for (unsigned int i = 0; i < element.subEntities(dim-1); ++i) {
            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;
Praetorius, Simon's avatar
Praetorius, Simon committed
72
      }
Stenger, Florian's avatar
Stenger, Florian committed
73
74
75
76
77

      /// Return number of grid cells
      std::uint64_t numCellsImpl () const
      {
        return gridView_.size(0);
Praetorius, Simon's avatar
Praetorius, Simon committed
78
      }
Stenger, Florian's avatar
Stenger, Florian committed
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

      /// \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();
        for (auto const& c : elements(gridView_, partition)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
95
          Vtk::CellType cellType(c.type(), Vtk::CellType::QUADRATIC);
Stenger, Florian's avatar
Stenger, Florian committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
          for (unsigned int j = 0; j < c.subEntities(dim); ++j) {
            int k = cellType.permutation(j);
            std::int64_t point_idx = indexSet.subIndex(c,k,dim);
            cells.connectivity.push_back(point_idx);
          }
          for (unsigned int j = 0; j < c.subEntities(dim-1); ++j) {
            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;
Praetorius, Simon's avatar
Praetorius, Simon committed
110
      }
Stenger, Florian's avatar
Stenger, Florian committed
111
112
113
114
115

      /// 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
      {
116
        std::vector<T> data(this->numPoints() * fct.numComponents());
Stenger, Florian's avatar
Stenger, Florian committed
117
118
119
120
        auto const& indexSet = gridView_.indexSet();
        auto localFct = localFunction(fct);
        for (auto const& e : elements(gridView_, partition)) {
          localFct.bind(e);
Praetorius, Simon's avatar
Praetorius, Simon committed
121
          Vtk::CellType cellType{e.type(), Vtk::CellType::QUADRATIC};
Stenger, Florian's avatar
Stenger, Florian committed
122
123
124
          auto refElem = referenceElement(e.geometry());
          for (unsigned int j = 0; j < e.subEntities(dim); ++j) {
            int k = cellType.permutation(j);
125
126
            std::size_t idx = fct.numComponents() * indexSet.subIndex(e, k, dim);
            for (int comp = 0; comp < fct.numComponents(); ++comp)
Stenger, Florian's avatar
Stenger, Florian committed
127
128
129
130
              data[idx + comp] = T(localFct.evaluate(comp, refElem.position(k, dim)));
          }
          for (unsigned int j = 0; j < e.subEntities(dim-1); ++j) {
            int k = cellType.permutation(e.subEntities(dim) + j);
131
132
            std::size_t idx = fct.numComponents() * (indexSet.subIndex(e, k, dim-1) + gridView_.size(dim));
            for (int comp = 0; comp < fct.numComponents(); ++comp)
Stenger, Florian's avatar
Stenger, Florian committed
133
134
135
136
137
              data[idx + comp] = T(localFct.evaluate(comp, refElem.position(k, dim-1)));
          }
          localFct.unbind();
        }
        return data;
Praetorius, Simon's avatar
Praetorius, Simon committed
138
      }
139

Stenger, Florian's avatar
Stenger, Florian committed
140
141
142
    protected:
      using Super::gridView_;
    };
Praetorius, Simon's avatar
Praetorius, Simon committed
143

Stenger, Florian's avatar
Stenger, Florian committed
144
  } // end namespace Vtk
145
} // end namespace Dune