dune-fem  2.8-git
codimindexset.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_CODIMINDEXSET_HH
2 #define DUNE_FEM_CODIMINDEXSET_HH
3 
4 #include <algorithm>
5 #include <set>
6 
7 #include <dune/grid/utility/persistentcontainer.hh>
8 #include <dune/grid/utility/persistentcontainervector.hh>
9 #include <dune/grid/utility/persistentcontainerwrapper.hh>
10 #include <dune/grid/utility/persistentcontainermap.hh>
11 
14 
15 namespace Dune
16 {
17 
18  namespace Fem
19  {
20 
21  //***********************************************************************
22  //
23  // Index Set for one codimension
24  // --CodimIndexSet
25  //
26  //***********************************************************************
27  template <class GridImp>
29  {
30  protected:
31  typedef GridImp GridType;
33 
34  private:
35  enum INDEXSTATE { UNUSED = 0, // unused indices
36  USED = 1, // used indices
37  NEW = 2 };// new indices
38 
39  // reference to grid
40  const GridType& grid_;
41 
42  public:
43  // type of exported index
44  typedef int IndexType ;
45 
46  // indices in this status have not been initialized
47  static IndexType invalidIndex() { return -1; }
48 
49  protected:
50  // array type for indices
53 
54  // use the imporved PersistentContainer
55  typedef PersistentContainer< GridType, IndexType > IndexContainerType;
56 
57  // the mapping of the global to leaf index
60 
61  // stack for holes
63 
64  // Array that only remeber the occuring
65  // holes (for compress of data)
68 
69  // last size of set before compress (needed in parallel runs)
71 
72  // codim for which index is provided
73  const int myCodim_;
74 
75  // actual number of holes
77 
78  public:
79 
81  CodimIndexSet (const GridType& grid,
82  const int codim,
83  const double memoryFactor = 1.1)
84  : grid_( grid )
85  , leafIndex_( grid, codim, invalidIndex() )
86  , indexState_( 0 )
87  , holes_(0)
88  , oldIdx_(0)
89  , newIdx_(0)
90  , lastSize_ (0)
91  , myCodim_( codim )
92  , numberHoles_(0)
93  {
94  setMemoryFactor(memoryFactor);
95  }
96 
98  void setMemoryFactor(const double memoryFactor)
99  {
100  indexState_.setMemoryFactor( memoryFactor );
101  holes_.setMemoryFactor(memoryFactor);
102  oldIdx_.setMemoryFactor(memoryFactor);
103  newIdx_.setMemoryFactor(memoryFactor);
104  }
105 
107  void resize () { leafIndex_.resize( invalidIndex() ); }
108 
110  void prepareCompress () {}
111 
112  public:
114  void clear()
115  {
116  // set all values to invalidIndex
117  leafIndex_.fill( invalidIndex() );
118  // free all indices
119  indexState_.resize( 0 );
120  }
121 
123  void resetUsed ()
124  {
125  std::fill( indexState_.begin(), indexState_.end(), UNUSED );
126  }
127 
128  bool consecutive ()
129  {
130  std::set< int > found ;
131  // Something's wrong here: This method _must_ always return true
132  typedef typename IndexContainerType::Iterator Iterator;
133  bool consecutive = true;
134  const Iterator end = leafIndex_.end();
135  for( Iterator it = leafIndex_.begin(); it != end; ++it )
136  {
137  if( *it != invalidIndex() )
138  {
139  if( found.find( *it ) != found.end() )
140  {
141  std::cout << "index " << *it << " exists twice " << std::endl;
142  }
143  assert( found.find( *it ) == found.end() );
144  found.insert( *it );
145  }
146  consecutive &= (*it < IndexType( indexState_.size() ));
147  }
148  return consecutive;
149  }
150 
152  void checkConsecutive () { assert( consecutive() ); }
153 
155  void clearHoles()
156  {
157  // set number of holes to zero
158  numberHoles_ = 0;
159  // remember actual size
161  }
162 
165  bool compress ()
166  {
167  const int sizeOfVecs = indexState_.size();
168  holes_.resize( sizeOfVecs );
169 
170  // true if a least one dof must be copied
171  bool haveToCopy = false;
172 
173  // mark holes
174  int actHole = 0;
175  for( int index = 0; index < sizeOfVecs; ++index )
176  {
177  // create vector with all holes
178  if( indexState_[ index ] == UNUSED )
179  holes_[ actHole++ ] = index;
180  }
181 
182  // the new size is the actual size minus the holes
183  const int actSize = sizeOfVecs - actHole;
184 
185  // resize hole storing vectors
186  oldIdx_.resize(actHole);
187  newIdx_.resize(actHole);
188 
189  // only compress if number of holes > 0
190  if(actHole > 0)
191  {
192  // close holes
193  int holes = 0; // number of real holes
194  typedef typename IndexContainerType::Iterator Iterator;
195  const Iterator end = leafIndex_.end();
196  for( Iterator it = leafIndex_.begin(); it != end; ++it )
197  {
198  IndexType& index = *it;
199  if( index == invalidIndex() )
200  {
201  continue ;
202  }
203  else if( indexState_[ index ] == UNUSED )
204  {
205  index = invalidIndex();
206  continue ;
207  }
208 
209  // a index that is used but larger then actual size
210  // has to move to a hole
211  // if used index lies behind size, then index has to move
212  // to one of the holes
213  if( index >= actSize )
214  {
215  //std::cout << "Check index " << index << std::endl;
216  // serach next hole that is smaler than actual size
217  --actHole;
218  // if actHole < 0 then error, because we have index larger then
219  // actual size
220  assert(actHole >= 0);
221  while ( holes_[actHole] >= actSize )
222  {
223  --actHole;
224  if(actHole < 0) break;
225  }
226 
227  assert(actHole >= 0);
228 
229 #if HAVE_MPI
230  // only for none-ghost elements hole storage is applied
231  // this is because ghost indices might have in introduced
232  // after the resize was done.
233  if( indexState_[ index ] == USED )
234 #endif
235  {
236  // remember old and new index
237  oldIdx_[holes] = index;
238  newIdx_[holes] = holes_[actHole];
239  ++holes;
240  }
241 
242  index = holes_[actHole];
243 
244  // means that dof manager has to copy the mem
245  haveToCopy = true;
246  }
247  }
248 
249  // this call only sets the size of the vectors
250  oldIdx_.resize(holes);
251  newIdx_.resize(holes);
252 
253  // mark holes as new
254  // note: This needs to be done after reassignment, so that their
255  // original entry will still see them as UNUSED.
256  for( int hole = 0; hole < holes; ++hole )
257  indexState_[ newIdx_[ hole ] ] = NEW;
258 
259  } // end if actHole > 0
260 
261  // store number of actual holes
263 
264  // adjust size of container to correct value
265  leafIndex_.resize( invalidIndex() );
266 
267  // resize vector of index states
268  indexState_.resize( actSize );
269 
270 #ifndef NDEBUG
271  for( int i=0; i<actSize; ++i )
272  assert( indexState_[ i ] == USED ||
273  indexState_[ i ] == UNUSED ||
274  indexState_[ i ] == NEW );
275 
277 #endif
278  return haveToCopy;
279  }
280 
283 
285  IndexType size () const { return indexState_.size(); }
286 
289  {
290  return leafIndex_.size();
291  }
292 
294  //- --index
295  template <class EntityType>
296  IndexType index ( const EntityType& entity ) const
297  {
298  assert( myCodim_ == EntityType :: codimension );
299  assert( checkValidIndex( leafIndex_[ entity ] ) );
300  return leafIndex_[ entity ];
301  }
302 
304  template <class EntityType>
305  IndexType subIndex ( const EntityType& entity,
306  const int subNumber ) const
307  {
308  return subIndex( entity, subNumber,
309  std::integral_constant<bool, EntityType::codimension == 0 > () );
310  }
311 
313  template <class EntityType>
314  IndexType subIndex ( const EntityType& entity,
315  const int subNumber,
316  const std::integral_constant<bool,false> codim0 ) const
317  {
318  DUNE_THROW(NotImplemented,"CodimIndexSet::subIndex: not implemented for entities with codim > 0" );
319  return IndexType( -1 );
320  }
321 
323  template <class EntityType>
324  IndexType subIndex ( const EntityType& entity,
325  const int subNumber,
326  const std::integral_constant<bool,true> codim0 ) const
327  {
328  assert( checkValidIndex( leafIndex_( entity, subNumber ) ) );
329  return leafIndex_( entity, subNumber );
330  }
331 
333  template <class EntityType>
334  bool exists ( const EntityType& entity ) const
335  {
336  assert( myCodim_ == EntityType :: codimension );
337  const IndexType &index = leafIndex_[ entity ];
338  // if index is invalid (-1) it does not exist
339  if (index==invalidIndex()) return false;
340  assert( index < IndexType( indexState_.size() ) );
341  return (indexState_[ index ] != UNUSED);
342  }
343 
344  template <class EntityType>
345  bool exists ( const EntityType& entity ,
346  const int subNumber ) const
347  {
348  assert( 0 == EntityType :: codimension );
349  const IndexType &index = leafIndex_( entity, subNumber );
350  // if index is invalid (-1) it does not exist
351  if (index==invalidIndex()) return false;
352  assert( index < IndexType( indexState_.size() ) );
353  return (indexState_[ index ] != UNUSED);
354  }
355 
358  {
359  return numberHoles_;
360  }
361 
363  IndexType oldIndex (int elNum ) const
364  {
365  assert( numberHoles_ == IndexType(oldIdx_.size()) );
366  return oldIdx_[elNum];
367  }
368 
370  IndexType newIndex (int elNum) const
371  {
372  assert( numberHoles_ == IndexType(newIdx_.size()) );
373  return newIdx_[elNum];
374  }
375 
376  // insert element and create index for element number
377  template <class EntityType>
378  void insert (const EntityType& entity )
379  {
380  assert( myCodim_ == EntityType :: codimension );
381  insertIdx( leafIndex_[ entity ] );
382  }
383 
384  // insert element and create index for element number
385  template <class EntityType>
386  void insertSubEntity (const EntityType& entity,
387  const int subNumber)
388  {
389  assert( 0 == EntityType :: codimension );
390  insertIdx( leafIndex_( entity, subNumber ) );
391  }
392 
393  // insert element as ghost and create index for element number
394  template <class EntityType>
395  void insertGhost (const EntityType& entity )
396  {
397  assert( myCodim_ == EntityType :: codimension );
398  // insert index
399  IndexType &leafIdx = leafIndex_[ entity ];
400  insertIdx( leafIdx );
401 
402  // if index is also larger than lastSize
403  // mark as new to skip old-new index lists
404  if( leafIdx >= lastSize_ )
405  indexState_[ leafIdx ] = NEW;
406  }
407 
408  // insert element and create index for element number
409  template <class EntityType>
410  void markForRemoval( const EntityType& entity )
411  {
412  assert( myCodim_ == EntityType :: codimension );
413  const IndexType &index = leafIndex_[ entity ];
414  if( index != invalidIndex() )
415  indexState_[ index ] = UNUSED;
416  }
417 
418  // insert element as ghost and create index for element number
419  template <class EntityType>
420  bool validIndex (const EntityType& entity ) const
421  {
422  assert( myCodim_ == EntityType :: codimension );
423  return (leafIndex_[ entity ] >= 0);
424  }
425 
426  void print( std::ostream& out ) const
427  {
428  typedef typename IndexContainerType::ConstIterator Iterator;
429  const Iterator end = leafIndex_.end();
430  for( Iterator it = leafIndex_.begin(); it != end; ++it )
431  {
432  const IndexType &leafIdx = *it;
433  if( leafIdx < indexState_.size() )
434  {
435  out << "idx: " << leafIdx << " stat: " << indexState_[ leafIdx ] << std::endl;
436  }
437  }
438  }
439 
440  protected:
441  // return true if the index idx is valid
442  bool checkValidIndex( const IndexType& idx ) const
443  {
444  assert( idx != invalidIndex() );
445  assert( idx < size() );
446  return (idx != invalidIndex() ) && ( idx < size() );
447  }
448 
449  // insert element and create index for element number
451  {
452  if( index == invalidIndex() )
453  {
454  index = indexState_.size();
455  indexState_.resize( index+1 );
456  }
457  assert( index < IndexType( indexState_.size() ) );
458  indexState_[ index ] = USED;
459  }
460 
461  public:
462  // write to stream
463  template <class StreamTraits>
465  {
466  // store current index set size
467  // don't write something like out << indexState_.size()
468  // since on read you then don't know exactly what
469  // type has been written, it must be the same types
470  const uint32_t indexSize = indexState_.size();
471  out << indexSize;
472 
473  // for consistency checking, write size as 64bit integer
474  const uint64_t mysize = leafIndex_.size();
475  out << mysize ;
476 
477  // backup indices
478  typedef typename IndexContainerType::ConstIterator ConstIterator;
479  const ConstIterator end = leafIndex_.end();
480  for( ConstIterator it = leafIndex_.begin(); it != end; ++it )
481  out << *it;
482 
483  return true;
484  }
485 
486  // read from stream
487  template <class StreamTraits>
489  {
490  // read current index set size
491  uint32_t size = 0;
492  in >> size;
493 
494  // mark all indices used
496  std::fill( indexState_.begin(), indexState_.end(), USED );
497 
498  // for consistency checking
499  uint64_t storedSize = 0;
500  in >> storedSize ;
501 
502  uint64_t leafsize = leafIndex_.size();
503  // the stored size can be larger (visualization of parallel grids in serial)
504  if( storedSize < leafsize )
505  {
506  DUNE_THROW(InvalidStateException,"CodimIndexSet: size consistency check failed during restore!");
507  }
508 
509  // restore indices
510  typedef typename IndexContainerType::Iterator Iterator;
511  const Iterator end = leafIndex_.end();
512  uint64_t count = 0 ;
513  for( Iterator it = leafIndex_.begin(); it != end; ++it, ++count )
514  in >> *it;
515 
516  // also read indices that were stored but are not needed on read
517  if( count < storedSize )
518  {
519  IndexType value ;
520  const uint64_t leftOver = storedSize - count ;
521  for( uint64_t i = 0; i < leftOver; ++i )
522  in >> value ;
523  }
524 
525  return true;
526  }
527 
528  }; // end of CodimIndexSet
529 
530  } // namespace Fem
531 
532 } // namespace Dune
533 
534 #endif // #ifndef DUNE_FEM_CODIMINDEXSET_HH
Definition: bindguard.hh:11
Definition: codimindexset.hh:29
DynamicArray< INDEXSTATE > IndexStateArrayType
Definition: codimindexset.hh:52
void setMemoryFactor(const double memoryFactor)
set memory overestimation factor
Definition: codimindexset.hh:98
bool exists(const EntityType &entity, const int subNumber) const
Definition: codimindexset.hh:345
IndexType additionalSizeEstimate() const
return how much extra memory is needed for restriction
Definition: codimindexset.hh:282
DynamicArray< IndexType > IndexArrayType
Definition: codimindexset.hh:51
IndexType subIndex(const EntityType &entity, const int subNumber, const std::integral_constant< bool, false > codim0) const
return leaf index for given entity
Definition: codimindexset.hh:314
void insert(const EntityType &entity)
Definition: codimindexset.hh:378
IndexType oldIndex(int elNum) const
return old index, for dof manager only
Definition: codimindexset.hh:363
void markForRemoval(const EntityType &entity)
Definition: codimindexset.hh:410
IndexType numberOfHoles() const
return number of holes
Definition: codimindexset.hh:357
void insertGhost(const EntityType &entity)
Definition: codimindexset.hh:395
bool read(InStreamInterface< StreamTraits > &in)
Definition: codimindexset.hh:488
const int myCodim_
Definition: codimindexset.hh:73
bool write(OutStreamInterface< StreamTraits > &out) const
Definition: codimindexset.hh:464
void insertSubEntity(const EntityType &entity, const int subNumber)
Definition: codimindexset.hh:386
IndexType size() const
return size of grid entities per level and codim
Definition: codimindexset.hh:285
void resetUsed()
set all entries to unused
Definition: codimindexset.hh:123
IndexArrayType holes_
Definition: codimindexset.hh:62
CodimIndexSet< GridType > ThisType
Definition: codimindexset.hh:32
IndexStateArrayType indexState_
Definition: codimindexset.hh:59
GridImp GridType
Definition: codimindexset.hh:31
CodimIndexSet(const GridType &grid, const int codim, const double memoryFactor=1.1)
Constructor taking memory factor (default = 1.1)
Definition: codimindexset.hh:81
IndexType subIndex(const EntityType &entity, const int subNumber, const std::integral_constant< bool, true > codim0) const
return leaf index for given entity
Definition: codimindexset.hh:324
bool compress()
Definition: codimindexset.hh:165
bool exists(const EntityType &entity) const
return state of index for given hierarchic number
Definition: codimindexset.hh:334
IndexArrayType newIdx_
Definition: codimindexset.hh:67
int IndexType
Definition: codimindexset.hh:44
bool checkValidIndex(const IndexType &idx) const
Definition: codimindexset.hh:442
void insertIdx(IndexType &index)
Definition: codimindexset.hh:450
IndexType numberHoles_
Definition: codimindexset.hh:76
bool consecutive()
Definition: codimindexset.hh:128
IndexArrayType oldIdx_
Definition: codimindexset.hh:66
void print(std::ostream &out) const
Definition: codimindexset.hh:426
IndexContainerType leafIndex_
Definition: codimindexset.hh:58
IndexType realSize() const
return size of grid entities per level and codim
Definition: codimindexset.hh:288
void prepareCompress()
prepare for setup (nothing to do here)
Definition: codimindexset.hh:110
IndexType subIndex(const EntityType &entity, const int subNumber) const
return leaf index for given entity
Definition: codimindexset.hh:305
void clear()
clear set
Definition: codimindexset.hh:114
void resize()
reallocate the vectors
Definition: codimindexset.hh:107
void checkConsecutive()
set all entries to unused
Definition: codimindexset.hh:152
void clearHoles()
clear holes, i.e. set number of holes to zero
Definition: codimindexset.hh:155
bool validIndex(const EntityType &entity) const
Definition: codimindexset.hh:420
PersistentContainer< GridType, IndexType > IndexContainerType
Definition: codimindexset.hh:55
IndexType index(const EntityType &entity) const
return leaf index for given entity
Definition: codimindexset.hh:296
IndexType newIndex(int elNum) const
return new index, for dof manager only returns index
Definition: codimindexset.hh:370
IndexType lastSize_
Definition: codimindexset.hh:70
static IndexType invalidIndex()
Definition: codimindexset.hh:47
abstract interface for an output stream
Definition: streams.hh:46
abstract interface for an input stream
Definition: streams.hh:179
size_type size() const
return size of array
Definition: dynamicarray.hh:170
void setMemoryFactor(double memFactor)
set memory factor
Definition: dynamicarray.hh:296
void resize(size_type nsize)
Definition: dynamicarray.hh:334