dune-functions  2.7.1
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 #include <dune/common/hybridutilities.hh>
12 
13 #include <dune/typetree/childextraction.hh>
14 
15 #include <dune/localfunctions/common/virtualinterface.hh>
16 
20 
26 
27 #include <dune/typetree/traversal.hh>
28 #include <dune/typetree/visitor.hh>
29 
30 namespace Dune {
31 namespace Functions {
32 
33 namespace Imp {
34 
35 struct AllTrueBitSetVector
36 {
37  struct AllTrueBitSet
38  {
39  bool test(int i) const { return true; }
40  } allTrue_;
41 
42  operator bool() const
43  {
44  return true;
45  }
46 
47  template<class I>
48  const AllTrueBitSetVector& operator[](const I&) const
49  {
50  return *this;
51  }
52 
53  template<class SP>
54  void resize(const SP&) const
55  {}
56 
57 };
58 
59 
60 
61 template <class B, class T, class NTRE, class HV, class LF, class HBV>
62 class LocalInterpolateVisitor
63  : public TypeTree::TreeVisitor
64  , public TypeTree::DynamicTraversal
65 {
66 
67 public:
68 
69  using Basis = B;
70  using LocalView = typename B::LocalView;
71  using MultiIndex = typename LocalView::MultiIndex;
72 
73  using LocalFunction = LF;
74 
75  using Tree = T;
76 
77  using VectorBackend = HV;
78  using BitVectorBackend = HBV;
79 
80  using NodeToRangeEntry = NTRE;
81 
82  using GridView = typename Basis::GridView;
83  using Element = typename GridView::template Codim<0>::Entity;
84 
85  using LocalDomain = typename Element::Geometry::LocalCoordinate;
86 
87  using GlobalDomain = typename Element::Geometry::GlobalCoordinate;
88 
89  LocalInterpolateVisitor(const B& basis, HV& coeff, const HBV& bitVector, const LF& localF, const LocalView& localView, const NodeToRangeEntry& nodeToRangeEntry) :
90  vector_(coeff),
91  localF_(localF),
92  bitVector_(bitVector),
93  localView_(localView),
94  nodeToRangeEntry_(nodeToRangeEntry)
95  {
96  static_assert(Dune::Functions::Concept::isCallable<LocalFunction, LocalDomain>(), "Function passed to LocalInterpolateVisitor does not model the Callable<LocalCoordinate> concept");
97  }
98 
99  template<typename Node, typename TreePath>
100  void pre(Node& node, TreePath treePath)
101  {}
102 
103  template<typename Node, typename TreePath>
104  void post(Node& node, TreePath treePath)
105  {}
106 
107  template<typename Node, typename TreePath>
108  void leaf(Node& node, TreePath treePath)
109  {
110  using FiniteElement = typename Node::FiniteElement;
111  using FiniteElementRange = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
112  using FiniteElementRangeField = typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
113  using FunctionBaseClass = typename Dune::LocalFiniteElementFunctionBase<FiniteElement>::type;
114 
115  // Note that we capture j by reference. Hence we can switch
116  // the selected component later on by modifying j. Maybe we
117  // should avoid this naughty statefull lambda hack in favor
118  // of a separate helper class.
119  std::size_t j=0;
120  auto localFj = [&](const LocalDomain& x){
121  const auto& y = localF_(x);
122  return flatVectorView(nodeToRangeEntry_(node, treePath, y))[j];
123  };
124 
125  using FunctionFromCallable = typename Dune::Functions::FunctionFromCallable<FiniteElementRange(LocalDomain), decltype(localFj), FunctionBaseClass>;
126 
127  auto interpolationCoefficients = std::vector<FiniteElementRangeField>();
128 
129  auto&& fe = node.finiteElement();
130 
131  // We loop over j defined above and thus over the components of the
132  // range type of localF_.
133 
134  auto blockSize = flatVectorView(vector_[localView_.index(0)]).size();
135 
136  for(j=0; j<blockSize; ++j)
137  {
138 
139  // We cannot use localFj directly because interpolate requires a Dune::VirtualFunction like interface
140  fe.localInterpolation().interpolate(FunctionFromCallable(localFj), interpolationCoefficients);
141  for (size_t i=0; i<fe.localBasis().size(); ++i)
142  {
143  auto multiIndex = localView_.index(node.localIndex(i));
144  auto bitVectorBlock = flatVectorView(bitVector_[multiIndex]);
145  if (bitVectorBlock[j])
146  {
147  auto vectorBlock = flatVectorView(vector_[multiIndex]);
148  vectorBlock[j] = interpolationCoefficients[i];
149  }
150  }
151  }
152  }
153 
154 
155 protected:
156 
157  VectorBackend& vector_;
158  const LocalFunction& localF_;
159  const BitVectorBackend& bitVector_;
160  const LocalView& localView_;
161  const NodeToRangeEntry& nodeToRangeEntry_;
162 };
163 
164 
165 } // namespace Imp
166 
167 
168 
169 
187 template <class B, class C, class F, class BV, class NTRE>
188 void interpolate(const B& basis, C&& coeff, const F& f, const BV& bv, const NTRE& nodeToRangeEntry)
189 {
190  using GridView = typename B::GridView;
191  using Element = typename GridView::template Codim<0>::Entity;
192 
193  using Tree = typename B::LocalView::Tree;
194 
195  using GlobalDomain = typename Element::Geometry::GlobalCoordinate;
196 
197  static_assert(Dune::Functions::Concept::isCallable<F, GlobalDomain>(), "Function passed to interpolate does not model the Callable<GlobalCoordinate> concept");
198 
199  auto&& gridView = basis.gridView();
200 
201  // Small helper functions to wrap vectors using istlVectorBackend
202  // if they do not already satisfy the VectorBackend interface.
203  auto toVectorBackend = [&](auto& v) -> decltype(auto) {
204  return Dune::Hybrid::ifElse(models<Concept::VectorBackend<B>, decltype(v)>(),
205  [&](auto id) -> decltype(auto) {
206  return id(v);
207  }, [&](auto id) -> decltype(auto) {
208  return istlVectorBackend(id(v));
209  });
210  };
211  auto toConstVectorBackend = [&](auto& v) -> decltype(auto) {
212  return Dune::Hybrid::ifElse(models<Concept::ConstVectorBackend<B>, decltype(v)>(),
213  [&](auto id) -> decltype(auto) {
214  return id(v);
215  }, [&](auto id) -> decltype(auto) {
216  return istlVectorBackend(id(v));
217  });
218  };
219 
220  auto&& bitVector = toConstVectorBackend(bv);
221  auto&& vector = toVectorBackend(coeff);
222  vector.resize(sizeInfo(basis));
223 
224 
225 
226  // Make a grid function supporting local evaluation out of f
227  auto gf = makeGridViewFunction(f, gridView);
228 
229  // Obtain a local view of f
230  auto localF = localFunction(gf);
231 
232  auto localView = basis.localView();
233 
234  for (const auto& e : elements(gridView))
235  {
236  localView.bind(e);
237  localF.bind(e);
238 
239  Imp::LocalInterpolateVisitor<B, Tree, NTRE, decltype(vector), decltype(localF), decltype(bitVector)> localInterpolateVisitor(basis, vector, bitVector, localF, localView, nodeToRangeEntry);
240  TypeTree::applyToTree(localView.tree(),localInterpolateVisitor);
241  }
242 }
243 
260 template <class B, class C, class F, class BV>
261 void interpolate(const B& basis, C&& coeff, const F& f, const BV& bitVector)
262 {
263  interpolate(basis, coeff, f, bitVector, HierarchicNodeToRangeMap());
264 }
265 
280 template <class B, class C, class F>
281 void interpolate(const B& basis, C&& coeff, const F& f)
282 {
283  interpolate (basis, coeff, f, Imp::AllTrueBitSetVector(), HierarchicNodeToRangeMap());
284 }
285 
286 } // namespace Functions
287 } // namespace Dune
288 
289 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
void localFunction(DiscreteGlobalBasisFunction< TT... > &&t)=delete
Construction of local functions from a temporary DiscreteGlobalBasisFunction (forbidden)
auto istlVectorBackend(Vector &v)
Return a vector backend wrapping non-const ISTL like containers.
Definition: istlvectorbackend.hh:345
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:188
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:159
SizeInfo< Basis > sizeInfo(const Basis &basis)
Definition: sizeinfo.hh:69
Definition: backends/concepts.hh:21
Definition: backends/concepts.hh:31
Definition: functionfromcallable.hh:20
A simple node to range map using the nested tree indices.
Definition: hierarchicnodetorangemap.hh:30