DiscreteLocalFunction.inc.hpp 11 KB
Newer Older
1
2
#pragma once

Praetorius, Simon's avatar
Praetorius, Simon committed
3
#include <amdis/Output.hpp>
4
5
#include <amdis/common/DerivativeTraits.hpp>
#include <amdis/common/FieldMatVec.hpp>
6
#include <amdis/common/Math.hpp>
7
8
#include <amdis/common/RecursiveForEach.hpp>
#include <amdis/functions/HierarchicNodeToRangeMap.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
9
#include <amdis/functions/NodeIndices.hpp>
10
#include <amdis/functions/Order.hpp>
11
#include <amdis/typetree/FiniteElementType.hpp>
12
13
14
15
16
#include <amdis/utility/LocalToGlobalAdapter.hpp>

#include <dune/common/ftraits.hh>

namespace AMDiS {
Praetorius, Simon's avatar
Praetorius, Simon committed
17
18
namespace Impl {

19
// specialization for containers with gather method
Praetorius, Simon's avatar
Praetorius, Simon committed
20
21
22
23
24
25
26
template <class Coeff, class LocalView, class LocalCoeff,
  class = decltype(std::declval<Coeff>().gather(std::declval<LocalView>(), std::declval<LocalCoeff&>()))>
void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoeff, Dune::PriorityTag<2>)
{
  coeff.gather(localView, localCoeff);
}

27
28
29
// implementation for containers with multi-index access
template <class Coeff, class LocalView, class LocalCoeff,
  class = decltype(std::declval<Coeff>()[std::declval<typename LocalView::MultiIndex>()])>
Praetorius, Simon's avatar
Praetorius, Simon committed
30
31
32
33
34
35
36
37
void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoeff, Dune::PriorityTag<1>)
{
  localCoeff.resize(localView.size());
  auto it = localCoeff.begin();
  for (auto const& idx : nodeIndices(localView))
    *it++ = coeff[idx];
}

38
39
40
41
42
43
44
45
46
47
48
49
// fallback implementation
template <class Coeff, class LocalView, class LocalCoeff,
  class = decltype(std::declval<Coeff>()[std::integral_constant<std::size_t,0>{}])>
void gather(Coeff const& coeff, LocalView const& localView, LocalCoeff& localCoeff, Dune::PriorityTag<0>)
{
  auto backend = Dune::Functions::istlVectorBackend(coeff);
  localCoeff.resize(localView.size());
  auto it = localCoeff.begin();
  for (auto const& idx : nodeIndices(localView))
    *it++ = backend[idx];
}

Praetorius, Simon's avatar
Praetorius, Simon committed
50
51
} // end namespace Impl

52

53
template <class C, class GB, class TP, class R>
54
  template <class Type>
55
class DiscreteFunction<C const,GB,TP,R>::LocalFunction
56
{
57
58
59
60
61
62
  using DomainType = typename DiscreteFunction::Domain;
  using RangeType = typename DiscreteFunction::Range;

  template <class T>
  using DerivativeRange = typename DerivativeTraits<RangeType(DomainType), T>::Range;

63
64
public:
  using Domain = typename EntitySet::LocalCoordinate;
65
  using Range = DerivativeRange<Type>;
66

67
  enum { hasDerivative = std::is_same<Type, tag::value>::value };
68
69
70
71
72

private:
  using LocalView = typename GlobalBasis::LocalView;
  using Element = typename EntitySet::Element;
  using Geometry = typename Element::Geometry;
73
  using NodeToRangeMap = HierarchicNodeToRangeMap;
74
75
76

public:
  /// Constructor. Stores a copy of the DiscreteFunction.
77
  template <class LV>
78
  LocalFunction(std::shared_ptr<LV> localView, TP const& treePath, C const& coefficients, Type type)
79
80
81
82
    : localView_(std::move(localView))
    , treePath_(treePath)
    , coefficients_(coefficients)
    , type_(type)
83
84
85
86
87
  {}

  /// \brief Bind the LocalView to the element
  void bind(Element const& element)
  {
88
89
90
    localView_->bind(element);
    geometry_.emplace(element.geometry());
    Impl::gather(coefficients_, *localView_, localCoefficients_, Dune::PriorityTag<4>{});
91
92
93
94
95
96
    bound_ = true;
  }

  /// \brief Unbind the LocalView from the element
  void unbind()
  {
97
98
    localView_->unbind();
    geometry_.reset();
99
100
101
102
    bound_ = false;
  }

  /// \brief Evaluate LocalFunction at bound element in local coordinates
103
  Range operator() (const Domain& local) const
104
  {
105
106
107
108
109
110
111
112
113
114
115
116
117
118
    assert(bound_);
    if constexpr (std::is_same_v<Type, tag::value>)
      return evaluateFunction(local);
    else if constexpr (std::is_same_v<Type, tag::gradient>)
      return evaluateJacobian(local);
    else if constexpr (std::is_same_v<Type, tag::partial>)
      return evaluatePartialDerivative(local);
    else if constexpr (std::is_same_v<Type, tag::divergence>) {
      static constexpr bool isVectorNode
        = SubTree::isPower && SubTree::degree() == GridView::dimensionworld;
      static_assert(isVectorNode);
      if (isVectorNode)
        return evaluateDivergence(local);
    }
119

120
    return Range{};
121
122
  }

123
124
125
  /// \brief Create a LocalFunction representing the derivative.
  template <class DerType>
  LocalFunction<DerType> makeDerivative(DerType type) const
126
  {
127
128
    static_assert(std::is_same_v<Type, tag::value>);
    return {localView_, treePath_, coefficients_, type};
129
130
  }

131
  /// \brief The polynomial degree of the LocalFunctions basis-functions
132
133
  auto order() const
    -> decltype(AMDiS::order(std::declval<SubTree>()))
134
  {
135
    assert(bound_);
136
    return AMDiS::order(subTree());
137
138
139
140
141
  }

  /// \brief Return the bound element
  Element const& localContext() const
  {
142
143
    assert(bound_);
    return localView_->element();
144
145
146
147
  }

private:

148
149
  template <class LocalCoordinate>
  auto evaluateFunction (const LocalCoordinate& local) const;
150

151
152
  template <class LocalCoordinate>
  auto evaluateJacobian (const LocalCoordinate& local) const;
153

154
155
  template <class LocalCoordinate>
  auto evaluateDivergence (const LocalCoordinate& local) const;
156

157
158
  template <class LocalCoordinate>
  auto evaluatePartialDerivative (const LocalCoordinate& local) const;
159

160
  // get the geometry of the bound element, that is constructed in the bind() method
161
  Geometry const& geometry() const
162
  {
163
164
    assert(bound_);
    return *geometry_;
165
166
  }

167
  // return the sub-tree of the basis-tree associated to the tree-path
168
  SubTree const& subTree() const
169
  {
170
    return Dune::TypeTree::child(localView_->tree(), treePath_);
171
172
  }

173
174
175
176
177
  // return the sub-tree of the basis-tree-cache associated to the tree-path
  // NOTE: this is only return in case the local-view does provide a treeCache() function
  template <class LV,
    class = decltype(std::declval<LV>().treeCache())>
  auto const& subTreeCache(Dune::PriorityTag<1>) const
178
  {
179
    return Dune::TypeTree::child(localView_->treeCache(), treePath_);
180
181
  }

182
183
184
185
186
187
188
189
190
191
192
193
  // construct a new tree-cache that stores a reference to the sub-tree
  template <class>
  auto subTreeCache(Dune::PriorityTag<0>) const
  {
    return makeNodeCache(subTree());
  }

  decltype(auto) subTreeCache() const
  {
    return subTreeCache<LocalView>(Dune::PriorityTag<2>());
  }

194
195
196
private:
  std::shared_ptr<LocalView> localView_;
  TP treePath_;
197
  Coefficients const& coefficients_;
198
  Type type_;
199
200
  NodeToRangeMap nodeToRangeMap_;

201
  std::optional<Geometry> geometry_;
202
  std::vector<Coefficient> localCoefficients_;
203
204
205
206
  bool bound_ = false;
};


207
template <class C, class GB, class TP, class R>
208
209
  template <class Type>
    template <class LocalCoordinate>
210
auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
211
::evaluateFunction(LocalCoordinate const& local) const
212
{
213
214
215
  Range y;
  Recursive::forEach(y, [](auto& v) { v = 0; });

216
  Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
217
  {
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
    auto const& shapeFunctionValues = node.localBasisValuesAt(local);
    std::size_t size = node.finiteElement().size();

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

    for (std::size_t i = 0; i < size; ++i) {
      // Get coefficient associated to i-th shape function
      auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);

      // Get value of i-th shape function
      auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);

      std::size_t dimC = c.size();
      std::size_t dimV = v.size();
      assert(dimC*dimV == std::size_t(re.size()));
      for(std::size_t j = 0; j < dimC; ++j) {
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*v[k];
238
      }
239
240
    }
  });
241

242
243
  return y;
}
244
245


246
template <class C, class GB, class TP, class R>
247
248
  template <class Type>
    template <class LocalCoordinate>
249
auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
250
::evaluateJacobian(LocalCoordinate const& local) const
251
{
252
253
254
  Range dy;
  Recursive::forEach(dy, [](auto& v) { v = 0; });

255
  Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
256
257
258
  {
    LocalToGlobalBasisAdapter localBasis(node, this->geometry());
    auto const& gradients = localBasis.gradientsAt(local);
259

260
261
    // Get range entry associated to this node
    auto re = Dune::Functions::flatVectorView(nodeToRangeMap_(node, tp, dy));
262

263
264
265
    for (std::size_t i = 0; i < localBasis.size(); ++i) {
      // Get coefficient associated to i-th shape function
      auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
266

267
268
      // Get value of i-th transformed reference gradient
      auto grad = Dune::Functions::flatVectorView(gradients[i]);
269

270
271
272
273
274
275
276
277
278
279
      std::size_t dimC = c.size();
      std::size_t dimV = grad.size();
      assert(dimC*dimV == std::size_t(re.size()));
      for(std::size_t j = 0; j < dimC; ++j) {
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*grad[k];
      }
    }
  });
280

281
282
  return dy;
}
283
284


285
template <class C, class GB, class TP, class R>
286
287
  template <class Type>
    template <class LocalCoordinate>
288
auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
289
290
291
::evaluateDivergence(LocalCoordinate const& local) const
{
  static_assert(static_size_v<Range> == 1, "Divergence implemented for scalar coefficients only.");
292

293
294
  Range dy;
  Recursive::forEach(dy, [](auto& v) { v = 0; });
295

296
  auto&& node = this->subTreeCache();
297

298
299
  LocalToGlobalBasisAdapter localBasis(node.child(0), this->geometry());
  auto const& gradients = localBasis.gradientsAt(local);
300

301
302
303
304
305
306
307
308
  auto re = Dune::Functions::flatVectorView(dy);
  assert(int(re.size()) == 1);
  for (std::size_t i = 0; i < localBasis.size(); ++i) {
    auto grad = Dune::Functions::flatVectorView(gradients[i]);

    assert(int(grad.size()) == GridView::dimensionworld);
    for (std::size_t j = 0; j < GridView::dimensionworld; ++j)
      re[0] += localCoefficients_[node.child(j).localIndex(i)] * grad[j];
309
310
  }

311
312
  return dy;
}
313
314


315
template <class C, class GB, class TP, class R>
316
317
  template <class Type>
    template <class LocalCoordinate>
318
auto DiscreteFunction<C const,GB,TP,R>::LocalFunction<Type>
319
::evaluatePartialDerivative(LocalCoordinate const& local) const
320
{
321
  std::size_t comp = this->type_.comp;
322

323
324
325
  Range dy;
  Recursive::forEach(dy, [](auto& v) { v = 0; });

326
  Traversal::forEachLeafNode(this->subTreeCache(), [&](auto const& node, auto const& tp)
327
328
329
  {
    LocalToGlobalBasisAdapter localBasis(node, this->geometry());
    auto const& partial = localBasis.partialsAt(local, comp);
330

331
332
    // Get range entry associated to this node
    auto re = Dune::Functions::flatVectorView(nodeToRangeMap_(node, tp, dy));
333

334
335
336
337
338
339
340
341
342
343
344
345
346
347
    for (std::size_t i = 0; i < localBasis.size(); ++i) {
      // Get coefficient associated to i-th shape function
      auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);

      // Get value of i-th transformed reference partial_derivative
      auto d_comp = Dune::Functions::flatVectorView(partial[i]);

      std::size_t dimC = c.size();
      std::size_t dimV = d_comp.size();
      assert(dimC*dimV == std::size_t(re.size()));
      for(std::size_t j = 0; j < dimC; ++j) {
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*d_comp[k];
348
      }
349
350
    }
  });
351

352
353
  return dy;
}
354
355
356


} // end namespace AMDiS