dune-fem  2.8-git
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>
13 #include <dune/fem/misc/functor.hh>
19 
20 namespace 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_;
93  char used_[ numOrders ];
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 
538  const CompiledLocalKeyType&
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  {
568  if( discontinuousMapper )
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) );
629  if( discontinuousMapper )
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  {
672  if( discontinuousMapper )
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  {
806  insertAllUsed();
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 
847  size_t insertAllUsed()
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 
901  void setUnused()
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 
1041  GenericAdaptiveDofMapper ( const ThisType& ) = delete;
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:1289
void addIndexSet(const IndexSetType &iset)
add index set to dof manager's list of index sets
Definition: dofmanager.hh:1259
int sequence() const
return number of sequence, if dofmanagers memory was changed by calling some method like resize,...
Definition: dofmanager.hh:981
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 constexpr T max(T a)
Definition: utility.hh:77
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: bartonnackmaninterface.hh:17
Definition: misc/functor.hh:31
Definition: dofmanager.hh:762
Definition: space/mapper/capabilities.hh:22
static const bool v
Definition: space/mapper/capabilities.hh:23
Interface for calculating the size of a function space for a grid on a specified level....
Definition: mapper/dofmapper.hh:43
Traits::ElementType ElementType
type of codimension 0 entities
Definition: mapper/dofmapper.hh:54
Extended interface for adaptive DoF mappers.
Definition: mapper/dofmapper.hh:219
static const Implementation & asImp(const ThisType &other)
Definition: bartonnackmaninterface.hh:27
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
Traits ::GlobalKeyType GlobalKeyType
type of global key
Definition: genericadaptivedofmapper.hh:60
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
static const Implementation & asImp(const ThisType &other)
Definition: bartonnackmaninterface.hh:27
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
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
void restore()
Definition: genericadaptivedofmapper.hh:1030
ThisType & operator=(const ThisType &)=delete
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
Traits ::CompiledLocalKeyVectorType CompiledLocalKeyVectorType
type of vector containing compiled local keys
Definition: genericadaptivedofmapper.hh:75
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
const CompiledLocalKeyType & compiledLocalKey(const int polOrd, const GeometryType type) const
Definition: genericadaptivedofmapper.hh:539
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
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
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
void mapEach(const ElementType &element, Functor f) const
Definition: genericadaptivedofmapper.hh:558
CompiledLocalKeyVectorType ::value_type ::value_type CompiledLocalKeyType
compiled local key type
Definition: genericadaptivedofmapper.hh:78
int numBlocks() const
Definition: genericadaptivedofmapper.hh:721
DofContainerType & dofContainer(const std::size_t codim) const
Definition: genericadaptivedofmapper.hh:531
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
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
EntityDofStorage & operator=(const EntityDofStorage &other)
Definition: genericadaptivedofmapper.hh:128
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()
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