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

3
#include <amdis/common/DerivativeTraits.hpp>
4
#include <amdis/common/FieldMatVec.hpp>
5
#include <amdis/utility/LocalBasisCache.hpp>
6
#include <amdis/utility/LocalToGlobalAdapter.hpp>
7

8
#include <dune/common/ftraits.hh>
9
10
#include <dune/grid/utility/hierarchicsearch.hh>

11
12
namespace AMDiS {

13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
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:
28
  /// Constructor. Stores a copy of the DiscreteFunction.
29
  LocalFunction(DiscreteFunction const& globalFunction)
30
    : globalFunction_(globalFunction)
31
    , localView_(globalFunction_.basis()->localView())
32
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
33
34
  {}

35
  /// Copy constructor.
36
37
  LocalFunction(LocalFunction const& other)
    : globalFunction_(other.globalFunction_)
38
    , localView_(globalFunction_.basis()->localView())
39
40
41
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
  {}

42
43
44
45
  /// \brief Bind the LocalView to the element
  void bind(Element const& element)
  {
    localView_.bind(element);
46
47

    globalFunction_.coefficients().gather(localView_, localCoefficients_);
48
49
50
51
52
53
54
55
56
57
58
59
60
61
    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
62
  GradientLocalFunction makeDerivative(tag::gradient type) const
63
  {
64
65
66
67
68
69
70
71
72
73
74
    return GradientLocalFunction{globalFunction_, type};
  }

  DivergenceLocalFunction makeDerivative(tag::divergence type) const
  {
    return DivergenceLocalFunction{globalFunction_, type};
  }

  PartialLocalFunction makeDerivative(tag::partial type) const
  {
    return PartialLocalFunction{globalFunction_, type};
75
76
77
  }

  /// \brief The \ref polynomialDegree() of the LocalFunctions
78
  int order() const
79
  {
80
81
    assert( bound_ );
    return polynomialDegree(*subTree_);
82
83
84
85
86
87
88
89
90
91
  }

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

private:
92
  DiscreteFunction globalFunction_;
93
94
  LocalView localView_;
  SubTree const* subTree_;
95
96

  std::vector<VT> localCoefficients_;
97
98
99
100
101
  bool bound_ = false;
};


template <class GB, class VT, class TP>
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
LocalFunction::operator()(Domain const& x) const
{
  assert( bound_ );
  Range y(0);

  auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
  for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
  {
    auto localBasisCache = makeNodeCache(node);
    auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(), x);
    std::size_t size = node.finiteElement().size();

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

    for (std::size_t i = 0; i < size; ++i) {
      // Get coefficient associated to i-th shape function
120
      auto c = Dune::Functions::flatVectorView(localCoefficients_[node.localIndex(i)]);
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

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

  return y;
}


template <class GB, class VT, class TP>
  template <class Type>
class DiscreteFunction<GB,VT,TP>::DerivativeLocalFunctionBase
143
144
145
146
{
  using R = typename DiscreteFunction::Range;
  using D = typename DiscreteFunction::Domain;
  using RawSignature = typename Dune::Functions::SignatureTraits<R(D)>::RawSignature;
147
  using Traits = DerivativeTraits<RawSignature, Type>;
148
149
150

public:
  using Domain = typename EntitySet::LocalCoordinate;
151
  using Range = typename Traits::Range;
152
153
154
155
156
157
158
159
160

  enum { hasDerivative = false };

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

public:
161
  /// Constructor. Stores a copy of the DiscreteFunction.
162
  DerivativeLocalFunctionBase(DiscreteFunction const& globalFunction, Type const& type)
163
    : globalFunction_(globalFunction)
164
    , type_(type)
165
    , localView_(globalFunction_.basis()->localView())
166
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
167
168
  {}

169
  /// Copy constructor.
170
  DerivativeLocalFunctionBase(DerivativeLocalFunctionBase const& other)
171
    : globalFunction_(other.globalFunction_)
172
    , type_(other.type_)
173
    , localView_(globalFunction_.basis()->localView())
174
175
176
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
  {}

177
178
179
180
  void bind(Element const& element)
  {
    localView_.bind(element);
    geometry_.emplace(element.geometry());
181
182

    globalFunction_.coefficients().gather(localView_, localCoefficients_);
183
184
185
186
187
188
189
190
191
192
    bound_ = true;
  }

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

193
  int order() const
194
  {
195
196
    assert( bound_ );
    return std::max(0, polynomialDegree(*subTree_)-1);
197
198
199
200
201
202
203
204
205
  }

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

206
207
208
209
210
211
212
213
214
215
216
  Geometry const& geometry() const
  {
    assert( bound_ );
    return geometry_.value();
  }

  LocalView const& localView() const
  {
    return localView_;
  }

217
protected:
218
  DiscreteFunction globalFunction_;
219
  Type type_;
220
221
222
  LocalView localView_;
  SubTree const* subTree_;
  Dune::Std::optional<Geometry> geometry_;
223
  std::vector<VT> localCoefficients_;
224
225
226
227
228
  bool bound_ = false;
};


template <class GB, class VT, class TP>
229
230
class DiscreteFunction<GB,VT,TP>::GradientLocalFunction
    : public DerivativeLocalFunctionBase<tag::gradient>
231
{
232
233
234
235
  using Super = DerivativeLocalFunctionBase<tag::gradient>;
public:
  using Range = typename Super::Range;
  using Domain = typename Super::Domain;
236

237
238
  // adopt constructor from base class
  using DerivativeLocalFunctionBase<tag::gradient>::DerivativeLocalFunctionBase;
239
240
241

  /// Evaluate Gradient at bound element in local coordinates
  Range operator()(Domain const& x) const
242
  {
243
    assert( this->bound_ );
244
245
    Range dy(0);

246
    auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_;
247
    for_each_leaf_node(*this->subTree_, [&](auto const& node, auto const& tp)
248
    {
249
      auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry());
250
251
252
253
254
255
256
      auto const& gradients = localBasis.gradientsAt(x);

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

      for (std::size_t i = 0; i < localBasis.size(); ++i) {
        // Get coefficient associated to i-th shape function
257
        auto c = Dune::Functions::flatVectorView(this->localCoefficients_[node.localIndex(i)]);
258
259
260
261
262
263
264
265
266
267
268
269
270
271

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

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

273
274
    return dy;
  }
275
276

  using Super::localCoefficients_;
277
};
278
279


280
281
282
283
284
template <class GB, class VT, class TP>
class DiscreteFunction<GB,VT,TP>::DivergenceLocalFunction
    : public DerivativeLocalFunctionBase<tag::divergence>
{
  using Super = DerivativeLocalFunctionBase<tag::divergence>;
285

286
287
288
public:
  using Range = typename Super::Range;
  using Domain = typename Super::Domain;
289

290
291
  // adopt constructor from base class
  using DerivativeLocalFunctionBase<tag::divergence>::DerivativeLocalFunctionBase;
292

293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
  /// Evaluate divergence at bound element in local coordinates
  Range operator()(Domain const& x) const
  {
    return evaluate_node(x, bool_t<SubTree::isPower && SubTree::degree() == GridView::dimensionworld>{});
  }

private:
  Range evaluate_node(Domain const& x, std::false_type) const
  {
    error_exit("Cannot calculate divergence(node) where node is not a power node.");
    return Range(0);
  }

  Range evaluate_node(Domain const& x, std::true_type) const
  {
308
    static_assert(static_size_v<Range> == 1, "Implemented for scalar coefficients only.");
309

310
    assert( this->bound_ );
311
312
    Range dy(0);

313
    auto&& node = *this->subTree_;
314

315
    auto localBasis = makeLocalToGlobalBasisAdapter(node.child(0), this->geometry());
316
317
318
319
320
321
322
323
    auto const& gradients = localBasis.gradientsAt(x);

    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);
324
325
      for (std::size_t j = 0; j < GridView::dimensionworld; ++j)
        re[0] += localCoefficients_[node.child(j).localIndex(i)] * grad[j];
326
327
    }

328
329
    return dy;
  }
330
331

  using Super::localCoefficients_;
332
};
333

Praetorius, Simon's avatar
Praetorius, Simon committed
334

335
template <class GB, class VT, class TP>
336
337
class DiscreteFunction<GB,VT,TP>::PartialLocalFunction
    : public DerivativeLocalFunctionBase<tag::partial>
338
{
339
  using Super = DerivativeLocalFunctionBase<tag::partial>;
340

341
342
343
public:
  using Range = typename Super::Range;
  using Domain = typename Super::Domain;
344

345
  using DerivativeLocalFunctionBase<tag::partial>::DerivativeLocalFunctionBase;
346

347
348
349
  /// Evaluate partial derivative at bound element in local coordinates
  Range operator()(Domain const& x) const
  {
350
    assert( this->bound_ );
351
352
    Range dy(0);

353
    std::size_t comp = this->type_.comp;
354

355
    auto&& nodeToRangeEntry = this->globalFunction_.nodeToRangeEntry_;
356
    for_each_leaf_node(*this->subTree_, [&](auto const& node, auto const& tp)
357
    {
358
      auto localBasis = makeLocalToGlobalBasisAdapter(node, this->geometry());
359
      auto const& partial = localBasis.partialsAt(x, comp);
360
361
362
363
364
365

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

      for (std::size_t i = 0; i < localBasis.size(); ++i) {
        // Get coefficient associated to i-th shape function
366
        auto c = Dune::Functions::flatVectorView(this->localCoefficients_[node.localIndex(i)]);
367
368
369
370
371
372
373
374
375
376
377
378
379
380

        // 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];
        }
      }
    });
381

382
383
    return dy;
  }
384
385

  using Super::localCoefficients_;
386
};
387

388

389
390
template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
391
392
operator()(Domain const& x) const
{
393
394
  using Grid = typename GlobalBasis::GridView::Grid;
  using IS = typename GlobalBasis::GridView::IndexSet;
395

396
  auto const& gv = this->basis()->gridView();
397
398
399
400
401
402
403
404
405
  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));
}

406
} // end namespace AMDiS