Dune Core Modules (2.6.0)

indexsets.hh
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>
10#include <dune/common/std/utility.hh>
11
14
16#include <dune/grid/albertagrid/misc.hh>
17#include <dune/grid/albertagrid/dofadmin.hh>
18#include <dune/grid/albertagrid/dofvector.hh>
21
22#if HAVE_ALBERTA
23
24namespace Dune
25{
26
27 namespace Alberta
28 {
29 typedef Dune::IndexStack< int, 100000 > IndexStack;
30
31 extern IndexStack *currentIndexStack;
32 }
33
34
35
36 // AlbertaGridHierarchicIndexSet
37 // -----------------------------
38
39 template< int dim, int dimworld >
40 class AlbertaGridHierarchicIndexSet
41 : public IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim,dimworld >, int, std::array< GeometryType, 1 > >
42 {
43 typedef AlbertaGridHierarchicIndexSet< dim, dimworld > This;
44 typedef IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim,dimworld >, int, std::array< GeometryType, 1 > > Base;
45
46 friend class AlbertaGrid< dim, dimworld >;
47
48 public:
50 typedef AlbertaGridFamily< dim, dimworld > GridFamily;
51
52 typedef typename Base::IndexType IndexType;
53
54 typedef typename Base::Types Types;
55
56 static const int dimension = GridFamily::dimension;
57
58 typedef Alberta::ElementInfo< dimension > ElementInfo;
59 typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering;
60
61 private:
62 typedef typename GridFamily::Traits Traits;
63
64 typedef Alberta::DofVectorPointer< IndexType > IndexVectorPointer;
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:
80 typedef Alberta::IndexStack IndexStack;
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 {
96 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
97 const EntityImp &entityImp = Grid::getRealImplementation( entity );
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 {
105 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
106 const EntityImp &entityImp = Grid::getRealImplementation( entity );
107
108 int k = i;
109 if( cc > 0 )
110 {
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
168 if( !IndexVectorPointer::supportsAdaptationData )
169 {
170 assert( Alberta::currentIndexStack == 0 );
171 Alberta::currentIndexStack = indexStack_;
172 }
173 }
174
175 void postAdapt ()
176 {
177 // remove global pointer to index stack
178 if( !IndexVectorPointer::supportsAdaptationData )
179 Alberta::currentIndexStack = 0;
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;
197 if( IndexVectorPointer::supportsAdaptationData )
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 >
224 class AlbertaGridHierarchicIndexSet< dim, dimworld >::InitEntityNumber
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,
251 AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
252
253 static void apply ( const std::string &filename,
254 const Alberta::MeshPointer< dimension > &mesh,
255 AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
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 >
331 class AlbertaGridIndexSet
332 : public IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > >
333 {
334 typedef AlbertaGridIndexSet< dim, dimworld > This;
335 typedef IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > > Base;
336
337 public:
338 typedef AlbertaGrid< dim, dimworld > Grid;
339
340 typedef typename Base::IndexType IndexType;
341
342 typedef typename Base::Types Types;
343
344 static const int dimension = Grid::dimension;
345
346 typedef Alberta::ElementInfo< dimension > ElementInfo;
347 typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering;
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
366 ~AlbertaGridIndexSet ()
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
377 const AlbertaGridEntity< codim, dim, const Grid > &entityImp
378 = Grid::getRealImplementation( entity );
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 {
394 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
395 const EntityImp &entityImp = Grid::getRealImplementation( entity );
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 {
403 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
404 const EntityImp &entityImp = Grid::getRealImplementation( entity );
405
406 int k = i;
407 if( cc > 0 )
408 {
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 {
457 const AlbertaGridEntity< 0, dim, const Grid > &entityImp
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 {
534 typedef IdSet< AlbertaGrid< dim, dimworld >, This, unsigned int > Base;
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
547 typedef typename Grid::HierarchicIndexSet HierarchicIndexSet;
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 the GridFamily for AlbertaGrid
hierarchic index set of AlbertaGrid
Definition: 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: indexsets.hh:565
IdType id(const Entity &e) const
Definition: indexsets.hh:557
Base::IdType IdType
export type of id
Definition: 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: indexsets.hh:573
Traits::HierarchicIndexSet HierarchicIndexSet
type of hierarchic index set
Definition: agrid.hh:187
Wrapper class for entities.
Definition: entity.hh:64
@ codimension
Know your own codimension.
Definition: entity.hh:105
static std::conditional< std::is_reference< InterfaceType >::value, typenamestd::add_lvalue_reference< typenameReturnImplementationType< typenamestd::remove_reference< InterfaceType >::type >::ImplementationType >::type, typenamestd::remove_const< typenameReturnImplementationType< typenamestd::remove_reference< InterfaceType >::type >::ImplementationType >::type >::type getRealImplementation(InterfaceType &&i)
return real implementation of interface class
Definition: grid.hh:1026
Id Set Interface.
Definition: indexidset.hh:441
IdTypeImp IdType
Type used to represent an id.
Definition: indexidset.hh:444
IndexType subIndex(const Entity &e, int i, unsigned int codim) const
Map a subentity to an index.
Definition: indexidset.hh:180
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
IndexType size(GeometryType type) const
Return total number of entities of given geometry type in entity set .
Definition: indexidset.hh:220
IndexType index(const Entity &e) const
Map entity to index. Easier to use than the above because codimension template parameter need not be ...
Definition: indexidset.hh:127
Definition: indexstack.hh:24
Different resources needed by all grid implementations.
provides a wrapper for ALBERTA's el_info structure
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition: apply.hh:58
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:696
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:308
Provides base classes for index and id sets.
Provides an index stack that supplies indices for element numbering for a grid (i....
Dune namespace.
Definition: alignedallocator.hh:10
Standard Dune debug streams.
Static tag representing a codimension.
Definition: dimension.hh:22
static const ReferenceElement & simplex()
get simplex reference elements
Definition: referenceelements.hh:202
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Oct 13, 22:30, 2024)