DUNE PDELab (2.8)

lfsindexcache.hh
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_PDELAB_GRIDFUNCTIONSPACE_LFSINDEXCACHE_HH
4#define DUNE_PDELAB_GRIDFUNCTIONSPACE_LFSINDEXCACHE_HH
5
6#include <vector>
7#include <stack>
8#include <algorithm>
9#include <unordered_map>
10
11#include <dune/common/reservedvector.hh>
12#include <dune/common/exceptions.hh>
13#include <dune/common/hash.hh>
14#include <dune/common/iteratorfacades.hh>
15
16#include <dune/typetree/typetree.hh>
17
18#include <dune/pdelab/constraints/common/constraintstransformation.hh>
19#include <dune/pdelab/gridfunctionspace/tags.hh>
20
21namespace Dune {
22 namespace PDELab {
23
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
29 >
30 {
31
32 friend class RandomAccessIteratorFacade<
33 DOFIndexViewIterator,
34 const typename std::iterator_traits<Iterator>::value_type::View,
35 const typename std::iterator_traits<Iterator>::value_type::View
36 >;
37
38 typedef typename std::iterator_traits<Iterator>::value_type::View View;
39
40 public:
41
42 // Add support for returning non-references from iterator.
43 // We need a little bit of magic to make operator->() work for this iterator
44 // because we return a temporary object from dereference(), and the standard
45 // implementation of operator->() in the facade tries to take the address of
46 // that temporary, which the compiler will vehemently object to... ;-)
47 //
48 // So I borrowed the following neat little trick from Boost's iterator library:
49 // The proxy object stores a copy of the temporary View object, and operator()->
50 // returns the proxy object to the caller. As mandated by the standard, the compiler
51 // will then attempt to repeat the operator->() on the returned object and get the
52 // address of the copy stored in the (temporary) proxy object. That proxy object
53 // is guaranteed to live until the next sequence point, and that is precisely as
54 // long as we have to guarantee the validity of the pointer to our View object.
55 // Problem solved - and another example of how difficult it is to get this low-level
56 // stuff implemented on the same level as Boost...
57 struct proxy
58 {
59
60 explicit proxy(const View& v)
61 : _tmp(v)
62 {}
63
64 View* operator->()
65 {
66 return &_tmp;
67 }
68
69 View _tmp;
70 };
71
72 // The proxy object will stand in as a pointer
73 typedef proxy pointer;
74
75 DOFIndexViewIterator()
76 : _iterator()
77 , _tail_length(0)
78 {}
79
80 explicit DOFIndexViewIterator(Iterator it, std::size_t tail_length = 0)
81 : _iterator(it)
82 , _tail_length(tail_length)
83 {}
84
85 void cut_back()
86 {
87 ++_tail_length;
88 }
89
90 void restore_back()
91 {
92 --_tail_length;
93 }
94
95 const typename std::iterator_traits<Iterator>::reference raw_index() const
96 {
97 return *_iterator;
98 }
99
100 bool equals(const DOFIndexViewIterator& other) const
101 {
102 return _iterator == other._iterator;
103 }
104
105 void increment()
106 {
107 ++_iterator;
108 }
109
110 void decrement()
111 {
112 --_iterator;
113 }
114
115 void advance(int n)
116 {
117 _iterator += n;
118 }
119
120 std::ptrdiff_t distanceTo(DOFIndexViewIterator& other) const
121 {
122 return other._iterator - _iterator;
123 }
124
125 const View dereference() const
126 {
127 return _iterator->view(_iterator->treeIndex().size() - _tail_length);
128 }
129
130 pointer operator->() const
131 {
132 return pointer(dereference());
133 }
134
135 private:
136
137 Iterator _iterator;
138 std::size_t _tail_length;
139
140 };
141
142
143 template<typename Iterator>
144 struct extract_lfs_leaf_size_visitor
145 : public TypeTree::TreeVisitor
146 , public TypeTree::DynamicTraversal
147 {
148
149 template<typename LeafLFS, typename TreePath>
150 void leaf(const LeafLFS& leaf_lfs, TreePath tp)
151 {
152 (*it) = leaf_lfs.size();
153 ++it;
154 }
155
156 extract_lfs_leaf_size_visitor(Iterator leaf_size_container_iterator)
157 : it(leaf_size_container_iterator)
158 {}
159
160 Iterator it;
161
162 };
163
164 template<typename LFS, typename Iterator>
165 Iterator extract_lfs_leaf_sizes(const LFS& lfs, Iterator it)
166 {
167 extract_lfs_leaf_size_visitor<Iterator> visitor(it);
168 TypeTree::applyToTree(lfs,visitor);
169 return visitor.it;
170 }
171
172
173 template<typename DOFIterator,
174 typename ContainerIterator,
175 typename LeafSizeIterator,
176 std::size_t tree_depth,
177 bool fast>
178 struct map_dof_indices_to_container_indices
179 : public TypeTree::TreeVisitor
180 , public TypeTree::DynamicTraversal
181 {
182
183 template<typename Ordering, typename TreePath>
184 void leaf(const Ordering& ordering, TreePath tp)
185 {
186 std::size_t leaf_size = *(leaf_size_pos++);
187 if (fast)
188 dof_end += 1;
189 else
190 dof_end += leaf_size;
191 ordering.map_lfs_indices(dof_pos,dof_end,container_pos);
192 dof_pos = dof_end;
193 container_pos += leaf_size;
194 }
195
196 template<typename Ordering, typename TreePath>
197 void post(const Ordering& ordering, TreePath tp)
198 {
199 if (Ordering::consume_tree_index)
200 {
201 dof_pos.restore_back();
202 dof_end.restore_back();
203 }
204 ordering.map_lfs_indices(dof_stack.top(),dof_end,container_stack.top());
205 dof_stack.pop();
206 container_stack.pop();
207 }
208
209 template<typename Ordering, typename TreePath>
210 void pre(const Ordering& ordering, TreePath tp)
211 {
212 dof_stack.push(dof_pos);
213 container_stack.push(container_pos);
214 if (Ordering::consume_tree_index)
215 {
216 dof_pos.cut_back();
217 dof_end.cut_back();
218 }
219 }
220
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)
229 {}
230
231
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;
238
239 };
240
241
242
243 template<typename LFS, typename C, typename CacheTag, bool fast>
244 class LFSIndexCacheBase
245 {
246
247 enum DOFFlags
248 {
249 DOF_NONCONSTRAINED = 0,
250 DOF_CONSTRAINED = 1<<0,
251 DOF_DIRICHLET = 1<<1
252 };
253
254 public:
255
256 typedef LFS LocalFunctionSpace;
257
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;
263 typedef DOFIndex DI;
264 typedef std::size_t size_type;
265
266 typedef std::vector<CI> CIVector;
267 typedef std::unordered_map<DI,CI> CIMap;
268
269 typedef std::unordered_map<const CI*,std::pair<size_type,bool> > InverseMap;
270
271 struct ConstraintsEntry
272 : public std::pair<const CI*,typename C::mapped_type::mapped_type>
273 {
274 typedef CI ContainerIndex;
275 typedef typename C::mapped_type::mapped_type Weight;
276
277 const ContainerIndex& containerIndex() const
278 {
279 return *(this->first);
280 }
281
282 const Weight& weight() const
283 {
284 return this->second;
285 }
286 };
287
288 //typedef std::pair<CI,typename C::mapped_type::mapped_type> ConstraintsEntry;
289
290 typedef std::vector<ConstraintsEntry> ConstraintsVector;
291 typedef typename ConstraintsVector::const_iterator ConstraintsIterator;
292
293 LFSIndexCacheBase(const LFS& lfs, const C& constraints, bool enable_constraints_caching)
294 : _lfs(lfs)
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)
300 , _gfs_constraints(constraints)
301 {
302 }
303
304 void update()
305 {
306 if(fast) {
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));
310 }
311 else {
312
313 // clear out existing state
314 _container_index_map.clear();
315 for (typename CIVector::iterator it = _container_indices.begin(); it != _container_indices.end(); ++it)
316 it->clear();
317
318 _inverse_map.clear();
319 _inverse_cache_built = false;
320
321 // extract size for all leaf spaces (into a flat list)
322 typedef ReservedVector<size_type,TypeTree::TreeInfo<LFS>::leafCount> LeafSizeVector;
323 LeafSizeVector leaf_sizes;
324 leaf_sizes.resize(TypeTree::TreeInfo<LFS>::leafCount);
325 extract_lfs_leaf_sizes(_lfs,leaf_sizes.begin());
326
327 // perform the actual mapping
328 map_dof_indices_to_container_indices<
329 typename LFS::Traits::DOFIndexContainer::const_iterator,
330 typename CIVector::iterator,
331 typename LeafSizeVector::const_iterator,
332 TypeTree::TreeInfo<Ordering>::depth,
333 fast
334 > index_mapper(_lfs._dof_indices->begin(),_container_indices.begin(),leaf_sizes.begin(),_lfs.subSpaceDepth());
335 TypeTree::applyToTree(_lfs.gridFunctionSpace().ordering(),index_mapper);
336
337 if (_enable_constraints_caching)
338 {
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)
344 {
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())
348 {
349 _dof_flags[i] = DOF_NONCONSTRAINED;
350 continue;
351 }
352
353 if (cit->second.size() == 0)
354 {
355 _dof_flags[i] = DOF_CONSTRAINED | DOF_DIRICHLET;
356 _constraints_iterators[i] = make_pair(_constraints.end(),_constraints.end());
357 }
358 else
359 {
360 _dof_flags[i] = DOF_CONSTRAINED;
361 constraint_entry_count += cit->second.size();
362 non_dirichlet_constrained_dofs.push_back(make_pair(i,cit));
363 }
364 }
365
366 if (constraint_entry_count > 0)
367 {
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();
372 ++it)
373 {
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)
376 {
377 eit->first = &(cit->first);
378 eit->second = cit->second;
379 }
380 _constraints_iterators[it->first].second = eit;
381 }
382 }
383 }
384 }
385 }
386
387 const DI& dofIndex(size_type i) const
388 {
389 return _lfs.dofIndex(i);
390 }
391
392 const CI& containerIndex(size_type i) const
393 {
394 return _container_indices[i];
395 }
396
397 const CI& containerIndex(const DI& i) const
398 {
399 // look up DOFIndex i
400 std::pair<typename CIMap::iterator,bool> r = _container_index_map.insert(std::make_pair(std::ref(i),CI()));
401
402 // i did not exist in the cache, map it into the newly inserted container index
403 if (r.second)
404 _lfs.gridFunctionSpace().ordering().mapIndex(i.view(),r.first->second);
405
406 // return cached container index
407 return r.first->second;
408 }
409
410 bool isConstrained(size_type i) const
411 {
412 return _dof_flags[i] & DOF_CONSTRAINED;
413 }
414
415 bool isDirichletConstraint(size_type i) const
416 {
417 return _dof_flags[i] & DOF_DIRICHLET;
418 }
419
420 ConstraintsIterator constraintsBegin(size_type i) const
421 {
422 assert(isConstrained(i));
423 return _constraints_iterators[i].first;
424 }
425
426 ConstraintsIterator constraintsEnd(size_type i) const
427 {
428 assert(isConstrained(i));
429 return _constraints_iterators[i].second;
430 }
431
432 const LocalFunctionSpace& localFunctionSpace() const
433 {
434 return _lfs;
435 }
436
437 size_type size() const
438 {
439 return _lfs.size();
440 }
441
442 std::pair<size_type,bool> localIndex(const ContainerIndex& ci) const
443 {
444 if (!_inverse_cache_built)
445 build_inverse_cache();
446 return _inverse_map[ci];
447 }
448
449 size_type offset(size_type i) const
450 {
451 if (!_inverse_cache_built)
452 build_inverse_cache();
453 return _offsets[i];
454 }
455
456 size_type extendedOffset(size_type i) const
457 {
458 if (!_inverse_cache_built)
459 build_inverse_cache();
460 return _extended_offsets[i];
461 }
462
463 bool constraintsCachingEnabled() const
464 {
465 return _enable_constraints_caching;
466 }
467
468 private:
469
470 struct sort_container_indices
471 {
472 template<typename T>
473 bool operator()(const T* a, const T* b) const
474 {
475 return std::lexicographical_compare(reversed_iterator(a->end()),reversed_iterator(a->begin()),
476 reversed_iterator(b->end()),reversed_iterator(b->begin())
477 );
478 }
479 };
480
481
482 void build_inverse_cache()
483 {
484 size_type i = 0;
485 size_type child = 0;
486 _offsets[0] = 0;
487 for (typename CIVector::const_iterator it = _container_indices.begin(),
488 endit = _container_indices.end();
489 it != endit;
490 ++it, ++i
491 )
492 {
493 _inverse_map.insert(std::make_pair(&(*it),std::make_pair(i,false)));
494 if (it->back() != child)
495 {
496 _offsets[child+1] = i;
497 ++child;
498 }
499 }
500
501 std::vector<const ContainerIndex*> extended_cis;
502 extended_cis.reserve(_constraints.size());
503
504 for (typename ConstraintsVector::const_iterator it = _constraints.begin(),
505 endit = _constraints.end();
506 it != endit;
507 ++it
508 )
509 {
510 if (_inverse_map.count(it->first) == 0)
511 extended_cis.push_back(it->first);
512 }
513
514 std::sort(extended_cis.begin(),extended_cis.end(),sort_container_indices());
515
516 typename std::vector<const ContainerIndex*>::const_iterator endit = std::unique(extended_cis.begin(),extended_cis.end());
517
518 i = 0;
519 child = 0;
520 for (typename std::vector<const ContainerIndex*>::const_iterator it = extended_cis.begin(); it != endit; ++it, ++i)
521 {
522 _inverse_map.insert(std::make_pair(&(*it),std::make_pair(i,true)));
523 if (it->back() != child)
524 {
525 _extended_offsets[child+1] = i;
526 ++child;
527 }
528 }
529
530 _inverse_cache_built = true;
531
532 }
533
534 const LFS& _lfs;
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;
545
546 const C& _gfs_constraints;
547
548 };
549
550
551 template<typename LFS, typename CacheTag, bool fast>
552 class LFSIndexCacheBase<LFS,EmptyTransformation,CacheTag,fast>
553 {
554
555 public:
556
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;
563 typedef DOFIndex DI;
564 typedef std::size_t size_type;
565
566 typedef std::vector<CI> CIVector;
567 typedef std::unordered_map<DI,CI> CIMap;
568
569 struct ConstraintsEntry
570 : public std::pair<const CI*,double>
571 {
572 typedef CI ContainerIndex;
573 typedef double Weight;
574
575 const ContainerIndex& containerIndex() const
576 {
577 return *(this->first);
578 }
579
580 const Weight& weight() const
581 {
582 return this->second;
583 }
584 };
585
586 typedef std::vector<ConstraintsEntry> ConstraintsVector;
587 typedef typename ConstraintsVector::const_iterator ConstraintsIterator;
588
589 explicit LFSIndexCacheBase(const LFS& lfs)
590 : _lfs(lfs)
591 , _container_indices(lfs.maxSize())
592 {
593 }
594
595 template<typename C>
596 LFSIndexCacheBase(const LFS& lfs, const C& c, bool enable_constraints_caching)
597 : _lfs(lfs)
598 , _container_indices(lfs.maxSize())
599 {
600 }
601
602
603 void update()
604 {
605 if(fast) {
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));
609 }
610 else {
611
612 // clear out existing state
613 _container_index_map.clear();
614 for (typename CIVector::iterator it = _container_indices.begin(); it != _container_indices.end(); ++it)
615 it->clear();
616
617 // extract size for all leaf spaces (into a flat list)
618 typedef ReservedVector<size_type,TypeTree::TreeInfo<LFS>::leafCount> LeafSizeVector;
619 LeafSizeVector leaf_sizes;
620 leaf_sizes.resize(TypeTree::TreeInfo<LFS>::leafCount);
621 extract_lfs_leaf_sizes(_lfs,leaf_sizes.begin());
622
623 // perform the actual mapping
624 map_dof_indices_to_container_indices<
625 typename LFS::Traits::DOFIndexContainer::const_iterator,
626 typename CIVector::iterator,
627 typename LeafSizeVector::const_iterator,
628 TypeTree::TreeInfo<Ordering>::depth,
629 fast
630 > index_mapper(_lfs._dof_indices->begin(),_container_indices.begin(),leaf_sizes.begin(),_lfs.subSpaceDepth());
631 TypeTree::applyToTree(_lfs.gridFunctionSpace().ordering(),index_mapper);
632 }
633 }
634
635 const DI& dofIndex(size_type i) const
636 {
637 return _lfs.dofIndex(i);
638 }
639
640 const CI& containerIndex(size_type i) const
641 {
642 return _container_indices[i];
643 }
644
645 const CI& containerIndex(const DI& i) const
646 {
647 // look up DOFIndex i
648 std::pair<typename CIMap::iterator,bool> r = _container_index_map.insert(std::make_pair(std::ref(i),CI()));
649
650 // i did not exist in the cache, map it into the newly inserted container index
651 if (r.second)
652 _lfs.gridFunctionSpace().ordering().mapIndex(i.view(),r.first->second);
653
654 // return cached container index
655 return r.first->second;
656 }
657
658 bool isConstrained(size_type i) const
659 {
660 return false;
661 }
662
663 bool isDirichletConstraint(size_type i) const
664 {
665 return false;
666 }
667
668 ConstraintsIterator constraintsBegin(size_type i) const
669 {
670 return _constraints.begin();
671 }
672
673 ConstraintsIterator constraintsEnd(size_type i) const
674 {
675 return _constraints.end();
676 }
677
678 const LocalFunctionSpace& localFunctionSpace() const
679 {
680 return _lfs;
681 }
682
683 size_type size() const
684 {
685 return _lfs.size();
686 }
687
688 bool constraintsCachingEnabled() const
689 {
690 return false;
691 }
692
693 private:
694
695 const LFS& _lfs;
696 CIVector _container_indices;
697 mutable CIMap _container_index_map;
698 const ConstraintsVector _constraints;
699
700 };
701
702
703
704 template<typename LFS, typename C, bool fast>
705 class LFSIndexCacheBase<LFS,C,SimpleLFSCacheTag,fast>
706 {
707
708 enum DOFFlags
709 {
710 DOF_NONCONSTRAINED = 0,
711 DOF_CONSTRAINED = 1<<0,
712 DOF_DIRICHLET = 1<<1
713 };
714
715 public:
716
717 typedef LFS LocalFunctionSpace;
718
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;
724
725 typedef std::vector<CI> CIVector;
726 typedef std::unordered_map<DI,CI> CIMap;
727
728 struct ConstraintsEntry
729 : public std::pair<CI,typename C::mapped_type::mapped_type>
730 {
731 typedef CI ContainerIndex;
732 typedef typename C::mapped_type::mapped_type Weight;
733
734 const ContainerIndex& containerIndex() const
735 {
736 return this->first;
737 }
738
739 const Weight& weight() const
740 {
741 return this->second;
742 }
743 };
744
745 typedef std::vector<ConstraintsEntry> ConstraintsVector;
746 typedef typename ConstraintsVector::const_iterator ConstraintsIterator;
747
748 LFSIndexCacheBase(const LFS& lfs, const C& constraints)
749 : _lfs(lfs)
750 , _dof_flags(lfs.maxSize())
751 , _constraints_iterators(lfs.maxSize())
752 , _gfs_constraints(constraints)
753 {
754 }
755
756 void update()
757 {
758 if(fast) {
759 DUNE_THROW(Dune::Exception, "This function shouldn't be called in fast mode");
760 }
761 else {
762
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)
767 {
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())
771 {
772 _dof_flags[i] = DOF_NONCONSTRAINED;
773 continue;
774 }
775
776 if (cit->second.size() == 0)
777 {
778 _dof_flags[i] = DOF_CONSTRAINED | DOF_DIRICHLET;
779 _constraints_iterators[i] = make_pair(_constraints.end(),_constraints.end());
780 }
781 else
782 {
783 _dof_flags[i] = DOF_CONSTRAINED;
784 constraint_entry_count += cit->second.size();
785 non_dirichlet_constrained_dofs.push_back(make_pair(i,cit));
786 }
787 }
788
789 if (constraint_entry_count > 0)
790 {
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();
795 ++it)
796 {
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)
799 {
800 eit->first = cit->first;
801 eit->second = cit->second;
802 }
803 _constraints_iterators[it->first].second = eit;
804 }
805 }
806 }
807 }
808
809 const DI& dofIndex(size_type i) const
810 {
811 return _lfs.dofIndex(i);
812 }
813
814 CI containerIndex(size_type i) const
815 {
816 return CI(_lfs.dofIndex(i)[0]);
817 }
818
819 const CI& containerIndex(const DI& i) const
820 {
821 return CI(i[0]);
822 }
823
824 bool isConstrained(size_type i) const
825 {
826 return _dof_flags[i] & DOF_CONSTRAINED;
827 }
828
829 bool isDirichletConstraint(size_type i) const
830 {
831 return _dof_flags[i] & DOF_DIRICHLET;
832 }
833
834 ConstraintsIterator constraintsBegin(size_type i) const
835 {
836 assert(isConstrained(i));
837 return _constraints_iterators[i].first;
838 }
839
840 ConstraintsIterator constraintsEnd(size_type i) const
841 {
842 assert(isConstrained(i));
843 return _constraints_iterators[i].second;
844 }
845
846 const LocalFunctionSpace& localFunctionSpace() const
847 {
848 return _lfs;
849 }
850
851 size_type size() const
852 {
853 return _lfs.size();
854 }
855
856 private:
857
858 const LFS& _lfs;
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;
864
865 const C& _gfs_constraints;
866
867 };
868
869
870 template<typename LFS, bool fast>
871 class LFSIndexCacheBase<LFS,EmptyTransformation,SimpleLFSCacheTag,fast>
872 {
873
874 public:
875
876 typedef LFS LocalFunctionSpace;
877 typedef typename LFS::Traits::GridFunctionSpace GFS;
878 typedef typename GFS::Ordering Ordering;
879 private:
880 typedef typename Ordering::Traits::ContainerIndex CI;
881 typedef typename Ordering::Traits::DOFIndex DI;
882 public:
883 typedef CI ContainerIndex;
884 typedef DI DOFIndex;
885 typedef std::size_t size_type;
886
887 typedef std::vector<CI> CIVector;
888 typedef std::unordered_map<DI,CI> CIMap;
889
890 struct ConstraintsEntry
891 : public std::pair<const CI*,double>
892 {
893 typedef CI ContainerIndex;
894 typedef double Weight;
895
896 const ContainerIndex& containerIndex() const
897 {
898 return *(this->first);
899 }
900
901 const Weight& weight() const
902 {
903 return this->second;
904 }
905 };
906
907 typedef std::vector<ConstraintsEntry> ConstraintsVector;
908 typedef typename ConstraintsVector::const_iterator ConstraintsIterator;
909
910 explicit LFSIndexCacheBase(const LFS& lfs)
911 : _lfs(lfs)
912 {
913 }
914
915 template<typename C>
916 LFSIndexCacheBase(const LFS& lfs, const C& c)
917 : _lfs(lfs)
918 {
919 }
920
921
922 void update()
923 {
924 // there's nothing to do here...
925 }
926
927 CI containerIndex(size_type i) const
928 {
929 return CI(_lfs.dofIndex(i)[0]);
930 }
931
932 CI containerIndex(const DI& i) const
933 {
934 return CI(i[0]);
935 }
936
937 bool isConstrained(size_type i) const
938 {
939 return false;
940 }
941
942 bool isDirichletConstraint(size_type i) const
943 {
944 return false;
945 }
946
947 ConstraintsIterator constraintsBegin(size_type i) const
948 {
949 return _constraints.begin();
950 }
951
952 ConstraintsIterator constraintsEnd(size_type i) const
953 {
954 return _constraints.end();
955 }
956
957 const LocalFunctionSpace& localFunctionSpace() const
958 {
959 return _lfs;
960 }
961
962 size_type size() const
963 {
964 return _lfs.size();
965 }
966
967 private:
968
969 const LFS& _lfs;
970 mutable CIMap _container_index_map;
971 const ConstraintsVector _constraints;
972
973 };
974
975
976 template<typename LFS, typename C = EmptyTransformation, bool fast = false>
977 class LFSIndexCache
978 : public LFSIndexCacheBase<LFS,C,typename LFS::Traits::GridFunctionSpace::Ordering::CacheTag,fast>
979 {
980
981 public:
982
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)
986 {
987 }
988
989 explicit LFSIndexCache(const LFS& lfs)
990 : LFSIndexCacheBase<LFS,C,typename LFS::Traits::GridFunctionSpace::Ordering::CacheTag,fast>(lfs)
991 {
992 }
993
994 };
995
996
997 } // namespace PDELab
998} // namespace Dune
999
1000#endif // DUNE_PDELAB_GRIDFUNCTIONSPACE_LFSINDEXCACHE_HH
Dune::Exception
Base class for Dune-Exceptions.
Definition: exceptions.hh:94
exceptions.hh
A few common exception classes.
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::PDELab::constraints
void constraints(const GFS &gfs, CG &cg, const bool verbose=false)
construct constraints
Definition: constraints.hh:749
Dune::Hybrid::equals
constexpr auto equals(T1 &&t1, T2 &&t2)
Equality comparison.
Definition: hybridutilities.hh:400
Dune::TypeTree::applyToTree
void applyToTree(Tree &&tree, Visitor &&visitor)
Apply visitor to TypeTree.
Definition: traversal.hh:237
Dune::TypeTree::child
ImplementationDefined child(Node &&node, Indices... indices)
Extracts the child of a node given by a sequence of compile-time and run-time indices.
Definition: childextraction.hh:126
hash.hh
Support for calculating hash values of objects.
iteratorfacades.hh
This file implements iterator facade classes for writing stl conformant iterators.
Dune
Dune namespace.
Definition: alignedallocator.hh:11
std
STL namespace.
reservedvector.hh
An stl-compliant random-access container which stores everything on the stack.
Dune::TypeTree::DefaultVisitor::pre
void pre(T &&t, TreePath treePath) const
Method for prefix tree traversal.
Definition: visitor.hh:58
Dune::TypeTree::DefaultVisitor::post
void post(T &&t, TreePath treePath) const
Method for postfix tree traversal.
Definition: visitor.hh:81
Dune::TypeTree::DefaultVisitor::leaf
void leaf(T &&t, TreePath treePath) const
Method for leaf traversal.
Definition: visitor.hh:91
Dune::TypeTree::TreeInfo::leafCount
static const std::size_t leafCount
The number of leaf nodes in the TypeTree.
Definition: utility.hh:183
Dune::TypeTree::TreeInfo::depth
static const std::size_t depth
The depth of the TypeTree.
Definition: utility.hh:177
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)