discontinuousgridcreator.hh 3.1 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
3
4
5
6
7
8
9
10
11
#pragma once

#include <cassert>
#include <cstdint>
#include <limits>
#include <vector>

#include <dune/common/exceptions.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/grid/common/gridfactory.hh>

Stenger, Florian's avatar
Stenger, Florian committed
12
#include <dune/vtk/types.hh>
Praetorius, Simon's avatar
Praetorius, Simon committed
13
14
15
#include <dune/vtk/gridcreatorinterface.hh>
namespace Dune
{
Stenger, Florian's avatar
Stenger, Florian committed
16
  namespace Vtk
Praetorius, Simon's avatar
Praetorius, Simon committed
17
  {
Stenger, Florian's avatar
Stenger, Florian committed
18
19
20
21
22
    // Create a grid where the input points are not connected and the connectivity
    // describes separated elements.
    template <class Grid>
    struct DiscontinuousGridCreator
        : public GridCreatorInterface<Grid, DiscontinuousGridCreator<Grid>>
Praetorius, Simon's avatar
Praetorius, Simon committed
23
    {
Stenger, Florian's avatar
Stenger, Florian committed
24
25
26
27
28
      using Self = DiscontinuousGridCreator;
      using Super = GridCreatorInterface<Grid, DiscontinuousGridCreator>;
      using GlobalCoordinate = typename Super::GlobalCoordinate;

      struct CoordLess
Praetorius, Simon's avatar
Praetorius, Simon committed
29
      {
Stenger, Florian's avatar
Stenger, Florian committed
30
31
32
33
34
35
36
37
38
        template <class T, int N>
        bool operator() (FieldVector<T,N> const& lhs, FieldVector<T,N> const& rhs) const
        {
          for (int i = 0; i < N; ++i) {
            if (std::abs(lhs[i] - rhs[i]) < std::numeric_limits<T>::epsilon())
              continue;
            return lhs[i] < rhs[i];
          }
          return false;
Praetorius, Simon's avatar
Praetorius, Simon committed
39
        }
Stenger, Florian's avatar
Stenger, Florian committed
40
      };
Praetorius, Simon's avatar
Praetorius, Simon committed
41

Stenger, Florian's avatar
Stenger, Florian committed
42
43
44
      DiscontinuousGridCreator (GridFactory<Grid>& factory)
        : Super(factory)
      {}
Praetorius, Simon's avatar
Praetorius, Simon committed
45

Stenger, Florian's avatar
Stenger, Florian committed
46
47
48
49
50
51
52
      using Super::factory;
      void insertVerticesImpl (std::vector<GlobalCoordinate> const& points,
                              std::vector<std::uint64_t> const& /*point_ids*/)
      {
        points_ = &points;
        uniquePoints_.clear();
        std::size_t idx = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
53

Stenger, Florian's avatar
Stenger, Florian committed
54
55
56
57
58
59
        for (auto const& p : points) {
          auto b = uniquePoints_.emplace(std::make_pair(p,idx));
          if (b.second) {
            factory().insertVertex(p);
            ++idx;
          }
Praetorius, Simon's avatar
Praetorius, Simon committed
60
61
62
        }
      }

Stenger, Florian's avatar
Stenger, Florian committed
63
64
65
66
67
68
69
70
71
      void insertElementsImpl (std::vector<std::uint8_t> const& types,
                              std::vector<std::int64_t> const& offsets,
                              std::vector<std::int64_t> const& connectivity)
      {
        assert(points_ != nullptr);
        std::size_t idx = 0;
        for (std::size_t i = 0; i < types.size(); ++i) {
          auto type = Vtk::to_geometry(types[i]);
          Vtk::CellType cellType{type};
Praetorius, Simon's avatar
Praetorius, Simon committed
72

Stenger, Florian's avatar
Stenger, Florian committed
73
74
75
76
77
78
79
80
          int nNodes = offsets[i] - (i == 0 ? 0 : offsets[i-1]);
          assert(nNodes > 0);
          std::vector<unsigned int> vtk_cell; vtk_cell.reserve(nNodes);
          for (int j = 0; j < nNodes; ++j) {
            std::size_t v_j = connectivity[idx++];
            std::size_t new_idx = uniquePoints_[(*points_)[v_j]];
            vtk_cell.push_back(new_idx);
          }
Praetorius, Simon's avatar
Praetorius, Simon committed
81

Stenger, Florian's avatar
Stenger, Florian committed
82
83
84
85
86
87
88
89
90
          if (cellType.noPermutation()) {
            factory().insertElement(type,vtk_cell);
          } else {
            // apply index permutation
            std::vector<unsigned int> cell(nNodes);
            for (int j = 0; j < nNodes; ++j)
              cell[j] = vtk_cell[cellType.permutation(j)];
            factory().insertElement(type,cell);
          }
Praetorius, Simon's avatar
Praetorius, Simon committed
91
92
93
        }
      }

Stenger, Florian's avatar
Stenger, Florian committed
94
95
96
97
    private:
      std::vector<GlobalCoordinate> const* points_ = nullptr;
      std::map<GlobalCoordinate, std::size_t, CoordLess> uniquePoints_;
    };
Praetorius, Simon's avatar
Praetorius, Simon committed
98

Stenger, Florian's avatar
Stenger, Florian committed
99
  } // end namespace Vtk
Praetorius, Simon's avatar
Praetorius, Simon committed
100
} // end namespace Dune