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

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

6
7
namespace AMDiS {

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

  auto&& coefficients = *globalFunction_->dofVector_;
  auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
17
  forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
18
19
20
21
  {
    auto&& fe = node.finiteElement();
    auto&& localBasis = fe.localBasis();

22
23
    LocalBasisCache<std::decay_t<decltype(localBasis)>> localBasisCache(localBasis);
    auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element(),x);
24
25

    // Get range entry associated to this node
26
    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));
27
28

    for (std::size_t i = 0; i < localBasis.size(); ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
29
      auto&& multiIndex = localView_.index(node.localIndex(i));
30
31

      // Get coefficient associated to i-th shape function
Praetorius, Simon's avatar
Praetorius, Simon committed
32
      auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
33
34

      // Get value of i-th shape function
Praetorius, Simon's avatar
Praetorius, Simon committed
35
      auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);
36

Praetorius, Simon's avatar
Praetorius, Simon committed
37
38
39
      std::size_t dimC = c.size();
      std::size_t dimV = v.size();
      assert(dimC*dimV == std::size_t(re.size()));
40
      for(std::size_t j = 0; j < dimC; ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
41
42
43
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*v[k];
44
45
46
47
48
49
50
      }
    }
  });

  return y;
}

Praetorius, Simon's avatar
Praetorius, Simon committed
51

52
53
template <class GB, class RT, class TP>
typename DOFVectorConstView<GB,RT,TP>::DerivativeRange DOFVectorConstView<GB,RT,TP>::
54
GradientLocalFunction::operator()(LocalDomain const& x) const
55
{
56
57
58
59
60
61
62
  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_;
63
  forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
64
65
66
67
68
69
70
71
  {
    // 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();

72
73
74
    LocalBasisCache<std::decay_t<decltype(localBasis)>> localBasisCache(localBasis);
    auto const& referenceGradients = localBasisCache.evaluateJacobian(localView_.element(),x);

75
76
77
78
79
80
81
82
83
    // The transposed inverse Jacobian of the map from the reference element to the element
    auto&& jacobian = geometry_.value().jacobianInverseTransposed(x);

    // 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
84
    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
85
86

    for (std::size_t i = 0; i < localBasis.size(); ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
87
      auto&& multiIndex = localView_.index(node.localIndex(i));
88
89

      // Get coefficient associated to i-th shape function
Praetorius, Simon's avatar
Praetorius, Simon committed
90
      auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);
91
92

      // Get value of i-th transformed reference gradient
Praetorius, Simon's avatar
Praetorius, Simon committed
93
      auto grad = Dune::Functions::flatVectorView(gradients[i]);
94

Praetorius, Simon's avatar
Praetorius, Simon committed
95
96
97
      std::size_t dimC = c.size();
      std::size_t dimV = grad.size();
      assert(dimC*dimV == std::size_t(re.size()));
98
      for(std::size_t j = 0; j < dimC; ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
99
100
101
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*grad[k];
102
      }
103
104
    }
  });
105

106
107
  return dy;
}
108
109

} // end namespace AMDiS