Dune Core Modules (2.5.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
9
12
14#include <dune/grid/albertagrid/misc.hh>
15#include <dune/grid/albertagrid/dofadmin.hh>
16#include <dune/grid/albertagrid/dofvector.hh>
19
20#if HAVE_ALBERTA
21
22namespace Dune
23{
24
25 namespace Alberta
26 {
27 typedef Dune::IndexStack< int, 100000 > IndexStack;
28
29 extern IndexStack *currentIndexStack;
30 }
31
32
33
34 // AlbertaGridHierarchicIndexSet
35 // -----------------------------
36
37 template< int dim, int dimworld >
38 class AlbertaGridHierarchicIndexSet
39 : public IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim,dimworld >, int, std::array< GeometryType, 1 > >
40 {
41 typedef AlbertaGridHierarchicIndexSet< dim, dimworld > This;
42 typedef IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim,dimworld >, int, std::array< GeometryType, 1 > > Base;
43
44 friend class AlbertaGrid< dim, dimworld >;
45
46 public:
48 typedef AlbertaGridFamily< dim, dimworld > GridFamily;
49
50 typedef typename Base::IndexType IndexType;
51
52 typedef typename Base::Types Types;
53
54 static const int dimension = GridFamily::dimension;
55
56 typedef Alberta::ElementInfo< dimension > ElementInfo;
57 typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering;
58
59 private:
60 typedef typename GridFamily::Traits Traits;
61
62 typedef Alberta::DofVectorPointer< IndexType > IndexVectorPointer;
63
64 class InitEntityNumber;
65
66 template< int codim >
67 struct CreateEntityNumbers;
68
69 template< int codim >
70 struct RefineNumbering;
71
72 template< int codim >
73 struct CoarsenNumbering;
74
75 explicit AlbertaGridHierarchicIndexSet ( const DofNumbering &dofNumbering );
76
77 public:
78 typedef Alberta::IndexStack IndexStack;
79
81 template< class Entity >
82 bool contains ( const Entity & ) const
83 {
84 return true;
85 }
86
87 using Base::index;
88 using Base::subIndex;
89
91 template< int cc >
92 IndexType index ( const typename Traits::template Codim< cc >::Entity &entity ) const
93 {
94 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
95 const EntityImp &entityImp = Grid::getRealImplementation( entity );
96 return subIndex( entityImp.elementInfo(), entityImp.subEntity(), cc );
97 }
98
100 template< int cc >
101 IndexType subIndex ( const typename Traits::template Codim< cc >::Entity &entity, int i, unsigned int codim ) const
102 {
103 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
104 const EntityImp &entityImp = Grid::getRealImplementation( entity );
105
106 int k = i;
107 if( cc > 0 )
108 {
109 const ReferenceElement< Alberta::Real, dimension > &refElement
111 k = refElement.subEntity( entityImp.subEntity(), cc, i, codim );
112 }
113
114 const int j = entityImp.grid().generic2alberta( codim, k );
115 return subIndex( entityImp.elementInfo(), j, codim );
116 }
117
119 IndexType size ( const GeometryType &type ) const
120 {
121 return (type.isSimplex() ? size( dimension - type.dim() ) : 0);
122 }
123
125 IndexType size ( int codim ) const
126 {
127 assert( (codim >= 0) && (codim <= dimension) );
128 return indexStack_[ codim ].size();
129 }
130
131 Types types ( int codim ) const
132 {
133 assert( (codim >= 0) && (codim <= dimension) );
134 std::array< GeometryType, 1 > types;
135 types[ 0 ] = GeometryType( GeometryType::simplex, dimension - codim );
136 return types;
137 }
138
140 const std::vector< GeometryType > &geomTypes( int codim ) const
141 {
142 assert( (codim >= 0) && (codim <= dimension) );
143 return geomTypes_[ codim ];
144 }
145
146 IndexType subIndex ( const ElementInfo &elementInfo, int i, unsigned int codim ) const
147 {
148 assert( !elementInfo == false );
149 return subIndex( elementInfo.element(), i, codim );
150 }
151
158 IndexType subIndex ( const Alberta::Element *element, int i, unsigned int codim ) const
159 {
160 IndexType *array = (IndexType *)entityNumbers_[ codim ];
161 const IndexType subIndex = array[ dofNumbering_( element, codim, i ) ];
162 assert( (subIndex >= 0) && (subIndex < size( codim )) );
163 return subIndex;
164 }
165
166 void preAdapt ()
167 {
168 // set global pointer to index stack
169 if( !IndexVectorPointer::supportsAdaptationData )
170 {
171 assert( Alberta::currentIndexStack == 0 );
172 Alberta::currentIndexStack = indexStack_;
173 }
174 }
175
176 void postAdapt ()
177 {
178 // remove global pointer to index stack
179 if( !IndexVectorPointer::supportsAdaptationData )
180 Alberta::currentIndexStack = 0;
181 }
182
183 void create ();
184 void read ( const std::string &filename );
185 bool write ( const std::string &filename ) const;
186
187 void release ()
188 {
189 for( int i = 0; i <= dimension; ++i )
190 entityNumbers_[ i ].release();
191 }
192
193 private:
194 template< int codim >
195 static IndexStack &getIndexStack ( const IndexVectorPointer &dofVector )
196 {
197 IndexStack *indexStack;
198 if( IndexVectorPointer::supportsAdaptationData )
199 indexStack = dofVector.template getAdaptationData< IndexStack >();
200 else
201 indexStack = &Alberta::currentIndexStack[ codim ];
202 assert( indexStack != 0 );
203 return *indexStack;
204 }
205
206 // access to the dof vectors
207 const DofNumbering &dofNumbering_;
208
209 // index stacks providing new numbers during adaptation
210 IndexStack indexStack_[ dimension+1 ];
211
212 // dof vectors storing the (persistent) numbering
213 IndexVectorPointer entityNumbers_[ dimension+1 ];
214
215 // all geometry types contained in the grid
216 std::vector< GeometryType > geomTypes_[ dimension+1 ];
217 };
218
219
220
221 // AlbertaGridHierarchicIndexSet::InitEntityNumber
222 // -----------------------------------------------
223
224 template< int dim, int dimworld >
225 class AlbertaGridHierarchicIndexSet< dim, dimworld >::InitEntityNumber
226 {
227 IndexStack &indexStack_;
228
229 public:
230 InitEntityNumber ( IndexStack &indexStack )
231 : indexStack_( indexStack )
232 {}
233
234 void operator() ( int &dof )
235 {
236 dof = indexStack_.getIndex();
237 }
238 };
239
240
241
242 // AlbertaGridHierarchicIndexSet::CreateEntityNumbers
243 // --------------------------------------------------
244
245 template< int dim, int dimworld >
246 template< int codim >
247 struct AlbertaGridHierarchicIndexSet< dim, dimworld >::CreateEntityNumbers
248 {
249 static void setup ( AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
250
251 static void apply ( const Alberta::HierarchyDofNumbering< dimension > &dofNumbering,
252 AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
253
254 static void apply ( const std::string &filename,
255 const Alberta::MeshPointer< dimension > &mesh,
256 AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
257 };
258
259
260
261 // AlbertaGridHierarchicIndexSet::RefineNumbering
262 // ----------------------------------------------
263
264 template< int dim, int dimworld >
265 template< int codim >
266 struct AlbertaGridHierarchicIndexSet< dim, dimworld >::RefineNumbering
267 {
268 static const int dimension = dim;
269 static const int codimension = codim;
270
271 private:
272 typedef Alberta::DofAccess< dimension, codimension > DofAccess;
273
274 explicit RefineNumbering ( const IndexVectorPointer &dofVector )
275 : indexStack_( getIndexStack< codimension >( dofVector ) ),
276 dofVector_( dofVector ),
277 dofAccess_( dofVector.dofSpace() )
278 {}
279
280 public:
281 void operator() ( const Alberta::Element *child, int subEntity );
282
283 typedef Alberta::Patch< dimension > Patch;
284 static void interpolateVector ( const IndexVectorPointer &dofVector,
285 const Patch &patch );
286
287 private:
288 IndexStack &indexStack_;
289 IndexVectorPointer dofVector_;
290 DofAccess dofAccess_;
291 };
292
293
294
295 // AlbertaGridHierarchicIndexSet::CoarsenNumbering
296 // -----------------------------------------------
297
298 template< int dim, int dimworld >
299 template< int codim >
300 struct AlbertaGridHierarchicIndexSet< dim, dimworld >::CoarsenNumbering
301 {
302 static const int dimension = dim;
303 static const int codimension = codim;
304
305 private:
306 typedef Alberta::DofAccess< dimension, codimension > DofAccess;
307
308 explicit CoarsenNumbering ( const IndexVectorPointer &dofVector )
309 : indexStack_( getIndexStack< codimension >( dofVector ) ),
310 dofVector_( dofVector ),
311 dofAccess_( dofVector.dofSpace() )
312 {}
313
314 public:
315 void operator() ( const Alberta::Element *child, int subEntity );
316
317 typedef Alberta::Patch< dimension > Patch;
318 static void restrictVector ( const IndexVectorPointer &dofVector,
319 const Patch &patch );
320 private:
321 IndexStack &indexStack_;
322 IndexVectorPointer dofVector_;
323 DofAccess dofAccess_;
324 };
325
326
327
328 // AlbertaGridIndexSet
329 // -------------------
330
331 template< int dim, int dimworld >
332 class AlbertaGridIndexSet
333 : public IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > >
334 {
335 typedef AlbertaGridIndexSet< dim, dimworld > This;
336 typedef IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > > Base;
337
338 public:
339 typedef AlbertaGrid< dim, dimworld > Grid;
340
341 typedef typename Base::IndexType IndexType;
342
343 typedef typename Base::Types Types;
344
345 static const int dimension = Grid::dimension;
346
347 typedef Alberta::ElementInfo< dimension > ElementInfo;
348 typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering;
349
350 private:
351 typedef typename Grid::Traits Traits;
352
353 template< int codim >
354 struct Insert;
355
356 public:
357 explicit AlbertaGridIndexSet ( const DofNumbering &dofNumbering )
358 : dofNumbering_( dofNumbering )
359 {
360 for( int codim = 0; codim <= dimension; ++codim )
361 {
362 indices_[ codim ] = 0;
363
364 const GeometryType type( GeometryType::simplex, dimension - codim );
365 geomTypes_[ codim ].push_back( type );
366 }
367 }
368
369 ~AlbertaGridIndexSet ()
370 {
371 for( int codim = 0; codim <= dimension; ++codim )
372 delete[] indices_[ codim ];
373 }
374
375 template< class Entity >
376 bool contains ( const Entity &entity ) const
377 {
378 const int codim = Entity::codimension;
379
380 const AlbertaGridEntity< codim, dim, const Grid > &entityImp
381 = Grid::getRealImplementation( entity );
382 const Alberta::Element *element = entityImp.elementInfo().el();
383
384 const IndexType *const array = indices_[ codim ];
385 const IndexType subIndex = array[ dofNumbering_( element, codim, entityImp.subEntity() ) ];
386
387 return (subIndex >= 0);
388 }
389
390 using Base::index;
391 using Base::subIndex;
392
394 template< int cc >
395 IndexType index ( const typename Traits::template Codim< cc >::Entity &entity ) const
396 {
397 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
398 const EntityImp &entityImp = Grid::getRealImplementation( entity );
399 return subIndex( entityImp.elementInfo(), entityImp.subEntity(), cc );
400 }
401
403 template< int cc >
404 IndexType subIndex ( const typename Traits::template Codim< cc >::Entity &entity, int i, unsigned int codim ) const
405 {
406 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
407 const EntityImp &entityImp = Grid::getRealImplementation( entity );
408
409 int k = i;
410 if( cc > 0 )
411 {
412 const ReferenceElement< Alberta::Real, dimension > &refElement
414 k = refElement.subEntity( entityImp.subEntity(), cc, i, codim );
415 }
416
417 const int j = entityImp.grid().generic2alberta( codim, k );
418 return subIndex( entityImp.elementInfo(), j, codim );
419 }
420
421 IndexType size ( const GeometryType &type ) const
422 {
423 return (type.isSimplex() ? size( dimension - type.dim() ) : 0);
424 }
425
426 IndexType size ( int codim ) const
427 {
428 assert( (codim >= 0) && (codim <= dimension) );
429 return size_[ codim ];
430 }
431
432 Types types ( int codim ) const
433 {
434 assert( (codim >= 0) && (codim <= dimension) );
435 std::array< GeometryType, 1 > types;
436 types[ 0 ] = GeometryType( GeometryType::simplex, dimension - codim );
437 return types;
438 }
439
440 const std::vector< GeometryType > &geomTypes( int codim ) const
441 {
442 assert( (codim >= 0) && (codim <= dimension) );
443 return geomTypes_[ codim ];
444 }
445
446 template< class Iterator >
447 void update ( const Iterator &begin, const Iterator &end )
448 {
449 for( int codim = 0; codim <= dimension; ++codim )
450 {
451 delete[] indices_[ codim ];
452
453 const unsigned int dofSize = dofNumbering_.size( codim );
454 indices_[ codim ] = new IndexType[ dofSize ];
455 for( unsigned int i = 0; i < dofSize; ++i )
456 indices_[ codim ][ i ] = -1;
457
458 size_[ codim ] = 0;
459 }
460
461 for( Iterator it = begin; it != end; ++it )
462 {
463 const AlbertaGridEntity< 0, dim, const Grid > &entityImp
465 const Alberta::Element *element = entityImp.elementInfo().el();
466 ForLoop< Insert, 0, dimension >::apply( element, *this );
467 }
468 }
469
470 private:
471 IndexType subIndex ( const ElementInfo &elementInfo, int i, unsigned int codim ) const
472 {
473 assert( !elementInfo == false );
474 return subIndex( elementInfo.element(), i, codim );
475 }
476
483 IndexType subIndex ( const Alberta::Element *element, int i, unsigned int codim ) const
484 {
485 const IndexType *const array = indices_[ codim ];
486 const IndexType subIndex = array[ dofNumbering_( element, codim, i ) ];
487 assert( (subIndex >= 0) && (subIndex < size( codim )) );
488 return subIndex;
489 }
490
491 // access to the dof vectors
492 const DofNumbering &dofNumbering_;
493
494 // an array of indices for each codimension
495 IndexType *indices_[ dimension+1 ];
496
497 // the size of each codimension
498 IndexType size_[ dimension+1 ];
499
500 // all geometry types contained in the grid
501 std::vector< GeometryType > geomTypes_[ dimension+1 ];
502 };
503
504
505
506 // AlbertaGridIndexSet::Insert
507 // ---------------------------
508
509 template< int dim, int dimworld >
510 template< int codim >
511 struct AlbertaGridIndexSet< dim, dimworld >::Insert
512 {
513 static void apply ( const Alberta::Element *const element,
514 AlbertaGridIndexSet< dim, dimworld > &indexSet )
515 {
516 int *const array = indexSet.indices_[ codim ];
517 IndexType &size = indexSet.size_[ codim ];
518
519 for( int i = 0; i < Alberta::NumSubEntities< dim, codim >::value; ++i )
520 {
521 int &index = array[ indexSet.dofNumbering_( element, codim, i ) ];
522 if( index < 0 )
523 index = size++;
524 }
525 }
526 };
527
528
529
530 // AlbertaGridIdSet
531 // ----------------
532
534 template< int dim, int dimworld >
536 : public IdSet< AlbertaGrid< dim, dimworld >, AlbertaGridIdSet< dim, dimworld >, unsigned int >
537 {
539 typedef IdSet< AlbertaGrid< dim, dimworld >, This, unsigned int > Base;
540
541 friend class AlbertaGrid< dim, dimworld >;
542
543 public:
545 typedef typename Base::IdType IdType;
546
547 private:
549
550 static const int dimension = Grid::dimension;
551
552 typedef typename Grid::HierarchicIndexSet HierarchicIndexSet;
553
554 // create id set, only allowed for AlbertaGrid
555 AlbertaGridIdSet ( const HierarchicIndexSet &hIndexSet )
556 : hIndexSet_( hIndexSet )
557 {}
558
559 public:
561 template< class Entity >
562 IdType id ( const Entity &e ) const
563 {
564 const int codim = Entity::codimension;
565 return id< codim >( e );
566 }
567
569 template< int codim >
570 IdType id ( const typename Grid::template Codim< codim >::Entity &e ) const
571 {
572 assert( (codim >= 0) && (codim <= dimension) );
573 const IdType index = hIndexSet_.index( e );
574 return (index << 2) | IdType( codim );
575 }
576
578 IdType subId ( const typename Grid::template Codim< 0 >::Entity &e, int i, unsigned int subcodim ) const
579 {
580 assert( int( subcodim ) <= dimension );
581 const IdType index = hIndexSet_.subIndex( e, i, subcodim );
582 return (index << 2) | IdType( subcodim );
583 }
584
585 template< int codim >
586 IdType subId ( const typename Grid::template Codim< codim >::Entity &e, int i, unsigned int subcodim ) const
587 {
588 assert( (codim >= 0) && (codim <= dimension) && (int( codim + subcodim ) <= dimension) );
589 const IdType index = hIndexSet_.subIndex( e, i, subcodim );
590 return (index << 2) | IdType( codim + subcodim );
591 }
592
593 template< class Entity >
594 IdType subId ( const Entity &e, int i, unsigned int subcodim ) const
595 {
596 return subId< Entity::codimension >( e, i, subcodim );
597 }
598
599 private:
600 // prohibit copying
601 AlbertaGridIdSet ( const This & );
602
603 const HierarchicIndexSet &hIndexSet_;
604 };
605
606} // namespace Dune
607
608#endif // #if HAVE_ALBERTA
609
610#endif // #ifndef DUNE_ALBERTAGRIDINDEXSETS_HH
provides the GridFamily for AlbertaGrid
hierarchic index set of AlbertaGrid
Definition: indexsets.hh:537
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 explicitely.
Definition: indexsets.hh:570
IdType id(const Entity &e) const
Definition: indexsets.hh:562
Base::IdType IdType
export type of id
Definition: indexsets.hh:545
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:578
Traits::HierarchicIndexSet HierarchicIndexSet
type of hierarchic index set
Definition: agrid.hh:189
Wrapper class for entities.
Definition: entity.hh:65
@ codimension
Know your own codimension.
Definition: entity.hh:107
@ simplex
Simplicial element in any nonnegative dimension.
Definition: type.hh:273
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:1115
Id Set Interface.
Definition: indexidset.hh:402
IdTypeImp IdType
Type used to represent an id.
Definition: indexidset.hh:405
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.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:178
provides a wrapper for ALBERTA's el_info structure
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: alignment.hh:11
Standard Dune debug streams.
Static tag representing a codimension.
Definition: dimension.hh:22
static const ReferenceElement< ctype, dim > & simplex()
get simplex reference elements
Definition: referenceelements.hh:763
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)