DOFVectorView.inc.hpp 2.03 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
#pragma once

namespace AMDiS
{

      {
        auto y = Range(0);

        auto const& coefficients = globalFunction_.coefficients();
        auto const& nodeToRangeEntry = globalFunction_.nodeToRangeEntry();
        forEachLeafNode(subTree_, [&,this](auto const& node, auto)
        {
          using Node = std::decay_t<decltype(node)>;
          using LocalBasisRange = typename LocalFunction::template LocalBasisRange<Node>;
          using MultiIndex = typename LocalIndexSet::MultiIndex;
          using CoefficientBlock = Range;
          using RangeBlock = std::decay_t<decltype(nodeToRangeEntry(node, y))>;

          auto&& fe = node.finiteElement();
          auto&& localBasis = fe.localBasis();

          auto&& shapeFunctionValues = shapeFunctionValueContainer_[node];
          localBasis.evaluateFunction(x, shapeFunctionValues);

          // Get range entry associated to this node
          auto&& re = nodeToRangeEntry(node, y);

          for (std::size_t i = 0; i < localBasis.size(); ++i) {
            auto&& multiIndex = localIndexSet_.index(node.localIndex(i));

            // Get coefficient associated to i-th shape function
            auto&& c = coefficients[multiIndex];

            // Get value of i-th shape function
            auto&& v = shapeFunctionValues[i];

            // Notice that the range entry re, the coefficient c, and the shape functions
            // value v may all be scalar, vector, matrix, or general container valued.
            auto dimC = Flat<CoefficientBlock>::size(c);
            auto dimV = Flat<LocalBasisRange>::size(v);
            assert(dimC*dimV == Flat<RangeBlock>::size(re));
            for(std::size_t j = 0; j < dimC; ++j) {
              auto&& c_j = Flat<CoefficientBlock>::getEntry(c, j);
              for(std::size_t k = 0; k < dimV; ++k) {
                auto&& v_k = Flat<LocalBasisRange>::getEntry(v, k);
                Flat<RangeBlock>::getEntry(re, j*dimV + k) += c_j*v_k;
              }
            }
          }
        });

        return y;
      }


} // end namespace AMDiS