dune-grid  2.7.1
albertagrid/indexsets.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ALBERTAGRIDINDEXSETS_HH
4 #define DUNE_ALBERTAGRIDINDEXSETS_HH
5 
6 #include <array>
7 
8 #include <dune/common/hybridutilities.hh>
9 #include <dune/common/stdstreams.hh>
10 #include <dune/common/std/utility.hh>
11 
12 #include <dune/grid/common/grid.hh>
14 
21 
22 #if HAVE_ALBERTA
23 
24 namespace Dune
25 {
26 
27  namespace Alberta
28  {
30 
32  }
33 
34 
35 
36  // AlbertaGridHierarchicIndexSet
37  // -----------------------------
38 
39  template< int dim, int dimworld >
41  : public IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim,dimworld >, int, std::array< GeometryType, 1 > >
42  {
45 
46  friend class AlbertaGrid< dim, dimworld >;
47 
48  public:
51 
52  typedef typename Base::IndexType IndexType;
53 
54  typedef typename Base::Types Types;
55 
56  static const int dimension = GridFamily::dimension;
57 
60 
61  private:
62  typedef typename GridFamily::Traits Traits;
63 
65 
66  class InitEntityNumber;
67 
68  template< int codim >
69  struct CreateEntityNumbers;
70 
71  template< int codim >
72  struct RefineNumbering;
73 
74  template< int codim >
75  struct CoarsenNumbering;
76 
77  explicit AlbertaGridHierarchicIndexSet ( const DofNumbering &dofNumbering );
78 
79  public:
81 
83  template< class Entity >
84  bool contains ( const Entity & ) const
85  {
86  return true;
87  }
88 
89  using Base::index;
90  using Base::subIndex;
91 
93  template< int cc >
94  IndexType index ( const typename Traits::template Codim< cc >::Entity &entity ) const
95  {
97  const EntityImp &entityImp = entity.impl();
98  return subIndex( entityImp.elementInfo(), entityImp.subEntity(), cc );
99  }
100 
102  template< int cc >
103  IndexType subIndex ( const typename Traits::template Codim< cc >::Entity &entity, int i, unsigned int codim ) const
104  {
106  const EntityImp &entityImp = entity.impl();
107 
108  int k = i;
109  if( cc > 0 )
110  {
111  auto refElement = ReferenceElements< Alberta::Real, dimension >::simplex();
112  k = refElement.subEntity( entityImp.subEntity(), cc, i, codim );
113  }
114 
115  const int j = entityImp.grid().generic2alberta( codim, k );
116  return subIndex( entityImp.elementInfo(), j, codim );
117  }
118 
120  IndexType size ( const GeometryType &type ) const
121  {
122  return (type.isSimplex() ? size( dimension - type.dim() ) : 0);
123  }
124 
126  IndexType size ( int codim ) const
127  {
128  assert( (codim >= 0) && (codim <= dimension) );
129  return indexStack_[ codim ].size();
130  }
131 
132  Types types ( int codim ) const
133  {
134  assert( (codim >= 0) && (codim <= dimension) );
135  return {{ GeometryTypes::simplex( dimension - codim ) }};
136  }
137 
139  const std::vector< GeometryType > &geomTypes( int codim ) const
140  {
141  assert( (codim >= 0) && (codim <= dimension) );
142  return geomTypes_[ codim ];
143  }
144 
145  IndexType subIndex ( const ElementInfo &elementInfo, int i, unsigned int codim ) const
146  {
147  assert( !elementInfo == false );
148  return subIndex( elementInfo.element(), i, codim );
149  }
150 
157  IndexType subIndex ( const Alberta::Element *element, int i, unsigned int codim ) const
158  {
159  IndexType *array = (IndexType *)entityNumbers_[ codim ];
160  const IndexType subIndex = array[ dofNumbering_( element, codim, i ) ];
161  assert( (subIndex >= 0) && (subIndex < size( codim )) );
162  return subIndex;
163  }
164 
165  void preAdapt ()
166  {
167  // set global pointer to index stack
169  {
170  assert( Alberta::currentIndexStack == 0 );
171  Alberta::currentIndexStack = indexStack_;
172  }
173  }
174 
175  void postAdapt ()
176  {
177  // remove global pointer to index stack
180  }
181 
182  void create ();
183  void read ( const std::string &filename );
184  bool write ( const std::string &filename ) const;
185 
186  void release ()
187  {
188  for( int i = 0; i <= dimension; ++i )
189  entityNumbers_[ i ].release();
190  }
191 
192  private:
193  template< int codim >
194  static IndexStack &getIndexStack ( const IndexVectorPointer &dofVector )
195  {
196  IndexStack *indexStack;
198  indexStack = dofVector.template getAdaptationData< IndexStack >();
199  else
200  indexStack = &Alberta::currentIndexStack[ codim ];
201  assert( indexStack != 0 );
202  return *indexStack;
203  }
204 
205  // access to the dof vectors
206  const DofNumbering &dofNumbering_;
207 
208  // index stacks providing new numbers during adaptation
209  IndexStack indexStack_[ dimension+1 ];
210 
211  // dof vectors storing the (persistent) numbering
212  IndexVectorPointer entityNumbers_[ dimension+1 ];
213 
214  // all geometry types contained in the grid
215  std::vector< GeometryType > geomTypes_[ dimension+1 ];
216  };
217 
218 
219 
220  // AlbertaGridHierarchicIndexSet::InitEntityNumber
221  // -----------------------------------------------
222 
223  template< int dim, int dimworld >
225  {
226  IndexStack &indexStack_;
227 
228  public:
229  InitEntityNumber ( IndexStack &indexStack )
230  : indexStack_( indexStack )
231  {}
232 
233  void operator() ( int &dof )
234  {
235  dof = indexStack_.getIndex();
236  }
237  };
238 
239 
240 
241  // AlbertaGridHierarchicIndexSet::CreateEntityNumbers
242  // --------------------------------------------------
243 
244  template< int dim, int dimworld >
245  template< int codim >
246  struct AlbertaGridHierarchicIndexSet< dim, dimworld >::CreateEntityNumbers
247  {
248  static void setup ( AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
249 
250  static void apply ( const Alberta::HierarchyDofNumbering< dimension > &dofNumbering,
252 
253  static void apply ( const std::string &filename,
256  };
257 
258 
259 
260  // AlbertaGridHierarchicIndexSet::RefineNumbering
261  // ----------------------------------------------
262 
263  template< int dim, int dimworld >
264  template< int codim >
265  struct AlbertaGridHierarchicIndexSet< dim, dimworld >::RefineNumbering
266  {
267  static const int dimension = dim;
268  static const int codimension = codim;
269 
270  private:
271  typedef Alberta::DofAccess< dimension, codimension > DofAccess;
272 
273  explicit RefineNumbering ( const IndexVectorPointer &dofVector )
274  : indexStack_( getIndexStack< codimension >( dofVector ) ),
275  dofVector_( dofVector ),
276  dofAccess_( dofVector.dofSpace() )
277  {}
278 
279  public:
280  void operator() ( const Alberta::Element *child, int subEntity );
281 
282  typedef Alberta::Patch< dimension > Patch;
283  static void interpolateVector ( const IndexVectorPointer &dofVector,
284  const Patch &patch );
285 
286  private:
287  IndexStack &indexStack_;
288  IndexVectorPointer dofVector_;
289  DofAccess dofAccess_;
290  };
291 
292 
293 
294  // AlbertaGridHierarchicIndexSet::CoarsenNumbering
295  // -----------------------------------------------
296 
297  template< int dim, int dimworld >
298  template< int codim >
299  struct AlbertaGridHierarchicIndexSet< dim, dimworld >::CoarsenNumbering
300  {
301  static const int dimension = dim;
302  static const int codimension = codim;
303 
304  private:
305  typedef Alberta::DofAccess< dimension, codimension > DofAccess;
306 
307  explicit CoarsenNumbering ( const IndexVectorPointer &dofVector )
308  : indexStack_( getIndexStack< codimension >( dofVector ) ),
309  dofVector_( dofVector ),
310  dofAccess_( dofVector.dofSpace() )
311  {}
312 
313  public:
314  void operator() ( const Alberta::Element *child, int subEntity );
315 
316  typedef Alberta::Patch< dimension > Patch;
317  static void restrictVector ( const IndexVectorPointer &dofVector,
318  const Patch &patch );
319  private:
320  IndexStack &indexStack_;
321  IndexVectorPointer dofVector_;
322  DofAccess dofAccess_;
323  };
324 
325 
326 
327  // AlbertaGridIndexSet
328  // -------------------
329 
330  template< int dim, int dimworld >
332  : public IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > >
333  {
335  typedef IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > > Base;
336 
337  public:
339 
340  typedef typename Base::IndexType IndexType;
341 
342  typedef typename Base::Types Types;
343 
344  static const int dimension = Grid::dimension;
345 
348 
349  private:
350  typedef typename Grid::Traits Traits;
351 
352  template< int codim >
353  struct Insert;
354 
355  public:
356  explicit AlbertaGridIndexSet ( const DofNumbering &dofNumbering )
357  : dofNumbering_( dofNumbering )
358  {
359  for( int codim = 0; codim <= dimension; ++codim )
360  {
361  indices_[ codim ] = 0;
362  geomTypes_[ codim ].push_back( GeometryTypes::simplex( dimension - codim ) );
363  }
364  }
365 
367  {
368  for( int codim = 0; codim <= dimension; ++codim )
369  delete[] indices_[ codim ];
370  }
371 
372  template< class Entity >
373  bool contains ( const Entity &entity ) const
374  {
375  const int codim = Entity::codimension;
376 
378  = entity.impl();
379  const Alberta::Element *element = entityImp.elementInfo().el();
380 
381  const IndexType *const array = indices_[ codim ];
382  const IndexType subIndex = array[ dofNumbering_( element, codim, entityImp.subEntity() ) ];
383 
384  return (subIndex >= 0);
385  }
386 
387  using Base::index;
388  using Base::subIndex;
389 
391  template< int cc >
392  IndexType index ( const typename Traits::template Codim< cc >::Entity &entity ) const
393  {
395  const EntityImp &entityImp = entity.impl();
396  return subIndex( entityImp.elementInfo(), entityImp.subEntity(), cc );
397  }
398 
400  template< int cc >
401  IndexType subIndex ( const typename Traits::template Codim< cc >::Entity &entity, int i, unsigned int codim ) const
402  {
404  const EntityImp &entityImp = entity.impl();
405 
406  int k = i;
407  if( cc > 0 )
408  {
409  auto refElement = ReferenceElements< Alberta::Real, dimension >::simplex();
410  k = refElement.subEntity( entityImp.subEntity(), cc, i, codim );
411  }
412 
413  const int j = entityImp.grid().generic2alberta( codim, k );
414  return subIndex( entityImp.elementInfo(), j, codim );
415  }
416 
417  IndexType size ( const GeometryType &type ) const
418  {
419  return (type.isSimplex() ? size( dimension - type.dim() ) : 0);
420  }
421 
422  IndexType size ( int codim ) const
423  {
424  assert( (codim >= 0) && (codim <= dimension) );
425  return size_[ codim ];
426  }
427 
428  Types types ( int codim ) const
429  {
430  assert( (codim >= 0) && (codim <= dimension) );
431  return {{ GeometryTypes::simplex( dimension - codim ) }};
432  }
433 
434  const std::vector< GeometryType > &geomTypes( int codim ) const
435  {
436  assert( (codim >= 0) && (codim <= dimension) );
437  return geomTypes_[ codim ];
438  }
439 
440  template< class Iterator >
441  void update ( const Iterator &begin, const Iterator &end )
442  {
443  for( int codim = 0; codim <= dimension; ++codim )
444  {
445  delete[] indices_[ codim ];
446 
447  const unsigned int dofSize = dofNumbering_.size( codim );
448  indices_[ codim ] = new IndexType[ dofSize ];
449  for( unsigned int i = 0; i < dofSize; ++i )
450  indices_[ codim ][ i ] = -1;
451 
452  size_[ codim ] = 0;
453  }
454 
455  for( Iterator it = begin; it != end; ++it )
456  {
458  = it->impl();
459  const Alberta::Element *element = entityImp.elementInfo().el();
460  Hybrid::forEach( Std::make_index_sequence< dimension+1 >{},
461  [ & ]( auto i ){ Insert< i >::apply( element, *this ); } );
462  }
463  }
464 
465  private:
466  IndexType subIndex ( const ElementInfo &elementInfo, int i, unsigned int codim ) const
467  {
468  assert( !elementInfo == false );
469  return subIndex( elementInfo.element(), i, codim );
470  }
471 
478  IndexType subIndex ( const Alberta::Element *element, int i, unsigned int codim ) const
479  {
480  const IndexType *const array = indices_[ codim ];
481  const IndexType subIndex = array[ dofNumbering_( element, codim, i ) ];
482  assert( (subIndex >= 0) && (subIndex < size( codim )) );
483  return subIndex;
484  }
485 
486  // access to the dof vectors
487  const DofNumbering &dofNumbering_;
488 
489  // an array of indices for each codimension
490  IndexType *indices_[ dimension+1 ];
491 
492  // the size of each codimension
493  IndexType size_[ dimension+1 ];
494 
495  // all geometry types contained in the grid
496  std::vector< GeometryType > geomTypes_[ dimension+1 ];
497  };
498 
499 
500 
501  // AlbertaGridIndexSet::Insert
502  // ---------------------------
503 
504  template< int dim, int dimworld >
505  template< int codim >
506  struct AlbertaGridIndexSet< dim, dimworld >::Insert
507  {
508  static void apply ( const Alberta::Element *const element,
509  AlbertaGridIndexSet< dim, dimworld > &indexSet )
510  {
511  int *const array = indexSet.indices_[ codim ];
512  IndexType &size = indexSet.size_[ codim ];
513 
514  for( int i = 0; i < Alberta::NumSubEntities< dim, codim >::value; ++i )
515  {
516  int &index = array[ indexSet.dofNumbering_( element, codim, i ) ];
517  if( index < 0 )
518  index = size++;
519  }
520  }
521  };
522 
523 
524 
525  // AlbertaGridIdSet
526  // ----------------
527 
529  template< int dim, int dimworld >
531  : public IdSet< AlbertaGrid< dim, dimworld >, AlbertaGridIdSet< dim, dimworld >, unsigned int >
532  {
535 
536  friend class AlbertaGrid< dim, dimworld >;
537 
538  public:
540  typedef typename Base::IdType IdType;
541 
542  private:
544 
545  static const int dimension = Grid::dimension;
546 
548 
549  // create id set, only allowed for AlbertaGrid
550  AlbertaGridIdSet ( const HierarchicIndexSet &hIndexSet )
551  : hIndexSet_( hIndexSet )
552  {}
553 
554  public:
556  template< class Entity >
557  IdType id ( const Entity &e ) const
558  {
559  const int codim = Entity::codimension;
560  return id< codim >( e );
561  }
562 
564  template< int codim >
565  IdType id ( const typename Grid::template Codim< codim >::Entity &e ) const
566  {
567  assert( (codim >= 0) && (codim <= dimension) );
568  const IdType index = hIndexSet_.index( e );
569  return (index << 2) | IdType( codim );
570  }
571 
573  IdType subId ( const typename Grid::template Codim< 0 >::Entity &e, int i, unsigned int subcodim ) const
574  {
575  assert( int( subcodim ) <= dimension );
576  const IdType index = hIndexSet_.subIndex( e, i, subcodim );
577  return (index << 2) | IdType( subcodim );
578  }
579 
580  template< int codim >
581  IdType subId ( const typename Grid::template Codim< codim >::Entity &e, int i, unsigned int subcodim ) const
582  {
583  assert( (codim >= 0) && (codim <= dimension) && (int( codim + subcodim ) <= dimension) );
584  const IdType index = hIndexSet_.subIndex( e, i, subcodim );
585  return (index << 2) | IdType( codim + subcodim );
586  }
587 
588  template< class Entity >
589  IdType subId ( const Entity &e, int i, unsigned int subcodim ) const
590  {
591  return subId< Entity::codimension >( e, i, subcodim );
592  }
593 
594  private:
595  // prohibit copying
596  AlbertaGridIdSet ( const This & );
597 
598  const HierarchicIndexSet &hIndexSet_;
599  };
600 
601 } // namespace Dune
602 
603 #endif // #if HAVE_ALBERTA
604 
605 #endif // #ifndef DUNE_ALBERTAGRIDINDEXSETS_HH
Provides base classes for index and id sets.
provides a wrapper for ALBERTA's el_info structure
Provides an index stack that supplies indices for element numbering for a grid (i....
Include standard header files.
Definition: agrid.hh:59
Dune::IndexStack< int, 100000 > IndexStack
Definition: albertagrid/indexsets.hh:29
ALBERTA EL Element
Definition: misc.hh:52
IndexStack * currentIndexStack
Definition: indexsets.cc:17
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:180
[ provides Dune::Grid ]
Definition: agrid.hh:139
static const int dimension
Definition: agrid.hh:175
int size(int codim) const
Definition: dofadmin.hh:160
static const bool supportsAdaptationData
Definition: dofvector.hh:185
Element * el() const
Definition: elementinfo.hh:735
const Element * element() const
Definition: elementinfo.hh:719
Definition: albertagrid/entity.hh:44
const ElementInfo & elementInfo() const
Definition: albertagrid/entity.hh:128
int subEntity() const
obtain number of the subentity within the element (in ALBERTA numbering)
Definition: albertagrid/entity.hh:146
Definition: albertagrid/indexsets.hh:42
IndexType subIndex(const typename Traits::template Codim< cc >::Entity &entity, int i, unsigned int codim) const
return subIndex of given enitiy's sub entity
Definition: albertagrid/indexsets.hh:103
Base::IndexType IndexType
Definition: albertagrid/indexsets.hh:52
void read(const std::string &filename)
Definition: indexsets.cc:151
IndexType size(int codim) const
return size of set
Definition: albertagrid/indexsets.hh:126
void release()
Definition: albertagrid/indexsets.hh:186
Alberta::IndexStack IndexStack
Definition: albertagrid/indexsets.hh:80
AlbertaGrid< dim, dimworld > Grid
Definition: albertagrid/indexsets.hh:49
IndexType index(const typename Traits::template Codim< cc >::Entity &entity) const
return hierarchic index of given entity
Definition: albertagrid/indexsets.hh:94
bool contains(const Entity &) const
return true if entity is contained in set
Definition: albertagrid/indexsets.hh:84
IndexType subIndex(const Alberta::Element *element, int i, unsigned int codim) const
obtain hierarchic subindex
Definition: albertagrid/indexsets.hh:157
void create()
Definition: indexsets.cc:143
Base::Types Types
Definition: albertagrid/indexsets.hh:54
bool write(const std::string &filename) const
Definition: indexsets.cc:159
Alberta::ElementInfo< dimension > ElementInfo
Definition: albertagrid/indexsets.hh:58
void postAdapt()
Definition: albertagrid/indexsets.hh:175
IndexType subIndex(const ElementInfo &elementInfo, int i, unsigned int codim) const
Definition: albertagrid/indexsets.hh:145
Alberta::HierarchyDofNumbering< dimension > DofNumbering
Definition: albertagrid/indexsets.hh:59
Types types(int codim) const
Definition: albertagrid/indexsets.hh:132
const std::vector< GeometryType > & geomTypes(int codim) const
return geometry types this set has indices for
Definition: albertagrid/indexsets.hh:139
AlbertaGridFamily< dim, dimworld > GridFamily
Definition: albertagrid/indexsets.hh:50
IndexType size(const GeometryType &type) const
return size of set for given GeometryType
Definition: albertagrid/indexsets.hh:120
static const int dimension
Definition: albertagrid/indexsets.hh:56
void preAdapt()
Definition: albertagrid/indexsets.hh:165
hierarchic index set of AlbertaGrid
Definition: albertagrid/indexsets.hh:532
IdType id(const typename Grid::template Codim< codim >::Entity &e) const
Get id of an entity of codim cc. Unhandy because template parameter must be supplied explicitly.
Definition: albertagrid/indexsets.hh:565
IdType id(const Entity &e) const
Definition: albertagrid/indexsets.hh:557
Base::IdType IdType
export type of id
Definition: albertagrid/indexsets.hh:540
IdType subId(const typename Grid::template Codim< 0 >::Entity &e, int i, unsigned int subcodim) const
Get id of subentity i of co-dimension codim of a co-dimension 0 entity.
Definition: albertagrid/indexsets.hh:573
IdType subId(const typename Grid::template Codim< codim >::Entity &e, int i, unsigned int subcodim) const
Definition: albertagrid/indexsets.hh:581
IdType subId(const Entity &e, int i, unsigned int subcodim) const
Definition: albertagrid/indexsets.hh:589
Definition: albertagrid/indexsets.hh:333
Base::Types Types
Definition: albertagrid/indexsets.hh:342
Alberta::ElementInfo< dimension > ElementInfo
Definition: albertagrid/indexsets.hh:346
IndexType size(const GeometryType &type) const
Definition: albertagrid/indexsets.hh:417
IndexType index(const typename Traits::template Codim< cc >::Entity &entity) const
return hierarchic index of given entity
Definition: albertagrid/indexsets.hh:392
AlbertaGridIndexSet(const DofNumbering &dofNumbering)
Definition: albertagrid/indexsets.hh:356
Alberta::HierarchyDofNumbering< dimension > DofNumbering
Definition: albertagrid/indexsets.hh:347
bool contains(const Entity &entity) const
Definition: albertagrid/indexsets.hh:373
IndexType subIndex(const typename Traits::template Codim< cc >::Entity &entity, int i, unsigned int codim) const
return subIndex of given enitiy's sub entity
Definition: albertagrid/indexsets.hh:401
Base::IndexType IndexType
Definition: albertagrid/indexsets.hh:340
Types types(int codim) const
Definition: albertagrid/indexsets.hh:428
IndexType size(int codim) const
Definition: albertagrid/indexsets.hh:422
~AlbertaGridIndexSet()
Definition: albertagrid/indexsets.hh:366
void update(const Iterator &begin, const Iterator &end)
Definition: albertagrid/indexsets.hh:441
static const int dimension
Definition: albertagrid/indexsets.hh:344
const std::vector< GeometryType > & geomTypes(int codim) const
Definition: albertagrid/indexsets.hh:434
AlbertaGrid< dim, dimworld > Grid
Definition: albertagrid/indexsets.hh:338
Definition: albertagrid/gridfamily.hh:81
static const int dimension
Definition: albertagrid/gridfamily.hh:86
Definition: albertagrid/gridfamily.hh:96
Definition: albertagrid/indexsets.hh:225
InitEntityNumber(IndexStack &indexStack)
Definition: albertagrid/indexsets.hh:229
Definition: indexstack.hh:24
int size() const
return maxIndex which is also the
Definition: indexstack.hh:77
Wrapper class for entities.
Definition: common/entity.hh:64
Implementation & impl()
access to the underlying implementation
Definition: common/entity.hh:78
@ codimension
Know your own codimension.
Definition: common/entity.hh:105
Index Set Interface base class.
Definition: indexidset.hh:76
IndexType subIndex(const typename Traits::template Codim< cc >::Entity &e, int i, unsigned int codim) const
Map a subentity to an index.
Definition: indexidset.hh:151
std::array< GeometryType, 1 > Types
iterator range for geometry types in domain
Definition: indexidset.hh:93
IndexType index(const typename Traits::template Codim< cc >::Entity &e) const
Map entity to index. The result of calling this method with an entity that is not in the index set is...
Definition: indexidset.hh:111
Id Set Interface.
Definition: indexidset.hh:441
IdTypeImp IdType
Type used to represent an id.
Definition: indexidset.hh:444
Different resources needed by all grid implementations.
provides the GridFamily for AlbertaGrid