3#ifndef DUNE_PDELAB_GRIDFUNCTIONSPACE_LFSINDEXCACHE_HH
4#define DUNE_PDELAB_GRIDFUNCTIONSPACE_LFSINDEXCACHE_HH
9#include <unordered_map>
16#include <dune/typetree/typetree.hh>
18#include <dune/pdelab/constraints/common/constraintstransformation.hh>
19#include <dune/pdelab/gridfunctionspace/tags.hh>
24 template<
typename Iterator>
25 class DOFIndexViewIterator
26 :
public RandomAccessIteratorFacade<DOFIndexViewIterator<Iterator>,
27 const typename std::iterator_traits<Iterator>::value_type::View,
28 const typename std::iterator_traits<Iterator>::value_type::View
32 friend class RandomAccessIteratorFacade<
34 const typename
std::iterator_traits<Iterator>::value_type::View,
35 const typename std::iterator_traits<Iterator>::value_type::View
38 typedef typename std::iterator_traits<Iterator>::value_type::View View;
60 explicit proxy(
const View& v)
73 typedef proxy pointer;
75 DOFIndexViewIterator()
80 explicit DOFIndexViewIterator(Iterator it, std::size_t tail_length = 0)
82 , _tail_length(tail_length)
95 const typename std::iterator_traits<Iterator>::reference raw_index()
const
100 bool equals(
const DOFIndexViewIterator& other)
const
102 return _iterator == other._iterator;
120 std::ptrdiff_t distanceTo(DOFIndexViewIterator& other)
const
122 return other._iterator - _iterator;
125 const View dereference()
const
127 return _iterator->view(_iterator->treeIndex().size() - _tail_length);
130 pointer operator->()
const
132 return pointer(dereference());
138 std::size_t _tail_length;
143 template<
typename Iterator>
144 struct extract_lfs_leaf_size_visitor
145 :
public TypeTree::TreeVisitor
146 ,
public TypeTree::DynamicTraversal
149 template<
typename LeafLFS,
typename TreePath>
150 void leaf(
const LeafLFS& leaf_lfs, TreePath tp)
152 (*it) = leaf_lfs.size();
156 extract_lfs_leaf_size_visitor(Iterator leaf_size_container_iterator)
157 : it(leaf_size_container_iterator)
164 template<
typename LFS,
typename Iterator>
165 Iterator extract_lfs_leaf_sizes(
const LFS& lfs, Iterator it)
167 extract_lfs_leaf_size_visitor<Iterator> visitor(it);
173 template<
typename DOFIterator,
174 typename ContainerIterator,
175 typename LeafSizeIterator,
176 std::size_t tree_depth,
178 struct map_dof_indices_to_container_indices
179 :
public TypeTree::TreeVisitor
180 ,
public TypeTree::DynamicTraversal
183 template<
typename Ordering,
typename TreePath>
184 void leaf(
const Ordering& ordering, TreePath tp)
186 std::size_t leaf_size = *(leaf_size_pos++);
190 dof_end += leaf_size;
191 ordering.map_lfs_indices(dof_pos,dof_end,container_pos);
193 container_pos += leaf_size;
196 template<
typename Ordering,
typename TreePath>
197 void post(
const Ordering& ordering, TreePath tp)
199 if (Ordering::consume_tree_index)
201 dof_pos.restore_back();
202 dof_end.restore_back();
204 ordering.map_lfs_indices(dof_stack.top(),dof_end,container_stack.top());
206 container_stack.pop();
209 template<
typename Ordering,
typename TreePath>
210 void pre(
const Ordering& ordering, TreePath tp)
212 dof_stack.push(dof_pos);
213 container_stack.push(container_pos);
214 if (Ordering::consume_tree_index)
221 map_dof_indices_to_container_indices(DOFIterator dof_begin,
222 ContainerIterator container_begin,
223 LeafSizeIterator leaf_size_begin,
224 std::size_t dof_index_tail_length = 0)
225 : dof_pos(dof_begin,dof_index_tail_length)
226 , dof_end(dof_begin,dof_index_tail_length)
227 , container_pos(container_begin)
228 , leaf_size_pos(leaf_size_begin)
232 DOFIndexViewIterator<DOFIterator> dof_pos;
233 DOFIndexViewIterator<DOFIterator> dof_end;
234 ContainerIterator container_pos;
235 LeafSizeIterator leaf_size_pos;
236 std::stack<DOFIndexViewIterator<DOFIterator>,ReservedVector<DOFIndexViewIterator<DOFIterator>,tree_depth> > dof_stack;
237 std::stack<ContainerIterator,ReservedVector<ContainerIterator,tree_depth> > container_stack;
243 template<
typename LFS,
typename C,
typename CacheTag,
bool fast>
244 class LFSIndexCacheBase
249 DOF_NONCONSTRAINED = 0,
250 DOF_CONSTRAINED = 1<<0,
256 typedef LFS LocalFunctionSpace;
258 typedef typename LFS::Traits::GridFunctionSpace GFS;
259 typedef typename GFS::Ordering Ordering;
260 typedef typename Ordering::Traits::ContainerIndex ContainerIndex;
261 typedef ContainerIndex CI;
262 typedef typename Ordering::Traits::DOFIndex DOFIndex;
264 typedef std::size_t size_type;
266 typedef std::vector<CI> CIVector;
267 typedef std::unordered_map<DI,CI> CIMap;
269 typedef std::unordered_map<const CI*,std::pair<size_type,bool> > InverseMap;
271 struct ConstraintsEntry
272 :
public std::pair<const CI*,typename C::mapped_type::mapped_type>
274 typedef CI ContainerIndex;
275 typedef typename C::mapped_type::mapped_type Weight;
277 const ContainerIndex& containerIndex()
const
279 return *(this->first);
282 const Weight& weight()
const
290 typedef std::vector<ConstraintsEntry> ConstraintsVector;
291 typedef typename ConstraintsVector::const_iterator ConstraintsIterator;
293 LFSIndexCacheBase(
const LFS& lfs,
const C&
constraints,
bool enable_constraints_caching)
295 , _enable_constraints_caching(enable_constraints_caching)
296 , _container_indices(lfs.maxSize())
297 , _dof_flags(lfs.maxSize(),0)
298 , _constraints_iterators(lfs.maxSize())
299 , _inverse_cache_built(false)
307 _container_indices[0].resize(2);
308 _container_indices[0][0] = 0;
309 _container_indices[0][1] = LFS::Traits::GridFunctionSpace::Ordering::Traits::DOFIndexAccessor::entityIndex(_lfs.dofIndex(0));
314 _container_index_map.clear();
315 for (
typename CIVector::iterator it = _container_indices.begin(); it != _container_indices.end(); ++it)
318 _inverse_map.clear();
319 _inverse_cache_built =
false;
322 typedef ReservedVector<size_type,TypeTree::TreeInfo<LFS>::leafCount> LeafSizeVector;
323 LeafSizeVector leaf_sizes;
325 extract_lfs_leaf_sizes(_lfs,leaf_sizes.begin());
328 map_dof_indices_to_container_indices<
329 typename LFS::Traits::DOFIndexContainer::const_iterator,
330 typename CIVector::iterator,
331 typename LeafSizeVector::const_iterator,
334 > index_mapper(_lfs._dof_indices->begin(),_container_indices.begin(),leaf_sizes.begin(),_lfs.subSpaceDepth());
337 if (_enable_constraints_caching)
339 _constraints.resize(0);
340 std::vector<std::pair<size_type,typename C::const_iterator> > non_dirichlet_constrained_dofs;
341 size_type constraint_entry_count = 0;
342 size_type end = fast ? 1 : _lfs.size();
343 for (size_type i = 0; i < end; ++i)
345 const CI& container_index = _container_indices[i];
346 const typename C::const_iterator cit = _gfs_constraints.find(container_index);
347 if (cit == _gfs_constraints.end())
349 _dof_flags[i] = DOF_NONCONSTRAINED;
353 if (cit->second.size() == 0)
355 _dof_flags[i] = DOF_CONSTRAINED | DOF_DIRICHLET;
356 _constraints_iterators[i] = make_pair(_constraints.end(),_constraints.end());
360 _dof_flags[i] = DOF_CONSTRAINED;
361 constraint_entry_count += cit->second.size();
362 non_dirichlet_constrained_dofs.push_back(make_pair(i,cit));
366 if (constraint_entry_count > 0)
368 _constraints.resize(constraint_entry_count);
369 typename ConstraintsVector::iterator eit = _constraints.begin();
370 for (
typename std::vector<std::pair<size_type,typename C::const_iterator> >::const_iterator it = non_dirichlet_constrained_dofs.begin();
371 it != non_dirichlet_constrained_dofs.end();
374 _constraints_iterators[it->first].first = eit;
375 for (
typename C::mapped_type::const_iterator cit = it->second->second.begin(); cit != it->second->second.end(); ++cit, ++eit)
377 eit->first = &(cit->first);
378 eit->second = cit->second;
380 _constraints_iterators[it->first].second = eit;
387 const DI& dofIndex(size_type i)
const
389 return _lfs.dofIndex(i);
392 const CI& containerIndex(size_type i)
const
394 return _container_indices[i];
397 const CI& containerIndex(
const DI& i)
const
400 std::pair<typename CIMap::iterator,bool> r = _container_index_map.insert(std::make_pair(std::ref(i),CI()));
404 _lfs.gridFunctionSpace().ordering().mapIndex(i.view(),r.first->second);
407 return r.first->second;
410 bool isConstrained(size_type i)
const
412 return _dof_flags[i] & DOF_CONSTRAINED;
415 bool isDirichletConstraint(size_type i)
const
417 return _dof_flags[i] & DOF_DIRICHLET;
420 ConstraintsIterator constraintsBegin(size_type i)
const
422 assert(isConstrained(i));
423 return _constraints_iterators[i].first;
426 ConstraintsIterator constraintsEnd(size_type i)
const
428 assert(isConstrained(i));
429 return _constraints_iterators[i].second;
432 const LocalFunctionSpace& localFunctionSpace()
const
437 size_type
size()
const
442 std::pair<size_type,bool> localIndex(
const ContainerIndex& ci)
const
444 if (!_inverse_cache_built)
445 build_inverse_cache();
446 return _inverse_map[ci];
449 size_type offset(size_type i)
const
451 if (!_inverse_cache_built)
452 build_inverse_cache();
456 size_type extendedOffset(size_type i)
const
458 if (!_inverse_cache_built)
459 build_inverse_cache();
460 return _extended_offsets[i];
463 bool constraintsCachingEnabled()
const
465 return _enable_constraints_caching;
470 struct sort_container_indices
473 bool operator()(
const T* a,
const T* b)
const
475 return std::lexicographical_compare(reversed_iterator(a->end()),reversed_iterator(a->begin()),
476 reversed_iterator(b->end()),reversed_iterator(b->begin())
482 void build_inverse_cache()
487 for (
typename CIVector::const_iterator it = _container_indices.begin(),
488 endit = _container_indices.end();
493 _inverse_map.insert(std::make_pair(&(*it),std::make_pair(i,
false)));
494 if (it->back() !=
child)
496 _offsets[
child+1] = i;
501 std::vector<const ContainerIndex*> extended_cis;
502 extended_cis.reserve(_constraints.size());
504 for (
typename ConstraintsVector::const_iterator it = _constraints.begin(),
505 endit = _constraints.end();
510 if (_inverse_map.count(it->first) == 0)
511 extended_cis.push_back(it->first);
514 std::sort(extended_cis.begin(),extended_cis.end(),sort_container_indices());
516 typename std::vector<const ContainerIndex*>::const_iterator endit = std::unique(extended_cis.begin(),extended_cis.end());
520 for (
typename std::vector<const ContainerIndex*>::const_iterator it = extended_cis.begin(); it != endit; ++it, ++i)
522 _inverse_map.insert(std::make_pair(&(*it),std::make_pair(i,
true)));
523 if (it->back() !=
child)
525 _extended_offsets[
child+1] = i;
530 _inverse_cache_built =
true;
535 const bool _enable_constraints_caching;
536 CIVector _container_indices;
537 std::vector<unsigned char> _dof_flags;
538 std::vector<std::pair<ConstraintsIterator,ConstraintsIterator> > _constraints_iterators;
539 mutable CIMap _container_index_map;
540 ConstraintsVector _constraints;
541 mutable std::array<size_type,TypeTree::StaticDegree<LFS>::value> _offsets;
542 mutable std::array<size_type,TypeTree::StaticDegree<LFS>::value> _extended_offsets;
543 mutable bool _inverse_cache_built;
544 mutable InverseMap _inverse_map;
546 const C& _gfs_constraints;
551 template<
typename LFS,
typename CacheTag,
bool fast>
552 class LFSIndexCacheBase<LFS,EmptyTransformation,CacheTag,fast>
557 typedef LFS LocalFunctionSpace;
558 typedef typename LFS::Traits::GridFunctionSpace GFS;
559 typedef typename GFS::Ordering Ordering;
560 typedef typename Ordering::Traits::ContainerIndex ContainerIndex;
561 typedef ContainerIndex CI;
562 typedef typename Ordering::Traits::DOFIndex DOFIndex;
564 typedef std::size_t size_type;
566 typedef std::vector<CI> CIVector;
567 typedef std::unordered_map<DI,CI> CIMap;
569 struct ConstraintsEntry
570 :
public std::pair<const CI*,double>
572 typedef CI ContainerIndex;
573 typedef double Weight;
575 const ContainerIndex& containerIndex()
const
577 return *(this->first);
580 const Weight& weight()
const
586 typedef std::vector<ConstraintsEntry> ConstraintsVector;
587 typedef typename ConstraintsVector::const_iterator ConstraintsIterator;
589 explicit LFSIndexCacheBase(
const LFS& lfs)
591 , _container_indices(lfs.maxSize())
596 LFSIndexCacheBase(
const LFS& lfs,
const C& c,
bool enable_constraints_caching)
598 , _container_indices(lfs.maxSize())
606 _container_indices[0].resize(2);
607 _container_indices[0][0] = 0;
608 _container_indices[0][1] = LFS::Traits::GridFunctionSpace::Ordering::Traits::DOFIndexAccessor::entityIndex(_lfs.dofIndex(0));
613 _container_index_map.clear();
614 for (
typename CIVector::iterator it = _container_indices.begin(); it != _container_indices.end(); ++it)
618 typedef ReservedVector<size_type,TypeTree::TreeInfo<LFS>::leafCount> LeafSizeVector;
619 LeafSizeVector leaf_sizes;
621 extract_lfs_leaf_sizes(_lfs,leaf_sizes.begin());
624 map_dof_indices_to_container_indices<
625 typename LFS::Traits::DOFIndexContainer::const_iterator,
626 typename CIVector::iterator,
627 typename LeafSizeVector::const_iterator,
630 > index_mapper(_lfs._dof_indices->begin(),_container_indices.begin(),leaf_sizes.begin(),_lfs.subSpaceDepth());
635 const DI& dofIndex(size_type i)
const
637 return _lfs.dofIndex(i);
640 const CI& containerIndex(size_type i)
const
642 return _container_indices[i];
645 const CI& containerIndex(
const DI& i)
const
648 std::pair<typename CIMap::iterator,bool> r = _container_index_map.insert(std::make_pair(std::ref(i),CI()));
652 _lfs.gridFunctionSpace().ordering().mapIndex(i.view(),r.first->second);
655 return r.first->second;
658 bool isConstrained(size_type i)
const
663 bool isDirichletConstraint(size_type i)
const
668 ConstraintsIterator constraintsBegin(size_type i)
const
670 return _constraints.begin();
673 ConstraintsIterator constraintsEnd(size_type i)
const
675 return _constraints.end();
678 const LocalFunctionSpace& localFunctionSpace()
const
683 size_type
size()
const
688 bool constraintsCachingEnabled()
const
696 CIVector _container_indices;
697 mutable CIMap _container_index_map;
698 const ConstraintsVector _constraints;
704 template<
typename LFS,
typename C,
bool fast>
705 class LFSIndexCacheBase<LFS,C,SimpleLFSCacheTag,fast>
710 DOF_NONCONSTRAINED = 0,
711 DOF_CONSTRAINED = 1<<0,
717 typedef LFS LocalFunctionSpace;
719 typedef typename LFS::Traits::GridFunctionSpace GFS;
720 typedef typename GFS::Ordering Ordering;
721 typedef typename Ordering::Traits::ContainerIndex CI;
722 typedef typename Ordering::Traits::DOFIndex DI;
723 typedef std::size_t size_type;
725 typedef std::vector<CI> CIVector;
726 typedef std::unordered_map<DI,CI> CIMap;
728 struct ConstraintsEntry
729 :
public std::pair<CI,typename C::mapped_type::mapped_type>
731 typedef CI ContainerIndex;
732 typedef typename C::mapped_type::mapped_type Weight;
734 const ContainerIndex& containerIndex()
const
739 const Weight& weight()
const
745 typedef std::vector<ConstraintsEntry> ConstraintsVector;
746 typedef typename ConstraintsVector::const_iterator ConstraintsIterator;
748 LFSIndexCacheBase(
const LFS& lfs,
const C&
constraints)
750 , _dof_flags(lfs.maxSize())
751 , _constraints_iterators(lfs.maxSize())
763 _constraints.resize(0);
764 std::vector<std::pair<size_type,typename C::const_iterator> > non_dirichlet_constrained_dofs;
765 size_type constraint_entry_count = 0;
766 for (size_type i = 0; i < _lfs.size(); ++i)
768 const DI& dof_index = _lfs.dofIndex(i);
769 const typename C::const_iterator cit = _gfs_constraints.find(dof_index);
770 if (cit == _gfs_constraints.end())
772 _dof_flags[i] = DOF_NONCONSTRAINED;
776 if (cit->second.size() == 0)
778 _dof_flags[i] = DOF_CONSTRAINED | DOF_DIRICHLET;
779 _constraints_iterators[i] = make_pair(_constraints.end(),_constraints.end());
783 _dof_flags[i] = DOF_CONSTRAINED;
784 constraint_entry_count += cit->second.size();
785 non_dirichlet_constrained_dofs.push_back(make_pair(i,cit));
789 if (constraint_entry_count > 0)
791 _constraints.resize(constraint_entry_count);
792 typename ConstraintsVector::iterator eit = _constraints.begin();
793 for (
typename std::vector<std::pair<size_type,typename C::const_iterator> >::const_iterator it = non_dirichlet_constrained_dofs.begin();
794 it != non_dirichlet_constrained_dofs.end();
797 _constraints_iterators[it->first].first = eit;
798 for (
typename C::mapped_type::const_iterator cit = it->second->second.begin(); cit != it->second->second.end(); ++cit, ++eit)
800 eit->first = cit->first;
801 eit->second = cit->second;
803 _constraints_iterators[it->first].second = eit;
809 const DI& dofIndex(size_type i)
const
811 return _lfs.dofIndex(i);
814 CI containerIndex(size_type i)
const
816 return CI(_lfs.dofIndex(i)[0]);
819 const CI& containerIndex(
const DI& i)
const
824 bool isConstrained(size_type i)
const
826 return _dof_flags[i] & DOF_CONSTRAINED;
829 bool isDirichletConstraint(size_type i)
const
831 return _dof_flags[i] & DOF_DIRICHLET;
834 ConstraintsIterator constraintsBegin(size_type i)
const
836 assert(isConstrained(i));
837 return _constraints_iterators[i].first;
840 ConstraintsIterator constraintsEnd(size_type i)
const
842 assert(isConstrained(i));
843 return _constraints_iterators[i].second;
846 const LocalFunctionSpace& localFunctionSpace()
const
851 size_type
size()
const
859 CIVector _container_indices;
860 std::vector<unsigned char> _dof_flags;
861 std::vector<std::pair<ConstraintsIterator,ConstraintsIterator> > _constraints_iterators;
862 mutable CIMap _container_index_map;
863 ConstraintsVector _constraints;
865 const C& _gfs_constraints;
870 template<
typename LFS,
bool fast>
871 class LFSIndexCacheBase<LFS,EmptyTransformation,SimpleLFSCacheTag,fast>
876 typedef LFS LocalFunctionSpace;
877 typedef typename LFS::Traits::GridFunctionSpace GFS;
878 typedef typename GFS::Ordering Ordering;
880 typedef typename Ordering::Traits::ContainerIndex CI;
881 typedef typename Ordering::Traits::DOFIndex DI;
883 typedef CI ContainerIndex;
885 typedef std::size_t size_type;
887 typedef std::vector<CI> CIVector;
888 typedef std::unordered_map<DI,CI> CIMap;
890 struct ConstraintsEntry
891 :
public std::pair<const CI*,double>
893 typedef CI ContainerIndex;
894 typedef double Weight;
896 const ContainerIndex& containerIndex()
const
898 return *(this->first);
901 const Weight& weight()
const
907 typedef std::vector<ConstraintsEntry> ConstraintsVector;
908 typedef typename ConstraintsVector::const_iterator ConstraintsIterator;
910 explicit LFSIndexCacheBase(
const LFS& lfs)
916 LFSIndexCacheBase(
const LFS& lfs,
const C& c)
927 CI containerIndex(size_type i)
const
929 return CI(_lfs.dofIndex(i)[0]);
932 CI containerIndex(
const DI& i)
const
937 bool isConstrained(size_type i)
const
942 bool isDirichletConstraint(size_type i)
const
947 ConstraintsIterator constraintsBegin(size_type i)
const
949 return _constraints.begin();
952 ConstraintsIterator constraintsEnd(size_type i)
const
954 return _constraints.end();
957 const LocalFunctionSpace& localFunctionSpace()
const
962 size_type
size()
const
970 mutable CIMap _container_index_map;
971 const ConstraintsVector _constraints;
976 template<
typename LFS,
typename C = EmptyTransformation,
bool fast = false>
978 :
public LFSIndexCacheBase<LFS,C,typename LFS::Traits::GridFunctionSpace::Ordering::CacheTag,fast>
983 template<
typename CC>
984 LFSIndexCache(
const LFS& lfs,
const CC& c,
bool enable_constraints_caching = !std::is_same<C,EmptyTransformation>::value)
985 : LFSIndexCacheBase<LFS,C,typename LFS::Traits::GridFunctionSpace::Ordering::CacheTag,fast>(lfs,c,enable_constraints_caching)
989 explicit LFSIndexCache(
const LFS& lfs)
990 : LFSIndexCacheBase<LFS,C,typename LFS::Traits::GridFunctionSpace::Ordering::CacheTag,fast>(lfs)
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
void constraints(const GFS &gfs, CG &cg, const bool verbose=false)
construct constraints
Definition: constraints.hh:749
constexpr auto equals(T1 &&t1, T2 &&t2)
Equality comparison.
Definition: hybridutilities.hh:587
void applyToTree(Tree &&tree, Visitor &&visitor)
Apply visitor to TypeTree.
Definition: traversal.hh:239
Support for calculating hash values of objects.
This file implements iterator facade classes for writing stl conformant iterators.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
An stl-compliant random-access container which stores everything on the stack.
void post(T &&, TreePath) const
Method for postfix tree traversal.
Definition: visitor.hh:83
void leaf(T &&, TreePath) const
Method for leaf traversal.
Definition: visitor.hh:93
void pre(T &&, TreePath) const
Method for prefix tree traversal.
Definition: visitor.hh:60
static const std::size_t leafCount
The number of leaf nodes in the TypeTree.
Definition: utility.hh:136
static const std::size_t depth
The depth of the TypeTree.
Definition: utility.hh:130