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

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

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

8
9
namespace AMDiS {

10
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:
  LocalFunction(DiscreteFunction const& globalFunction)
26
27
28
    : globalFunction_(globalFunction)
    , localView_(globalFunction_.basis().localView())
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
  {}

  /// \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)
  {
51
    return GradientLocalFunction{localFunction.globalFunction_};
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
  }

  /// \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:
69
  DiscreteFunction globalFunction_;
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
  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)
97
98
99
    : globalFunction_(globalFunction)
    , localView_(globalFunction_.basis().localView())
    , subTree_(&child(localView_.tree(), globalFunction_.treePath()))
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
  {}

  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:
133
  DiscreteFunction globalFunction_;
134
135
136
137
138
139
140
141
142
143
  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
144
145
146
147
{
  assert( bound_ );
  auto y = Range(0);

148
149
  auto&& coefficients = *globalFunction_.dofVector_;
  auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
150
  for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
151
152
153
154
  {
    auto&& fe = node.finiteElement();
    auto&& localBasis = fe.localBasis();

155
    NodeCache<TYPEOF(node)> localBasisCache(localBasis);
156
    auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element().type(),x);
157
158

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

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
170
171
172
      std::size_t dimC = c.size();
      std::size_t dimV = v.size();
      assert(dimC*dimV == std::size_t(re.size()));
173
      for(std::size_t j = 0; j < dimC; ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
174
175
176
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*v[k];
177
178
179
180
181
182
183
      }
    }
  });

  return y;
}

Praetorius, Simon's avatar
Praetorius, Simon committed
184

185
186
187
template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::GradientLocalFunction::Range DiscreteFunction<GB,VT,TP>::
GradientLocalFunction::operator()(Domain const& x) const
188
{
189
190
191
192
193
  assert( bound_ );
  Range dy;
  for (std::size_t j = 0; j < dy.size(); ++j)
    dy[j] = 0;

194
195
  auto&& coefficients = *globalFunction_.dofVector_;
  auto&& nodeToRangeEntry = globalFunction_.nodeToRangeEntry_;
196
  for_each_leaf_node(*subTree_, [&,this](auto const& node, auto const& tp)
197
198
  {
    // TODO: may DOFVectorView::Range to FieldVector type if necessary
199
200
    using LocalDerivativeTraits
      = Dune::Functions::DefaultDerivativeTraits<Dune::FieldVector<double,1>(Domain)>;
201
202
203
204
205
    using GradientBlock = typename LocalDerivativeTraits::Range;

    auto&& fe = node.finiteElement();
    auto&& localBasis = fe.localBasis();

206
    NodeCache<TYPEOF(node)> localBasisCache(localBasis);
207
    auto const& referenceGradients = localBasisCache.evaluateJacobian(localView_.element().type(),x);
208

209
210
211
212
213
214
215
216
217
    // 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());
    for (std::size_t i = 0; i < gradients.size(); ++i)
      multiplies_ABt(referenceGradients[i], jacobian, gradients[i]); // D[phi] * J^(-1) -> grad

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

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
229
230
231
      std::size_t dimC = c.size();
      std::size_t dimV = grad.size();
      assert(dimC*dimV == std::size_t(re.size()));
232
      for(std::size_t j = 0; j < dimC; ++j) {
Praetorius, Simon's avatar
Praetorius, Simon committed
233
234
235
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*grad[k];
236
      }
237
238
    }
  });
239

240
241
  return dy;
}
242

243

244
245
template <class GB, class VT, class TP>
typename DiscreteFunction<GB,VT,TP>::Range DiscreteFunction<GB,VT,TP>::
246
247
operator()(Domain const& x) const
{
248
249
  using Grid = typename GlobalBasis::GridView::Grid;
  using IS = typename GlobalBasis::GridView::IndexSet;
250
251
252
253
254
255
256
257
258
259
260

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

261
} // end namespace AMDiS