DUNE PDELab (git)

interpolate.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
5// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
6
7#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
8#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
9
10#include <memory>
11#include <vector>
12
15#include <dune/common/referencehelper.hh>
16
17#include <dune/typetree/traversal.hh>
18
19#include <dune/functions/gridfunctions/gridviewfunction.hh>
20#include <dune/functions/common/functionconcepts.hh>
21
22#include <dune/functions/backends/concepts.hh>
23#include <dune/functions/backends/istlvectorbackend.hh>
24#include <dune/functions/functionspacebases/flatvectorview.hh>
25#include <dune/functions/functionspacebases/hierarchicnodetorangemap.hh>
26
27namespace Dune {
28namespace Functions {
29
30namespace Imp {
31
32struct AllTrueBitSetVector
33{
34 struct AllTrueBitSet
35 {
36 bool test(int) const { return true; }
37 } allTrue_;
38
39 operator bool() const
40 {
41 return true;
42 }
43
44 template<class I>
45 const AllTrueBitSetVector& operator[](const I&) const
46 {
47 return *this;
48 }
49
50 template<class SP>
51 void resize(const SP&) const
52 {}
53
54};
55
56
57
58// This helper function implements the restriction of some given function of type F.
59// The restriction is a simple callback that is applied to the values of the
60// function and the values of its derivative.
61template<class F, class Restriction>
62class ComponentFunction
63{
64public:
65
66 ComponentFunction(F f, Restriction restriction) :
67 f_(std::move(f)),
68 restriction_(std::move(restriction))
69 {}
70
71 template<class Domain>
72 auto operator()(const Domain& x) const
73 {
74 return restriction_(f_(x));
75 }
76
77 friend auto derivative(const ComponentFunction& cf)
78 {
79 // This provides support for capturing the derivative of the function by reference
80 // using forwardCapture for perfect forwarding capture. If the function caches its
81 // derivative, this saves a potentially costly copy.
82 auto&& df = derivative(Dune::resolveRef(cf.f_));
83 return [&, df=forwardCapture(std::forward<decltype(df)>(df))](auto&& x) {
84 return cf.restriction_(df.forward()(x));
85 };
86 }
87
88private:
89 F f_;
90 Restriction restriction_;
91};
92
93
94
95
96// This helper function implements caching of the derivative for local functions.
97// When using an algorithm that gets a LocalFunction and calls its derivative
98// on each element, this leads to a costly call of derivative(f). E.g. for a
99// DiscreteGlobalBasisFunction, this will allocate several buffer.
100// To avoid this, this helper function caches the derivative and hands
101// out the cached derivative by reference. To ensure that the handed out
102// derivative is appropriately bound, binding the function will automatically
103// bind the cached derivative.
104//
105// Notice that we cannot simply create the derivative in the constructor,
106// because this may throw for functions that do not implement the derivative.
107template<class F>
108class CachedDerivativeLocalFunction
109{
110 using Derivative = std::decay_t<decltype(derivative(Dune::resolveRef(std::declval<const F&>())))>;
111
112public:
113
114 CachedDerivativeLocalFunction(F f) :
115 f_(f)
116 {}
117
118 template<class Element>
119 void bind(const Element& element)
120 {
121 Dune::resolveRef(f_).bind(element);
122 if (derivative_)
123 derivative_.value().bind(element);
124 }
125
126 template<class X>
127 auto operator()(const X& x) const
128 {
129 return f_(x);
130 }
131
132 friend const Derivative& derivative(const CachedDerivativeLocalFunction& cdlf)
133 {
134 if (not cdlf.derivative_)
135 {
136 auto&& lf = Dune::resolveRef(cdlf.f_);
137 cdlf.derivative_ = derivative(lf);
138 if (lf.bound())
139 cdlf.derivative_.value().bind(lf.localContext());
140 }
141 return cdlf.derivative_.value();
142 }
143
144private:
145 F f_;
146 mutable std::optional<Derivative> derivative_;
147};
148
149
150
151template<class VectorBackend, class BitVectorBackend, class LocalFunction, class LocalView, class NodeToRangeEntry>
152void interpolateLocal(VectorBackend& vector, const BitVectorBackend& bitVector, const LocalFunction& localF, const LocalView& localView, const NodeToRangeEntry& nodeToRangeEntry)
153{
154 Dune::TypeTree::forEachLeafNode(localView.tree(), [&](auto&& node, auto&& treePath) {
155 using Node = std::decay_t<decltype(node)>;
156 using FiniteElement = typename Node::FiniteElement;
157 using FiniteElementRangeField = typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
158
159 auto interpolationCoefficients = std::vector<FiniteElementRangeField>();
160 auto&& fe = node.finiteElement();
161 auto localF_RE = ComponentFunction(std::cref(localF), [&](auto&& y) { return nodeToRangeEntry(node, treePath, y); });
162
163 fe.localInterpolation().interpolate(localF_RE, interpolationCoefficients);
164 for (size_t i=0; i<fe.localBasis().size(); ++i)
165 {
166 auto multiIndex = localView.index(node.localIndex(i));
167 if ( bitVector[multiIndex] )
168 vector[multiIndex] = interpolationCoefficients[i];
169 }
170 });
171}
172
173
174struct HasDerivative
175{
176 template<class F>
177 auto require(F&& f) -> decltype(derivative(f));
178};
179
180} // namespace Imp
181
182
183
184
202template <class B, class C, class F, class BV, class NTRE>
203void interpolate(const B& basis, C&& coeff, const F& f, const BV& bv, const NTRE& nodeToRangeEntry)
204{
205 using GridView = typename B::GridView;
206 using Element = typename GridView::template Codim<0>::Entity;
207 using GlobalDomain = typename Element::Geometry::GlobalCoordinate;
208
209 static_assert(Dune::Functions::Concept::isCallable<F, GlobalDomain>(), "Function passed to interpolate does not model the Callable<GlobalCoordinate> concept");
210
211 auto&& gridView = basis.gridView();
212
213 // Small helper functions to wrap vectors using istlVectorBackend
214 // if they do not already satisfy the VectorBackend interface.
215 auto toVectorBackend = [&](auto& v) -> decltype(auto) {
216 if constexpr (models<Concept::VectorBackend<B>, decltype(v)>()) {
217 return v;
218 } else {
219 return istlVectorBackend(v);
220 }
221 };
222
223 auto toConstVectorBackend = [&](auto& v) -> decltype(auto) {
224 if constexpr (models<Concept::ConstVectorBackend<B>, decltype(v)>()) {
225 return v;
226 } else {
227 return istlVectorBackend(v);
228 }
229 };
230
231 auto&& bitVector = toConstVectorBackend(bv);
232 auto&& vector = toVectorBackend(coeff);
233 vector.resize(basis);
234
235 // Make a grid function supporting local evaluation out of f
236 auto gf = makeGridViewFunction(f, gridView);
237
238 // Obtain a local view of f
239 // To avoid costly reconstruction of the derivative on each element,
240 // we use the CachedDerivativeLocalFunction wrapper if the function
241 // is differentiable. This wrapper will handout
242 // a reference to a single cached derivative object.
243 auto localF = [&](){
244 if constexpr (models<Imp::HasDerivative, decltype(localFunction(gf))>())
245 return Imp::CachedDerivativeLocalFunction(localFunction(gf));
246 else
247 return localFunction(gf);
248 }();
249
250 auto localView = basis.localView();
251
252 for (const auto& e : elements(gridView))
253 {
254 localView.bind(e);
255 localF.bind(e);
256 Imp::interpolateLocal(vector, bitVector, localF, localView, nodeToRangeEntry);
257 }
258}
259
276template <class B, class C, class F, class BV>
277void interpolate(const B& basis, C&& coeff, const F& f, const BV& bitVector)
278{
279 interpolate(basis, coeff, f, bitVector, HierarchicNodeToRangeMap());
280}
281
296template <class B, class C, class F>
297void interpolate(const B& basis, C&& coeff, const F& f)
298{
299 interpolate (basis, coeff, f, Imp::AllTrueBitSetVector(), HierarchicNodeToRangeMap());
300}
301
302} // namespace Functions
303} // namespace Dune
304
305#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
Efficient implementation of a dynamic array of static arrays of booleans.
A few common exception classes.
constexpr T & resolveRef(T &gf) noexcept
Helper function to resolve std::reference_wrapper.
Definition: referencehelper.hh:47
TrigonometricFunction< K, -cosFactor, sinFactor > derivative(const TrigonometricFunction< K, sinFactor, cosFactor > &f)
Obtain derivative of TrigonometricFunction function.
Definition: trigonometricfunction.hh:43
auto istlVectorBackend(Vector &v)
Return a vector backend wrapping non-const ISTL like containers.
Definition: istlvectorbackend.hh:350
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:177
constexpr auto treePath(const T &... t)
Constructs a new HybridTreePath from the given indices.
Definition: treepath.hh:326
void forEachLeafNode(Tree &&tree, LeafFunc &&leafFunc)
Traverse tree and visit each leaf node.
Definition: traversal.hh:269
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)