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

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

5
6
7
namespace AMDiS {

template <class GlobalBasis, class TreePath>
Praetorius, Simon's avatar
Praetorius, Simon committed
8
typename DOFVectorConstView<GlobalBasis, TreePath>::Range DOFVectorConstView<GlobalBasis, TreePath>::
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
59
LocalFunction::operator()(LocalDomain const& x) const
{
  assert( bound_ );
  auto y = Range(0);

  auto&& coefficients = *globalFunction_->dofVector_;
  auto&& 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 = typename std::decay<decltype(std::declval<Vector>()[std::declval<MultiIndex>()])>::type;
    using RangeBlock = std::decay_t<decltype(std::declval<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.
      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
60

61
template <class GlobalBasis, class TreePath>
Praetorius, Simon's avatar
Praetorius, Simon committed
62
typename DOFVectorConstView<GlobalBasis, TreePath>::DerivativeRange DOFVectorConstView<GlobalBasis, TreePath>::
63
GradientLocalFunction::operator()(LocalDomain const& x) const
64
{
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
  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_;
  forEachLeafNode(*subTree_, [&,this](auto const& node, auto)
  {
    // TODO: may DOFVectorView::Range to FieldVector type if necessary
    using LocalDerivativeTraits = Dune::Functions::DefaultDerivativeTraits<Dune::FieldVector<double,1>(Domain)>;
    using GradientBlock = typename LocalDerivativeTraits::Range;
    using MultiIndex = typename LocalIndexSet::MultiIndex;
    using CoefficientBlock = typename std::decay<decltype(std::declval<Vector>()[std::declval<MultiIndex>()])>::type;
    using RangeBlock = std::decay_t<decltype(std::declval<NodeToRangeEntry>()(node, dy))>;

    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);

    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 transformed reference gradient
      auto&& grad = gradients[i];
106

107
108
109
110
      // 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);
111
      assert(dimC*dimV == std::size_t(Flat<RangeBlock>::size(re)));
112
113
114
115
116
117
      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;
        }
118
      }
119
120
    }
  });
121

122
123
  return dy;
}
124
125

} // end namespace AMDiS