dune-grid  2.2.0
albertagrid/indexsets.hh
Go to the documentation of this file.
1 #ifndef DUNE_ALBERTAGRIDINDEXSETS_HH
2 #define DUNE_ALBERTAGRIDINDEXSETS_HH
3 
4 #include <dune/common/stdstreams.hh>
5 
8 
15 
16 #if HAVE_ALBERTA
17 
18 namespace Dune
19 {
20 
21  namespace Alberta
22  {
24 
26  }
27 
28 
29 
30  // AlbertaGridHierarchicIndexSet
31  // -----------------------------
32 
33  template< int dim, int dimworld >
35  : public IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim,dimworld >, int >
36  {
39 
40  friend class AlbertaGrid< dim, dimworld >;
41 
42  public:
45 
46  typedef typename Base::IndexType IndexType;
47 
48  static const int dimension = GridFamily::dimension;
49 
52 
53  private:
54  typedef typename GridFamily::Traits Traits;
55 
57 
58  class InitEntityNumber;
59 
60  template< int codim >
61  struct CreateEntityNumbers;
62 
63  template< int codim >
64  class RefineNumbering;
65 
66  template< int codim >
67  class CoarsenNumbering;
68 
69  explicit AlbertaGridHierarchicIndexSet ( const DofNumbering &dofNumbering );
70 
71  public:
73 
75  template< class Entity >
76  bool contains ( const Entity & ) const
77  {
78  return true;
79  }
80 
81  using Base::index;
82  using Base::subIndex;
83 
85  template< int cc >
86  IndexType index ( const typename Traits::template Codim< cc >::Entity &entity ) const
87  {
89  const EntityImp &entityImp = Grid::getRealImplementation( entity );
90  return subIndex( entityImp.elementInfo(), entityImp.subEntity(), cc );
91  }
92 
94  template< int cc >
95  IndexType subIndex ( const typename Traits::template Codim< cc >::Entity &entity, int i, unsigned int codim ) const
96  {
98  const EntityImp &entityImp = Grid::getRealImplementation( entity );
99 
100  int k = i;
101  if( cc > 0 )
102  {
103  const GenericReferenceElement< Alberta::Real, dimension > &refElement
105  k = refElement.subEntity( entityImp.subEntity(), cc, i, codim );
106  }
107 
108  const int j = entityImp.grid().generic2alberta( codim, k );
109  return subIndex( entityImp.elementInfo(), j, codim );
110  }
111 
113  IndexType size ( const GeometryType &type ) const
114  {
115  return (type.isSimplex() ? size( dimension - type.dim() ) : 0);
116  }
117 
119  IndexType size ( int codim ) const
120  {
121  assert( (codim >= 0) && (codim <= dimension) );
122  return indexStack_[ codim ].size();
123  }
124 
126  const std::vector< GeometryType > &geomTypes( int codim ) const
127  {
128  assert( (codim >= 0) && (codim <= dimension) );
129  return geomTypes_[ codim ];
130  }
131 
132  IndexType subIndex ( const ElementInfo &elementInfo, int i, unsigned int codim ) const
133  {
134  assert( !elementInfo == false );
135  return subIndex( elementInfo.element(), i, codim );
136  }
137 
144  IndexType subIndex ( const Alberta::Element *element, int i, unsigned int codim ) const
145  {
146  IndexType *array = (IndexType *)entityNumbers_[ codim ];
147  const IndexType subIndex = array[ dofNumbering_( element, codim, i ) ];
148  assert( (subIndex >= 0) && (subIndex < size( codim )) );
149  return subIndex;
150  }
151 
152  void preAdapt ()
153  {
154  // set global pointer to index stack
156  {
157  assert( Alberta::currentIndexStack == 0 );
158  Alberta::currentIndexStack = indexStack_;
159  }
160  }
161 
162  void postAdapt ()
163  {
164  // remove global pointer to index stack
167  }
168 
169  void create ();
170  void read ( const std::string &filename );
171  bool write ( const std::string &filename ) const;
172 
173  void release ()
174  {
175  for( int i = 0; i <= dimension; ++i )
176  entityNumbers_[ i ].release();
177  }
178 
179  private:
180  template< int codim >
181  static IndexStack &getIndexStack ( const IndexVectorPointer &dofVector )
182  {
183  IndexStack *indexStack;
185  indexStack = dofVector.template getAdaptationData< IndexStack >();
186  else
187  indexStack = &Alberta::currentIndexStack[ codim ];
188  assert( indexStack != 0 );
189  return *indexStack;
190  }
191 
192  // access to the dof vectors
193  const DofNumbering &dofNumbering_;
194 
195  // index stacks providing new numbers during adaptation
196  IndexStack indexStack_[ dimension+1 ];
197 
198  // dof vectors storing the (persistent) numbering
199  IndexVectorPointer entityNumbers_[ dimension+1 ];
200 
201  // all geometry types contained in the grid
202  std::vector< GeometryType > geomTypes_[ dimension+1 ];
203  };
204 
205 
206 
207  // AlbertaGridHierarchicIndexSet::InitEntityNumber
208  // -----------------------------------------------
209 
210  template< int dim, int dimworld >
212  {
213  IndexStack &indexStack_;
214 
215  public:
216  InitEntityNumber ( IndexStack &indexStack )
217  : indexStack_( indexStack )
218  {}
219 
220  void operator() ( int &dof )
221  {
222  dof = indexStack_.getIndex();
223  }
224  };
225 
226 
227 
228  // AlbertaGridHierarchicIndexSet::CreateEntityNumbers
229  // --------------------------------------------------
230 
231  template< int dim, int dimworld >
232  template< int codim >
234  {
235  static void setup ( AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
236 
237  static void apply ( const Alberta::HierarchyDofNumbering< dimension > &dofNumbering,
239 
240  static void apply ( const std::string &filename,
243  };
244 
245 
246 
247  // AlbertaGridHierarchicIndexSet::RefineNumbering
248  // ----------------------------------------------
249 
250  template< int dim, int dimworld >
251  template< int codim >
253  {
254  static const int dimension = dim;
255  static const int codimension = codim;
256 
257  private:
259 
260  explicit RefineNumbering ( const IndexVectorPointer &dofVector )
261  : indexStack_( getIndexStack< codimension >( dofVector ) ),
262  dofVector_( dofVector ),
263  dofAccess_( dofVector.dofSpace() )
264  {}
265 
266  public:
267  void operator() ( const Alberta::Element *child, int subEntity );
268 
270  static void interpolateVector ( const IndexVectorPointer &dofVector,
271  const Patch &patch );
272 
273  private:
274  IndexStack &indexStack_;
275  IndexVectorPointer dofVector_;
276  DofAccess dofAccess_;
277  };
278 
279 
280 
281  // AlbertaGridHierarchicIndexSet::CoarsenNumbering
282  // -----------------------------------------------
283 
284  template< int dim, int dimworld >
285  template< int codim >
287  {
288  static const int dimension = dim;
289  static const int codimension = codim;
290 
291  private:
293 
294  explicit CoarsenNumbering ( const IndexVectorPointer &dofVector )
295  : indexStack_( getIndexStack< codimension >( dofVector ) ),
296  dofVector_( dofVector ),
297  dofAccess_( dofVector.dofSpace() )
298  {}
299 
300  public:
301  void operator() ( const Alberta::Element *child, int subEntity );
302 
304  static void restrictVector ( const IndexVectorPointer &dofVector,
305  const Patch &patch );
306  private:
307  IndexStack &indexStack_;
308  IndexVectorPointer dofVector_;
309  DofAccess dofAccess_;
310  };
311 
312 
313 
314  // AlbertaGridIndexSet
315  // -------------------
316 
317  template< int dim, int dimworld >
319  : public IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int >
320  {
323 
324  public:
326 
327  typedef typename Base::IndexType IndexType;
328 
329  static const int dimension = Grid::dimension;
330 
333 
334  private:
335  typedef typename Grid::Traits Traits;
336 
337  template< int codim >
338  struct Insert;
339 
340  public:
341  explicit AlbertaGridIndexSet ( const DofNumbering &dofNumbering )
342  : dofNumbering_( dofNumbering )
343  {
344  for( int codim = 0; codim <= dimension; ++codim )
345  {
346  indices_[ codim ] = 0;
347 
348  const GeometryType type( GeometryType::simplex, dimension - codim );
349  geomTypes_[ codim ].push_back( type );
350  }
351  }
352 
354  {
355  for( int codim = 0; codim <= dimension; ++codim )
356  delete[] indices_[ codim ];
357  }
358 
359  template< class Entity >
360  bool contains ( const Entity &entity ) const
361  {
362  const int codim = Entity::codimension;
363 
365  = Grid::getRealImplementation( entity );
366  const Alberta::Element *element = entityImp.elementInfo().el();
367 
368  const IndexType *const array = indices_[ codim ];
369  const IndexType subIndex = array[ dofNumbering_( element, codim, entityImp.subEntity() ) ];
370 
371  return (subIndex >= 0);
372  }
373 
374  using Base::index;
375  using Base::subIndex;
376 
378  template< int cc >
379  IndexType index ( const typename Traits::template Codim< cc >::Entity &entity ) const
380  {
382  const EntityImp &entityImp = Grid::getRealImplementation( entity );
383  return subIndex( entityImp.elementInfo(), entityImp.subEntity(), cc );
384  }
385 
387  template< int cc >
388  IndexType subIndex ( const typename Traits::template Codim< cc >::Entity &entity, int i, unsigned int codim ) const
389  {
391  const EntityImp &entityImp = Grid::getRealImplementation( entity );
392 
393  int k = i;
394  if( cc > 0 )
395  {
396  const GenericReferenceElement< Alberta::Real, dimension > &refElement
398  k = refElement.subEntity( entityImp.subEntity(), cc, i, codim );
399  }
400 
401  const int j = entityImp.grid().generic2alberta( codim, k );
402  return subIndex( entityImp.elementInfo(), j, codim );
403  }
404 
405  IndexType size ( const GeometryType &type ) const
406  {
407  return (type.isSimplex() ? size( dimension - type.dim() ) : 0);
408  }
409 
410  IndexType size ( int codim ) const
411  {
412  assert( (codim >= 0) && (codim <= dimension) );
413  return size_[ codim ];
414  }
415 
416  const std::vector< GeometryType > &geomTypes( int codim ) const
417  {
418  assert( (codim >= 0) && (codim <= dimension) );
419  return geomTypes_[ codim ];
420  }
421 
422  template< class Iterator >
423  void update ( const Iterator &begin, const Iterator &end )
424  {
425  for( int codim = 0; codim <= dimension; ++codim )
426  {
427  delete[] indices_[ codim ];
428 
429  const unsigned int dofSize = dofNumbering_.size( codim );
430  indices_[ codim ] = new IndexType[ dofSize ];
431  for( unsigned int i = 0; i < dofSize; ++i )
432  indices_[ codim ][ i ] = -1;
433 
434  size_[ codim ] = 0;
435  }
436 
437  for( Iterator it = begin; it != end; ++it )
438  {
441  const Alberta::Element *element = entityImp.elementInfo().el();
442  ForLoop< Insert, 0, dimension >::apply( element, *this );
443  }
444  }
445 
446  private:
447  IndexType subIndex ( const ElementInfo &elementInfo, int i, unsigned int codim ) const
448  {
449  assert( !elementInfo == false );
450  return subIndex( elementInfo.element(), i, codim );
451  }
452 
459  IndexType subIndex ( const Alberta::Element *element, int i, unsigned int codim ) const
460  {
461  const IndexType *const array = indices_[ codim ];
462  const IndexType subIndex = array[ dofNumbering_( element, codim, i ) ];
463  assert( (subIndex >= 0) && (subIndex < size( codim )) );
464  return subIndex;
465  }
466 
467  // access to the dof vectors
468  const DofNumbering &dofNumbering_;
469 
470  // an array of indices for each codimension
471  IndexType *indices_[ dimension+1 ];
472 
473  // the size of each codimension
474  IndexType size_[ dimension+1 ];
475 
476  // all geometry types contained in the grid
477  std::vector< GeometryType > geomTypes_[ dimension+1 ];
478  };
479 
480 
481 
482  // AlbertaGridIndexSet::Insert
483  // ---------------------------
484 
485  template< int dim, int dimworld >
486  template< int codim >
487  struct AlbertaGridIndexSet< dim, dimworld >::Insert
488  {
489  static void apply ( const Alberta::Element *const element,
491  {
492  int *const array = indexSet.indices_[ codim ];
493  IndexType &size = indexSet.size_[ codim ];
494 
495  for( int i = 0; i < Alberta::NumSubEntities< dim, codim >::value; ++i )
496  {
497  int &index = array[ indexSet.dofNumbering_( element, codim, i ) ];
498  if( index < 0 )
499  index = size++;
500  }
501  }
502  };
503 
504 
505 
506  // AlbertaGridIdSet
507  // ----------------
508 
510  template< int dim, int dimworld >
512  : public IdSet< AlbertaGrid< dim, dimworld >, AlbertaGridIdSet< dim, dimworld >, unsigned int >
513  {
515  typedef IdSet< AlbertaGrid< dim, dimworld >, This, unsigned int > Base;
516 
517  friend class AlbertaGrid< dim, dimworld >;
518 
519  public:
521  typedef typename Base::IdType IdType;
522 
523  private:
525 
526  static const int dimension = Grid::dimension;
527 
529 
530  // create id set, only allowed for AlbertaGrid
531  AlbertaGridIdSet ( const HierarchicIndexSet &hIndexSet )
532  : hIndexSet_( hIndexSet )
533  {}
534 
535  public:
537  template< class Entity >
538  IdType id ( const Entity &e ) const
539  {
540  const int codim = Entity::codimension;
541  return id< codim >( e );
542  }
543 
545  template< int codim >
546  IdType id ( const typename Grid::template Codim< codim >::Entity &e ) const
547  {
548  assert( (codim >= 0) && (codim <= dimension) );
549  const IdType index = hIndexSet_.index( e );
550  return (index << 2) | IdType( codim );
551  }
552 
554  IdType subId ( const typename Grid::template Codim< 0 >::Entity &e, int i, unsigned int subcodim ) const
555  {
556  assert( int( subcodim ) <= dimension );
557  const IdType index = hIndexSet_.subIndex( e, i, subcodim );
558  return (index << 2) | IdType( subcodim );
559  }
560 
561  template< int codim >
562  IdType subId ( const typename Grid::template Codim< codim >::Entity &e, int i, unsigned int subcodim ) const
563  {
564  assert( (codim >= 0) && (codim <= dimension) && (int( codim + subcodim ) <= dimension) );
565  const IdType index = hIndexSet_.subIndex( e, i, subcodim );
566  return (index << 2) | IdType( codim + subcodim );
567  }
568 
569  template< class Entity >
570  IdType subId ( const Entity &e, int i, unsigned int subcodim ) const
571  {
572  return subId< Entity::codimension >( e, i, subcodim );
573  }
574 
575  private:
576  // prohibit copying
577  AlbertaGridIdSet ( const This & );
578 
579  const HierarchicIndexSet &hIndexSet_;
580  };
581 
582 } // namespace Dune
583 
584 #endif // #if HAVE_ALBERTA
585 
586 #endif // #ifndef DUNE_ALBERTAGRIDINDEXSETS_HH