lagrangegridfunction.hh 9.82 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
51
52
53
54
55
56
57
58
#pragma once

#include <optional>
#include <vector>

#include <dune/common/dynvector.hh>
#include <dune/localfunctions/lagrange.hh>
#include <dune/vtk/gridfunctions/common.hh>
#include <dune/vtk/gridfunctions/continuousgridfunction.hh>
#include <dune/vtk/utility/errors.hh>
#include <dune/vtk/utility/lagrangepoints.hh>

namespace Dune
{
  namespace Vtk
  {
    template <class GridType>
    class LagrangeGridCreator;

    /// \brief Grid-function representing values from a VTK file with local Lagrange
    /// interpolation of the values stored on the Lagrange nodes.
    template <class GridType, class FieldType, class Context>
    class LagrangeGridFunction
    {
      using Grid = GridType;
      using Field = FieldType;

      using Factory = GridFactory<Grid>;
      using GridCreator = LagrangeGridCreator<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
      {
        using LFE = LagrangeLocalFiniteElement<LagrangePointSet, LC::dimension, FieldType, FieldType>;
        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:
        /// Constructor. Stores references to the passed data.
59
60
61
62
        PointDataLocalFunction (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)
63
64
65
66
67
68
69
70
          : creator_(creator)
          , values_(values)
          , comp_(comp)
          , types_(types)
          , offsets_(offsets)
          , connectivity_(connectivity)
        {}

71
72
        PointDataLocalFunction () = default;

73
74
75
76
77
78
79
80
81
        /// Binding the local-function to an element.
        /**
         * Constructs a new local finite-element with a polynomial order given
         * by the number of Lagrange nodes on the element. Extracts values on all
         * Lagrange nodes stored in file and stores these local coefficients in
         * a class member variable.
         **/
        void bind (LocalContext const& element)
        {
82
          unsigned int insertionIndex = creator_->factory().insertionIndex(element);
83

84
85
          std::int64_t shift = (insertionIndex == 0 ? 0 : (*offsets_)[insertionIndex-1]);
          std::int64_t numNodes = (*offsets_)[insertionIndex] - shift;
86
87
88
          [[maybe_unused]] std::int64_t maxNumNodes = numLagrangePoints(element.type().id(), element.type().dim(), 20);
          VTK_ASSERT(numNodes > 0 && numNodes < maxNumNodes);

89
          int order = creator_->order(element.type(), numNodes);
90
91
92
93
94
          VTK_ASSERT(order > 0 && order < 20);

          // construct a local finite-element with the corresponding order and Lagrange points
          // as stored in the file
          lfe_.emplace(LFE{element.type(), (unsigned int)(order)});
95
          lgeo_.emplace(creator_->localGeometry(element));
96
97
98

          // collect values on lagrange nodes
          localValues_.resize(numNodes);
99
100
          for (std::int64_t i = shift, i0 = 0; i < (*offsets_)[insertionIndex]; ++i, ++i0) {
            std::int64_t idx = (*connectivity_)[i];
101
102
103
            DynamicVector<Field>& v = localValues_[i0];
            v.resize(comp_);
            for (unsigned int j = 0; j < comp_; ++j)
104
              v[j] = (*values_)[comp_*idx + j];
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
          }
        }

        /// Unbind the local-function and the local finite-element from the element.
        void unbind ()
        {
          lfe_.reset();
          lgeo_.reset();
        }

        /// Evaluation in element local coordinates. Essentially, a local Lagrange
        /// interpolation with coefficients extracted in \ref bind().
        // NOTE: do we need to transform the local coordinates?
        Range operator() (Domain const& local) const
        {
          assert(!!lfe_);
          auto const& lb = lfe_->localBasis();
          lb.evaluateFunction(lgeo_->global(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:
133
134
        GridCreator const* creator_ = nullptr;
        std::vector<Field> const* values_ = nullptr;
135
        unsigned int comp_;
136
137
138
        std::vector<std::uint8_t> const* types_ = nullptr;
        std::vector<std::int64_t> const* offsets_ = nullptr;
        std::vector<std::int64_t> const* connectivity_ = nullptr;
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160

        // Local Finite-Element
        std::optional<LFE> lfe_ = std::nullopt;
        std::optional<typename GridCreator::LocalGeometry> lgeo_ = std::nullopt;

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

      /// Evaluation of data on the cells of the grid
      template <class LC>
      class CellDataLocalFunction
      {
      public:
        using LocalContext = LC;
        using Domain = typename LC::Geometry::LocalCoordinate;
        using Range = DynamicVector<Field>;
        using Signature = Range(Domain);

      public:
        /// Constructor. Stores references to the passed data.
161
162
163
164
165
        CellDataLocalFunction (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*/)
          : creator_(creator)
166
167
168
169
          , values_(values)
          , comp_(comp)
        {}

170
171
        CellDataLocalFunction () = default;

172
173
174
175
        /// Binding the local-function to an element extract the cell-value from the vector
        /// of data.
        void bind (LocalContext const& element)
        {
176
          unsigned int idx = creator_->factory().insertionIndex(element);
177
178
179
180
181
182

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

          for (unsigned int j = 0; j < comp_; ++j)
183
            v[j] = (*values_)[comp_*idx + j];
184
185
186
187
188
189
190
191
192
193
194
195
196
197
        }

        /// Unbinds from the bound element. Does nothing
        void unbind ()
        {}

        /// Evaluation in local element coordinates. Returns the constant value
        /// extracted in \ref bind().
        Range operator() (Domain const& local) const
        {
          return localValue_;
        }

      private:
198
199
        GridCreator const* creator_ = nullptr;
        std::vector<Field> const* values_ = nullptr;
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
        unsigned int comp_;

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

      /// Type switch for local-functions depending on the Context
      template <class LC>
      using LocalFunction = std::conditional_t< std::is_same_v<Context,PointContext>,
        PointDataLocalFunction<LC>,
        CellDataLocalFunction<LC>>;

    public:
      /// Construct a grid-function. Passed in data is stroed by reference, thus must have
      /// a life-time greater than that of the grid-function and corresponding local-function.
215
216
      LagrangeGridFunction (GridCreator const& creator, std::vector<Field> const& values,
                            std::string name ,unsigned int ncomps, Vtk::DataTypes dataType,
217
218
219
                            std::vector<std::uint8_t> const& types,
                            std::vector<std::int64_t> const& offsets,
                            std::vector<std::int64_t> const& connectivity)
220
221
        : creator_(&creator)
        , values_(&values)
222
223
224
        , name_(std::move(name))
        , ncomps_(ncomps)
        , dataType_(dataType)
225
226
227
        , types_(&types)
        , offsets_(&offsets)
        , connectivity_(&connectivity)
228
229
      {}

230
231
      LagrangeGridFunction () = default;

232
233
234
235
      /// Global evaluation. Not supported!
      Range operator() (Domain const& global) const
      {
        DUNE_THROW(Dune::NotImplemented, "Evaluation in global coordinates not implemented.");
236
        return Range(ncomps_, 0);
237
238
239
240
241
242
243
244
      }

      /// Return a type that defines the element that can be iterated.
      EntitySet const& entitySet () const
      {
        return entitySet_;
      }

245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
      std::string const& name () const
      {
        return name_;
      }

      int numComponents () const
      {
        return ncomps_;
      }

      Vtk::DataTypes dataType () const
      {
        return dataType_;
      }

260
261
262
263
      /// Construct a local-function depending on the Context type either PointDataLocalFunction
      /// or CellDataLocalFunction
      friend LocalFunction<typename EntitySet::Element> localFunction (LagrangeGridFunction const& gf)
      {
264
        return {gf.creator_, gf.values_, gf.ncomps_, gf.types_, gf.offsets_, gf.connectivity_};
265
266
267
      }

    private:
268
269
      GridCreator const* creator_ = nullptr;
      std::vector<Field> const* values_ = nullptr;
270
271
272
      std::string name_ = "GridFunction";
      unsigned int ncomps_ = 0;
      Vtk::DataTypes dataType_ = Vtk::DataTypes::UNKNOWN;
273
274
275
      std::vector<std::uint8_t> const* types_ = nullptr;
      std::vector<std::int64_t> const* offsets_ = nullptr;
      std::vector<std::int64_t> const* connectivity_ = nullptr;
276
277
278
279
280
281

      EntitySet entitySet_;
    };

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