dune-functions  2.9.0
interpolate.hh
Go to the documentation of this file.
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 
9 #include <dune/common/exceptions.hh>
10 #include <dune/common/bitsetvector.hh>
11 
12 #include <dune/typetree/childextraction.hh>
13 
16 
22 
23 #include <dune/typetree/traversal.hh>
24 #include <dune/typetree/visitor.hh>
25 
26 namespace Dune {
27 namespace Functions {
28 
29 namespace Imp {
30 
31 struct 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 
57 template <class B, class T, class NTRE, class HV, class LF, class HBV>
58 class LocalInterpolateVisitor
59  : public TypeTree::TreeVisitor
60  , public TypeTree::DynamicTraversal
61 {
62 
63 public:
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 
169 protected:
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 
201 template <class B, class C, class F, class BV, class NTRE>
202 void 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 
273 template <class B, class C, class F, class BV>
274 void interpolate(const B& basis, C&& coeff, const F& f, const BV& bitVector)
275 {
276  interpolate(basis, coeff, f, bitVector, HierarchicNodeToRangeMap());
277 }
278 
293 template <class B, class C, class F>
294 void 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
auto istlVectorBackend(Vector &v)
Return a vector backend wrapping non-const ISTL like containers.
Definition: istlvectorbackend.hh:346
Definition: polynomial.hh:10
void interpolate(const B &basis, C &&coeff, const F &f, const BV &bv, const NTRE &nodeToRangeEntry)
Interpolate given function in discrete function space.
Definition: interpolate.hh:202
std::decay< F >::type makeGridViewFunction(F &&f, const GridView &gridView)
Construct a function modeling GridViewFunction from function and grid view.
Definition: gridviewfunction.hh:68
auto flatVectorView(T &t)
Create flat vector view of passed mutable container.
Definition: flatvectorview.hh:179
SizeInfo< Basis > sizeInfo(const Basis &basis)
Definition: sizeinfo.hh:69
Definition: backends/concepts.hh:21
Definition: backends/concepts.hh:31
A simple node to range map using the nested tree indices.
Definition: hierarchicnodetorangemap.hh:30