dune-functions  2.6-dev
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 
14 #include <dune/localfunctions/common/virtualinterface.hh>
15 
19 
24 
25 #include <dune/typetree/traversal.hh>
26 #include <dune/typetree/visitor.hh>
27 
28 namespace Dune {
29 namespace Functions {
30 
31 namespace Imp {
32 
33 struct AllTrueBitSetVector
34 {
35  struct AllTrueBitSet
36  {
37  bool test(int i) const { return true; }
38  } allTrue_;
39 
40  operator bool() const
41  {
42  return true;
43  }
44 
45  template<class I>
46  const AllTrueBitSetVector& operator[](const I&) const
47  {
48  return *this;
49  }
50 
51  template<class SP>
52  void resize(const SP&) const
53  {}
54 
55 };
56 
57 
58 
59 template <class B, class T, class NTRE, class HV, class LF, class HBV>
60 class LocalInterpolateVisitor
61  : public TypeTree::TreeVisitor
62  , public TypeTree::DynamicTraversal
63 {
64 
65 public:
66 
67  using Basis = B;
68  using LocalIndexSet = typename B::LocalIndexSet;
69  using MultiIndex = typename LocalIndexSet::MultiIndex;
70 
71  using LocalFunction = LF;
72 
73  using Tree = T;
74 
75  using HierarchicVector = HV;
76  using HierarchicBitVector = HBV;
77 
78  using NodeToRangeEntry = NTRE;
79 
80  using GridView = typename Basis::GridView;
81  using Element = typename GridView::template Codim<0>::Entity;
82 
83  using LocalDomain = typename Element::Geometry::LocalCoordinate;
84 
85  using GlobalDomain = typename Element::Geometry::GlobalCoordinate;
86 
87  using CoefficientBlock = typename std::decay<decltype(std::declval<HierarchicVector>()[std::declval<MultiIndex>()])>::type;
88  using BitVectorBlock = typename std::decay<decltype(std::declval<HierarchicBitVector>()[std::declval<MultiIndex>()])>::type;
89 
90  LocalInterpolateVisitor(const B& basis, HV& coeff, const HBV& bitVector, const LF& localF, const LocalIndexSet& localIndexSet, const NodeToRangeEntry& nodeToRangeEntry) :
91  vector_(coeff),
92  localF_(localF),
93  bitVector_(bitVector),
94  localIndexSet_(localIndexSet),
95  nodeToRangeEntry_(nodeToRangeEntry)
96  {
97  static_assert(Dune::Functions::Concept::isCallable<LocalFunction, LocalDomain>(), "Function passed to LocalInterpolateVisitor does not model the Callable<LocalCoordinate> concept");
98  }
99 
100  template<typename Node, typename TreePath>
101  void pre(Node& node, TreePath treePath)
102  {}
103 
104  template<typename Node, typename TreePath>
105  void post(Node& node, TreePath treePath)
106  {}
107 
108  template<typename Node, typename TreePath>
109  void leaf(Node& node, TreePath treePath)
110  {
111  using FiniteElement = typename Node::FiniteElement;
112  using FiniteElementRange = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
113  using FiniteElementRangeField = typename FiniteElement::Traits::LocalBasisType::Traits::RangeFieldType;
114  using FunctionBaseClass = typename Dune::LocalFiniteElementFunctionBase<FiniteElement>::type;
115 
116  // Note that we capture j by reference. Hence we can switch
117  // the selected component later on by modifying j. Maybe we
118  // should avoid this naughty statefull lambda hack in favor
119  // of a separate helper class.
120  std::size_t j=0;
121  auto localFj = [&](const LocalDomain& x){
122  const auto& y = localF_(x);
123  const auto& y_node = nodeToRangeEntry_(node, y);
124  using FunctionRange = typename std::decay<decltype(y_node)>::type;
126  };
127 
128  using FunctionFromCallable = typename Dune::Functions::FunctionFromCallable<FiniteElementRange(LocalDomain), decltype(localFj), FunctionBaseClass>;
129 
130  auto interpolationCoefficients = std::vector<FiniteElementRangeField>();
131 
132  auto&& fe = node.finiteElement();
133 
134  // We loop over j defined above and thus over the components of the
135  // range type of localF_.
136 
137  auto blockSize = FlatVectorBackend<CoefficientBlock>::size(vector_[localIndexSet_.index(0)]);
138 
139  for(j=0; j<blockSize; ++j)
140  {
141 
142  // We cannot use localFj directly because interpolate requires a Dune::VirtualFunction like interface
143  fe.localInterpolation().interpolate(FunctionFromCallable(localFj), interpolationCoefficients);
144  for (size_t i=0; i<fe.localBasis().size(); ++i)
145  {
146  auto multiIndex = localIndexSet_.index(node.localIndex(i));
147  const auto& bitVectorBlock = bitVector_[multiIndex];
148  const auto& interpolateHere = FlatVectorBackend<BitVectorBlock>::getEntry(bitVectorBlock,j);
149 
150  if (interpolateHere)
151  {
152  auto&& vectorBlock = vector_[multiIndex];
153  FlatVectorBackend<CoefficientBlock>::getEntry(vectorBlock, j) = interpolationCoefficients[i];
154  }
155  }
156  }
157  }
158 
159 
160 protected:
161 
162  HierarchicVector& vector_;
163  const LocalFunction& localF_;
164  const HierarchicBitVector& bitVector_;
165  const LocalIndexSet& localIndexSet_;
166  const NodeToRangeEntry& nodeToRangeEntry_;
167 };
168 
169 
170 } // namespace Imp
171 
172 
173 
174 
192 template <class B, class... TreeIndices, class NTRE, class C, class F, class BV>
193 void interpolateTreeSubset(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const NTRE& nodeToRangeEntry, const BV& bv)
194 {
195  using GridView = typename B::GridView;
196  using Element = typename GridView::template Codim<0>::Entity;
197 
198  using Tree = typename std::decay<decltype(TypeTree::child(basis.localView().tree(),treePath))>::type;
199 
200  using GlobalDomain = typename Element::Geometry::GlobalCoordinate;
201 
202  static_assert(Dune::Functions::Concept::isCallable<F, GlobalDomain>(), "Function passed to interpolate does not model the Callable<GlobalCoordinate> concept");
203 
204  auto&& gridView = basis.gridView();
205 
206  auto&& bitVector = makeHierarchicVectorForMultiIndex<typename B::MultiIndex>(bv);
207  auto&& vector = makeHierarchicVectorForMultiIndex<typename B::MultiIndex>(coeff);
208  vector.resize(sizeInfo(basis));
209 
210 
211 
212  // Make a grid function supporting local evaluation out of f
213  auto gf = makeGridViewFunction(f, gridView);
214 
215  // Obtain a local view of f
216  auto localF = localFunction(gf);
217 
218  auto localView = basis.localView();
219  auto localIndexSet = basis.localIndexSet();
220 
221  for (const auto& e : elements(gridView))
222  {
223  localView.bind(e);
224  localIndexSet.bind(localView);
225  localF.bind(e);
226 
227  auto&& subTree = TypeTree::child(localView.tree(),treePath);
228 
229  Imp::LocalInterpolateVisitor<B, Tree, NTRE, decltype(vector), decltype(localF), decltype(bitVector)> localInterpolateVisitor(basis, vector, bitVector, localF, localIndexSet, nodeToRangeEntry);
230  TypeTree::applyToTree(subTree,localInterpolateVisitor);
231  }
232 }
233 
234 
235 template <class B, class... TreeIndices, class C, class F, class BV>
236 void interpolateTreeSubset(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const BV& bitVector)
237 {
238  interpolateTreeSubset(basis, treePath, coeff, f, makeDefaultNodeToRangeMap(basis, treePath), bitVector);
239 }
240 
241 
242 template <class B, class... TreeIndices, class NTRE, class C, class F>
243 void interpolateTree(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const NTRE& nodeToRangeEntry)
244 {
245  interpolateTreeSubset(basis, treePath, coeff, f, nodeToRangeEntry, Imp::AllTrueBitSetVector());
246 }
247 
248 
249 template <class B, class... TreeIndices, class C, class F>
250 void interpolateTree(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f)
251 {
252  interpolateTreeSubset(basis, treePath, coeff, f, makeDefaultNodeToRangeMap(basis, treePath), Imp::AllTrueBitSetVector());
253 }
254 
255 
256 
273 template <class B, class... TreeIndices, class C, class F, class BV>
274 void interpolate(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const BV& bitVector)
275 {
276  interpolateTreeSubset(basis, treePath, coeff, f, makeDefaultNodeToRangeMap(basis, treePath), bitVector);
277 }
278 
279 namespace Imp {
280 
281  template<class... T>
282  constexpr std::true_type isHybridTreePath(const TypeTree::HybridTreePath<T...>&)
283  { return {}; }
284 
285  template<class T>
286  constexpr std::false_type isHybridTreePath(const T&)
287  { return {}; }
288 
289  template<class T>
290  constexpr auto isHybridTreePath() -> decltype(isHybridTreePath(std::declval<std::decay_t<T>>()))
291  { return {}; }
292 
293 }
294 
311 template <class B, class C, class F, class BV,
312  std::enable_if_t<not Imp::isHybridTreePath<C>(), int> = 0>
313 void interpolate(const B& basis, C&& coeff, const F& f, const BV& bitVector)
314 {
315  auto root = Dune::TypeTree::hybridTreePath();
316  interpolateTreeSubset(basis, root, coeff, f, makeDefaultNodeToRangeMap(basis, root), bitVector);
317 }
318 
319 
320 
335 template <class B, class C, class F>
336 void interpolate(const B& basis, C&& coeff, const F& f)
337 {
338  interpolate (basis, Dune::TypeTree::hybridTreePath(), coeff, f, Imp::AllTrueBitSetVector());
339 }
340 
356 template <class B, class... TreeIndices, class C, class F>
357 void interpolate(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f)
358 {
359  interpolate (basis, treePath, coeff, f, Imp::AllTrueBitSetVector());
360 }
361 
362 } // namespace Functions
363 } // namespace Dune
364 
365 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
gridviewfunction.hh
hierarchicvectorwrapper.hh
Dune::Functions::localFunction
void localFunction(DiscreteGlobalBasisFunction< TT... > &&t)=delete
Construction of local functions from a temporary DiscreteGlobalBasisFunction (forbidden)
Dune
Definition: polynomial.hh:7
Dune::Functions::interpolateTree
void interpolateTree(const B &basis, const TypeTree::HybridTreePath< TreeIndices... > &treePath, C &&coeff, const F &f, const NTRE &nodeToRangeEntry)
Definition: interpolate.hh:243
Dune::Functions::makeDefaultNodeToRangeMap
DefaultNodeToRangeMap< Tree > makeDefaultNodeToRangeMap(const Tree &tree)
Definition: defaultnodetorangemap.hh:106
Dune::Functions::interpolateTreeSubset
void interpolateTreeSubset(const B &basis, const TypeTree::HybridTreePath< TreeIndices... > &treePath, C &&coeff, const F &f, const NTRE &nodeToRangeEntry, const BV &bv)
Interpolate given function in discrete function space.
Definition: interpolate.hh:193
sizeinfo.hh
Dune::Functions::FlatVectorBackend::getEntry
static auto getEntry(VV &&v, const Index &i) -> decltype(v[i])
Definition: flatvectorbackend.hh:25
functionfromcallable.hh
Dune::Functions::FunctionFromCallable
Definition: functionfromcallable.hh:20
Dune::Functions::interpolate
void interpolate(const B &basis, const TypeTree::HybridTreePath< TreeIndices... > &treePath, C &&coeff, const F &f, const BV &bitVector)
Interpolate given function in discrete function space.
Definition: interpolate.hh:274
flatvectorbackend.hh
Dune::Functions::sizeInfo
SizeInfo< Basis > sizeInfo(const Basis &basis)
Definition: sizeinfo.hh:69
Dune::Functions::FlatVectorBackend::size
static auto size(VV &&v) -> decltype(v.size())
Definition: flatvectorbackend.hh:41
functionconcepts.hh
Dune::Functions::makeGridViewFunction
std::decay< F >::type makeGridViewFunction(F &&f, const GridView &gridView)
Construct a function modeling GridViewFunction from function and grid view.
Definition: gridviewfunction.hh:68
defaultnodetorangemap.hh