VectorBase.hpp 15.1 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
#pragma once

#include <memory>
#include <string>
#include <type_traits>
#include <utility>

#include <dune/common/classname.hh>
#include <dune/common/deprecated.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/std/type_traits.hh>
#include <dune/functions/functionspacebases/concepts.hh>
#include <dune/functions/functionspacebases/sizeinfo.hh>

#include <amdis/Output.hpp>
#include <amdis/common/Concepts.hpp>
#include <amdis/common/FakeContainer.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/functions/NodeIndices.hpp>
#include <amdis/operations/Assigner.hpp>
#include <amdis/typetree/MultiIndex.hpp>

namespace AMDiS
{
  /// \brief The basic container that stores a base vector and a corresponding basis
  /**
29
30
   * A vector storing all the assembled Operators indexed with DOF indices. The
   * vector data is associated to a global basis.
31
32
33
34
35
36
37
38
39
40
41
42
43
44
   *
   * \tparam GB  Basis of the vector
   * \tparam B   A linear algebra backend implementing the storage and operations.
   **/
  template <class GB, class B>
  class VectorBase
  {
    using Self = VectorBase;

  public:
    /// The type of the functionspace basis associated to this linearform
    using GlobalBasis = GB;
    using LocalView = typename GlobalBasis::LocalView;

45
    /// The Linear-Algebra backend used to store the assembled coefficients
46
47
48
    using Backend = B;

  private:
49
    // proxy class redirecting the mutable element access to the insertValue method
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
    template <class Index>
    struct AccessProxy;

    enum class VectorState
    {
      unknown = 0,
      synchronized = 1,
      insert_values = 2,
      add_values = 3
    };

    template <class A>
    using VectorStateOf_t = std::conditional_t<std::is_same<A,Assigner::plus_assign>::value,
      std::integral_constant<VectorState, VectorState::add_values>,
      std::integral_constant<VectorState, VectorState::insert_values>>;

  public:
67
68
    /// (1) Constructor. Stores the shared_ptr of the basis.
    VectorBase(std::shared_ptr<GB> basis)
69
      : basis_(std::move(basis))
70
      , backend_(*basis_)
71
72
73
74
    {
      resizeZero();
    }

75
    /// (2) Constructor. Forwards to (1) by wrapping reference into non-destroying shared_ptr.
76
77
    template <class GB_,
      REQUIRES(Concepts::Similar<GB_,GB>)>
78
    VectorBase(GB_&& basis)
79
80
81
      : VectorBase(Dune::wrap_or_move(FWD(basis)))
    {}

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
110
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
138
139
140
141
142
143
144
145
146
147
148
149
150
    /// Return the basis \ref basis_ associated to the vector
    std::shared_ptr<GlobalBasis> const& basis() const { return basis_; }
    std::shared_ptr<GlobalBasis> const& basis()       { return basis_; }

    /// Return the underlying linear algebra backend
    Backend const& backend() const { return backend_; }
    Backend&       backend()       { return backend_; }


    template <class B_>
    using HasLocalSize = decltype(std::declval<B_>().localSize());

    template <class B_>
    using HasGlobalSize = decltype(std::declval<B_>().globalSize());

    /// Return the number of entries in the local part of the vector
    std::size_t localSize() const
    {
      return Dune::Hybrid::ifElse(Dune::Std::is_detected<HasLocalSize,Backend>{},
        [&](auto id) { return id(backend_).localSize(); },
        [&](auto id) { return id(backend_).size(); });
    }

    /// Return the number of entries in the global vector
    std::size_t globalSize() const
    {
      return Dune::Hybrid::ifElse(Dune::Std::is_detected<HasGlobalSize,Backend>{},
        [&](auto id) { return id(backend_).globalSize(); },
        [&](auto id) { return id(backend_).size(); });
    }

    /// Resize the \ref vector to the size of the \ref basis
    void resize()
    {
      backend_.init(sizeInfo(*basis_), false);
      state_ = VectorState::unknown;
    }

    /// Resize the \ref vector to the size of the \ref basis and set to zero
    void resizeZero()
    {
      backend_.init(sizeInfo(*basis_), true);
      state_ = VectorState::synchronized;
    }

    /// Prepare the vector for insertion of values, finish the insertion with
    /// \ref finish().
    void init(bool clear)
    {
      backend_.init(sizeInfo(*basis_), clear);
      state_ = clear ? VectorState::synchronized : VectorState::unknown;
    }

    /// Finish the insertion of values, see \ref init()
    void finish()
    {
      backend_.finish();
      state_ = VectorState::unknown;
    }

    /// Access the entry \p i of the \ref vector with read-only-access.
    template <class Index>
    typename Backend::value_type DUNE_DEPRECATED_MSG("Do not use element-wise vector access")
    operator[](Index const& idx) const
    {
      assert(state_ == VectorState::synchronized || state_ == VectorState::unknown);
      return at(idx);
    }

151
152
    /// Access the entry \p i of the \ref vector with write-access. Uses a proxy
    /// class the redirects to \ref insertValue.
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
    template <class Index>
    AccessProxy<Index> DUNE_DEPRECATED_MSG("Do not use element-wise vector access")
    operator[](Index const& idx)
    {
      return AccessProxy<Index>{this, idx};
    }

    /// Return the value of the vector at the local index \ref idx
    template <class Index>
    typename Backend::value_type at(Index const& idx) const
    {
      const_cast<Self*>(this)->synchronize();

      return backend_.at(flatMultiIndex(idx));
    }

    /// \brief Insert a single value into the matrix (add or overwrite to existing value)
    /**
171
172
173
174
     * Inserts or adds a value into a certain location \p dof (given as dof multi-index)
     * of a vector. The insertion mode is determined by the \p assign functor. Use
     * \ref Assigner::plus_assign for adding values (default) or \ref Assigner::assign
     * for overwriting (setting) values. Different insertion modes can not be mixed!
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
     *
     * Insertion must be closed with a call to \ref finish().
     *
     * [[not collective]]
     */
    template <class Index, class Assign = Assigner::plus_assign>
    void insert(Index const& dof, typename Backend::value_type const& value, Assign assign = {})
    {
      test_exit_dbg(state_ == VectorStateOf_t<Assign>::value ||
                    state_ == VectorState::unknown ||
                    state_ == VectorState::synchronized,
        "Vector in invalid state {} for insertion by {}.", to_string(state_), Dune::className<Assign>());

      backend_.insert(flatMultiIndex(dof), value, assign);

      // set the state to insert_values or add_values
      state_ = VectorStateOf_t<Assign>::value;
    }

    /// See \ref insert for assignment operation \ref Assigner::assign
    template <class Index>
    void set(Index const& dof, typename Backend::value_type const& value)
    {
      insert(dof, value, Assigner::assign{});
    }

    /// See \ref insert for assignment operation \ref Assigner::plus_assign
    template <class Index>
    void add(Index const& dof, typename Backend::value_type const& value)
    {
      insert(dof, value, Assigner::plus_assign{});
    }


209
210
    /// \brief Extract values from the vector referring to the given local indices
    /// and store it into a buffer
211
    /**
212
213
214
     * Collect value of indices and store them into a buffer. The buffer must be
     * a vector-like container with `buffer.resize()` and `buffer.begin()`. The
     * indices must be stored in an iterable container.
215
     *
216
217
     * If the vector is not in synchronized state, collects all necessary values
     * possibly from neighbouring processors.
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
     *
     * [[expects: localView is bound to an element]]
     * [[expects: node is in localView.tree()]]
     * [[possibly neighbor-wise collective operation]]
     */
    template <class Node, class Buffer,
      REQUIRES(Dune::models<Dune::Functions::Concept::BasisNode, Node>())>
    void gather(LocalView const& localView, Node const& node, Buffer& buffer) const
    {
      test_exit(state_ == VectorState::unknown ||
                state_ == VectorState::synchronized,
        "Vector in invalid state {} for gather operations.", to_string(state_));

      const_cast<Self*>(this)->synchronize();

      buffer.resize(node.size());
      backend_.gather(nodeIndices(localView, node), buffer.begin());
    }

    // [[expects: localView is bound to an element]]
    template <class Buffer>
    void gather(LocalView const& localView, Buffer& buffer) const
    {
      gather(localView, localView.tree(), buffer);
    }

    /// Insert a block of values into the vector (add or overwrite to existing values)
    /**
246
247
248
249
250
251
     * Inserts or adds values into certain locations of a vector. Insertion indices
     * are extracted from the given \p localView. The insertion mode is determined
     * by the \p assign functor. Use \ref Assigner::plus_assign for adding values
     * (default) or \ref Assigner::assign for overwriting (setting) values. Different
     * insertion modes can not be mixed! The \p localVector is assumed to be a continuous
     * memory container with a `data()` method to get a pointer to the beginning.
252
     *
253
254
255
     * The \p mask models a boolean range with at least a `begin()` method. Must
     * be forward iterable for at least `localVector.size()` elements. Does not
     * need an `end()` method. See, e.g. \ref FakeContainer.
256
     *
257
258
     * Insertion must be closed with a call to \ref finish(). It is not allowed
     * to switch insertion mode before calling `finish()`.
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
     *
     * [[expects: localView is bound to an element]]
     * [[expects: node is in localView.tree()]]
     * [[not collective]]
     */
    template <class Node, class NodeVector, class MaskRange, class Assign,
      REQUIRES(Dune::models<Dune::Functions::Concept::BasisNode, Node>())>
    void scatter(LocalView const& localView, Node const& node, NodeVector const& localVector,
                 MaskRange const& mask, Assign assign)
    {
      test_exit(state_ == VectorStateOf_t<Assign>::value ||
                state_ == VectorState::unknown ||
                state_ == VectorState::synchronized,
        "Vector in invalid state {} for insertion by {}.", to_string(state_), Dune::className<Assign>());

      assert(localVector.size() == node.size());

      // create a range of DOF indices on the node
      backend_.scatter(nodeIndices(localView, node), localVector, mask, assign);

      // set the state to insert_values or add_values
      state_ = VectorStateOf_t<Assign>::value;
    }

    /// Call \ref scatter with `MaskRange` set to \ref FakeContainer.
    // [[expects: localView is bound to an element]]
    // [[expects: node is in localView.tree()]]
    template <class Node, class NodeVector, class Assign,
      REQUIRES(Dune::models<Dune::Functions::Concept::BasisNode, Node>())>
    void scatter(LocalView const& localView, Node const& node, NodeVector const& localVector, Assign assign)
    {
      scatter(localView, node, localVector, FakeContainer<bool,true>{}, assign);
    }

    /// Call \ref scatter with `Node` given by the tree of the \p localView.
    // [[expects: localView is bound to an element]]
    template <class LocalVector, class Assign>
    void scatter(LocalView const& localView, LocalVector const& localVector, Assign assign)
    {
      scatter(localView, localView.tree(), localVector, assign);
    }

    /// Copies a block of values from \p vector to this
    /**
     * The values to copy might be masked with the range \p mask.
     * The assignment mode \p assign must be one of \ref Assigner::assign, or
     * \ref Assigner::plus_Assign.
     *
     * Requires the (local) size of the vector to be equal to the basis dimension.
     *
     * [[not collective]]
     **/
    template <class Vector, class MaskRange, class Assign>
    void copy(Vector const& vector, MaskRange const& mask, Assign assign)
    {
      test_exit(state_ == VectorStateOf_t<Assign>::value ||
                state_ == VectorState::unknown ||
                state_ == VectorState::synchronized,
        "Vector in invalid state {} for insertion by {}.", to_string(state_), Dune::className<Assign>());

      using size_type = typename Backend::size_type;
      size_type size = sizeInfo(*basis_);
      backend_.scatter(Dune::range(size), vector, mask, assign);

      // set the state to insert_values or add_values
      state_ = VectorStateOf_t<Assign>::value;
    }

    /// Call \ref copy with `MaskRange` set to \ref FakeContainer.
    template <class Vector, class Assign>
    void copy(Vector const& vector, Assign assign)
    {
      copy(vector, FakeContainer<bool,true>{}, assign);
    }

    /// Apply \p func to each value at given indices \p localInd
    /**
336
337
     * First, synchronizes the values of the vector, then applies the functor to
     * each local value associated to the local indices \p localInd.
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
     *
     * [[mutable]]
     **/
    template <class LocalInd, class Func>
    void forEach(LocalInd const& localInd, Func&& func)
    {
      synchronize();
      backend_.forEach(localInd, FWD(func));
    }

    /// Call \ref forEach for all entries in the vector.
    template <class Func>
    void forEach(Func&& func)
    {
      using size_type = typename Backend::size_type;
      size_type size = sizeInfo(*basis_);
      forEach(Dune::range(size), FWD(func));
    }

    /// Apply \p func to each value at given indices \p localInd
    /**
359
360
     * First, synchronizes the values of the vector, then applies the functor to
     * each local value associated to the local indices \p localInd.
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
     *
     * [[const]]
     **/
    template <class LocalInd, class Func>
    void forEach(LocalInd const& localInd, Func&& func) const
    {
      const_cast<Self*>(this)->synchronize();
      backend_.forEach(localInd, FWD(func));
    }

    /// Call \ref forEach for all entries in the vector.
    template <class Func>
    void forEach(Func&& func) const
    {
      using size_type = typename Backend::size_type;
      size_type size = sizeInfo(*basis_);
      forEach(Dune::range(size), FWD(func));
    }


  private:
    // synchronizes ghost values. Does not touch owned values
    void synchronize()
    {
      if (state_ != VectorState::synchronized)
        backend_.synchronize();

      state_ = VectorState::synchronized;
    }

    // print the vector state to string for debugging purposes
    static std::string to_string(VectorState state)
    {
      switch (state) {
        case VectorState::synchronized: return "synchronized";
        case VectorState::insert_values: return "insert_values";
        case VectorState::add_values: return "add_values";
        default: return "unknown";
      }
    }

  private:
    /// The finite element space / basis associated with the data vector
    std::shared_ptr<GlobalBasis> basis_;

    /// Data backend
    Backend backend_;

409
410
    /// The current state of the vector, one of {synchronized, insert_values,
    /// add_values, unknown}
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
    VectorState state_ = VectorState::unknown;
  };


  template <class GB, class B>
    template <class Index>
  struct VectorBase<GB,B>::AccessProxy
  {
    using value_type = typename Backend::value_type;

    void operator=(value_type const& value)
    {
      self->set(i, value);
    }

    void operator+=(value_type const& value)
    {
      self->add(i, value);
    }

    void operator-=(value_type const& value)
    {
      self->add(i, -value);
    }

    operator value_type() const
    {
      return self->at(i);
    }

    VectorBase* self;
    Index i;
  };

} // end namespace AMDiS