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