Dune Core Modules (2.9.0)

indexsets.hh
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>
13
16
18#include <dune/grid/albertagrid/misc.hh>
19#include <dune/grid/albertagrid/dofadmin.hh>
20#include <dune/grid/albertagrid/dofvector.hh>
23
24#if HAVE_ALBERTA
25
26namespace Dune
27{
28
29 namespace Alberta
30 {
31 typedef Dune::IndexStack< int, 100000 > IndexStack;
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 static Alberta::IndexStack *currentIndexStack;
80
81 public:
82 typedef Alberta::IndexStack IndexStack;
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 {
98 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
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 {
107 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
108 const EntityImp &entityImp = entity.impl();
109
110 int k = i;
111 if( cc > 0 )
112 {
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
170 if( !IndexVectorPointer::supportsAdaptationData )
171 {
172 assert( currentIndexStack == nullptr );
173 currentIndexStack = indexStack_;
174 }
175 }
176
177 void postAdapt ()
178 {
179 // remove global pointer to index stack
180 if( !IndexVectorPointer::supportsAdaptationData )
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;
199 if( IndexVectorPointer::supportsAdaptationData )
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 >
226 class AlbertaGridHierarchicIndexSet< dim, dimworld >::InitEntityNumber
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,
253 AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
254
255 static void apply ( const std::string &filename,
256 const Alberta::MeshPointer< dimension > &mesh,
257 AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
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 >
333 class AlbertaGridIndexSet
334 : public IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > >
335 {
336 typedef AlbertaGridIndexSet< dim, dimworld > This;
337 typedef IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > > Base;
338
339 public:
340 typedef AlbertaGrid< dim, dimworld > Grid;
341
342 typedef typename Base::IndexType IndexType;
343
344 typedef typename Base::Types Types;
345
346 static const int dimension = Grid::dimension;
347
348 typedef Alberta::ElementInfo< dimension > ElementInfo;
349 typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering;
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
368 ~AlbertaGridIndexSet ()
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
379 const AlbertaGridEntity< codim, dim, const Grid > &entityImp
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 {
396 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
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 {
405 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
406 const EntityImp &entityImp = entity.impl();
407
408 int k = i;
409 if( cc > 0 )
410 {
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 {
459 const AlbertaGridEntity< 0, dim, const Grid > &entityImp
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 {
536 typedef IdSet< AlbertaGrid< dim, dimworld >, This, unsigned int > Base;
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
549 typedef typename Grid::HierarchicIndexSet HierarchicIndexSet;
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 the GridFamily for AlbertaGrid
hierarchic index set of AlbertaGrid
Definition: indexsets.hh:534
IdType id(const typename Grid::template Codim< codim >::Entity &e) const
Definition: indexsets.hh:567
IdType id(const Entity &e) const
Definition: indexsets.hh:559
Base::IdType IdType
export type of id
Definition: indexsets.hh:542
IdType subId(const typename Grid::template Codim< 0 >::Entity &e, int i, unsigned int subcodim) const
Definition: indexsets.hh:575
Traits::HierarchicIndexSet HierarchicIndexSet
type of hierarchic index set
Definition: agrid.hh:157
Wrapper class for entities.
Definition: entity.hh:66
static constexpr int codimension
Know your own codimension.
Definition: entity.hh:106
Id Set Interface.
Definition: indexidset.hh:452
IdTypeImp IdType
Type used to represent an id.
Definition: indexidset.hh:458
auto size(GeometryType type) const
Return total number of entities of given geometry type in entity set .
Definition: indexidset.hh:223
IndexType subIndex(const Entity &e, int i, unsigned int codim) const
Map a subentity to an index.
Definition: indexidset.hh:182
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
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:129
Definition: indexstack.hh:26
Different resources needed by all grid implementations.
provides a wrapper for ALBERTA's el_info structure
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:463
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:268
ImplementationDefined child(Node &&node, Indices... indices)
Extracts the child of a node given by a sequence of compile-time and run-time indices.
Definition: childextraction.hh:126
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:13
Standard Dune debug streams.
Static tag representing a codimension.
Definition: dimension.hh:24
static const ReferenceElement & simplex()
get simplex reference elements
Definition: referenceelements.hh:204
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)