dune-fem 2.8.0
Loading...
Searching...
No Matches
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
15namespace 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;
32 typedef CodimIndexSet < GridType > ThisType;
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
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
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
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 {
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 {
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