DiscreteFunction.inc.hpp 4.36 KB
Newer Older
1
2
#pragma once

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

6
7
#include <dune/grid/utility/hierarchicsearch.hh>

8
9
namespace AMDiS {

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

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

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

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

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

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

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

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

  return y;
}

Praetorius, Simon's avatar
Praetorius, Simon committed
53

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

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

77
78
79
80
81
82
83
84
85
    // 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
86
    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
87
88

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

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

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

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

108
109
  return dy;
}
110

111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

template <class GB, class RT, class TP>
typename DiscreteFunction<GB,RT,TP>::Range DiscreteFunction<GB,RT,TP>::
operator()(Domain const& x) const
{
  using Grid = typename GB::GridView::Grid;
  using IS = typename GB::GridView::IndexSet;

  auto const& gv = this->basis().gridView();
  Dune::HierarchicSearch<Grid,IS> hsearch{gv.grid(), gv.indexSet()};

  auto element = hsearch.findEntity(x);
  auto geometry = element.geometry();
  auto localFct = localFunction(*this);
  localFct.bind(element);
  return localFct(geometry.local(x));
}

129
} // end namespace AMDiS