DUNE-FEM (unstable)

indexsets.hh
1// SPDX-FileCopyrightText: Copyright © 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>
15#include <dune/grid/common/indexidset.hh>
16
17#include <dune/grid/albertagrid/indexstack.hh>
18#include <dune/grid/albertagrid/misc.hh>
19#include <dune/grid/albertagrid/dofadmin.hh>
20#include <dune/grid/albertagrid/dofvector.hh>
21#include <dune/grid/albertagrid/elementinfo.hh>
22#include <dune/grid/albertagrid/gridfamily.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:
49 typedef AlbertaGrid< dim, dimworld > Grid;
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 {
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
135 Types types ( int codim ) const
136 {
137 assert( (codim >= 0) && (codim <= dimension) );
138 return {{ GeometryTypes::simplex( dimension - codim ) }};
139 }
140
141 IndexType subIndex ( const ElementInfo &elementInfo, int i, unsigned int codim ) const
142 {
143 assert( !elementInfo == false );
144 return subIndex( elementInfo.element(), i, codim );
145 }
146
153 IndexType subIndex ( const Alberta::Element *element, int i, unsigned int codim ) const
154 {
155 IndexType *array = (IndexType *)entityNumbers_[ codim ];
156 const IndexType subIndex = array[ dofNumbering_( element, codim, i ) ];
157 assert( (subIndex >= 0) && (subIndex < IndexType(size( codim ))) );
158 return subIndex;
159 }
160
161 void preAdapt ()
162 {
163 // set global pointer to index stack
164 if( !IndexVectorPointer::supportsAdaptationData )
165 {
166 assert( currentIndexStack == nullptr );
167 currentIndexStack = indexStack_;
168 }
169 }
170
171 void postAdapt ()
172 {
173 // remove global pointer to index stack
174 if( !IndexVectorPointer::supportsAdaptationData )
175 currentIndexStack = nullptr;
176 }
177
178 void create ();
179 void read ( const std::string &filename );
180 bool write ( const std::string &filename ) const;
181
182 void release ()
183 {
184 for( int i = 0; i <= dimension; ++i )
185 entityNumbers_[ i ].release();
186 }
187
188 private:
189 template< int codim >
190 static IndexStack &getIndexStack ( const IndexVectorPointer &dofVector )
191 {
192 IndexStack *indexStack;
193 if( IndexVectorPointer::supportsAdaptationData )
194 indexStack = dofVector.template getAdaptationData< IndexStack >();
195 else
196 indexStack = &currentIndexStack[ codim ];
197 assert( indexStack != 0 );
198 return *indexStack;
199 }
200
201 // access to the dof vectors
202 const DofNumbering &dofNumbering_;
203
204 // index stacks providing new numbers during adaptation
205 IndexStack indexStack_[ dimension+1 ];
206
207 // dof vectors storing the (persistent) numbering
208 IndexVectorPointer entityNumbers_[ dimension+1 ];
209 };
210
211 template< int dim, int dimworld >
212 Alberta::IndexStack* AlbertaGridHierarchicIndexSet<dim,dimworld>::currentIndexStack = nullptr;
213
214
215
216 // AlbertaGridHierarchicIndexSet::InitEntityNumber
217 // -----------------------------------------------
218
219 template< int dim, int dimworld >
220 class AlbertaGridHierarchicIndexSet< dim, dimworld >::InitEntityNumber
221 {
222 IndexStack &indexStack_;
223
224 public:
225 InitEntityNumber ( IndexStack &indexStack )
226 : indexStack_( indexStack )
227 {}
228
229 void operator() ( int &dof )
230 {
231 dof = indexStack_.getIndex();
232 }
233 };
234
235
236
237 // AlbertaGridHierarchicIndexSet::CreateEntityNumbers
238 // --------------------------------------------------
239
240 template< int dim, int dimworld >
241 template< int codim >
242 struct AlbertaGridHierarchicIndexSet< dim, dimworld >::CreateEntityNumbers
243 {
244 static void setup ( AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
245
246 static void apply ( const Alberta::HierarchyDofNumbering< dimension > &dofNumbering,
247 AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
248
249 static void apply ( const std::string &filename,
250 const Alberta::MeshPointer< dimension > &mesh,
251 AlbertaGridHierarchicIndexSet< dim, dimworld > &indexSet );
252 };
253
254
255
256 // AlbertaGridHierarchicIndexSet::RefineNumbering
257 // ----------------------------------------------
258
259 template< int dim, int dimworld >
260 template< int codim >
261 struct AlbertaGridHierarchicIndexSet< dim, dimworld >::RefineNumbering
262 {
263 static const int dimension = dim;
264 static const int codimension = codim;
265
266 private:
267 typedef Alberta::DofAccess< dimension, codimension > DofAccess;
268
269 explicit RefineNumbering ( const IndexVectorPointer &dofVector )
270 : indexStack_( getIndexStack< codimension >( dofVector ) ),
271 dofVector_( dofVector ),
272 dofAccess_( dofVector.dofSpace() )
273 {}
274
275 public:
276 void operator() ( const Alberta::Element *child, int subEntity );
277
278 typedef Alberta::Patch< dimension > Patch;
279 static void interpolateVector ( const IndexVectorPointer &dofVector,
280 const Patch &patch );
281
282 private:
283 IndexStack &indexStack_;
284 IndexVectorPointer dofVector_;
285 DofAccess dofAccess_;
286 };
287
288
289
290 // AlbertaGridHierarchicIndexSet::CoarsenNumbering
291 // -----------------------------------------------
292
293 template< int dim, int dimworld >
294 template< int codim >
295 struct AlbertaGridHierarchicIndexSet< dim, dimworld >::CoarsenNumbering
296 {
297 static const int dimension = dim;
298 static const int codimension = codim;
299
300 private:
301 typedef Alberta::DofAccess< dimension, codimension > DofAccess;
302
303 explicit CoarsenNumbering ( const IndexVectorPointer &dofVector )
304 : indexStack_( getIndexStack< codimension >( dofVector ) ),
305 dofVector_( dofVector ),
306 dofAccess_( dofVector.dofSpace() )
307 {}
308
309 public:
310 void operator() ( const Alberta::Element *child, int subEntity );
311
312 typedef Alberta::Patch< dimension > Patch;
313 static void restrictVector ( const IndexVectorPointer &dofVector,
314 const Patch &patch );
315 private:
316 IndexStack &indexStack_;
317 IndexVectorPointer dofVector_;
318 DofAccess dofAccess_;
319 };
320
321
322
323 // AlbertaGridIndexSet
324 // -------------------
325
326 template< int dim, int dimworld >
327 class AlbertaGridIndexSet
328 : public IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > >
329 {
330 typedef AlbertaGridIndexSet< dim, dimworld > This;
331 typedef IndexSet< AlbertaGrid< dim, dimworld >, AlbertaGridIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > > Base;
332
333 public:
334 typedef AlbertaGrid< dim, dimworld > Grid;
335
336 typedef typename Base::IndexType IndexType;
337
338 typedef typename Base::Types Types;
339
340 static const int dimension = Grid::dimension;
341
342 typedef Alberta::ElementInfo< dimension > ElementInfo;
343 typedef Alberta::HierarchyDofNumbering< dimension > DofNumbering;
344
345 private:
346 typedef typename Grid::Traits Traits;
347
348 template< int codim >
349 struct Insert;
350
351 public:
352 explicit AlbertaGridIndexSet ( const DofNumbering &dofNumbering )
353 : dofNumbering_( dofNumbering )
354 {
355 for( int codim = 0; codim <= dimension; ++codim )
356 {
357 indices_[ codim ] = 0;
358 }
359 }
360
361 ~AlbertaGridIndexSet ()
362 {
363 for( int codim = 0; codim <= dimension; ++codim )
364 delete[] indices_[ codim ];
365 }
366
367 template< class Entity >
368 bool contains ( const Entity &entity ) const
369 {
370 const int codim = Entity::codimension;
371
372 const AlbertaGridEntity< codim, dim, const Grid > &entityImp
373 = entity.impl();
374 const Alberta::Element *element = entityImp.elementInfo().el();
375
376 const IndexType *const array = indices_[ codim ];
377 const IndexType subIndex = array[ dofNumbering_( element, codim, entityImp.subEntity() ) ];
378
379 return (subIndex >= 0);
380 }
381
382 using Base::index;
383 using Base::subIndex;
384
386 template< int cc >
387 IndexType index ( const typename Traits::template Codim< cc >::Entity &entity ) const
388 {
389 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
390 const EntityImp &entityImp = entity.impl();
391 return subIndex( entityImp.elementInfo(), entityImp.subEntity(), cc );
392 }
393
395 template< int cc >
396 IndexType subIndex ( const typename Traits::template Codim< cc >::Entity &entity, int i, unsigned int codim ) const
397 {
398 typedef AlbertaGridEntity< cc, dim, const Grid > EntityImp;
399 const EntityImp &entityImp = entity.impl();
400
401 int k = i;
402 if( cc > 0 )
403 {
404 auto refElement = ReferenceElements< Alberta::Real, dimension >::simplex();
405 k = refElement.subEntity( entityImp.subEntity(), cc, i, codim );
406 }
407
408 const int j = entityImp.grid().generic2alberta( codim, k );
409 return subIndex( entityImp.elementInfo(), j, codim );
410 }
411
412 std::size_t size ( const GeometryType &type ) const
413 {
414 return (type.isSimplex() ? size( dimension - type.dim() ) : 0);
415 }
416
417 std::size_t size ( int codim ) const
418 {
419 assert( (codim >= 0) && (codim <= dimension) );
420 return size_[ codim ];
421 }
422
423 Types types ( int codim ) const
424 {
425 assert( (codim >= 0) && (codim <= dimension) );
426 return {{ GeometryTypes::simplex( dimension - codim ) }};
427 }
428
429 template< class Iterator >
430 void update ( const Iterator &begin, const Iterator &end )
431 {
432 for( int codim = 0; codim <= dimension; ++codim )
433 {
434 delete[] indices_[ codim ];
435
436 const unsigned int dofSize = dofNumbering_.size( codim );
437 indices_[ codim ] = new IndexType[ dofSize ];
438 for( unsigned int i = 0; i < dofSize; ++i )
439 indices_[ codim ][ i ] = -1;
440
441 size_[ codim ] = 0;
442 }
443
444 for( Iterator it = begin; it != end; ++it )
445 {
446 const AlbertaGridEntity< 0, dim, const Grid > &entityImp
447 = it->impl();
448 const Alberta::Element *element = entityImp.elementInfo().el();
449 Hybrid::forEach( std::make_index_sequence< dimension+1 >{},
450 [ & ]( auto i ){ Insert< i >::apply( element, *this ); } );
451 }
452 }
453
454 private:
455 IndexType subIndex ( const ElementInfo &elementInfo, int i, unsigned int codim ) const
456 {
457 assert( !elementInfo == false );
458 return subIndex( elementInfo.element(), i, codim );
459 }
460
467 IndexType subIndex ( const Alberta::Element *element, int i, unsigned int codim ) const
468 {
469 const IndexType *const array = indices_[ codim ];
470 const IndexType subIndex = array[ dofNumbering_( element, codim, i ) ];
471 assert( (subIndex >= 0) && (static_cast<unsigned int>(subIndex) < size( codim )) );
472 return subIndex;
473 }
474
475 // access to the dof vectors
476 const DofNumbering &dofNumbering_;
477
478 // an array of indices for each codimension
479 IndexType *indices_[ dimension+1 ];
480
481 // the size of each codimension
482 IndexType size_[ dimension+1 ];
483 };
484
485
486
487 // AlbertaGridIndexSet::Insert
488 // ---------------------------
489
490 template< int dim, int dimworld >
491 template< int codim >
492 struct AlbertaGridIndexSet< dim, dimworld >::Insert
493 {
494 static void apply ( const Alberta::Element *const element,
495 AlbertaGridIndexSet< dim, dimworld > &indexSet )
496 {
497 int *const array = indexSet.indices_[ codim ];
498 IndexType &size = indexSet.size_[ codim ];
499
500 for( int i = 0; i < Alberta::NumSubEntities< dim, codim >::value; ++i )
501 {
502 int &index = array[ indexSet.dofNumbering_( element, codim, i ) ];
503 if( index < 0 )
504 index = size++;
505 }
506 }
507 };
508
509
510
511 // AlbertaGridIdSet
512 // ----------------
513
515 template< int dim, int dimworld >
516 class AlbertaGridIdSet
517 : public IdSet< AlbertaGrid< dim, dimworld >, AlbertaGridIdSet< dim, dimworld >, unsigned int >
518 {
519 typedef AlbertaGridIdSet< dim, dimworld > This;
520 typedef IdSet< AlbertaGrid< dim, dimworld >, This, unsigned int > Base;
521
522 friend class AlbertaGrid< dim, dimworld >;
523
524 public:
526 typedef typename Base::IdType IdType;
527
528 private:
529 typedef AlbertaGrid< dim, dimworld > Grid;
530
531 static const int dimension = Grid::dimension;
532
533 typedef typename Grid::HierarchicIndexSet HierarchicIndexSet;
534
535 // create id set, only allowed for AlbertaGrid
536 AlbertaGridIdSet ( const HierarchicIndexSet &hIndexSet )
537 : hIndexSet_( hIndexSet )
538 {}
539
540 public:
542 template< class Entity >
543 IdType id ( const Entity &e ) const
544 {
545 const int codim = Entity::codimension;
546 return id< codim >( e );
547 }
548
550 template< int codim >
551 IdType id ( const typename Grid::template Codim< codim >::Entity &e ) const
552 {
553 assert( (codim >= 0) && (codim <= dimension) );
554 const IdType index = hIndexSet_.index( e );
555 return (index << 2) | IdType( codim );
556 }
557
559 IdType subId ( const typename Grid::template Codim< 0 >::Entity &e, int i, unsigned int subcodim ) const
560 {
561 assert( int( subcodim ) <= dimension );
562 const IdType index = hIndexSet_.subIndex( e, i, subcodim );
563 return (index << 2) | IdType( subcodim );
564 }
565
566 template< int codim >
567 IdType subId ( const typename Grid::template Codim< codim >::Entity &e, int i, unsigned int subcodim ) const
568 {
569 assert( (codim >= 0) && (codim <= dimension) && (int( codim + subcodim ) <= dimension) );
570 const IdType index = hIndexSet_.subIndex( e, i, subcodim );
571 return (index << 2) | IdType( codim + subcodim );
572 }
573
574 template< class Entity >
575 IdType subId ( const Entity &e, int i, unsigned int subcodim ) const
576 {
577 return subId< Entity::codimension >( e, i, subcodim );
578 }
579
580 private:
581 // prohibit copying
582 AlbertaGridIdSet ( const This & );
583
584 const HierarchicIndexSet &hIndexSet_;
585 };
586
587} // namespace Dune
588
589#endif // #if HAVE_ALBERTA
590
591#endif // #ifndef DUNE_ALBERTAGRIDINDEXSETS_HH
gridfamily.hh
provides the GridFamily for AlbertaGrid
Dune::AlbertaGridIdSet
hierarchic index set of AlbertaGrid
Definition: indexsets.hh:518
Dune::AlbertaGridIdSet::id
IdType id(const typename Grid::template Codim< codim >::Entity &e) const
Definition: indexsets.hh:551
Dune::AlbertaGridIdSet::id
IdType id(const Entity &e) const
Definition: indexsets.hh:543
Dune::AlbertaGridIdSet::IdType
Base::IdType IdType
export type of id
Definition: indexsets.hh:526
Dune::AlbertaGridIdSet::subId
IdType subId(const typename Grid::template Codim< 0 >::Entity &e, int i, unsigned int subcodim) const
Definition: indexsets.hh:559
Dune::AlbertaGrid< dim, dimworld >::HierarchicIndexSet
Traits::HierarchicIndexSet HierarchicIndexSet
type of hierarchic index set
Definition: agrid.hh:157
Dune::Entity
Wrapper class for entities.
Definition: entity.hh:66
Dune::Entity::codimension
static constexpr int codimension
Know your own codimension.
Definition: entity.hh:106
Dune::IdSet
Id Set Interface.
Definition: indexidset.hh:447
Dune::IdSet::IdType
IdTypeImp IdType
Type used to represent an id.
Definition: indexidset.hh:453
Dune::IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > >::subIndex
IndexType subIndex(const Entity &e, int i, unsigned int codim) const
Map a subentity to an index.
Definition: indexidset.hh:182
Dune::IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > >::subIndex
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
Dune::IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > >::Types
std::array< GeometryType, 1 > Types
iterator range for geometry types in domain
Definition: indexidset.hh:95
Dune::IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > >::index
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
Dune::IndexSet< AlbertaGridFamily< dim, dimworld >, AlbertaGridHierarchicIndexSet< dim, dimworld >, int, std::array< GeometryType, 1 > >::index
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
Dune::IndexStack
Definition: indexstack.hh:26
grid.hh
Different resources needed by all grid implementations.
indexidset.hh
Provides base classes for index and id sets.
elementinfo.hh
provides a wrapper for ALBERTA's el_info structure
Dune::GeometryTypes::simplex
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:453
Dune::Hybrid::forEach
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
indexstack.hh
Provides an index stack that supplies indices for element numbering for a grid (i....
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::size
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Dune::contains
constexpr std::bool_constant<((II==value)||...)> contains(std::integer_sequence< T, II... >, std::integral_constant< T, value >)
Checks whether or not a given sequence contains a value.
Definition: integersequence.hh:137
stdstreams.hh
Standard Dune debug streams.
Dune::Codim
Static tag representing a codimension.
Definition: dimension.hh:24
Dune::Geo::ReferenceElements::simplex
static const ReferenceElement & simplex()
get simplex reference elements
Definition: referenceelements.hh:162
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)