analyticdiscretefunction.hh 8.09 KB
Newer Older
1
2
#ifndef DUNE_CURVEDSURFACEGRID_ANALYTIC_DISCRETEFUNCTION_HH
#define DUNE_CURVEDSURFACEGRID_ANALYTIC_DISCRETEFUNCTION_HH
3
4
5
6
7

#include <optional>
#include <type_traits>
#include <utility>

8
9
10
11
#if !HAVE_DUNE_LOCALFUNCTIONS
#error "Need dune-localfunctions for the definition of the AnalyticDiscreteFunction"
#endif

12
13
14
15
16
#include <dune/common/typeutilities.hh>
#include <dune/curvedsurfacegrid/gridfunctions/discretegridviewfunction.hh>
#include <dune/curvedsurfacegrid/gridfunctions/gridentityset.hh>
#include <dune/localfunctions/lagrange/lfecache.hh>

17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
namespace Dune {

/// \brief LocalFunction associated to the \ref AnalyticDiscreteFunction
/**
 * \tparam LocalContext  Context this localFunction can be bound to, e.g. a grid element
 * \tparam F             Type of a function that can be evaluated in global coordinates.
 * \tparam derivativeOrder  Order of the derivative of the functor to return on evaluation
 **/
template <class LocalContext, class F, int derivativeOrder = 0>
class LocalAnalyticDiscreteFunction;

/// \brief Generator for \ref LocalAnalyticDiscreteFunction
/**
 * \param  ff            Function that can be evaluated at global coordinates of the \ref LocalContext
 * \param  order         Polynomial order of the approximation of ff
 * \tparam LocalContext  Context this localFunction can be bound to, e.g. a grid element
 **/
template <class LocalContext, class FF>
auto localAnalyticDiscreteFunction (FF&& ff, int order)
{
  using F = std::decay_t<FF>;
  return LocalAnalyticDiscreteFunction<LocalContext, F>{std::forward<FF>(ff), order};
}


template <class LC, class Functor, int derivativeOrder>
class LocalAnalyticDiscreteFunction
44
{
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
public:
  using LocalContext = LC;
  using Geometry = typename LocalContext::Geometry;
  using ctype = typename Geometry::ctype;

  using Domain = typename Geometry::GlobalCoordinate;
  using RangeType = std::result_of_t<Functor(Domain)>;

  template <int degree>
  using DerivativeRange = typename Impl::DerivativeRangeType<RangeType(Domain), degree, Functions::DefaultDerivativeTraits>::type;

  using LocalDomain = typename Geometry::LocalCoordinate;
  using Range = DerivativeRange<derivativeOrder>;
  using Signature = Range(LocalDomain);

public:
  /// \brief Constructor. Stores the functor f by value
  template <class FF>
  LocalAnalyticDiscreteFunction (FF&& ff, int order)
    : f_(std::forward<FF>(ff))
    , order_(order)
    , cache_(order)
  {}

  /// \brief bind the LocalFunction to the local context
70
  /**
71
   * Stores the localContext and its geometry in a cache variable
72
   **/
73
74
75
76
  void bind (const LocalContext& localContext)
  {
    localContext_.emplace(localContext);
    geometry_.emplace(localContext.geometry());
77

78
79
80
81
82
83
84
85
    lfe_ = &cache_.get(localContext_->type());
    lfe_->localInterpolation().interpolate([&](const auto& local)
    {
      return f_(geometry_->global(local));
    }, coefficients_);
  }

  /// \brief unbind the localContext from the localFunction
86
  /**
87
   * Reset the geometry
88
   **/
89
  void unbind ()
90
  {
91
92
93
    lfe_ = nullptr;
    geometry_.reset();
    localContext_.reset();
94
95
  }

96
97
98
99
100
101
  /// \brief evaluate the stored function in local coordinates
  /**
   * Transform the local coordinates to global coordinates first,
   * then evalute the stored functor.
   **/
  Range operator() (const LocalDomain& local) const
102
  {
103
    static_assert(derivativeOrder < 2, "Higher-order derivatives not implemented");
104

105
106
107
108
109
    if constexpr (derivativeOrder == 0)
      return evaluateFunction(local);
    else if constexpr (derivativeOrder == 1)
      return evaluateJacobian(local);
  }
110

111
112
113
114
  friend LocalAnalyticDiscreteFunction<LC,Functor,derivativeOrder+1> derivative (const LocalAnalyticDiscreteFunction& lf)
  {
    return {lf.f_, lf.order_};
  }
115

116
117
118
119
120
121
  /// \brief return the bound localContext.
  const LocalContext& localContext () const
  {
    assert(!!localContext_);
    return *localContext_;
  }
122

123
124
125
126
127
  /// \brief obtain the functor
  const Functor& f () const
  {
    return f_;
  }
128

129
130
131
132
private:
  DerivativeRange<0> evaluateFunction (const LocalDomain& local) const
  {
    assert(!!lfe_);
133

134
    const auto& localBasis = lfe_->localBasis();
135

136
137
    localBasis.evaluateFunction(local, shapeValues_);
    assert(coefficients_.size() == shapeValues_.size());
138

139
140
141
    DerivativeRange<0> nh(0);
    for (std::size_t i = 0; i < shapeValues_.size(); ++i)
      nh.axpy(shapeValues_[i], coefficients_[i]);
142

143
144
    return nh;
  }
145

146
147
148
149
  DerivativeRange<1> evaluateJacobian (const LocalDomain& local) const
  {
    assert(!!geometry_);
    assert(!!lfe_);
150

151
    const auto& localBasis = lfe_->localBasis();
152

153
154
155
    // evaluate basis functions in local coordinate
    localBasis.evaluateJacobian(local, shapeGradients_);
    assert(coefficients_.size() == shapeGradients_.size());
156

157
158
159
160
161
    // transform gradients to global coordinates
    auto jit = geometry_->jacobianInverseTransposed(local);
    gradients_.resize(shapeGradients_.size());
    for (std::size_t i = 0; i < shapeGradients_.size(); ++i)
      jit.mv(shapeGradients_[i][0], gradients_[i]);
162

163
164
165
166
    DerivativeRange<1> J(0);
    for (std::size_t i = 0; i < coefficients_.size(); ++i)
      for (std::size_t j = 0; j < J.N(); ++j)
        J[j].axpy(coefficients_[i][j], gradients_[i]);
167

168
169
    return J;
  }
170

171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
private:
  Functor f_;
  int order_;

  using LFECache = LagrangeLFECache<ctype,ctype,Geometry::mydimension>;
  LFECache cache_;

  using LFE = typename LFECache::FiniteElementType;
  LFE const* lfe_ = nullptr;

  std::vector<FieldVector<ctype,Geometry::coorddimension>> coefficients_;

  using LB = typename LFE::Traits::LocalBasisType;
  mutable std::vector<typename LB::Traits::RangeType> shapeValues_;
  mutable std::vector<typename LB::Traits::JacobianType> shapeGradients_;
  mutable std::vector<FieldVector<ctype,Geometry::coorddimension>> gradients_;

  // some caches
  std::optional<LocalContext> localContext_;
  std::optional<Geometry> geometry_;
};


/// \brief GridFunction associated to the mapping F
/**
 * \tparam Grid       The grid type with elements the corresponding LocalFunction can be bound to
 * \tparam F          Type of a function that can be evaluated in global coordinates.
 * \tparam ORDER      Polynomial order of the local lagrange basis functions (optional)
 **/
template <class Grid, class F, int ORDER = -1>
class AnalyticDiscreteFunction;

/// \brief Generator for \ref AnalyticDiscreteFunction
/**
 * \param  ff    Function that can be evaluated at global coordinates of the \ref Grid
 * \tparam Grid  The grid type with elements the corresponding LocalFunction can be bound to
 **/
template <class FF, class Grid>
auto analyticDiscreteFunction (FF&& ff, const Grid& /*grid*/, int order)
{
  using F = std::decay_t<FF>;
  return AnalyticDiscreteFunction<Grid, F>{std::forward<FF>(ff), order};
}
214

215
216
217
218
219
220
template <int ORDER, class FF, class Grid>
auto analyticDiscreteFunction (FF&& ff, const Grid& /*grid*/)
{
  using F = std::decay_t<FF>;
  return AnalyticDiscreteFunction<Grid, F, ORDER>{std::forward<FF>(ff)};
}
221
222


223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
template <class GridType, class Functor, int ORDER>
class AnalyticDiscreteFunction
{
public:
  using Grid = GridType;
  using EntitySet = GridEntitySet<Grid,0>;

  using Domain = typename EntitySet::GlobalCoordinate;
  using Range = std::result_of_t<Functor(Domain)>;
  using Signature = Range(Domain);

public:
  /// \brief Constructor. Stores the functor f by value
  template <class FF,
    disableCopyMove<AnalyticDiscreteFunction, FF> = 0>
  AnalyticDiscreteFunction (FF&& ff, int order = ORDER)
    : f_(std::forward<FF>(ff))
    , order_(order)
  {}

  /// \brief evaluate the stored function in global coordinates
  Range operator() (const Domain& x) const
245
  {
246
    return f_(x);
247
248
  }

249
250
251
  /// \brief construct the \ref LocalAnalyticDiscreteFunction
  using LocalFunction = LocalAnalyticDiscreteFunction<typename EntitySet::Element, Functor>;
  friend LocalFunction localFunction (const AnalyticDiscreteFunction& t)
252
  {
253
    return LocalFunction(t.f_, t.order_);
254
255
  }

256
257
  /// \brief obtain the stored \ref GridEntitySet
  const EntitySet& entitySet () const
258
  {
259
260
    return entitySet_;
  }
261

262
263
264
265
266
  /// \brief obtain the functor
  const Functor& f () const
  {
    return f_;
  }
267

268
269
270
private:
  Functor f_;
  int order_;
271

272
273
  EntitySet entitySet_;
};
274
275
276

} // end namespace Dune

277
#endif // DUNE_CURVEDSURFACEGRID_ANALYTIC_DISCRETEFUNCTION_HH