DiscreteFunction.inc.hpp 8.22 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
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:
26
  /// Constructor. Stores a copy of the DiscreteFunction.
27
  LocalFunction(DiscreteFunction const& globalFunction)
28
    : globalFunction_(globalFunction)
29
    , localView_(globalFunction_.basis()->localView())
30
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
31
32
  {}

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

40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
  /// \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)
  {
60
    return GradientLocalFunction{localFunction.globalFunction_};
61
62
63
  }

  /// \brief The \ref polynomialDegree() of the LocalFunctions
64
  int order() const
65
  {
66
67
    assert( bound_ );
    return polynomialDegree(*subTree_);
68
69
70
71
72
73
74
75
76
77
  }

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

private:
78
  DiscreteFunction globalFunction_;
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
104
  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:
105
  /// Constructor. Stores a copy of the DiscreteFunction.
106
  GradientLocalFunction(DiscreteFunction const& globalFunction)
107
    : globalFunction_(globalFunction)
108
    , localView_(globalFunction_.basis()->localView())
109
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
110
111
  {}

112
  /// Copy constructor.
113
114
  GradientLocalFunction(GradientLocalFunction const& other)
    : globalFunction_(other.globalFunction_)
115
    , localView_(globalFunction_.basis()->localView())
116
117
118
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
  {}

119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
  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;

136
  int order() const
137
  {
138
139
    assert( bound_ );
    return std::max(0, polynomialDegree(*subTree_)-1);
140
141
142
143
144
145
146
147
148
149
  }

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

private:
150
  DiscreteFunction globalFunction_;
151
152
153
154
155
156
157
158
159
160
  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
161
162
{
  assert( bound_ );
163
  Range y(0);
164

165
166
  auto&& coefficients = *globalFunction_.dofVector_;
  auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
167
  for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
168
169
170
171
  {
    auto&& fe = node.finiteElement();
    auto&& localBasis = fe.localBasis();

172
    NodeCache<TYPEOF(node)> localBasisCache(localBasis);
173
    auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(),x);
174
175

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

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

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

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

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

  return y;
}

Praetorius, Simon's avatar
Praetorius, Simon committed
201

202
203
204
template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::GradientLocalFunction::Range DiscreteFunction<GB,VT,TP>::
GradientLocalFunction::operator()(Domain const& x) const
205
{
206
  assert( bound_ );
207
  Range dy(0);
208

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

    auto&& fe = node.finiteElement();
    auto&& localBasis = fe.localBasis();
220
    std::size_t size = localBasis.size();
221

222
    NodeCache<TYPEOF(node)> localBasisCache(localBasis);
223
    auto const& referenceGradients = localBasisCache.evaluateJacobian(localView_.element().type(),x);
224

225
226
227
228
229
    // 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());
230
231
    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]);
232
233

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

236
    for (std::size_t i = 0; i < size; ++i) {
Praetorius, Simon's avatar
Praetorius, Simon committed
237
      auto&& multiIndex = localView_.index(node.localIndex(i));
238
239

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

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

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

256
257
  return dy;
}
258

259

260
261
template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
262
263
operator()(Domain const& x) const
{
264
265
  using Grid = typename GlobalBasis::GridView::Grid;
  using IS = typename GlobalBasis::GridView::IndexSet;
266

267
  auto const& gv = this->basis()->gridView();
268
269
270
271
272
273
274
275
276
  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));
}

277
} // end namespace AMDiS