DOFVectorView.inc.hpp 4.54 KB
Newer Older
1
2
#pragma once

3
#include <amdis/common/FieldMatVec.hpp>
4

5
6
namespace AMDiS {

7
8
template <class Traits, class TreePath>
typename DOFVectorConstView<Traits, TreePath>::Range DOFVectorConstView<Traits, TreePath>::
9
10
11
12
13
14
15
LocalFunction::operator()(LocalDomain const& x) const
{
  assert( bound_ );
  auto y = Range(0);

  auto&& coefficients = *globalFunction_->dofVector_;
  auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
16
  AMDiS::forEachLeafNode_(*subTree_, [&,this](auto const& node, auto)
17
18
19
20
21
22
23
24
25
26
27
28
  {
    using Node = std::decay_t<decltype(node)>;
    using LocalBasisRange = typename LocalFunction::template LocalBasisRange<Node>;

    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);
29
    using RangeBlock = std::decay_t<decltype(re)>;
30
31
32
33
34
35

    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];
36
      using CoefficientBlock = std::decay_t<decltype(c)>;
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

      // 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.
      std::size_t dimC = Flat<CoefficientBlock>::size(c);
      std::size_t 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;
}

Praetorius, Simon's avatar
Praetorius, Simon committed
59

60
61
template <class Traits, class TreePath>
typename DOFVectorConstView<Traits, TreePath>::DerivativeRange DOFVectorConstView<Traits, TreePath>::
62
GradientLocalFunction::operator()(LocalDomain const& x) const
63
{
64
65
66
67
68
69
70
  assert( bound_ );
  Range dy;
  for (std::size_t j = 0; j < dy.size(); ++j)
    dy[j] = 0;

  auto&& coefficients = *globalFunction_->dofVector_;
  auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
71
  AMDiS::forEachLeafNode_(*subTree_, [&,this](auto const& node, auto)
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
  {
    // TODO: may DOFVectorView::Range to FieldVector type if necessary
    using LocalDerivativeTraits = Dune::Functions::DefaultDerivativeTraits<Dune::FieldVector<double,1>(Domain)>;
    using GradientBlock = typename LocalDerivativeTraits::Range;

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

    // The transposed inverse Jacobian of the map from the reference element to the element
    auto&& jacobian = geometry_.value().jacobianInverseTransposed(x);

    auto&& referenceGradients = referenceGradientContainer_[node];
    localBasis.evaluateJacobian(x, referenceGradients);

    // Compute the shape function gradients on the real element
    std::vector<GradientBlock> gradients(referenceGradients.size());
    for (std::size_t i = 0; i < gradients.size(); ++i)
      multiplies_ABt(referenceGradients[i], jacobian, gradients[i]); // D[phi] * J^(-1) -> grad

    // Get range entry associated to this node
    auto&& re = nodeToRangeEntry(node, dy);
93
    using RangeBlock = std::decay_t<decltype(re)>;
94
95
96
97
98
99

    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];
100
      using CoefficientBlock = std::decay_t<decltype(c)>;
101
102
103

      // Get value of i-th transformed reference gradient
      auto&& grad = gradients[i];
104

105
106
107
108
      // Notice that the range entry re, the coefficient c, and the transformed
      // reference gradient grad may all be scalar, vector, matrix, or general container valued.
      std::size_t dimC = Flat<CoefficientBlock>::size(c);
      std::size_t dimV = Flat<GradientBlock>::size(grad);
109
      assert(dimC*dimV == std::size_t(Flat<RangeBlock>::size(re)));
110
111
112
113
114
115
      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<GradientBlock>::getEntry(grad, k);
          Flat<RangeBlock>::getEntry(re, j*dimV + k) += c_j*v_k;
        }
116
      }
117
118
    }
  });
119

120
121
  return dy;
}
122
123

} // end namespace AMDiS