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

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

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

8
9
namespace AMDiS {

10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP>::LocalFunction
{
public:
  using Domain = typename EntitySet::LocalCoordinate;
  using Range = typename DiscreteFunction::Range;

  enum { hasDerivative = true };

private:
  using LocalView = typename GlobalBasis::LocalView;
  using Element = typename EntitySet::Element;
  using Geometry = typename Element::Geometry;

public:
  LocalFunction(DiscreteFunction const& globalFunction)
26
27
28
    : globalFunction_(globalFunction)
    , localView_(globalFunction_.basis().localView())
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
  {}

  /// \brief Bind the LocalView to the element
  void bind(Element const& element)
  {
    localView_.bind(element);
    bound_ = true;
  }

  /// \brief Unbind the LocalView from the element
  void unbind()
  {
    localView_.unbind();
    bound_ = false;
  }

  /// \brief Evaluate LocalFunction at bound element in local coordinates
  Range operator()(Domain const& x) const;

  /// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction
  friend GradientLocalFunction derivative(LocalFunction const& localFunction)
  {
51
    return GradientLocalFunction{localFunction.globalFunction_};
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
  }

  /// \brief The \ref polynomialDegree() of the LocalFunctions
  friend int order(LocalFunction const& self)
  {
    assert( self.bound_ );
    return polynomialDegree(*self.subTree_);
  }

  /// \brief Return the bound element
  Element const& localContext() const
  {
    assert( bound_ );
    return localView_.element();
  }

private:
69
  DiscreteFunction globalFunction_;
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
  LocalView localView_;
  SubTree const* subTree_;
  bool bound_ = false;
};


template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP>::GradientLocalFunction
{
  using R = typename DiscreteFunction::Range;
  using D = typename DiscreteFunction::Domain;
  using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
  using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>;

public:
  using Domain = typename EntitySet::LocalCoordinate;
  using Range = typename DerivativeTraits::Range;

  enum { hasDerivative = false };

private:
  using LocalView = typename GlobalBasis::LocalView;
  using Element = typename EntitySet::Element;
  using Geometry = typename Element::Geometry;

public:
  GradientLocalFunction(DiscreteFunction const& globalFunction)
97
98
99
    : globalFunction_(globalFunction)
    , localView_(globalFunction_.basis().localView())
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
  {}

  void bind(Element const& element)
  {
    localView_.bind(element);
    geometry_.emplace(element.geometry());
    bound_ = true;
  }

  void unbind()
  {
    localView_.unbind();
    geometry_.reset();
    bound_ = false;
  }

  /// Evaluate Gradient at bound element in local coordinates
  Range operator()(Domain const& x) const;

  friend int order(GradientLocalFunction const& self)
  {
    assert( self.bound_ );
    return std::max(0, polynomialDegree(*self.subTree_)-1);
  }

  /// Return the bound element
  Element const& localContext() const
  {
    assert( bound_ );
    return localView_.element();
  }

private:
133
  DiscreteFunction globalFunction_;
134
135
136
137
138
139
140
141
142
143
  LocalView localView_;
  SubTree const* subTree_;
  Dune::Std::optional<Geometry> geometry_;
  bool bound_ = false;
};


template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
LocalFunction::operator()(Domain const& x) const
144
145
146
147
{
  assert( bound_ );
  auto y = Range(0);

148
149
  auto&& coefficients = *globalFunction_.dofVector_;
  auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
150
  forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
151
152
153
154
  {
    auto&& fe = node.finiteElement();
    auto&& localBasis = fe.localBasis();

155
156
157
    using Node = std::decay_t<decltype(node)>;
    NodeCache<Node> localBasisCache(localBasis);
    auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(),x);
158
159

    // Get range entry associated to this node
160
    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));
161
162

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
171
172
173
      std::size_t dimC = c.size();
      std::size_t dimV = v.size();
      assert(dimC*dimV == std::size_t(re.size()));
174
      for(std::size_t j = 0; j < dimC; ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
175
176
177
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*v[k];
178
179
180
181
182
183
184
      }
    }
  });

  return y;
}

Praetorius, Simon's avatar
Praetorius, Simon committed
185

186
187
188
template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::GradientLocalFunction::Range DiscreteFunction<GB,VT,TP>::
GradientLocalFunction::operator()(Domain const& x) const
189
{
190
191
192
193
194
  assert( bound_ );
  Range dy;
  for (std::size_t j = 0; j < dy.size(); ++j)
    dy[j] = 0;

195
196
  auto&& coefficients = *globalFunction_.dofVector_;
  auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
197
  forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
198
199
  {
    // TODO: may DOFVectorView::Range to FieldVector type if necessary
200
201
    using LocalDerivativeTraits
      = Dune::Functions::DefaultDerivativeTraits<Dune::FieldVector<double,1>(Domain)>;
202
203
204
205
206
    using GradientBlock = typename LocalDerivativeTraits::Range;

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

207
208
209
    using Node = std::decay_t<decltype(node)>;
    NodeCache<Node> localBasisCache(localBasis);
    auto const& referenceGradients = localBasisCache.evaluateJacobian(localView_.element().type(),x);
210

211
212
213
214
215
216
217
218
219
    // 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
220
    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));
221
222

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
231
232
233
      std::size_t dimC = c.size();
      std::size_t dimV = grad.size();
      assert(dimC*dimV == std::size_t(re.size()));
234
      for(std::size_t j = 0; j < dimC; ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
235
236
237
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*grad[k];
238
      }
239
240
    }
  });
241

242
243
  return dy;
}
244

245

246
247
template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
248
249
operator()(Domain const& x) const
{
250
251
  using Grid = typename GlobalBasis::GridView::Grid;
  using IS = typename GlobalBasis::GridView::IndexSet;
252
253
254
255
256
257
258
259
260
261
262

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

263
} // end namespace AMDiS