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

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

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

9
10
namespace AMDiS {

11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
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)
27
28
29
    : globalFunction_(globalFunction)
    , localView_(globalFunction_.basis().localView())
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
30
31
  {}

32
33
34
35
36
37
  LocalFunction(LocalFunction const& other)
    : globalFunction_(other.globalFunction_)
    , localView_(globalFunction_.basis().localView())
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
  {}

38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
  /// \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)
  {
58
    return GradientLocalFunction{localFunction.globalFunction_};
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
  }

  /// \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:
76
  DiscreteFunction globalFunction_;
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
  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)
104
105
106
    : globalFunction_(globalFunction)
    , localView_(globalFunction_.basis().localView())
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
107
108
  {}

109
110
111
112
113
114
  GradientLocalFunction(GradientLocalFunction const& other)
    : globalFunction_(other.globalFunction_)
    , localView_(globalFunction_.basis().localView())
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
  {}

115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
  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:
146
  DiscreteFunction globalFunction_;
147
148
149
150
151
152
153
154
155
156
  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
157
158
{
  assert( bound_ );
159
  Range y(0);
160

161
162
  auto&& coefficients = *globalFunction_.dofVector_;
  auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
163
  for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
164
165
166
167
  {
    auto&& fe = node.finiteElement();
    auto&& localBasis = fe.localBasis();

168
    NodeCache<TYPEOF(node)> localBasisCache(localBasis);
169
    auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(),x);
170
171

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

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
183
184
185
      std::size_t dimC = c.size();
      std::size_t dimV = v.size();
      assert(dimC*dimV == std::size_t(re.size()));
186
      for(std::size_t j = 0; j < dimC; ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
187
188
189
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*v[k];
190
191
192
193
194
195
196
      }
    }
  });

  return y;
}

Praetorius, Simon's avatar
Praetorius, Simon committed
197

198
199
200
template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::GradientLocalFunction::Range DiscreteFunction<GB,VT,TP>::
GradientLocalFunction::operator()(Domain const& x) const
201
{
202
  assert( bound_ );
203
  Range dy(0);
204

205
206
  auto&& coefficients = *globalFunction_.dofVector_;
  auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
207
  for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
208
209
  {
    // TODO: may DOFVectorView::Range to FieldVector type if necessary
210
    using LocalDerivativeTraits
211
      = Dune::Functions::DefaultDerivativeTraits<VT(Domain)>;
212
213
214
215
    using GradientBlock = typename LocalDerivativeTraits::Range;

    auto&& fe = node.finiteElement();
    auto&& localBasis = fe.localBasis();
216
    std::size_t size = localBasis.size();
217

218
    NodeCache<TYPEOF(node)> localBasisCache(localBasis);
219
    auto const& referenceGradients = localBasisCache.evaluateJacobian(localView_.element().type(),x);
220

221
222
223
224
225
    // 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());
226
227
    for (std::size_t i = 0; i < gradients.size(); ++i) // J^(-T) * D[phi] -> grad^T
      Dune::MatVec::as_vector(gradients[i]) = jacobian * Dune::MatVec::as_vector(referenceGradients[i]);
228
229

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

232
    for (std::size_t i = 0; i < size; ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
233
      auto&& multiIndex = localView_.index(node.localIndex(i));
234
235

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
241
242
243
      std::size_t dimC = c.size();
      std::size_t dimV = grad.size();
      assert(dimC*dimV == std::size_t(re.size()));
244
      for(std::size_t j = 0; j < dimC; ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
245
246
247
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*grad[k];
248
      }
249
250
    }
  });
251

252
253
  return dy;
}
254

255

256
257
template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
258
259
operator()(Domain const& x) const
{
260
261
  using Grid = typename GlobalBasis::GridView::Grid;
  using IS = typename GlobalBasis::GridView::IndexSet;
262
263
264
265
266
267
268
269
270
271
272

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

273
} // end namespace AMDiS