Dune Core Modules (2.9.0)

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#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
5
6#include <memory>
7#include <vector>
8
11
12#include <dune/typetree/childextraction.hh>
13
14#include <dune/functions/gridfunctions/gridviewfunction.hh>
15#include <dune/functions/common/functionconcepts.hh>
16
17#include <dune/functions/backends/concepts.hh>
18#include <dune/functions/backends/istlvectorbackend.hh>
19#include <dune/functions/functionspacebases/sizeinfo.hh>
20#include <dune/functions/functionspacebases/flatvectorview.hh>
21#include <dune/functions/functionspacebases/hierarchicnodetorangemap.hh>
22
23#include <dune/typetree/traversal.hh>
24#include <dune/typetree/visitor.hh>
25
26namespace Dune {
27namespace Functions {
28
29namespace Imp {
30
31struct AllTrueBitSetVector
32{
33 struct AllTrueBitSet
34 {
35 bool test(int) const { return true; }
36 } allTrue_;
37
38 operator bool() const
39 {
40 return true;
41 }
42
43 template<class I>
44 const AllTrueBitSetVector& operator[](const I&) const
45 {
46 return *this;
47 }
48
49 template<class SP>
50 void resize(const SP&) const
51 {}
52
53};
54
55
56
57template <class B, class T, class NTRE, class HV, class LF, class HBV>
58class LocalInterpolateVisitor
59 : public TypeTree::TreeVisitor
60 , public TypeTree::DynamicTraversal
61{
62
63public:
64
65 using Basis = B;
66 using LocalView = typename B::LocalView;
67 using MultiIndex = typename LocalView::MultiIndex;
68
69 using LocalFunction = LF;
70
71 using Tree = T;
72
73 using VectorBackend = HV;
74 using BitVectorBackend = HBV;
75
76 using NodeToRangeEntry = NTRE;
77
78 using GridView = typename Basis::GridView;
79 using Element = typename GridView::template Codim<0>::Entity;
80
81 using LocalDomain = typename Element::Geometry::LocalCoordinate;
82
83 using GlobalDomain = typename Element::Geometry::GlobalCoordinate;
84
85 LocalInterpolateVisitor(const B& /*basis*/, HV& coeff, const HBV& bitVector, const LF& localF, const LocalView& localView, const NodeToRangeEntry& nodeToRangeEntry) :
86 vector_(coeff),
87 localF_(localF),
88 bitVector_(bitVector),
89 localView_(localView),
90 nodeToRangeEntry_(nodeToRangeEntry)
91 {
92 static_assert(Dune::Functions::Concept::isCallable<LocalFunction, LocalDomain>(), "Function passed to LocalInterpolateVisitor does not model the Callable<LocalCoordinate> concept");
93 }
94
95 template<typename Node, typename TreePath>
96 void pre(Node&, TreePath)
97 {}
98
99 template<typename Node, typename TreePath>
100 void post(Node&, TreePath)
101 {}
102
103 template<typename Node, typename TreePath>
104 void leaf(Node& node, TreePath treePath)
105 {
106 using FiniteElement = typename Node::FiniteElement;
107 using FiniteElementRange = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
108 using FiniteElementRangeField = typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
109
110 auto interpolationCoefficients = std::vector<FiniteElementRangeField>();
111 auto&& fe = node.finiteElement();
112
113 // backward compatibility: for scalar basis functions and possibly vector valued coefficients
114 // (like used in dune-fufem for power bases) loop over the components
115 // of the coefficients
116 if constexpr ( FiniteElement::Traits::LocalBasisType::Traits::dimRange == 1 )
117 {
118 // Note that we capture j by reference. Hence we can switch
119 // the selected component later on by modifying j. Maybe we
120 // should avoid this naughty statefull lambda hack in favor
121 // of a separate helper class.
122 std::size_t j=0;
123 auto localFj = [&](const LocalDomain& x){
124 const auto& y = localF_(x);
125 return FiniteElementRange(flatVectorView(nodeToRangeEntry_(node, treePath, y))[j]);
126 };
127
128 // We loop over j defined above and thus over the components of the
129 // range type of localF_.
130
131 auto blockSize = flatVectorView(vector_[localView_.index(0)]).size();
132
133 for(j=0; j<blockSize; ++j)
134 {
135 fe.localInterpolation().interpolate(localFj, interpolationCoefficients);
136 for (size_t i=0; i<fe.localBasis().size(); ++i)
137 {
138 auto multiIndex = localView_.index(node.localIndex(i));
139 auto bitVectorBlock = flatVectorView(bitVector_[multiIndex]);
140 if (bitVectorBlock[j])
141 {
142 auto vectorBlock = flatVectorView(vector_[multiIndex]);
143 vectorBlock[j] = interpolationCoefficients[i];
144 }
145 }
146 }
147 }
148 else // ( FiniteElement::Traits::LocalBasisType::Traits::dimRange != 1 )
149 {
150 // for all other finite elements: use the FiniteElementRange directly for the interpolation
151 auto localF = [&](const LocalDomain& x){
152 const auto& y = localF_(x);
153 return FiniteElementRange(nodeToRangeEntry_(node, treePath, y));
154 };
155
156 fe.localInterpolation().interpolate(localF, interpolationCoefficients);
157 for (size_t i=0; i<fe.localBasis().size(); ++i)
158 {
159 auto multiIndex = localView_.index(node.localIndex(i));
160 if ( bitVector_[multiIndex] )
161 {
162 vector_[multiIndex] = interpolationCoefficients[i];
163 }
164 }
165 }
166 }
167
168
169protected:
170
171 VectorBackend& vector_;
172 const LocalFunction& localF_;
173 const BitVectorBackend& bitVector_;
174 const LocalView& localView_;
175 const NodeToRangeEntry& nodeToRangeEntry_;
176};
177
178
179} // namespace Imp
180
181
182
183
201template <class B, class C, class F, class BV, class NTRE>
202void interpolate(const B& basis, C&& coeff, const F& f, const BV& bv, const NTRE& nodeToRangeEntry)
203{
204 using GridView = typename B::GridView;
205 using Element = typename GridView::template Codim<0>::Entity;
206
207 using Tree = typename B::LocalView::Tree;
208
209 using GlobalDomain = typename Element::Geometry::GlobalCoordinate;
210
211 static_assert(Dune::Functions::Concept::isCallable<F, GlobalDomain>(), "Function passed to interpolate does not model the Callable<GlobalCoordinate> concept");
212
213 auto&& gridView = basis.gridView();
214
215 // Small helper functions to wrap vectors using istlVectorBackend
216 // if they do not already satisfy the VectorBackend interface.
217 auto toVectorBackend = [&](auto& v) -> decltype(auto) {
218 if constexpr (models<Concept::VectorBackend<B>, decltype(v)>()) {
219 return v;
220 } else {
221 return istlVectorBackend(v);
222 }
223 };
224
225 auto toConstVectorBackend = [&](auto& v) -> decltype(auto) {
226 if constexpr (models<Concept::ConstVectorBackend<B>, decltype(v)>()) {
227 return v;
228 } else {
229 return istlVectorBackend(v);
230 }
231 };
232
233 auto&& bitVector = toConstVectorBackend(bv);
234 auto&& vector = toVectorBackend(coeff);
235 vector.resize(sizeInfo(basis));
236
237
238
239 // Make a grid function supporting local evaluation out of f
240 auto gf = makeGridViewFunction(f, gridView);
241
242 // Obtain a local view of f
243 auto localF = localFunction(gf);
244
245 auto localView = basis.localView();
246
247 for (const auto& e : elements(gridView))
248 {
249 localView.bind(e);
250 localF.bind(e);
251
252 Imp::LocalInterpolateVisitor<B, Tree, NTRE, decltype(vector), decltype(localF), decltype(bitVector)> localInterpolateVisitor(basis, vector, bitVector, localF, localView, nodeToRangeEntry);
253 TypeTree::applyToTree(localView.tree(),localInterpolateVisitor);
254 }
255}
256
273template <class B, class C, class F, class BV>
274void interpolate(const B& basis, C&& coeff, const F& f, const BV& bitVector)
275{
276 interpolate(basis, coeff, f, bitVector, HierarchicNodeToRangeMap());
277}
278
293template <class B, class C, class F>
294void interpolate(const B& basis, C&& coeff, const F& f)
295{
296 interpolate (basis, coeff, f, Imp::AllTrueBitSetVector(), HierarchicNodeToRangeMap());
297}
298
299} // namespace Functions
300} // namespace Dune
301
302#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
Efficient implementation of a dynamic array of static arrays of booleans.
A few common exception classes.
auto istlVectorBackend(Vector &v)
Return a vector backend wrapping non-const ISTL like containers.
Definition: istlvectorbackend.hh:346
constexpr HybridTreePath< T... > treePath(const T &... t)
Constructs a new HybridTreePath from the given indices.
Definition: treepath.hh:191
void applyToTree(Tree &&tree, Visitor &&visitor)
Apply visitor to TypeTree.
Definition: traversal.hh:237
Dune namespace.
Definition: alignedallocator.hh:13
void post(T &&, TreePath) const
Method for postfix tree traversal.
Definition: visitor.hh:81
void leaf(T &&, TreePath) const
Method for leaf traversal.
Definition: visitor.hh:91
void pre(T &&, TreePath) const
Method for prefix tree traversal.
Definition: visitor.hh:58
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 21, 23:30, 2024)