dune-functions  2.5-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 
34 {
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>
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 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  const auto& y_node = nodeToRangeEntry_(node, y);
123  using FunctionRange = typename std::decay<decltype(y_node)>::type;
125  };
126 
128 
129  auto interpolationValues = std::vector<FiniteElementRange>();
130 
131  auto&& fe = node.finiteElement();
132 
133  // We loop over j defined above and thus over the components of the
134  // range type of localF_.
135 
136  auto blockSize = FlatVectorBackend<CoefficientBlock>::size(vector_[localIndexSet_.index(0)]);
137 
138  for(j=0; j<blockSize; ++j)
139  {
140 
141  // We cannot use localFj directly because interpolate requires a Dune::VirtualFunction like interface
142  fe.localInterpolation().interpolate(FunctionFromCallable(localFj), interpolationValues);
143  for (size_t i=0; i<fe.localBasis().size(); ++i)
144  {
145  auto multiIndex = localIndexSet_.index(node.localIndex(i));
146  const auto& bitVectorBlock = bitVector_[multiIndex];
147  const auto& interpolateHere = FlatVectorBackend<BitVectorBlock>::getEntry(bitVectorBlock,j);
148 
149  if (interpolateHere)
150  {
151  auto&& vectorBlock = vector_[multiIndex];
152  FlatVectorBackend<CoefficientBlock>::getEntry(vectorBlock, j) = interpolationValues[i];
153  }
154  }
155  }
156  }
157 
158 
159 protected:
160 
166 };
167 
168 
169 } // namespace Imp
170 
171 
172 
173 
191 template <class B, class... TreeIndices, class NTRE, class C, class F, class BV>
192 void interpolateTreeSubset(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const NTRE& nodeToRangeEntry, const BV& bv)
193 {
194  using GridView = typename B::GridView;
195  using Element = typename GridView::template Codim<0>::Entity;
196 
197  using Tree = typename std::decay<decltype(TypeTree::child(basis.localView().tree(),treePath))>::type;
198 
199  using GlobalDomain = typename Element::Geometry::GlobalCoordinate;
200 
201  static_assert(Dune::Functions::Concept::isCallable<F, GlobalDomain>(), "Function passed to interpolate does not model the Callable<GlobalCoordinate> concept");
202 
203  auto&& gridView = basis.gridView();
204 
205  auto&& bitVector = makeHierarchicVectorForMultiIndex<typename B::MultiIndex>(bv);
206  auto&& vector = makeHierarchicVectorForMultiIndex<typename B::MultiIndex>(coeff);
207  vector.resize(sizeInfo(basis));
208 
209 
210 
211  // Make a grid function supporting local evaluation out of f
212  auto gf = makeGridViewFunction(f, gridView);
213 
214  // Obtain a local view of f
215  auto localF = localFunction(gf);
216 
217  auto localView = basis.localView();
218  auto localIndexSet = basis.localIndexSet();
219 
220  for (const auto& e : elements(gridView))
221  {
222  localView.bind(e);
223  localIndexSet.bind(localView);
224  localF.bind(e);
225 
226  auto&& subTree = TypeTree::child(localView.tree(),treePath);
227 
228  Imp::LocalInterpolateVisitor<B, Tree, NTRE, decltype(vector), decltype(localF), decltype(bitVector)> localInterpolateVisitor(basis, vector, bitVector, localF, localIndexSet, nodeToRangeEntry);
229  TypeTree::applyToTree(subTree,localInterpolateVisitor);
230  }
231 }
232 
233 
234 template <class B, class... TreeIndices, class C, class F, class BV>
235 void interpolateTreeSubset(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const BV& bitVector)
236 {
237  interpolateTreeSubset(basis, treePath, coeff, f, makeDefaultNodeToRangeMap(basis, treePath), bitVector);
238 }
239 
240 
241 template <class B, class... TreeIndices, class NTRE, class C, class F>
242 void interpolateTree(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const NTRE& nodeToRangeEntry)
243 {
244  interpolateTreeSubset(basis, treePath, coeff, f, nodeToRangeEntry, Imp::AllTrueBitSetVector());
245 }
246 
247 
248 template <class B, class... TreeIndices, class C, class F>
249 void interpolateTree(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f)
250 {
251  interpolateTreeSubset(basis, treePath, coeff, f, makeDefaultNodeToRangeMap(basis, treePath), Imp::AllTrueBitSetVector());
252 }
253 
254 
255 
272 template <class B, class... TreeIndices, class C, class F, class BV>
273 void interpolate(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f, const BV& bitVector)
274 {
275  interpolateTreeSubset(basis, treePath, coeff, f, makeDefaultNodeToRangeMap(basis, treePath), bitVector);
276 }
277 
278 namespace Imp {
279 
280  template<class... T>
281  constexpr std::true_type isHybridTreePath(const TypeTree::HybridTreePath<T...>&)
282  { return {}; }
283 
284  template<class T>
285  constexpr std::false_type isHybridTreePath(const T&)
286  { return {}; }
287 
288  template<class T>
289  constexpr auto isHybridTreePath() -> decltype(isHybridTreePath(std::declval<std::decay_t<T>>()))
290  { return {}; }
291 
292 }
293 
310 template <class B, class C, class F, class BV,
311  std::enable_if_t<not Imp::isHybridTreePath<C>(), int> = 0>
312 void interpolate(const B& basis, C&& coeff, const F& f, const BV& bitVector)
313 {
314  auto root = Dune::TypeTree::hybridTreePath();
315  interpolateTreeSubset(basis, root, coeff, f, makeDefaultNodeToRangeMap(basis, root), bitVector);
316 }
317 
318 
319 
334 template <class B, class C, class F>
335 void interpolate(const B& basis, C&& coeff, const F& f)
336 {
337  interpolate (basis, Dune::TypeTree::hybridTreePath(), coeff, f, Imp::AllTrueBitSetVector());
338 }
339 
355 template <class B, class... TreeIndices, class C, class F>
356 void interpolate(const B& basis, const TypeTree::HybridTreePath<TreeIndices...>& treePath, C&& coeff, const F& f)
357 {
358  interpolate (basis, treePath, coeff, f, Imp::AllTrueBitSetVector());
359 }
360 
361 } // namespace Functions
362 } // namespace Dune
363 
364 #endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_INTERPOLATE_HH
typename LocalIndexSet::MultiIndex MultiIndex
Definition: interpolate.hh:69
void interpolateTree(const B &basis, const TypeTree::HybridTreePath< TreeIndices... > &treePath, C &&coeff, const F &f, const NTRE &nodeToRangeEntry)
Definition: interpolate.hh:242
void resize(const SP &) const
Definition: interpolate.hh:52
Definition: functionfromcallable.hh:20
NTRE NodeToRangeEntry
Definition: interpolate.hh:78
void leaf(Node &node, TreePath treePath)
Definition: interpolate.hh:109
typename Basis::GridView GridView
Definition: interpolate.hh:80
const NodeToRangeEntry & nodeToRangeEntry_
Definition: interpolate.hh:165
typename std::decay< decltype(std::declval< HierarchicBitVector >()[std::declval< MultiIndex >()])>::type BitVectorBlock
Definition: interpolate.hh:88
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:273
typename GridView::template Codim< 0 >::Entity Element
Definition: interpolate.hh:81
const AllTrueBitSetVector & operator[](const I &) const
Definition: interpolate.hh:46
Definition: polynomial.hh:7
struct Dune::Functions::Imp::AllTrueBitSetVector::AllTrueBitSet allTrue_
void post(Node &node, TreePath treePath)
Definition: interpolate.hh:105
const LocalFunction & localF_
Definition: interpolate.hh:162
typename Element::Geometry::LocalCoordinate LocalDomain
Definition: interpolate.hh:83
const LocalIndexSet & localIndexSet_
Definition: interpolate.hh:164
std::decay< F >::type makeGridViewFunction(F &&f, const GridView &gridView)
Construct a function modeling GridViewFunction from function and grid view.
Definition: gridviewfunction.hh:68
typename Element::Geometry::GlobalCoordinate GlobalDomain
Definition: interpolate.hh:85
static auto size(VV &&v) -> decltype(v.size())
Definition: flatvectorbackend.hh:41
HV HierarchicVector
Definition: interpolate.hh:75
typename std::decay< decltype(std::declval< HierarchicVector >()[std::declval< MultiIndex >()])>::type CoefficientBlock
Definition: interpolate.hh:87
constexpr std::true_type isHybridTreePath(const TypeTree::HybridTreePath< T... > &)
Definition: interpolate.hh:281
HBV HierarchicBitVector
Definition: interpolate.hh:76
typename B::LocalIndexSet LocalIndexSet
Definition: interpolate.hh:68
bool test(int i) const
Definition: interpolate.hh:37
HierarchicVector & vector_
Definition: interpolate.hh:161
DefaultNodeToRangeMap< Tree > makeDefaultNodeToRangeMap(const Tree &tree)
Definition: defaultnodetorangemap.hh:107
SizeInfo< Basis > sizeInfo(const Basis &basis)
Definition: sizeinfo.hh:69
void pre(Node &node, TreePath treePath)
Definition: interpolate.hh:101
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:192
T Tree
Definition: interpolate.hh:73
Definition: interpolate.hh:33
B Basis
Definition: interpolate.hh:67
LF LocalFunction
Definition: interpolate.hh:71
LocalInterpolateVisitor(const B &basis, HV &coeff, const HBV &bitVector, const LF &localF, const LocalIndexSet &localIndexSet, const NodeToRangeEntry &nodeToRangeEntry)
Definition: interpolate.hh:90
const HierarchicBitVector & bitVector_
Definition: interpolate.hh:163
static auto getEntry(VV &&v, const Index &i) -> decltype(v[i])
Definition: flatvectorbackend.hh:25