continuousgridfunction.hh 5.86 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#pragma once

#include <vector>

#include <dune/common/dynvector.hh>
#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
#include <dune/vtk/gridfunctions/common.hh>

namespace Dune
{
  namespace Vtk
  {
    /// \brief A GridFunction representing data stored on the grid vertices in a continuous manner.
    template <class GridType, class FieldType, class Context>
    class ContinuousGridFunction
    {
      using Grid = GridType;
      using Field = FieldType;

      using Factory = GridFactory<Grid>;

    public:
      struct EntitySet
      {
        using Grid = GridType;
        using Element = typename GridType::template Codim<0>::Entity;
        using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
        using GlobalCoordinate = typename Element::Geometry::GlobalCoordinate;
      };

      using Domain = typename EntitySet::GlobalCoordinate;
      using Range = DynamicVector<Field>;
      using Signature = Range(Domain);

    private:
      template <class LC>
      class PointDataLocalFunction
      {
        // NOTE: local finite-element fixed to Lagrange P1/Q1
        using LFECache = LagrangeLocalFiniteElementCache<Field,Field,LC::mydimension,1>;
        using LFE = typename LFECache::FiniteElementType;
        using LB = typename LFE::Traits::LocalBasisType;

      public:
        using LocalContext = LC;
        using Domain = typename LC::Geometry::LocalCoordinate;
        using Range = DynamicVector<Field>;
        using Signature = Range(Domain);

      public:
51
        PointDataLocalFunction (Factory const* factory, std::vector<Field> const* values, unsigned int comp)
52
53
54
55
56
          : factory_(factory)
          , values_(values)
          , comp_(comp)
        {}

57
58
        PointDataLocalFunction () = default;

59
60
61
62
63
64
65
66
        void bind (LocalContext const& element)
        {
          lfe_ = &cache_.get(element.type());

          // collect values on vertices
          // NOTE: assumes, that Lagrange nodes are ordered like element vertices
          localValues_.resize(element.subEntities(Grid::dimension));
          for (unsigned int i = 0; i < element.subEntities(Grid::dimension); ++i) {
67
            unsigned int idx = factory_->insertionIndex(element.template subEntity<Grid::dimension>(i));
68
69
70
            DynamicVector<Field>& v = localValues_[i];
            v.resize(comp_);
            for (unsigned int j = 0; j < comp_; ++j)
71
              v[j] = (*values_)[comp_*idx + j];
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
          }
        }

        void unbind ()
        {
          lfe_ = nullptr;
        }

        Range operator() (Domain const& local) const
        {
          assert(!!lfe_);
          auto const& lb = lfe_->localBasis();
          lb.evaluateFunction(local, shapeValues_);
          assert(shapeValues_.size() == localValues_.size());

          Range y(comp_, Field(0));
          for (std::size_t i = 0; i < shapeValues_.size(); ++i)
            y.axpy(shapeValues_[i], localValues_[i]);

          return y;
        }

      private:
95
96
        Factory const* factory_ = nullptr;
        std::vector<Field> const* values_ = nullptr;
97
98
99
        unsigned int comp_;

        // Local Finite-Element
100
        LFECache cache_;
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
        LFE const* lfe_ = nullptr;

        // cache of local values
        std::vector<DynamicVector<Field>> localValues_;
        mutable std::vector<typename LB::Traits::RangeType> shapeValues_;
      };

      template <class LC>
      class CellDataLocalFunction
      {
      public:
        using LocalContext = LC;
        using Domain = typename LC::Geometry::LocalCoordinate;
        using Range = DynamicVector<Field>;
        using Signature = Range(Domain);

      public:
118
        CellDataLocalFunction (Factory const* factory, std::vector<Field> const* values, unsigned int comp)
119
120
121
122
123
          : factory_(factory)
          , values_(values)
          , comp_(comp)
        {}

124
125
        CellDataLocalFunction () = default;

126
127
        void bind (LocalContext const& element)
        {
128
          unsigned int idx = factory_->insertionIndex(element);
129
130
131
132
133
134

          // collect values on cells
          DynamicVector<Field>& v = localValue_;
          v.resize(comp_);

          for (unsigned int j = 0; j < comp_; ++j)
135
            v[j] = (*values_)[comp_*idx + j];
136
137
138
139
140
141
142
143
144
145
146
        }

        void unbind ()
        {}

        Range operator() (Domain const& local) const
        {
          return localValue_;
        }

      private:
147
148
        Factory const* factory_ = nullptr;
        std::vector<Field> const* values_ = nullptr;
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
        unsigned int comp_;

        // cache of local values
        DynamicVector<Field> localValue_;
      };

      template <class LC>
      using LocalFunction = std::conditional_t< std::is_same_v<Context,PointContext>,
        PointDataLocalFunction<LC>,
        CellDataLocalFunction<LC>>;

    public:
      template <class GridCreator>
      ContinuousGridFunction (GridCreator const& creator, std::vector<Field> const& values, unsigned int comp,
                              std::vector<std::uint8_t> const& /*types*/,
                              std::vector<std::int64_t> const& /*offsets*/,
                              std::vector<std::int64_t> const& /*connectivity*/)
166
167
        : factory_(&creator.factory())
        , values_(&values)
168
169
170
        , comp_(comp)
      {}

171
172
      ContinuousGridFunction () = default;

173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
      Range operator() (Domain const& global) const
      {
        DUNE_THROW(Dune::NotImplemented, "Evaluation in global coordinates not implemented.");
        return Range(comp_, 0);
      }

      EntitySet const& entitySet () const
      {
        return entitySet_;
      }

      friend LocalFunction<typename EntitySet::Element> localFunction (ContinuousGridFunction const& gf)
      {
        return {gf.factory_, gf.values_, gf.comp_};
      }

    private:
190
191
192
      Factory const* factory_;
      std::vector<Field> const* values_ = nullptr;
      unsigned int comp_ = 0;
193
194
195
196
197
198

      EntitySet entitySet_;
    };

  } // end namespace Vtk
} // end namespace Dune