dune-functions  2.8.0
discreteglobalbasisfunction.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_GRIDFUNCTIONS_DISCRETEGLOBALBASISFUNCTIONS_HH
4 #define DUNE_FUNCTIONS_GRIDFUNCTIONS_DISCRETEGLOBALBASISFUNCTIONS_HH
5 
6 #include <memory>
7 
8 #include <dune/common/typetraits.hh>
9 
10 #include <dune/typetree/treecontainer.hh>
11 
18 
19 namespace Dune {
20 namespace Functions {
21 
22 
23 
66 template<typename B, typename V,
67  typename NTRE = HierarchicNodeToRangeMap,
68  typename R = typename V::value_type>
70 {
71 public:
72  using Basis = B;
73  using Vector = V;
74 
75  // In order to make the cache work for proxy-references
76  // we have to use AutonomousValue<T> instead of std::decay_t<T>
77  using Coefficient = Dune::AutonomousValue<decltype(std::declval<Vector>()[std::declval<typename Basis::MultiIndex>()])>;
78 
79  using GridView = typename Basis::GridView;
81  using Tree = typename Basis::LocalView::Tree;
82  using NodeToRangeEntry = NTRE;
83 
85  using Range = R;
86 
88  using Element = typename EntitySet::Element;
89 
90  using Traits = Imp::GridFunctionTraits<Range(Domain), EntitySet, DefaultDerivativeTraits, 16>;
91 
93  {
94  using LocalView = typename Basis::LocalView;
95  using size_type = typename Tree::size_type;
96 
97  template<class Node>
98  using LocalBasisRange = typename Node::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
99 
100  template<class Node>
101  using NodeData = typename std::vector<LocalBasisRange<Node>>;
102 
103  using PerNodeEvaluationBuffer = typename TypeTree::TreeContainer<NodeData,Tree>;
104 
105  public:
106 
111 
113  : globalFunction_(&globalFunction)
114  , localView_(globalFunction.basis().localView())
115  , evaluationBuffer_(localView_.tree())
116  , bound_(false)
117  {
118  localDoFs_.reserve(localView_.maxSize());
119  }
120 
122  : globalFunction_(other.globalFunction_)
123  , localView_(other.localView_)
124  , evaluationBuffer_(localView_.tree())
125  , bound_(other.bound_)
126  {
127  localDoFs_.reserve(localView_.maxSize());
128  if (bound())
129  localDoFs_ = other.localDoFs_;
130  }
131 
133  {
134  globalFunction_ = other.globalFunction_;
135  localView_ = other.localView_;
136  bound_ = other.bound_;
137  if (bound())
138  localDoFs_ = other.localDoFs_;
139  return *this;
140  }
141 
148  void bind(const Element& element)
149  {
150  localView_.bind(element);
151  bound_ = true;
152  // Use cache of full local view size. For a subspace basis,
153  // this may be larger than the number of local DOFs in the
154  // tree. In this case only cache entries associated to local
155  // DOFs in the subspace are filled. Cache entries associated
156  // to local DOFs which not contained in the subspace will not
157  // be touched.
158  //
159  // Alternatively one could use a cache that exactly fits
160  // the size of the tree. However, this would require to
161  // subtract an offset from localIndex(i) on each cache
162  // access in operator().
163  localDoFs_.resize(localView_.size());
164  for (size_type i = 0; i < localView_.tree().size(); ++i)
165  {
166  // For a subspace basis the index-within-tree i
167  // is not the same as the localIndex withn the
168  // full local view.
169  size_t localIndex = localView_.tree().localIndex(i);
170  localDoFs_[localIndex] = globalFunction_->dofs()[localView_.index(localIndex)];
171  }
172  }
173 
174  void unbind()
175  {
176  localView_.unbind();
177  bound_ = false;
178  }
179 
183  bool bound() const {
184  return bound_;
185  }
186 
196  Range operator()(const Domain& x) const
197  {
198  auto y = Range(0);
199 
200  TypeTree::forEachLeafNode(localView_.tree(), [&](auto&& node, auto&& treePath) {
201  const auto& nodeToRangeEntry = globalFunction_->nodeToRangeEntry();
202  const auto& fe = node.finiteElement();
203  const auto& localBasis = fe.localBasis();
204  auto& shapeFunctionValues = evaluationBuffer_[treePath];
205 
206  localBasis.evaluateFunction(x, shapeFunctionValues);
207 
208  // Get range entry associated to this node
209  auto re = flatVectorView(nodeToRangeEntry(node, treePath, y));
210 
211  for (size_type i = 0; i < localBasis.size(); ++i)
212  {
213  // Get coefficient associated to i-th shape function
214  auto c = flatVectorView(localDoFs_[node.localIndex(i)]);
215 
216  // Get value of i-th shape function
217  auto v = flatVectorView(shapeFunctionValues[i]);
218 
219  // Notice that the range entry re, the coefficient c, and the shape functions
220  // value v may all be scalar, vector, matrix, or general container valued.
221  // The matching of their entries is done via the multistage procedure described
222  // in the class documentation of DiscreteGlobalBasisFunction.
223  auto&& dimC = c.size();
224  auto dimV = v.size();
225  assert(dimC*dimV == re.size());
226  for(size_type j=0; j<dimC; ++j)
227  {
228  auto&& c_j = c[j];
229  for(size_type k=0; k<dimV; ++k)
230  re[j*dimV + k] += c_j*v[k];
231  }
232  }
233  });
234 
235  return y;
236  }
237 
238  const Element& localContext() const
239  {
240  return localView_.element();
241  }
242 
244  {
245  DUNE_THROW(NotImplemented,"not implemented");
246  }
247 
248  private:
249 
250  const DiscreteGlobalBasisFunction* globalFunction_;
251  LocalView localView_;
252  mutable PerNodeEvaluationBuffer evaluationBuffer_;
253  bool bound_ = false;
254  std::vector<Coefficient> localDoFs_;
255  };
256 
257  template<class B_T, class V_T, class NTRE_T>
258  DiscreteGlobalBasisFunction(B_T && basis, V_T && coefficients, NTRE_T&& nodeToRangeEntry) :
259  entitySet_(basis.gridView()),
260  basis_(wrap_or_move(std::forward<B_T>(basis))),
261  coefficients_(wrap_or_move(std::forward<V_T>(coefficients))),
262  nodeToRangeEntry_(wrap_or_move(std::forward<NTRE_T>(nodeToRangeEntry)))
263  {}
264 
265  DiscreteGlobalBasisFunction(std::shared_ptr<const Basis> basis, std::shared_ptr<const V> coefficients, std::shared_ptr<const NodeToRangeEntry> nodeToRangeEntry) :
266  entitySet_(basis->gridView()),
267  basis_(basis),
268  coefficients_(coefficients),
269  nodeToRangeEntry_(nodeToRangeEntry)
270  {}
271 
272  const Basis& basis() const
273  {
274  return *basis_;
275  }
276 
277  const V& dofs() const
278  {
279  return *coefficients_;
280  }
281 
283  {
284  return *nodeToRangeEntry_;
285  }
286 
287  // TODO: Implement this using hierarchic search
288  Range operator() (const Domain& x) const
289  {
290  DUNE_THROW(NotImplemented,"not implemented");
291  }
292 
294  {
295  DUNE_THROW(NotImplemented,"not implemented");
296  }
297 
314  {
315  return LocalFunction(t);
316  }
317 
321  const EntitySet& entitySet() const
322  {
323  return entitySet_;
324  }
325 
326 private:
327 
328  EntitySet entitySet_;
329  std::shared_ptr<const Basis> basis_;
330  std::shared_ptr<const V> coefficients_;
331  std::shared_ptr<const NodeToRangeEntry> nodeToRangeEntry_;
332 };
333 
334 
335 
346 template<typename... TT>
348 
349 
350 
351 template<typename R, typename B, typename V>
352 auto makeDiscreteGlobalBasisFunction(B&& basis, V&& vector)
353 {
354  using Basis = std::decay_t<B>;
355  using NTREM = HierarchicNodeToRangeMap;
356 
357  // Small helper functions to wrap vectors using istlVectorBackend
358  // if they do not already satisfy the VectorBackend interface.
359  auto toConstVectorBackend = [&](auto&& v) -> decltype(auto) {
360  if constexpr (models<Concept::ConstVectorBackend<Basis>, decltype(v)>()) {
361  return std::forward<decltype(v)>(v);
362  } else {
363  return istlVectorBackend(v);
364  }
365  };
366 
367  using Vector = std::decay_t<decltype(toConstVectorBackend(std::forward<V>(vector)))>;
369  std::forward<B>(basis),
370  toConstVectorBackend(std::forward<V>(vector)),
372 }
373 
374 
375 
376 } // namespace Functions
377 } // namespace Dune
378 
379 #endif // DUNE_FUNCTIONS_GRIDFUNCTIONS_DISCRETEGLOBALBASISFUNCTIONS_HH
friend LocalFunction localFunction(const DiscreteGlobalBasisFunction &t)
Construct local function from a DiscreteGlobalBasisFunction.
Definition: discreteglobalbasisfunction.hh:313
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
auto makeDiscreteGlobalBasisFunction(B &&basis, V &&vector)
Definition: discreteglobalbasisfunction.hh:352
Definition: backends/concepts.hh:21
Default implementation for derivative traits.
Definition: defaultderivativetraits.hh:37
A simple node to range map using the nested tree indices.
Definition: hierarchicnodetorangemap.hh:30
A grid function induced by a global basis and a coefficient vector.
Definition: discreteglobalbasisfunction.hh:70
Dune::AutonomousValue< decltype(std::declval< Vector >()[std::declval< typename Basis::MultiIndex >()])> Coefficient
Definition: discreteglobalbasisfunction.hh:77
Imp::GridFunctionTraits< Range(Domain), EntitySet, DefaultDerivativeTraits, 16 > Traits
Definition: discreteglobalbasisfunction.hh:90
DiscreteGlobalBasisFunction(B_T &&basis, V_T &&coefficients, NTRE_T &&nodeToRangeEntry)
Definition: discreteglobalbasisfunction.hh:258
friend Traits::DerivativeInterface derivative(const DiscreteGlobalBasisFunction &t)
Definition: discreteglobalbasisfunction.hh:293
const V & dofs() const
Definition: discreteglobalbasisfunction.hh:277
const Basis & basis() const
Definition: discreteglobalbasisfunction.hh:272
const EntitySet & entitySet() const
Get associated EntitySet.
Definition: discreteglobalbasisfunction.hh:321
DiscreteGlobalBasisFunction(std::shared_ptr< const Basis > basis, std::shared_ptr< const V > coefficients, std::shared_ptr< const NodeToRangeEntry > nodeToRangeEntry)
Definition: discreteglobalbasisfunction.hh:265
typename Basis::GridView GridView
Definition: discreteglobalbasisfunction.hh:79
GridViewEntitySet< GridView, 0 > EntitySet
Definition: discreteglobalbasisfunction.hh:80
typename EntitySet::LocalCoordinate LocalDomain
Definition: discreteglobalbasisfunction.hh:87
typename EntitySet::Element Element
Definition: discreteglobalbasisfunction.hh:88
B Basis
Definition: discreteglobalbasisfunction.hh:72
typename EntitySet::GlobalCoordinate Domain
Definition: discreteglobalbasisfunction.hh:84
R Range
Definition: discreteglobalbasisfunction.hh:85
V Vector
Definition: discreteglobalbasisfunction.hh:73
typename Basis::LocalView::Tree Tree
Definition: discreteglobalbasisfunction.hh:81
NTRE NodeToRangeEntry
Definition: discreteglobalbasisfunction.hh:82
const NodeToRangeEntry & nodeToRangeEntry() const
Definition: discreteglobalbasisfunction.hh:282
Definition: discreteglobalbasisfunction.hh:93
GlobalFunction::Element Element
Definition: discreteglobalbasisfunction.hh:110
void bind(const Element &element)
Bind LocalFunction to grid element.
Definition: discreteglobalbasisfunction.hh:148
bool bound() const
Check if LocalFunction is already bound to an element.
Definition: discreteglobalbasisfunction.hh:183
LocalDomain Domain
Definition: discreteglobalbasisfunction.hh:108
const Element & localContext() const
Definition: discreteglobalbasisfunction.hh:238
GlobalFunction::Range Range
Definition: discreteglobalbasisfunction.hh:109
friend Traits::LocalFunctionTraits::DerivativeInterface derivative(const LocalFunction &t)
Definition: discreteglobalbasisfunction.hh:243
void unbind()
Definition: discreteglobalbasisfunction.hh:174
LocalFunction & operator=(const LocalFunction &other)
Definition: discreteglobalbasisfunction.hh:132
LocalFunction(const DiscreteGlobalBasisFunction &globalFunction)
Definition: discreteglobalbasisfunction.hh:112
Range operator()(const Domain &x) const
Evaluate LocalFunction at bound element.
Definition: discreteglobalbasisfunction.hh:196
LocalFunction(const LocalFunction &other)
Definition: discreteglobalbasisfunction.hh:121
Definition: gridfunction.hh:32
GridView::template Codim< codim >::Entity Element
Type of Elements contained in this EntitySet.
Definition: gridviewentityset.hh:32
Element::Geometry::LocalCoordinate LocalCoordinate
Type of local coordinates with respect to the Element.
Definition: gridviewentityset.hh:35
Element::Geometry::GlobalCoordinate GlobalCoordinate
Definition: gridviewentityset.hh:36