dune-fem 2.8.0
Loading...
Searching...
No Matches
genericadaptivedofmapper.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_GENERICPADAPTIVEDOFMAPPER_HH
2#define DUNE_FEM_GENERICPADAPTIVEDOFMAPPER_HH
3
4#include <dune/common/exceptions.hh>
5
6#include <dune/geometry/type.hh>
7
8#include <dune/grid/utility/persistentcontainer.hh>
9
12#include <dune/fem/misc/metaprogramming.hh>
19
20namespace Dune
21{
22
23 namespace Fem
24 {
25
26 // GenericAdaptiveDofMapper
27 // ------------------------
28
29 template< class TraitsImp >
31 : public AdaptiveDofMapper< TraitsImp >
32 {
35
36 public:
37 typedef TraitsImp Traits;
38 typedef std::size_t SizeType;
39
40 // true if all dofs are associated with entities of codim 0
41 static const bool discontinuousMapper = Traits :: discontinuousMapper ;
42
43 protected:
44 using BaseType::asImp;
45
46 public:
48 typedef typename Traits::GridPartType GridPartType;
49
51 typedef typename Traits::ElementType ElementType;
52
54 typedef typename GridPartType::GridType GridType;
55
57 typedef typename GridPartType::IndexSetType IndexSetType;
58
60 typedef typename Traits :: GlobalKeyType GlobalKeyType ;
61
63 typedef typename GridType::ctype FieldType;
64
66 static const int dimension = GridType::dimension;
67
69 static const int highestDimension = ( discontinuousMapper ) ? 0 : dimension ;
70
72 static const int polynomialOrder = Traits::polynomialOrder;
73
75 typedef typename Traits :: CompiledLocalKeyVectorType CompiledLocalKeyVectorType;
76
78 typedef typename CompiledLocalKeyVectorType :: value_type :: value_type CompiledLocalKeyType;
79
82
83 enum { minOrder = 1 };
85 enum { numOrders = maxOrder - minOrder + 1 };
86
88 {
89 typedef std::vector< int > DofVectorType;
90 std::vector< DofVectorType > dofs_;
91
92 GeometryType type_;
94
97 type_()
98 {
99 // set used to zero
100 for( int i=0; i<numOrders; ++i )
101 used_[ i ] = 0;
102 }
103
104 void assign( const EntityDofStorage& other )
105 {
106 assert( dofs_.size() == other.dofs_.size() );
107 type_ = other.type_;
108
109 // set used to zero
110 for( int k=0; k<numOrders; ++k )
111 {
112 used_[ k ] = other.used_[ k ];
113 DofVectorType& dofs = dofs_[ k ];
114 const DofVectorType& otherDofs = other.dofs_[ k ];
115 const int dofSize = otherDofs.size();
116 dofs.resize( dofSize );
117 for( int d = 0; d<dofSize; ++d )
118 dofs[ d ] = otherDofs[ d ];
119 }
120 }
121
123 : dofs_( numOrders )
124 {
125 assign( other );
126 }
127
129 {
130 assign( other );
131 return *this;
132 }
133
134 bool exists( const int codim, const int polOrd ) const
135 {
136 const int entry = determineVectorEntry( codim, polOrd );
137 return dofs_[ entry ].size() > 0 ;
138 }
139
141 bool use( const int codim, const int polOrd )
142 {
143 const int entry = determineVectorEntry( codim, polOrd ) ;
144 ++used_[ entry ];
145 return ( used_[ entry ] == 1 );
146 }
147
148 void insert( const GeometryType type,
149 const int codim,
150 const int polOrd,
151 const int numDofs, const int startDof )
152 {
153 use( codim, polOrd );
154 assert( ! exists ( codim, polOrd ) );
155 {
156 type_ = type ;
157 DofVectorType& dofs = dofs_[ determineVectorEntry( codim, polOrd ) ];
158
159 dofs.resize( numDofs );
160 for(int i=0, dof=startDof ; i<numDofs; ++i, ++dof )
161 dofs[ i ] = dof;
162 }
163 }
164
165 int determineVectorEntry( const int codim, const int polOrd ) const
166 {
167 assert( codim >= 0 );
168 assert( codim <= highestDimension );
169 // also for codim == 0 we have more then
170 // one storage because of different number of
171 // dofs per polynmomial degree
172 return (codim < dimension) ? (polOrd-minOrder) : 0;
173 }
174
175 const GeometryType& type () const { return type_ ; }
176
177 void remove( const int codim, const int polOrd )
178 {
179 const int entry = determineVectorEntry( codim, polOrd );
180 if( used_[ entry ] > 0 )
181 --used_[ entry ] ;
182 }
183
184 void reset()
185 {
186 for( int k=0; k<numOrders; ++k )
187 used_[ k ] = 0;
188 }
189
190 int dof ( const int codim, const int polOrd, const size_t dofNumber ) const
191 {
192 const unsigned int entry = determineVectorEntry( codim, polOrd );
193 assert( entry < dofs_.size() );
194 assert( type_ != GeometryType() );
195 assert( dofNumber < dofs_[ entry ].size() );
196 return dofs_[ entry ][ dofNumber ];
197 }
198
199 int entityDof ( int dofNumber ) const
200 {
201 for( int k = 0; k<numOrders; ++k )
202 {
203 const int dofSize = dofs_[ k ].size();
204 if( dofNumber < dofSize )
205 return dofs_[ k ][ dofNumber ];
206 else
207 dofNumber -= dofSize;
208 }
209 // we should not get here
210 assert( false );
211 abort();
212 return -1;
213 }
214
215 int entityDofs () const
216 {
217 int dofSize = 0;
218 for( int k = 0; k<numOrders; ++k )
219 {
220 dofSize += dofs_[ k ].size();
221 }
222 return dofSize;
223 }
224
225 template <class VectorType>
226 void detectUnusedDofs( VectorType& isHole,
227 const int actSize )
228 {
229 for( int k=0; k<numOrders; ++k )
230 {
231 DofVectorType& dofs = dofs_[ k ];
232 const int dofSize = dofs.size();
233
234 if( dofSize > 0 )
235 {
236 if( used_[ k ] )
237 {
238 for( int d = 0; d<dofSize; ++d )
239 {
240 const int dof = dofs[ d ] ;
241 if( dof < actSize )
242 {
243 assert( dof < (int)isHole.size() );
244 isHole[ dof ] = false ;
245 }
246 }
247 }
248 else
249 {
250 dofs.resize( 0 );
251 }
252 }
253 }
254 }
255
256 void printDofs() const
257 {
258 for( int k = 0; k<numOrders; ++k )
259 {
260 const DofVectorType& dofs = dofs_[ k ];
261 const int dofSize = dofs.size();
262 for( int d = 0; d<dofSize; ++d )
263 std::cout << dofs[ d ] << " dofs " << std::endl;
264 }
265 }
266
267 template <class VectorType>
268 bool removeHoles( VectorType& oldIdx, VectorType& newIdx,
269 VectorType& holesVec, int& currentHole,
270 const int usedSize, int& holes )
271 {
272 bool haveToCopy = false ;
273 for( int k=0; k<numOrders; ++k )
274 {
275 DofVectorType& dofs = dofs_[ k ];
276 const int dofSize = dofs.size();
277 for( int dof = 0; dof<dofSize; ++dof )
278 {
279 assert( used_[ k ] );
280 // get global DoF number
281 int& currDof = dofs[ dof ] ;
282
283 // if dof >= usedSize it has to be migrated to a hole
284 if( currDof >= usedSize )
285 {
286 // serach next hole that is smaler than actual size
287 --currentHole;
288
289 // if currentHole < 0 then error, because we have index larger then
290 // actual size
291 assert(currentHole >= 0);
292 assert( holesVec[currentHole] < usedSize );
293
294 // remember old and new index
295 oldIdx[ holes ] = currDof;
296 currDof = holesVec[ currentHole ];
297 newIdx[ holes ] = currDof ;
298
299 // increase number of holes
300 ++holes;
301
302 haveToCopy = true;
303 }
304 }
305 }
306 return haveToCopy;
307 }
308 };
309
311
313 {
314 private:
315 // PolynomialOrderStorage() : k_( minOrder ), active_( -std::abs(k_) ) {}
316 signed char k_; // stores current polynomial order
317 signed char active_ ; // stores active/non-active and suggested pol order
318
319 public:
320 // PolynomialOrderStorage() : k_( minOrder ), active_( -std::abs(k_) ) {}
321 PolynomialOrderStorage( const int k ) : k_( k ), active_( -std::abs(k_) ) {}
322 int order () const { return k_; }
323 void suggest ( const int k )
324 {
325 if( active() )
326 active_ = std::abs( k );
327 else
328 active_ = -std::abs( k );
329 }
330 void set ( const int k ) { k_ = k; active_ = std::abs( k ) ; }
331 void activate() { active_ = std::abs( active_ ); }
332 bool active () const { return active_ > 0; }
333 bool deactivate ( int& k )
334 {
335 k = k_;
336 if( active() )
337 {
338 active_ = -active_;
339 return true ;
340 }
341 return false;
342 }
343
344 int suggested () const { return std::abs( active_ ); }
345 void update() { set( suggested() ); }
346 };
347
349
350 typedef PersistentContainer< GridType, EntityDofStorageType > DofContainerType ;
351
352 typedef PersistentContainer< GridType, PolynomialOrderStorageType > PolyOrderContainerType ;
353 protected:
354 template <int codim, bool dg>
355 struct NumDofs
356 {
357 static int numDofs( const ElementType& entity,
358 const CompiledLocalKeyType& clk,
359 const int subEntity )
360 {
361 return clk.numDofs( codim, subEntity );
362 }
363 };
364
365 template <int codim>
366 struct NumDofs<codim, true>
367 {
368 static int numDofs( const ElementType& entity,
369 const CompiledLocalKeyType& clk,
370 const int subEntity )
371 {
372 if( codim == 0 )
373 return clk.size();
374 else
375 return 0;
376 }
377 };
378
379 template < int codim >
381 {
382 static void insertDofs( const ElementType& entity,
383 const CompiledLocalKeyType& clk,
384 const int polOrd,
385 const int subEntity,
386 unsigned int& globalSize,
387 unsigned int& notAlreadyCounted,
388 EntityDofStorage& entityDofs )
389 {
390 const int numDofs = NumDofs<codim, discontinuousMapper>::numDofs( entity, clk, subEntity );
391
392 // only if dofs exists on this entity do something
393 if( numDofs > 0 )
394 {
395 bool notCountedYet = false ;
396 if( ! entityDofs.exists( codim, polOrd ) )
397 {
398 entityDofs.insert( entity.type(), codim, polOrd, numDofs, globalSize );
399 globalSize += numDofs;
400 notCountedYet = true ;
401 }
402 else
403 {
404 // if refcount is only 1 then true is returned
405 notCountedYet = entityDofs.use( codim, polOrd );
406 }
407
408 // if not counted yet, count !
409 if( notCountedYet )
410 {
411 notAlreadyCounted += numDofs;
412 }
413 }
414 }
415
416 static void apply( const ElementType& entity,
417 const CompiledLocalKeyType& clk,
418 const int polOrd,
419 unsigned int& globalSize,
420 unsigned int& notAlreadyCounted,
421 std::vector< DofContainerType* > dofContainers )
422 {
423 DofContainerType &dofContainer = *dofContainers[ codim ];
424 if( codim == 0 )
425 {
426 insertDofs( entity, clk, polOrd, 0, globalSize,
427 notAlreadyCounted, dofContainer[ entity ] );
428 }
429 else
430 {
431 const int count = entity.subEntities( codim );
432 for(int i=0; i<count; ++i )
433 {
434 insertDofs( entity, clk, polOrd, i, globalSize,
435 notAlreadyCounted, dofContainer( entity, i ) );
436 }
437 }
438 }
439 };
440
441 template < int codim >
443 {
444 static void apply( const ElementType& entity,
445 const int polOrd,
446 std::vector< DofContainerType* > dofContainers )
447 {
448 DofContainerType &dofContainer = *dofContainers[ codim ];
449 const int count = entity.subEntities( codim );
450 for(int i=0; i<count; ++i )
451 {
452 EntityDofStorage& entityDofs = dofContainer( entity, i );
453 entityDofs.remove( codim, polOrd );
454 }
455 }
456 };
457
458 public:
461 const int order,
462 CompiledLocalKeyVectorType &compiledLocalKeyVector )
463 : gridPart_( gridPart ),
464 dm_( DofManagerType :: instance(gridPart.grid()) ),
465 compiledLocalKeys_( compiledLocalKeyVector ),
466 order_( ( order >= minOrder && order <= maxOrder ) ? order : maxOrder ),
467 entityPolynomOrder_( gridPart.grid(), 0, PolynomialOrderStorageType( order_ ) ),
468 dofContainer_( dimension+1, nullptr ),
469 numberOfHoles_( 0 ),
470 oldIndex_(),
471 newIndex_(),
472 size_(0),
473 maxNumDofs_( 0 ),
474 sequence_( dm_.sequence() )
475 {
476 for( int codim = 0; codim <= highestDimension; ++codim )
477 dofContainer_[ codim ] = new DofContainerType( gridPart.grid(), codim );
478
479 for( size_t i=0; i<compiledLocalKeys_.size(); ++i )
480 {
481 maxNumDofs_ = std :: max( maxNumDofs_, compiledLocalKeys_[ i ].maxSize() );
482 }
483
484 resize();
485 // register to DofManager for adaptation
486 dm_.addIndexSet( asImp() );
487 }
488
491 const int order,
492 CompiledLocalKeyVectorType &compiledLocalKeyVector )
493 : gridPart_( other.gridPart_ ),
494 dm_( other.dm_ ),
495 compiledLocalKeys_( compiledLocalKeyVector ),
496 order_( ( order >= minOrder && order <= maxOrder ) ? order : maxOrder ),
497 entityPolynomOrder_( other.entityPolynomOrder_ ),
498 dofContainer_( dimension+1, nullptr ),
499 numberOfHoles_( other.numberOfHoles_ ),
500 oldIndex_( other.oldIndex_ ),
501 newIndex_( other.newIndex_ ),
502 size_( other.size_ ),
503 maxNumDofs_( other.maxNumDofs_ ),
504 sequence_( other.sequence_ )
505 {
506 for( int codim = 0; codim <= highestDimension; ++codim )
507 dofContainer_[ codim ] = new DofContainerType( *(other.dofContainer_[ codim ]) );
508
509 dm_.addIndexSet( asImp() );
510 }
511
512 int polynomOrder( const ElementType& entity ) const
513 {
514 return entityPolynomOrder_[ entity ].order();
515 }
516
517 int suggestedOrder( const ElementType& entity ) const
518 {
519 return entityPolynomOrder_[ entity ].suggestedOrder();
520 }
521
522 void suggestPolynomOrder( const ElementType& entity, const int polOrd )
523 {
524 // minOrder is static but order_ is dynamically set
525 if( polOrd < minOrder || polOrd > order_ )
526 return ;
527
528 entityPolynomOrder_[ entity ].suggest( polOrd );
529 }
530
531 DofContainerType &dofContainer ( const std::size_t codim ) const
532 {
533 assert( codim < dofContainer_.size() );
534 assert( dofContainer_[ codim ] );
535 return *(dofContainer_[ codim ]);
536 }
537
539 compiledLocalKey( const int polOrd, const GeometryType type ) const
540 {
541 // add polOrd here
542 return compiledLocalKeys_[ polOrd ][ type ];
543 }
544
547 {
548 dm_.removeIndexSet( asImp() );
549 }
550
552 int size () const
553 {
554 return size_;
555 }
556
557 template< class Functor >
558 void mapEach ( const ElementType &element, Functor f ) const
559 {
560 const int n = numDofs( element );
561 for( int i = 0; i < n; ++i )
562 f( i, mapToGlobal( element, i ) );
563 }
564
566 int mapToGlobal ( const ElementType &entity, const int localDof ) const
567 {
569 {
570 // get polynomial order
571 const int polOrd = polynomOrder( entity );
572
573 // return global dof
574 return dofContainer( 0 )[ entity ].dof( 0, polOrd, localDof );
575 }
576 else
577 {
578 // the continuous case
579 const int polOrd = polynomOrder( entity );
580
581 const CompiledLocalKeyType& compLocalKey = compiledLocalKey( polOrd, entity.type() );
582 // get dof info for entity and local dof
583 const Fem::LocalKey &dofInfo = compLocalKey.localKey( localDof );
584
585 const unsigned int codim = dofInfo.codim();
586 const unsigned int subEntity = dofInfo.subEntity();
587
588 unsigned int index = dofInfo.index() ;
589
590 // account for possible twists in the grid (only 2d)
591 if( dimension == 2 && codim == 1 )
592 {
593 auto refElem = referenceElement< FieldType, dimension >( entity.type() );
594
595#ifndef NDEBUG
596 const int vxSize = refElem.size( subEntity, codim, dimension );
597 // two vertices per edge in 2d
598 assert( vxSize == 2 );
599#endif
600 const int vx[ 2 ] = { refElem.subEntity ( subEntity, codim, 0, dimension ),
601 refElem.subEntity ( subEntity, codim, 1, dimension ) };
602
603 // flip index if face is twisted
604 if( gridPart_.grid().localIdSet().subId( entity, vx[ 0 ], dimension ) >
605 gridPart_.grid().localIdSet().subId( entity, vx[ 1 ], dimension ) )
606 {
607 const unsigned int numDofsSubEntity = compLocalKey.numDofs( codim, subEntity );
608 index = numDofsSubEntity - index - 1;
609 }
610 }
611
612 assert( index < compLocalKey.numDofs( codim, subEntity ) );
613 return dofContainer( codim )( entity, subEntity ).dof( codim, polOrd, index );
614 }
615 }
616
618 template< class Entity, class Functor >
619 void mapEachEntityDof ( const Entity &entity, Functor f ) const
620 {
621 const int n = numEntityDofs( entity );
622 for( int i = 0; i < n; ++i )
623 f( i, dofContainer( Entity::codimension )[ entity ].entityDof( i ) );
624 }
625
626 void onSubEntity ( const ElementType &element, int i, int c, std::vector< bool > &indices ) const
627 {
628 indices.resize( numDofs(element) );
630 {
631 if (c == 0)
632 std::fill(indices.begin(),indices.end(),true);
633 else
634 std::fill(indices.begin(),indices.end(),false);
635 }
636 else
637 {
638 DUNE_THROW( NotImplemented, "Method onSubEntity(...) not yet implemented for TupleMapper" );
639 }
640 }
641
642 void map ( const ElementType &element, std::vector< SizeType > &indices ) const
643 {
644 indices.resize( numDofs( element ) );
645 mapEach( element, AssignFunctor< std::vector< SizeType > >( indices ) );
646 }
647
648 template< class Entity >
649 void mapEntityDofs ( const Entity &entity, std::vector< SizeType > &indices ) const
650 {
651 indices.resize( numEntityDofs( entity ) );
652 mapEachEntityDof( entity, AssignFunctor< std::vector< SizeType > >( indices ) );
653 }
654
656 int maxNumDofs () const
657 {
658 return maxNumDofs_;
659 }
660
662 int numDofs ( const ElementType &entity ) const
663 {
664 const int polOrd = polynomOrder( entity );
665 return compiledLocalKey( polOrd, entity.type() ).size();
666 }
667
669 template< class Entity >
670 int numEntityDofs ( const Entity &entity ) const
671 {
673 {
674 if( Entity :: codimension == 0 )
675 return dofContainer( 0 )[ entity ].entityDofs();
676 else
677 return 0;
678 }
679 else
680 {
681 // !!! to be revised
682 // This implementation only works for nonhybrid grids (simplices or cubes)
683 return dofContainer( Entity :: codimension )[ entity ].entityDofs();
684 }
685 }
686
688 bool contains ( int codim ) const
689 {
690 return (discontinuousMapper) ? (codim == 0) : true;
691 }
692
694 bool fixedDataSize ( int codim ) const
695 {
696 return false;
697 }
698
700 int oldIndex ( const int hole, const int block ) const
701 {
702 assert( oldIndex_[ hole ] >= 0 );
703 return oldIndex_[ hole ];
704 }
705
707 int newIndex ( const int hole, const int block ) const
708 {
709 assert( newIndex_[ hole ] >= 0 );
710 return newIndex_[ hole ];
711 }
712
714 int numberOfHoles ( const int block ) const
715 {
716 return numberOfHoles_;
717 }
718
721 int numBlocks () const
722 {
723 return 1;
724 }
725
728 int oldOffSet ( const int block ) const
729 {
730 return 0;
731 }
732
735 int offSet ( const int block ) const
736 {
737 return 0;
738 }
739
741 bool consecutive () const
742 {
743 return true;
744 }
745
746 // Adaptation Methods (as for Index Sets)
748 {
749 entityPolynomOrder_.resize( PolynomialOrderStorageType( order_ ) );
750 entityPolynomOrder_.shrinkToFit();
751 for( int codim = 0; codim <= highestDimension; ++codim )
752 {
753 dofContainer( codim ).resize();
754 dofContainer( codim ).shrinkToFit();
755 }
756 }
757
758 void insertEntity ( const ElementType &entity )
759 {
761 insertEntityDofs( entity );
762 }
763
764 // return number of local dofs that were not visited yet
765 unsigned int insertEntityDofs( const ElementType &entity )
766 {
767 PolynomialOrderStorageType& polyStorage = entityPolynomOrder_[ entity ];
768 if( ! polyStorage.active() )
769 {
770 unsigned int notAlreadyCounted = 0;
771
772 const int polOrd = polyStorage.order();
773 // get lagrange point set
774 const CompiledLocalKeyType& clk = compiledLocalKey( polOrd, entity.type() );
775
776 //std::cout << "Insert Entity " << gridPart_.grid().localIdSet().id( entity ) << std::endl;
777
778 // activate dofs for this entity
779 polyStorage.activate();
780
781 // insert for all sub entities
783 apply( entity, clk, polOrd, size_, notAlreadyCounted, dofContainer_ );
784
785 //printEntityDofs( entity );
786 return notAlreadyCounted ;
787 }
788
789 return 0;
790 }
791
792 void removeEntity ( const ElementType &entity )
793 {
794 int polOrd;
795 // polOrd ist set on call of deactivate
796 if( entityPolynomOrder_[ entity ].deactivate( polOrd ) )
797 {
799 apply( entity, polOrd, dofContainer_ );
800 }
801 }
802
803 void resize ()
804 {
807 }
808
810 void adapt()
811 {
812 // set new polynomial orders for entities
813 for( auto& pol : entityPolynomOrder_ )
814 pol.update();
815
816 sequence_ = -1;
817 compress();
818 }
819
820 // insert father element when conforming refinement is enabled
821 // (due to refinement of more than one level)
822 unsigned int insertFather( const ElementType &entity )
823 {
824 if( entity.hasFather() )
825 {
826 // if father is a new element, insert it
827 ElementType dad = entity.father();
828 if( dad.isNew() )
829 {
830 unsigned int usedSize = insertEntityDofs( dad );
831 // also insert dad's fathers
832 usedSize += insertFather( dad );
833 return usedSize;
834 }
835 }
836 return 0;
837 }
838
840 bool considerHierarchy () const
841 {
842 return DGFGridInfo< GridType > :: refineStepsForHalf() > 1 ||
843 ! Dune::Capabilities::hasSingleGeometryType<GridType> :: v ;
844 }
845
848 {
849 // reset all dof entries
850 setUnused();
851
852 // count current size
853 size_t usedSize = 0;
854
855 const bool considerHierarchyOfElements = considerHierarchy();
856
857 typedef typename GridPartType :: template Codim< 0 > :: IteratorType IteratorType;
858 const IteratorType end = gridPart_.template end<0>();
859 for( IteratorType it = gridPart_.template begin<0>();
860 it != end ; ++it )
861 {
862 const ElementType &element = *it;
863 if( considerHierarchyOfElements )
864 {
865 // insert father elements (conforming and hybrid grids only)
866 usedSize += insertFather( element );
867 }
868
869 // number of new dofs (not already counted) is returned
870 usedSize += insertEntityDofs( element );
871 }
872
873 //std::cout << "Insert Size = " << size_ << std::endl;
874 //printDofs();
875 //std::cout << "Insert done! " << std::endl;
876 return usedSize ;
877 }
878
879 void printDofs() const
880 {
881 //std::cout << "Size " << size_ << std::endl;
882 typedef typename GridPartType :: template Codim< 0 > :: IteratorType IteratorType;
883 const IteratorType end = gridPart_.template end<0>();
884 for( IteratorType it = gridPart_.template begin<0>();
885 it != end ; ++it )
886 {
887 printEntityDofs( *it );
888 }
889 }
890
891 void printEntityDofs( const ElementType& entity ) const
892 {
893 std::cout << "Print entity " << gridPart_.grid().localIdSet().id( entity ) << " with " << std::endl;
894 for( int i = 0; i<numDofs( entity ); ++i )
895 {
896 std::cout << "en[ " << i << " ] = " << mapToGlobal( entity, i ) << std::endl;
897 }
898 }
899
902 {
903 // deactivate all entries in the polynomial order container
904 {
905 typedef typename PolyOrderContainerType :: Iterator Iterator;
906 const Iterator endit = entityPolynomOrder_.end();
907 for( Iterator it = entityPolynomOrder_.begin(); it != endit; ++it )
908 {
910 int pOrd;
911 p.deactivate( pOrd );
912 }
913 }
914
915 for( int codim = 0; codim <= highestDimension; ++codim )
916 {
917 DofContainerType& codimContainer = dofContainer( codim );
918 typedef typename DofContainerType :: Iterator Iterator;
919 const Iterator endit = codimContainer.end();
920 for( Iterator it = codimContainer.begin(); it != endit; ++it )
921 {
922 it->reset();
923 }
924 }
925 }
926
927 // --compress
928 bool compress ()
929 {
930 if( sequence_ == dm_.sequence() )
931 {
932 numberOfHoles_ = 0;
933 return false ;
934 }
935
936 // adjust size of containers
938
939 // make all dofs that are currently used
940 // (corresponding to GridPart's iterator)
941 const size_t usedSize = insertAllUsed();
942
943 //std::cout << "Size " << size_ << " holes " << numberOfHoles_ << std::endl;
944 //printDofs();
945
946 // reset number of holes
947 numberOfHoles_ = 0;
948
949 // true if a least one dof must be copied
950 bool haveToCopy = false;
951
952 std::vector<int> validHoles;
953
954 {
955 // default is true which means entry is hole
956 std::vector< bool > holeMarker( size_, true );
957
958 // mark holes
959 for( int codim = 0; codim <= highestDimension; ++codim )
960 {
961 DofContainerType& codimContainer = dofContainer( codim );
962 typedef typename DofContainerType :: Iterator Iterator;
963 const Iterator endit = codimContainer.end();
964 for( Iterator it = codimContainer.begin(); it != endit; ++it )
965 {
966 EntityDofStorageType& dof = *it;
967 // store holes if unused dofs exits, otherwise increase actSize
968 dof.detectUnusedDofs( holeMarker, usedSize );
969 }
970 }
971
972 // check for real holes
973 validHoles.reserve( usedSize );
974 validHoles.resize( 0 );
975
976 // check the vector of hole flags and store numbers
977 for( size_t i=0; i<usedSize; ++i )
978 {
979 // it's only a valid hole, if it was not marked otherwise
980 if( holeMarker[ i ] )
981 {
982 //std::cout << "Dof number " << i << " is not used" << std::endl;
983 validHoles.push_back( i );
984 }
985 }
986 }
987
988 if( validHoles.size() > 0 )
989 {
990 // counter for current hole to use
991 int currentHole = validHoles.size();
992
993 // resize old-new index storage to correct size
994 oldIndex_.resize( currentHole, -1) ;
995 newIndex_.resize( currentHole, -1) ;
996
997 for( int codim = 0; codim <= highestDimension; ++codim )
998 {
999 DofContainerType& codimContainer = dofContainer( codim );
1000 typedef typename DofContainerType :: Iterator Iterator;
1001 const Iterator endit = codimContainer.end();
1002 for( Iterator it = codimContainer.begin(); it != endit; ++it )
1003 {
1004 // get dof storage
1005 EntityDofStorageType& dof = *it;
1006
1007 // replace DoF larger than size with DoF from holes
1008 // set haveToCopy to true if true is returned
1009 haveToCopy |= dof.removeHoles( oldIndex_, newIndex_,
1010 validHoles, currentHole,
1011 usedSize, numberOfHoles_ );
1012 }
1013 }
1014 }
1015
1016 // store new size
1017 size_ = usedSize ;
1018
1019 //std::cout << "Size " << size_ << " holes " << numberOfHoles_ << std::endl;
1020 //printDofs();
1021
1022 sequence_ = dm_.sequence();
1023
1024 return haveToCopy;
1025 }
1026
1027 void backup () const
1028 {}
1029
1030 void restore ()
1031 {}
1032
1033 template< class InStream >
1034 void read ( InStream &in )
1035 {}
1036
1037 template< class OutStream >
1038 void write ( OutStream &out )
1039 {}
1040
1042 ThisType& operator=( const ThisType& ) = delete;
1043
1044 private:
1045
1046 const GridPartType& gridPart_;
1047 // reference to dof manager
1048 DofManagerType& dm_;
1049
1050 CompiledLocalKeyVectorType &compiledLocalKeys_;
1051
1052 // default polynomial order to be set to newly inserted elements
1053 const int order_;
1054
1055 PolyOrderContainerType entityPolynomOrder_;
1056
1057 mutable std::vector< DofContainerType* > dofContainer_;
1058
1059 int numberOfHoles_ ;
1060 std::vector< int > oldIndex_ ;
1061 std::vector< int > newIndex_ ;
1062
1063 mutable unsigned int size_;
1064 unsigned int maxNumDofs_;
1065 int sequence_ ;
1066 }; // class GenericAdaptiveDofMapper
1067
1068
1069 namespace Capabilities
1070 {
1071 // isConsecutiveIndexSet
1072 // ---------------------
1073
1074 template< class Traits >
1076 {
1077 static const bool v = true;
1078 };
1079
1080 template< class Traits >
1082 {
1083 static const bool v = true;
1084 };
1085
1086 } // namespace Capabilities
1087
1088 } // namespace Fem
1089
1090} // namespace Dune
1091
1092#endif // #ifndef DUNE_FEM_GENERICPADAPTIVEDOFMAPPER_HH
void removeIndexSet(const IndexSetType &iset)
removed index set from dof manager's list of index sets
Definition: dofmanager.hh:1291
void addIndexSet(const IndexSetType &iset)
add index set to dof manager's list of index sets
Definition: dofmanager.hh:1256
int sequence() const
return number of sequence, if dofmanagers memory was changed by calling some method like resize,...
Definition: dofmanager.hh:978
STL namespace.
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
Definition: bindguard.hh:11
Double abs(const Double &a)
Definition: double.hh:871
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
static DUNE_PRIVATE void apply(Args &&... args)
Definition: forloop.hh:23
specialize with true if index set implements the interface for consecutive index sets
Definition: common/indexset.hh:42
static const bool v
Definition: common/indexset.hh:49
Definition: misc/functor.hh:31
Definition: dofmanager.hh:761
Definition: space/mapper/capabilities.hh:22
static const bool v
Definition: space/mapper/capabilities.hh:23
Extended interface for adaptive DoF mappers.
Definition: mapper/dofmapper.hh:219
const Implementation & asImp() const
Definition: bartonnackmaninterface.hh:37
Definition: genericadaptivedofmapper.hh:32
int offSet(const int block) const
Definition: genericadaptivedofmapper.hh:735
GenericAdaptiveDofMapper(const ThisType &)=delete
@ minOrder
Definition: genericadaptivedofmapper.hh:83
void printEntityDofs(const ElementType &entity) const
Definition: genericadaptivedofmapper.hh:891
EntityDofStorage EntityDofStorageType
Definition: genericadaptivedofmapper.hh:310
unsigned int insertEntityDofs(const ElementType &entity)
Definition: genericadaptivedofmapper.hh:765
static const bool discontinuousMapper
Definition: genericadaptivedofmapper.hh:41
bool consecutive() const
Definition: genericadaptivedofmapper.hh:741
int numberOfHoles(const int block) const
Definition: genericadaptivedofmapper.hh:714
int newIndex(const int hole, const int block) const
Definition: genericadaptivedofmapper.hh:707
PolynomialOrderStorage PolynomialOrderStorageType
Definition: genericadaptivedofmapper.hh:348
void setUnused()
reset all used flags of all DoF entries
Definition: genericadaptivedofmapper.hh:901
void insertEntity(const ElementType &entity)
Definition: genericadaptivedofmapper.hh:758
void mapEntityDofs(const Entity &entity, std::vector< SizeType > &indices) const
Definition: genericadaptivedofmapper.hh:649
std::size_t SizeType
Definition: genericadaptivedofmapper.hh:38
DofManager< GridType > DofManagerType
type of the DoF manager
Definition: genericadaptivedofmapper.hh:81
@ maxOrder
Definition: genericadaptivedofmapper.hh:84
Traits::ElementType ElementType
type of entities (codim 0)
Definition: genericadaptivedofmapper.hh:51
GenericAdaptiveDofMapper(const GenericAdaptiveDofMapper &other, const int order, CompiledLocalKeyVectorType &compiledLocalKeyVector)
sort of copy constructor
Definition: genericadaptivedofmapper.hh:490
Traits::GridPartType GridPartType
type of the grid part
Definition: genericadaptivedofmapper.hh:48
GridPartType::GridType GridType
type of the underlying grid
Definition: genericadaptivedofmapper.hh:54
PersistentContainer< GridType, PolynomialOrderStorageType > PolyOrderContainerType
Definition: genericadaptivedofmapper.hh:352
void onSubEntity(const ElementType &element, int i, int c, std::vector< bool > &indices) const
Definition: genericadaptivedofmapper.hh:626
@ numOrders
Definition: genericadaptivedofmapper.hh:85
GridPartType::IndexSetType IndexSetType
type of the index set
Definition: genericadaptivedofmapper.hh:57
Traits::CompiledLocalKeyVectorType CompiledLocalKeyVectorType
type of vector containing compiled local keys
Definition: genericadaptivedofmapper.hh:75
int suggestedOrder(const ElementType &entity) const
Definition: genericadaptivedofmapper.hh:517
void mapEachEntityDof(const Entity &entity, Functor f) const
map each local DoF number to a global key
Definition: genericadaptivedofmapper.hh:619
ThisType & operator=(const ThisType &)=delete
void restore()
Definition: genericadaptivedofmapper.hh:1030
TraitsImp Traits
Definition: genericadaptivedofmapper.hh:37
int numDofs(const ElementType &entity) const
obtain number of DoFs on an entity
Definition: genericadaptivedofmapper.hh:662
void write(OutStream &out)
Definition: genericadaptivedofmapper.hh:1038
void read(InStream &in)
Definition: genericadaptivedofmapper.hh:1034
GenericAdaptiveDofMapper(const GridPartType &gridPart, const int order, CompiledLocalKeyVectorType &compiledLocalKeyVector)
constructor
Definition: genericadaptivedofmapper.hh:460
int numEntityDofs(const Entity &entity) const
obtain number of DoFs actually belonging to an entity
Definition: genericadaptivedofmapper.hh:670
static const int polynomialOrder
order of the Lagrange polynoms
Definition: genericadaptivedofmapper.hh:72
int mapToGlobal(const ElementType &entity, const int localDof) const
Definition: genericadaptivedofmapper.hh:566
Traits::GlobalKeyType GlobalKeyType
type of global key
Definition: genericadaptivedofmapper.hh:60
void suggestPolynomOrder(const ElementType &entity, const int polOrd)
Definition: genericadaptivedofmapper.hh:522
PersistentContainer< GridType, EntityDofStorageType > DofContainerType
Definition: genericadaptivedofmapper.hh:350
size_t insertAllUsed()
return number of DoFs currently used for space
Definition: genericadaptivedofmapper.hh:847
bool fixedDataSize(int codim) const
Check, whether the data in a codimension has fixed size.
Definition: genericadaptivedofmapper.hh:694
DofContainerType & dofContainer(const std::size_t codim) const
Definition: genericadaptivedofmapper.hh:531
int oldOffSet(const int block) const
Definition: genericadaptivedofmapper.hh:728
int size() const
return overall number of degrees of freedom
Definition: genericadaptivedofmapper.hh:552
void adapt()
adjust mapper to newly set polynomial orders
Definition: genericadaptivedofmapper.hh:810
static const int dimension
dimension of the grid
Definition: genericadaptivedofmapper.hh:66
void removeEntity(const ElementType &entity)
Definition: genericadaptivedofmapper.hh:792
unsigned int insertFather(const ElementType &entity)
Definition: genericadaptivedofmapper.hh:822
CompiledLocalKeyVectorType::value_type::value_type CompiledLocalKeyType
compiled local key type
Definition: genericadaptivedofmapper.hh:78
void mapEach(const ElementType &element, Functor f) const
Definition: genericadaptivedofmapper.hh:558
int numBlocks() const
Definition: genericadaptivedofmapper.hh:721
int polynomOrder(const ElementType &entity) const
Definition: genericadaptivedofmapper.hh:512
virtual ~GenericAdaptiveDofMapper()
destructor
Definition: genericadaptivedofmapper.hh:546
void backup() const
Definition: genericadaptivedofmapper.hh:1027
GridType::ctype FieldType
type of coordinates within the grid
Definition: genericadaptivedofmapper.hh:63
bool contains(int codim) const
Check, whether any DoFs are associated with a codimension.
Definition: genericadaptivedofmapper.hh:688
bool compress()
Definition: genericadaptivedofmapper.hh:928
static const int highestDimension
highest codimension used to attach dofs
Definition: genericadaptivedofmapper.hh:69
const CompiledLocalKeyType & compiledLocalKey(const int polOrd, const GeometryType type) const
Definition: genericadaptivedofmapper.hh:539
void map(const ElementType &element, std::vector< SizeType > &indices) const
Definition: genericadaptivedofmapper.hh:642
int maxNumDofs() const
obtain maximal number of DoFs on one entity
Definition: genericadaptivedofmapper.hh:656
void resize()
Definition: genericadaptivedofmapper.hh:803
void printDofs() const
Definition: genericadaptivedofmapper.hh:879
int oldIndex(const int hole, const int block) const
Definition: genericadaptivedofmapper.hh:700
void resizeContainers()
Definition: genericadaptivedofmapper.hh:747
bool considerHierarchy() const
return true if elements can be refined more than once during adaptation
Definition: genericadaptivedofmapper.hh:840
Definition: genericadaptivedofmapper.hh:88
bool use(const int codim, const int polOrd)
returns true if entry has a reference count of 1
Definition: genericadaptivedofmapper.hh:141
int entityDofs() const
Definition: genericadaptivedofmapper.hh:215
int dof(const int codim, const int polOrd, const size_t dofNumber) const
Definition: genericadaptivedofmapper.hh:190
bool removeHoles(VectorType &oldIdx, VectorType &newIdx, VectorType &holesVec, int &currentHole, const int usedSize, int &holes)
Definition: genericadaptivedofmapper.hh:268
std::vector< DofVectorType > dofs_
Definition: genericadaptivedofmapper.hh:90
void insert(const GeometryType type, const int codim, const int polOrd, const int numDofs, const int startDof)
Definition: genericadaptivedofmapper.hh:148
int entityDof(int dofNumber) const
Definition: genericadaptivedofmapper.hh:199
EntityDofStorage(const EntityDofStorage &other)
Definition: genericadaptivedofmapper.hh:122
char used_[numOrders]
Definition: genericadaptivedofmapper.hh:93
int determineVectorEntry(const int codim, const int polOrd) const
Definition: genericadaptivedofmapper.hh:165
void printDofs() const
Definition: genericadaptivedofmapper.hh:256
void detectUnusedDofs(VectorType &isHole, const int actSize)
Definition: genericadaptivedofmapper.hh:226
GeometryType type_
Definition: genericadaptivedofmapper.hh:92
void assign(const EntityDofStorage &other)
Definition: genericadaptivedofmapper.hh:104
EntityDofStorage & operator=(const EntityDofStorage &other)
Definition: genericadaptivedofmapper.hh:128
EntityDofStorage()
Definition: genericadaptivedofmapper.hh:95
const GeometryType & type() const
Definition: genericadaptivedofmapper.hh:175
void reset()
Definition: genericadaptivedofmapper.hh:184
void remove(const int codim, const int polOrd)
Definition: genericadaptivedofmapper.hh:177
bool exists(const int codim, const int polOrd) const
Definition: genericadaptivedofmapper.hh:134
std::vector< int > DofVectorType
Definition: genericadaptivedofmapper.hh:89
Definition: genericadaptivedofmapper.hh:313
void set(const int k)
Definition: genericadaptivedofmapper.hh:330
int order() const
Definition: genericadaptivedofmapper.hh:322
void suggest(const int k)
Definition: genericadaptivedofmapper.hh:323
void update()
Definition: genericadaptivedofmapper.hh:345
int suggested() const
Definition: genericadaptivedofmapper.hh:344
bool deactivate(int &k)
Definition: genericadaptivedofmapper.hh:333
void activate()
Definition: genericadaptivedofmapper.hh:331
PolynomialOrderStorage(const int k)
Definition: genericadaptivedofmapper.hh:321
bool active() const
Definition: genericadaptivedofmapper.hh:332
Definition: genericadaptivedofmapper.hh:356
static int numDofs(const ElementType &entity, const CompiledLocalKeyType &clk, const int subEntity)
Definition: genericadaptivedofmapper.hh:357
static int numDofs(const ElementType &entity, const CompiledLocalKeyType &clk, const int subEntity)
Definition: genericadaptivedofmapper.hh:368
Definition: genericadaptivedofmapper.hh:381
static void insertDofs(const ElementType &entity, const CompiledLocalKeyType &clk, const int polOrd, const int subEntity, unsigned int &globalSize, unsigned int &notAlreadyCounted, EntityDofStorage &entityDofs)
Definition: genericadaptivedofmapper.hh:382
static void apply(const ElementType &entity, const CompiledLocalKeyType &clk, const int polOrd, unsigned int &globalSize, unsigned int &notAlreadyCounted, std::vector< DofContainerType * > dofContainers)
Definition: genericadaptivedofmapper.hh:416
Definition: genericadaptivedofmapper.hh:443
static void apply(const ElementType &entity, const int polOrd, std::vector< DofContainerType * > dofContainers)
Definition: genericadaptivedofmapper.hh:444
Definition: localkey.hh:21
unsigned int subEntity() const
Definition: localkey.hh:26
unsigned int index() const
Definition: localkey.hh:28
unsigned int codim() const
Definition: localkey.hh:27