quadraticdatacollector.hh 5.44 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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109

      /// \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)) {
          Vtk::CellType cellType(c.type(), Vtk::QUADRATIC);
          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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

      /// 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);
        for (auto const& e : elements(gridView_, partition)) {
          localFct.bind(e);
          Vtk::CellType cellType{e.type(), Vtk::QUADRATIC};
          auto refElem = referenceElement(e.geometry());
          for (unsigned int j = 0; j < e.subEntities(dim); ++j) {
            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)));
          }
          for (unsigned int j = 0; j < e.subEntities(dim-1); ++j) {
            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;
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