DiscreteLocalFunction.inc.hpp 9.46 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>
Praetorius, Simon's avatar
Praetorius, Simon committed
7
#include <amdis/functions/NodeIndices.hpp>
8
#include <amdis/functions/Order.hpp>
9
#include <amdis/typetree/FiniteElementType.hpp>
10
11
12
13
14
#include <amdis/utility/LocalToGlobalAdapter.hpp>

#include <dune/common/ftraits.hh>

namespace AMDiS {
Praetorius, Simon's avatar
Praetorius, Simon committed
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
namespace Impl {

// specialization of Coeff has gather method
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);
}

// fallback implementation
template <class Coeff, class LocalView, class LocalCoeff>
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];
}

} // end namespace Impl

37

Praetorius, Simon's avatar
Praetorius, Simon committed
38
template <class Coeff, class GB, class TP>
39
  template <class Type>
Praetorius, Simon's avatar
Praetorius, Simon committed
40
class DiscreteFunction<Coeff const,GB,TP>::LocalFunction
41
{
42
43
44
45
46
47
  using DomainType = typename DiscreteFunction::Domain;
  using RangeType = typename DiscreteFunction::Range;

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

48
49
public:
  using Domain = typename EntitySet::LocalCoordinate;
50
  using Range = DerivativeRange<Type>;
51

52
  enum { hasDerivative = std::is_same<Type, tag::value>::value };
53
54
55
56
57

private:
  using LocalView = typename GlobalBasis::LocalView;
  using Element = typename EntitySet::Element;
  using Geometry = typename Element::Geometry;
58
  using NodeToRangeMap = Dune::Functions::DefaultNodeToRangeMap<SubTree>;
59
60
61

public:
  /// Constructor. Stores a copy of the DiscreteFunction.
62
63
64
65
66
67
68
  template <class LV>
  LocalFunction(std::shared_ptr<LV> localView, TP const& treePath, Coeff const& coefficients, Type type)
    : localView_(std::move(localView))
    , treePath_(treePath)
    , coefficients_(coefficients)
    , type_(type)
    , nodeToRangeMap_(subTree())
69
70
71
72
73
  {}

  /// \brief Bind the LocalView to the element
  void bind(Element const& element)
  {
74
75
76
    localView_->bind(element);
    geometry_.emplace(element.geometry());
    Impl::gather(coefficients_, *localView_, localCoefficients_, Dune::PriorityTag<4>{});
77
78
79
80
81
82
    bound_ = true;
  }

  /// \brief Unbind the LocalView from the element
  void unbind()
  {
83
84
    localView_->unbind();
    geometry_.reset();
85
86
87
88
    bound_ = false;
  }

  /// \brief Evaluate LocalFunction at bound element in local coordinates
89
  Range operator() (const Domain& local) const
90
  {
91
92
93
94
95
96
97
98
99
100
101
102
103
104
    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);
    }
105

106
    return Range(0);
107
108
  }

109
110
111
  /// \brief Create a LocalFunction representing the derivative.
  template <class DerType>
  LocalFunction<DerType> makeDerivative(DerType type) const
112
  {
113
114
    static_assert(std::is_same_v<Type, tag::value>);
    return {localView_, treePath_, coefficients_, type};
115
116
  }

117
  /// \brief The polynomial degree of the LocalFunctions basis-functions
118
119
  auto order() const
    -> decltype(AMDiS::order(std::declval<SubTree>()))
120
  {
121
122
    assert(bound_);
    return AMDiS::order(subTreeCache());
123
124
125
126
127
  }

  /// \brief Return the bound element
  Element const& localContext() const
  {
128
129
    assert(bound_);
    return localView_->element();
130
131
132
133
  }

private:

134
135
  template <class LocalCoordinate>
  auto evaluateFunction (const LocalCoordinate& local) const;
136

137
138
  template <class LocalCoordinate>
  auto evaluateJacobian (const LocalCoordinate& local) const;
139

140
141
  template <class LocalCoordinate>
  auto evaluateDivergence (const LocalCoordinate& local) const;
142

143
144
  template <class LocalCoordinate>
  auto evaluatePartialDerivative (const LocalCoordinate& local) const;
145

146
  Geometry const& geometry() const
147
  {
148
149
    assert(bound_);
    return *geometry_;
150
151
  }

152
  SubTree const& subTree() const
153
  {
154
    return Dune::TypeTree::child(localView_->tree(), treePath_);
155
156
  }

157
  SubTreeCache const& subTreeCache() const
158
  {
159
    return Dune::TypeTree::child(localView_->treeCache(), treePath_);
160
161
  }

162
163
164
165
private:
  std::shared_ptr<LocalView> localView_;
  TP treePath_;
  Coeff const& coefficients_;
166
  Type type_;
167
168
  NodeToRangeMap nodeToRangeMap_;

169
  std::optional<Geometry> geometry_;
Praetorius, Simon's avatar
Praetorius, Simon committed
170
  std::vector<ValueType> localCoefficients_;
171
172
173
174
  bool bound_ = false;
};


Praetorius, Simon's avatar
Praetorius, Simon committed
175
template <class Coeff, class GB, class TP>
176
177
178
179
  template <class Type>
    template <class LocalCoordinate>
auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
::evaluateFunction(LocalCoordinate const& local) const
180
{
181
182
  Range y(0);
  for_each_leaf_node(this->subTreeCache(), [&](auto const& node, auto const& tp)
183
  {
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
    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];
204
      }
205
206
    }
  });
207

208
209
  return y;
}
210
211


Praetorius, Simon's avatar
Praetorius, Simon committed
212
template <class Coeff, class GB, class TP>
213
214
215
216
  template <class Type>
    template <class LocalCoordinate>
auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
::evaluateJacobian(LocalCoordinate const& local) const
217
{
218
219
220
221
222
  Range dy(0);
  for_each_leaf_node(this->subTreeCache(), [&](auto const& node, auto const& tp)
  {
    LocalToGlobalBasisAdapter localBasis(node, this->geometry());
    auto const& gradients = localBasis.gradientsAt(local);
223

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

227
228
229
    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)]);
230

231
232
      // Get value of i-th transformed reference gradient
      auto grad = Dune::Functions::flatVectorView(gradients[i]);
233

234
235
236
237
238
239
240
241
242
243
      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];
      }
    }
  });
244

245
246
  return dy;
}
247
248


249
250
251
252
253
254
255
template <class Coeff, class GB, class TP>
  template <class Type>
    template <class LocalCoordinate>
auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
::evaluateDivergence(LocalCoordinate const& local) const
{
  static_assert(static_size_v<Range> == 1, "Divergence implemented for scalar coefficients only.");
256

257
  Range dy(0);
258

259
  auto&& node = this->subTreeCache();
260

261
262
  LocalToGlobalBasisAdapter localBasis(node.child(0), this->geometry());
  auto const& gradients = localBasis.gradientsAt(local);
263

264
265
266
267
268
269
270
271
  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];
272
273
  }

274
275
  return dy;
}
276
277


Praetorius, Simon's avatar
Praetorius, Simon committed
278
template <class Coeff, class GB, class TP>
279
280
281
282
  template <class Type>
    template <class LocalCoordinate>
auto DiscreteFunction<Coeff const,GB,TP>::LocalFunction<Type>
::evaluatePartialDerivative(LocalCoordinate const& local) const
283
{
284
  std::size_t comp = this->type_.comp;
285

286
287
288
289
290
  Range dy(0);
  for_each_leaf_node(this->subTreeCache(), [&](auto const& node, auto const& tp)
  {
    LocalToGlobalBasisAdapter localBasis(node, this->geometry());
    auto const& partial = localBasis.partialsAt(local, comp);
291

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

295
296
297
298
299
300
301
302
303
304
305
306
307
308
    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];
309
      }
310
311
    }
  });
312

313
314
  return dy;
}
315
316
317


} // end namespace AMDiS