dune-fem  2.6-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  assert( 0 == EntityType :: codimension );
309  assert( checkValidIndex( leafIndex_( entity, subNumber ) ) );
310  return leafIndex_( entity, subNumber );
311  }
312 
314  template <class EntityType>
315  bool exists ( const EntityType& entity ) const
316  {
317  assert( myCodim_ == EntityType :: codimension );
318  const IndexType &index = leafIndex_[ entity ];
319  // if index is invalid (-1) it does not exist
320  if (index==invalidIndex()) return false;
321  assert( index < IndexType( indexState_.size() ) );
322  return (indexState_[ index ] != UNUSED);
323  }
324 
325  template <class EntityType>
326  bool exists ( const EntityType& entity ,
327  const int subNumber ) const
328  {
329  assert( 0 == EntityType :: codimension );
330  const IndexType &index = leafIndex_( entity, subNumber );
331  // if index is invalid (-1) it does not exist
332  if (index==invalidIndex()) return false;
333  assert( index < IndexType( indexState_.size() ) );
334  return (indexState_[ index ] != UNUSED);
335  }
336 
339  {
340  return numberHoles_;
341  }
342 
344  IndexType oldIndex (int elNum ) const
345  {
346  assert( numberHoles_ == IndexType(oldIdx_.size()) );
347  return oldIdx_[elNum];
348  }
349 
351  IndexType newIndex (int elNum) const
352  {
353  assert( numberHoles_ == IndexType(newIdx_.size()) );
354  return newIdx_[elNum];
355  }
356 
357  // insert element and create index for element number
358  template <class EntityType>
359  void insert (const EntityType& entity )
360  {
361  assert( myCodim_ == EntityType :: codimension );
362  insertIdx( leafIndex_[ entity ] );
363  }
364 
365  // insert element and create index for element number
366  template <class EntityType>
367  void insertSubEntity (const EntityType& entity,
368  const int subNumber)
369  {
370  assert( 0 == EntityType :: codimension );
371  insertIdx( leafIndex_( entity, subNumber ) );
372  }
373 
374  // insert element as ghost and create index for element number
375  template <class EntityType>
376  void insertGhost (const EntityType& entity )
377  {
378  assert( myCodim_ == EntityType :: codimension );
379  // insert index
380  IndexType &leafIdx = leafIndex_[ entity ];
381  insertIdx( leafIdx );
382 
383  // if index is also larger than lastSize
384  // mark as new to skip old-new index lists
385  if( leafIdx >= lastSize_ )
386  indexState_[ leafIdx ] = NEW;
387  }
388 
389  // insert element and create index for element number
390  template <class EntityType>
391  void markForRemoval( const EntityType& entity )
392  {
393  assert( myCodim_ == EntityType :: codimension );
394  const IndexType &index = leafIndex_[ entity ];
395  if( index != invalidIndex() )
396  indexState_[ index ] = UNUSED;
397  }
398 
399  // insert element as ghost and create index for element number
400  template <class EntityType>
401  bool validIndex (const EntityType& entity ) const
402  {
403  assert( myCodim_ == EntityType :: codimension );
404  return (leafIndex_[ entity ] >= 0);
405  }
406 
407  void print( std::ostream& out ) const
408  {
409  typedef typename IndexContainerType::ConstIterator Iterator;
410  const Iterator end = leafIndex_.end();
411  for( Iterator it = leafIndex_.begin(); it != end; ++it )
412  {
413  const IndexType &leafIdx = *it;
414  out << "idx: " << leafIdx << " stat: " << indexState_[ leafIdx ] << std::endl;
415  }
416  }
417 
418  protected:
419  // return true if the index idx is valid
420  bool checkValidIndex( const IndexType& idx ) const
421  {
422  assert( idx != invalidIndex() );
423  assert( idx < size() );
424  return (idx != invalidIndex() ) && ( idx < size() );
425  }
426 
427  // insert element and create index for element number
429  {
430  if( index == invalidIndex() )
431  {
432  index = indexState_.size();
433  indexState_.resize( index+1 );
434  }
435  assert( index < IndexType( indexState_.size() ) );
436  indexState_[ index ] = USED;
437  }
438 
439  public:
440  // write to stream
441  template <class StreamTraits>
443  {
444  // store current index set size
445  // don't write something like out << indexState_.size()
446  // since on read you then don't know exactly what
447  // type has been written, it must be the same types
448  const uint32_t indexSize = indexState_.size();
449  out << indexSize;
450 
451  // for consistency checking, write size as 64bit integer
452  const uint64_t mysize = leafIndex_.size();
453  out << mysize ;
454 
455  // backup indices
456  typedef typename IndexContainerType::ConstIterator ConstIterator;
457  const ConstIterator end = leafIndex_.end();
458  for( ConstIterator it = leafIndex_.begin(); it != end; ++it )
459  out << *it;
460 
461  return true;
462  }
463 
464  // read from stream
465  template <class StreamTraits>
467  {
468  // read current index set size
469  uint32_t size = 0;
470  in >> size;
471 
472  // mark all indices used
474  std::fill( indexState_.begin(), indexState_.end(), USED );
475 
476  // for consistency checking
477  uint64_t storedSize = 0;
478  in >> storedSize ;
479 
480  uint64_t leafsize = leafIndex_.size();
481  // the stored size can be larger (visualization of parallel grids in serial)
482  if( storedSize < leafsize )
483  {
484  DUNE_THROW(InvalidStateException,"CodimIndexSet: size consistency check failed during restore!");
485  }
486 
487  // restore indices
488  typedef typename IndexContainerType::Iterator Iterator;
489  const Iterator end = leafIndex_.end();
490  uint64_t count = 0 ;
491  for( Iterator it = leafIndex_.begin(); it != end; ++it, ++count )
492  in >> *it;
493 
494  // also read indices that were stored but are not needed on read
495  if( count < storedSize )
496  {
497  IndexType value ;
498  const uint64_t leftOver = storedSize - count ;
499  for( uint64_t i = 0; i < leftOver; ++i )
500  in >> value ;
501  }
502 
503  return true;
504  }
505 
506  }; // end of CodimIndexSet
507 
508  } // namespace Fem
509 
510 } // namespace Dune
511 
512 #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:326
IndexType additionalSizeEstimate() const
return how much extra memory is needed for restriction
Definition: codimindexset.hh:282
DynamicArray< IndexType > IndexArrayType
Definition: codimindexset.hh:51
void insert(const EntityType &entity)
Definition: codimindexset.hh:359
IndexType oldIndex(int elNum) const
return old index, for dof manager only
Definition: codimindexset.hh:344
void markForRemoval(const EntityType &entity)
Definition: codimindexset.hh:391
IndexType numberOfHoles() const
return number of holes
Definition: codimindexset.hh:338
void insertGhost(const EntityType &entity)
Definition: codimindexset.hh:376
bool read(InStreamInterface< StreamTraits > &in)
Definition: codimindexset.hh:466
const int myCodim_
Definition: codimindexset.hh:73
bool write(OutStreamInterface< StreamTraits > &out) const
Definition: codimindexset.hh:442
void insertSubEntity(const EntityType &entity, const int subNumber)
Definition: codimindexset.hh:367
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
bool compress()
Definition: codimindexset.hh:165
bool exists(const EntityType &entity) const
return state of index for given hierarchic number
Definition: codimindexset.hh:315
IndexArrayType newIdx_
Definition: codimindexset.hh:67
int IndexType
Definition: codimindexset.hh:44
bool checkValidIndex(const IndexType &idx) const
Definition: codimindexset.hh:420
void insertIdx(IndexType &index)
Definition: codimindexset.hh:428
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:407
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:401
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:351
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:164
void setMemoryFactor(double memFactor)
set memory factor
Definition: dynamicarray.hh:290
void resize(size_type nsize)
Definition: dynamicarray.hh:328