DUNE-FEM (unstable)

tuplemapper.hh
1#ifndef DUNE_FEM_SPACE_COMBINEDSPACE_TUPLEMAPPER_HH
2#define DUNE_FEM_SPACE_COMBINEDSPACE_TUPLEMAPPER_HH
3
4#include <array>
5#include <tuple>
6#include <utility>
7
8#include <dune/common/hybridutilities.hh>
9
10#include <dune/fem/common/utility.hh>
11#include <dune/fem/space/common/adaptationmanager.hh>
12#include <dune/fem/space/common/dofmanager.hh>
13#include <dune/fem/space/mapper/dofmapper.hh>
14#include <dune/fem/space/mapper/nonblockmapper.hh>
15#include <dune/fem/misc/functor.hh>
16
17// save diagnostic state
18#pragma GCC diagnostic push
19// turn off the specific warning, caused by code in line 419.
20#pragma GCC diagnostic ignored "-Wattributes"
21
22namespace Dune
23{
24
25 namespace Fem
26 {
27
28 // Internal forward declaration
29 // ----------------------------
30 template< class GridPart, class ... Mapper >
31 class TupleMapper;
32
33
34#ifndef DOXYGEN
35
36 namespace __TupleMapper
37 {
38
39 // Traits
40 // ------
41
42 template< class GridPart, class ... Mapper >
43 struct Traits
44 {
45 static_assert( Std::are_all_same< typename Mapper::ElementType ... >::value,
46 "TupleMapper needs common ElementType" );
47
48 typedef typename std::tuple_element< 0, std::tuple< Mapper ... > >::type FirstMapperType;
49 typedef typename FirstMapperType::ElementType ElementType;
50 typedef typename FirstMapperType::SizeType SizeType;
51 typedef typename FirstMapperType::GlobalKeyType GlobalKeyType;
52
53 typedef TupleMapper< GridPart, Mapper ... > DofMapperType;
54 };
55
56 // CombinedIndex
57 // -------------
58
59 template< class Index, class Int, Int i >
60 struct CombinedIndex
61 {
62 constexpr CombinedIndex ( Index index, Index offset ) : index_( index ), offset_( offset ) {}
63
64 static constexpr Int component () { return i; }
65
66 constexpr operator Index () const { return index_ + offset_; }
67
68 constexpr Index index () const { return index_; }
69
70 constexpr Index offset () const { return offset_; }
71
72 private:
73 Index index_, offset_;
74 };
75
76
77 // DofMapper
78 // ---------
79
80 template< class T, template< class > class Base = Dune::Fem::DofMapper >
81 class DofMapper;
82
83 template< class GridPart, class ... Mapper, template< class > class Base >
84 class DofMapper< Traits< GridPart, Mapper ... >, Base >
85 : public Base< Traits< GridPart, Mapper ... > >
86 {
87 typedef Base< Traits< GridPart, Mapper ... > > BaseType;
88
89 // FunctorWrapper
90 // --------------
91
92 template< class Functor, int i >
93 struct FunctorWrapper
94 {
95 FunctorWrapper ( Functor functor, int localOffset, int globalOffset )
96 : functor_( functor ),
97 localOffset_( localOffset ),
98 globalOffset_( globalOffset )
99 {}
100
101 template< class GlobalKey >
102 void operator() ( int localDof, const GlobalKey &globalKey ) const
103 {
104 functor_( localDof + localOffset_, CombinedIndex< GlobalKey, int, i >( globalKey, globalOffset_ ) );
105 }
106
107 template< class GlobalKey >
108 void operator() ( const GlobalKey &globalKey ) const
109 {
110 functor_( CombinedIndex< GlobalKey, int, i >( globalKey, globalOffset_ ) );
111 }
112
113 private:
114 Functor functor_;
115 const int localOffset_;
116 const int globalOffset_;
117 };
118
119 protected:
120 // size of the Mapper Tuple
121 static const int mapperTupleSize = sizeof ... ( Mapper );
122
123 typedef std::array< typename BaseType::SizeType, mapperTupleSize + 1 > OffsetType;
124
125 public:
126 typedef typename BaseType::ElementType ElementType;
127 typedef typename BaseType::SizeType SizeType;
128 typedef typename BaseType::Traits::GlobalKeyType GlobalKeyType;
129
130 typedef GridPart GridPartType;
131
132 DofMapper ( GridPartType &gridPart, Mapper & ... mapper )
133 : gridPart_( gridPart ),
134 mapperTuple_( mapper ... )
135 {
136 computeOffset();
137 }
138
139 DofMapper ( GridPartType &gridPart, Mapper && ... mapper )
140 : gridPart_( gridPart ),
141 mapperTuple_( std::move( mapper ) ... )
142 {
143 computeOffset();
144 }
145
146 SizeType size () const { return size( std::index_sequence_for< Mapper ... >() ); }
147
148 bool contains ( const int codim ) const { return contains( codim, std::index_sequence_for< Mapper ... >() ); }
149
150 bool fixedDataSize ( int codim ) const { return fixedDataSize( codim, std::index_sequence_for< Mapper ... >() ); }
151
152 template< class Functor >
153 void mapEach ( const ElementType &element, Functor f ) const
154 {
155 OffsetType localOffset;
156 localOffset[ 0 ] = 0;
157 Hybrid::forEach( std::make_index_sequence< mapperTupleSize >{},
158 [ & ]( auto i )
159 {
160 FunctorWrapper< Functor, i > wrappedFunctor( f, localOffset[ i ], globalOffset_[ i ] );
161 std::get< i >( mapperTuple_ ).mapEach( element, wrappedFunctor );
162 localOffset[ i + 1 ] = localOffset[ i ] + std::get< i >( mapperTuple_ ).numDofs( element );
163 } );
164 }
165
166 template< class Entity, class Functor >
167 void mapEachEntityDof ( const Entity &entity, Functor f ) const
168 {
169 OffsetType localOffset;
170 localOffset[ 0 ] = 0;
171 Hybrid::forEach( std::make_index_sequence< mapperTupleSize >{},
172 [ & ]( auto i )
173 {
174 FunctorWrapper< Functor, i > wrappedFunctor( f, localOffset[ i ], globalOffset_[ i ] );
175 std::get< i >( mapperTuple_ ).mapEachEntityDof( entity, wrappedFunctor );
176 localOffset[ i + 1 ] = localOffset[ i ] + std::get< i >( mapperTuple_ ).numEntityDofs( entity );
177 } );
178 }
179
180 [[deprecated("Use onSubEntity method with char vector instead")]]
181 void onSubEntity ( const ElementType &element, int i, int c, std::vector< bool > &indices ) const
182 {
183 std::vector< char > _idx;
184 onSubEntity(element, i, c, _idx);
185 indices.resize( _idx.size() );
186 for (std::size_t i=0; i<_idx.size();++i)
187 _idx[i] = indices[i] > 0;
188 }
189 // this method returns which local dofs are attached to the given subentity.
190 // indices[locDofNr] =
191 // 0 : not attached, not equal to 0 : attached
192 // (so this method can still be used in the way the deprecated method was).
193 // New: In case the dof can be associated to a component of the
194 // space, the value returned is that component+1. In other
195 // cases (normal velocity for RT for example) the value is -1).
196 // So indices[i] is in [-1,dimRange+1]
197 void onSubEntity ( const ElementType &element, int i, int c, std::vector< char > &indices ) const
198 {
199 indices.resize( numDofs( element ) );
200 OffsetType localOffset;
201 localOffset[ 0 ] = 0;
202 int rangeOffset = 0;
203 Hybrid::forEach( std::make_index_sequence< mapperTupleSize >{},
204 [ & ]( auto t )
205 {
206 std::vector< char > subIndices;
207 std::get< t >( mapperTuple_ ).onSubEntity( element, i,c, subIndices );
208 auto localSize = std::get< t >( mapperTuple_ ).numDofs( element );
209 assert( subIndices.size() == localSize );
210 assert( localOffset[ t ] + localSize <= indices.size() );
211 for (std::size_t d=0;d<subIndices.size();++d)
212 {
213 indices[ localOffset[t] + d ] = subIndices[d]==0? 0 : subIndices[d] + rangeOffset;
214 }
215 // std::copy( subIndices.begin(), subIndices.end(), indices.begin() + localOffset[ t ] );
216 /*
217 std::cout << "\t space<" << t << ">:";
218 for (std::size_t d=0;d<localSize;++d)
219 std::cout << " " << subIndices[d];
220 std::cout << " ->";
221 for (std::size_t d=0;d<indices.size();++d)
222 std::cout << " " << indices[d];
223 std::cout << std::endl;
224 */
225 // FIXME: here we need 'dimRange' of the subspace to cover
226 // cases where a space is vector valued by has blockSize=1
227 // like RT. Otherwise we will get for RTxRT something like
228 // [1,0,0,1, 2,0,0,2] instead of [1,0,0,1, 3,0,0,3]
229 rangeOffset += std::get< t >( mapperTuple_ ).blockSize;
230 localOffset[ t + 1 ] = localOffset[ t ] + localSize;
231 } );
232 /*
233 std::cout << "element (" << i << "," << c << "):";
234 for (std::size_t d=0;d<indices.size();++d)
235 std::cout << " " << indices[d];
236 std::cout << std::endl;
237 */
238 // DUNE_THROW( NotImplemented, "Method onSubEntity(...) not yet implemented for TupleMapper" );
239 }
240
241 int maxNumDofs () const { return maxNumDofs( std::index_sequence_for< Mapper ... >() ); }
242
243 SizeType numDofs ( const ElementType &element ) const { return numDofs( element, std::index_sequence_for< Mapper ... >() ); }
244
245 template< class Entity >
246 SizeType numEntityDofs ( const Entity &entity ) const { return numEntityDofs( entity, std::index_sequence_for< Mapper ... >() ); }
247
248 void map ( const ElementType &element, std::vector< SizeType > &indices ) const
249 {
250 indices.resize( numDofs( element ) );
251 mapEach( element, AssignFunctor< std::vector< SizeType > >( indices ) );
252 }
253 template< class Entity >
254 void mapEntityDofs ( const Entity &entity, std::vector< SizeType > &indices ) const
255 {
256 indices.resize( numEntityDofs( entity ) );
257 mapEachEntityDof( entity, AssignFunctor< std::vector< SizeType > >( indices ) );
258 }
259
260 static constexpr bool consecutive () noexcept { return false; }
261
262 SizeType numBlocks () const
263 {
264 DUNE_THROW( NotImplemented, "Method numBlocks() called on non-adaptive block mapper" );
265 }
266
267 SizeType numberOfHoles ( int ) const
268 {
269 DUNE_THROW( NotImplemented, "Method numberOfHoles() called on non-adaptive block mapper" );
270 }
271
272 GlobalKeyType oldIndex ( int hole, int ) const
273 {
274 DUNE_THROW( NotImplemented, "Method oldIndex() called on non-adaptive block mapper" );
275 }
276
277 GlobalKeyType newIndex ( int hole, int ) const
278 {
279 DUNE_THROW( NotImplemented, "Method newIndex() called on non-adaptive block mapper" );
280 }
281
282 SizeType oldOffSet ( int ) const
283 {
284 DUNE_THROW( NotImplemented, "Method oldOffSet() called on non-adaptive block mapper" );
285 }
286
287 SizeType offSet ( int ) const
288 {
289 DUNE_THROW( NotImplemented, "Method offSet() called on non-adaptive block mapper" );
290 }
291
292 void update()
293 {
294 // compute update for each mapper (if any)
295 Hybrid::forEach( std::make_index_sequence< mapperTupleSize >{},
296 [ & ](auto i){ std::get< i >( mapperTuple_ ).update(); } );
297
298 computeOffset();
299 }
300
301 /*** NonInterface Methods ***/
302
303 SizeType offset ( int i ) const { return globalOffset_[ i ]; }
304
305 template< int i >
306 SizeType subSize () const { return std::get< i >( mapperTuple_ ).size(); }
307
308 protected:
309 template< std::size_t ... i >
310 SizeType size ( std::index_sequence< i ... > ) const
311 {
312 return Std::sum( std::get< i >( mapperTuple_ ).size() ... );
313 }
314
315 template< std::size_t ... i >
316 bool fixedDataSize ( const int codim, std::index_sequence< i ... > ) const
317 {
318 return Std::And( std::get< i >( mapperTuple_ ).fixedDataSize( codim ) ... );
319 }
320
321 template< std::size_t ... i >
322 bool contains ( const int codim, std::index_sequence< i ... > ) const
323 {
324 return Std::Or( std::get< i >( mapperTuple_ ).contains( codim ) ... );
325 }
326
327 template< std::size_t ... i >
328 int maxNumDofs ( std::index_sequence< i ... > ) const
329 {
330 return Std::sum( std::get< i >( mapperTuple_ ).maxNumDofs() ... );
331 }
332
333 template< std::size_t ... i >
334 SizeType numDofs ( const ElementType &element, std::index_sequence< i ... > ) const
335 {
336 return Std::sum( std::get< i >( mapperTuple_ ).numDofs( element ) ... );
337 }
338
339 template< class Entity, std::size_t ... i >
340 SizeType numEntityDofs ( const Entity &entity, std::index_sequence< i ... > ) const
341 {
342 return Std::sum( std::get< i >( mapperTuple_ ).numEntityDofs( entity ) ... );
343 }
344
345 void computeOffset ()
346 {
347 globalOffset_[ 0 ] = 0;
348 // compute new offsets
349 Hybrid::forEach( std::make_index_sequence< mapperTupleSize >{},
350 [ & ]( auto i ){ globalOffset_[ i + 1 ] = globalOffset_[ i ] + std::get< i >( mapperTuple_ ).size(); } );
351 }
352
353 GridPartType &gridPart_;
354 std::tuple< Mapper ... > mapperTuple_;
355 OffsetType globalOffset_;
356 };
357
358
359
360 // AdaptiveDofMapper
361 // -----------------
362
363 template< class T >
364 class AdaptiveDofMapper;
365
366 template< class GridPart, class ... Mapper >
367 class AdaptiveDofMapper< Traits< GridPart, Mapper ... > >
368 : public DofMapper< Traits< GridPart, Mapper ... >, Dune::Fem::AdaptiveDofMapper >
369 {
370 typedef DofMapper< Traits< GridPart, Mapper ... >, Dune::Fem::AdaptiveDofMapper > BaseType;
371
372 protected:
373 typedef typename GridPart :: GridType GridType;
374 typedef AdaptationMethod< GridType > AdaptationMethodType;
375
376 typedef typename BaseType::OffsetType OffsetType;
377
378 using BaseType::mapperTupleSize;
379 using BaseType::mapperTuple_;
380 using BaseType::gridPart_;
381 using BaseType::globalOffset_;
382
383 public:
384 typedef typename BaseType::ElementType ElementType;
385 typedef typename BaseType::SizeType SizeType;
386 typedef typename BaseType::GlobalKeyType GlobalKeyType;
387 typedef GridPart GridPartType;
388
389 AdaptiveDofMapper ( GridPartType &gridPart, Mapper & ... mapper )
390 : BaseType( gridPart, mapper ... ),
391 numBlocks_( computeNumBlocks() ),
392 isCallBackAdapt_( AdaptationMethodType( gridPart.grid() ).isCallBackAdaptation() ),
393 needsFullUpdate_(true)
394 {
395 oldGlobalOffset_ = globalOffset_;
397 }
398
399 AdaptiveDofMapper ( GridPartType &gridPart, Mapper && ... mapper )
400 : BaseType( gridPart, std::move( mapper ) ... ),
401 numBlocks_( computeNumBlocks() ),
402 isCallBackAdapt_( AdaptationMethodType( gridPart.grid() ).isCallBackAdaptation() ),
403 needsFullUpdate_(true)
404 {
405 oldGlobalOffset_ = globalOffset_;
407 }
408
409 ~AdaptiveDofMapper () { DofManager< typename GridPartType::GridType >::instance( gridPart_.grid() ).removeIndexSet( *this ); }
410
411 AdaptiveDofMapper ( const AdaptiveDofMapper & ) = delete;
412 AdaptiveDofMapper ( AdaptiveDofMapper && ) = delete;
413
414 static constexpr bool consecutive () noexcept { return true; }
415
416 SizeType numBlocks () const { return numBlocks( std::index_sequence_for< Mapper ... >() ); }
417
418 SizeType numberOfHoles ( const int block ) const
419 {
420 SizeType nHoles = 0;
421 Hybrid::forEach( std::make_index_sequence< mapperTupleSize >{},
422 [ & ]( auto i )
423 {
424 const int localBlock = computeBlock( i, block );
425 if( localBlock >= 0 )
426 {
427 nHoles = std::get< i >( this->mapperTuple_ ).numberOfHoles( localBlock );
428 return;
429 }
430 } );
431 return nHoles;
432 }
433
434 GlobalKeyType oldIndex ( const int hole, const int block ) const
435 {
436 SizeType oIndex = 0;
437 Hybrid::forEach( std::make_index_sequence< mapperTupleSize >{},
438 [ & ]( auto i )
439 {
440 const int localBlock = computeBlock( i, block );
441 if( localBlock >= 0 )
442 {
443 oIndex = std::get< i >( this->mapperTuple_ ).oldIndex( hole, localBlock ) + globalOffset_[ i ];
444 return;
445 }
446 } );
447 return oIndex;
448 }
449
450 GlobalKeyType newIndex ( const int hole, const int block ) const
451 {
452 SizeType nIndex = 0;
453 Hybrid::forEach( std::make_index_sequence< mapperTupleSize >{},
454 [ & ]( auto i )
455 {
456 const int localBlock = computeBlock( i, block );
457 if( localBlock >= 0 )
458 {
459 nIndex = std::get< i >( this->mapperTuple_ ).newIndex( hole, localBlock ) + globalOffset_[ i ];
460 return ;
461 }
462 } );
463 return nIndex;
464 }
465
466 SizeType oldOffSet ( int block ) const
467 {
468 SizeType oOffset = 0;
469 Hybrid::forEach( std::make_index_sequence< mapperTupleSize >{},
470 [ & ]( auto i )
471 {
472 const int localBlock = computeBlock( i, block );
473 if( localBlock >= 0 )
474 {
475 oOffset = std::get< i >( this->mapperTuple_ ).oldOffSet( localBlock ) + oldGlobalOffset_[ i ];
476 return ;
477 }
478 } );
479 return oOffset;
480 }
481
482 SizeType offSet ( int block ) const
483 {
484 SizeType offset = 0;
485 Hybrid::forEach( std::make_index_sequence< mapperTupleSize >{},
486 [ & ]( auto i )
487 {
488 const int localBlock = computeBlock( i, block );
489 if( localBlock >= 0 )
490 {
491 offset = std::get< i >( this->mapperTuple_ ).offSet( localBlock ) + globalOffset_[ i ];
492 return ;
493 }
494 } );
495 return offset;
496 }
497
498 void resize () { update(); }
499
500 bool compress ()
501 {
502 update();
503 // after this step for both methods we need the full update
504 needsFullUpdate_ = true;
505 return true;
506 }
507
508 void backup () const {}
509
510 void restore () { compress(); }
511
512 template< class IOStream >
513 void read ( IOStream &in ) { compress(); }
514
515 template< class IOStream >
516 void write ( IOStream &out ) const {}
517
518 template< class Entity >
519 void insertEntity ( const Entity & )
520 {
521 if( needsFullUpdate_ )
522 {
523 // update also offsets in this case
524 update();
525 }
526 else
527 {
528 // call BaseType::update to avoid changing oldGlobalOffSet_
529 // do not use this->update here!
530 BaseType::update();
531 }
532 }
533
534 template< class Entity >
535 void removeEntity ( const Entity & ) {}
536
537 void update ()
538 {
539 // store previous offset
540 oldGlobalOffset_ = globalOffset_;
541
542 // in callback we always need the full update, otherwise set to false here
543 needsFullUpdate_ = isCallBackAdapt_ ;
544
545 // update component mappers and compute offsets
546 BaseType::update();
547 }
548
549 protected:
550 int computeBlock( const int i, const unsigned int block ) const
551 {
552 if( block >= blocks_[ i ] && block < blocks_[ i + 1 ] )
553 return block - blocks_[ i ];
554 else
555 return -1;
556 }
557
558 void printOffSet( const OffsetType& offset, const std::string& name ) const
559 {
560 for( const auto& off : offset )
561 {
562 std::cout << name << " = " << off << std::endl;
563 }
564 }
565
566 template< std::size_t ... i >
567 SizeType numBlocks ( std::index_sequence< i ... > ) const
568 {
569 return Std::sum( std::get< i >( mapperTuple_ ).numBlocks() ... );
570 }
571
572 SizeType computeNumBlocks ()
573 {
574 // compute blocks (only needs to be done once)
575 // length of blocks_ is mapperTupleSize + 1
576 blocks_[ 0 ] = 0;
577 Hybrid::forEach( std::make_index_sequence< mapperTupleSize >{},
578 [ & ]( auto i ){ blocks_[ i + 1 ] = blocks_[ i ] + std::get< i >( mapperTuple_ ).numBlocks(); } );
579
580 // last block number determines overall block numbers
581 return blocks_.back();
582 }
583
584
585 OffsetType oldGlobalOffset_;
586 OffsetType blocks_;
587
588 const SizeType numBlocks_;
589 const bool isCallBackAdapt_;
590 bool needsFullUpdate_;
591 };
592
593
594
595 // Implementation
596 // --------------
597
598 template< class GridPart, class ... Mapper >
599 struct Implementation
600 {
601 typedef typename std::conditional<
602 Std::And( Capabilities::isAdaptiveDofMapper< Mapper >::v ... ),
603 AdaptiveDofMapper< Traits< GridPart, Mapper ... > >,
604 DofMapper< Traits< GridPart, Mapper ... > > >::type Type;
605 };
606
607
608 } // namespace __TupleMapper
609
610#endif // #ifndef DOXYGEN
611
612
613 // TupleMapper
614 // -----------
615
626 template< class GridPart, class ... Mapper >
628 : public __TupleMapper::template Implementation< GridPart, Mapper ... >::Type
629 {
630 typedef typename __TupleMapper::template Implementation< GridPart, Mapper ... >::Type BaseType;
631
632 public:
633 TupleMapper ( GridPart &gridPart, Mapper & ... mapper ) : BaseType( gridPart, mapper ... ) {}
634 TupleMapper ( GridPart &gridPart, Mapper && ... mapper ) : BaseType( gridPart, std::move( mapper ) ... ) {}
635 };
636
637 // Capabilities
638 // ------------
639
640 namespace Capabilities
641 {
642 template< class GridPart, class ... Mapper >
643 struct isAdaptiveDofMapper< TupleMapper< GridPart, Mapper ... > >
644 {
645 static const bool v = Std::And( isAdaptiveDofMapper< Mapper >::v ... );
646 };
647
648 template< class GridPart, class ... Mapper >
649 struct isConsecutiveIndexSet< __TupleMapper::AdaptiveDofMapper< __TupleMapper::Traits< GridPart, Mapper ... > > >
650 {
651 static const bool v = true;
652 };
653
654 } // namespace Capabilities
655
656 } // namespace Fem
657
658} // namespace Dune
659
660// turn the warnings back on
661#pragma GCC diagnostic pop
662
663#endif // #ifndef DUNE_FEM_SPACE_COMBINEDSPACE_TUPLEMAPPER_HH
Extended interface for adaptive DoF mappers.
Definition: dofmapper.hh:219
Interface for calculating the size of a function space for a grid on a specified level....
Definition: dofmapper.hh:43
mapper allocating one DoF per subentity of a given codimension
Definition: tuplemapper.hh:629
Mapper interface.
Definition: mapper.hh:110
void removeIndexSet(const IndexSetType &iset)
removed index set from dof manager's list of index sets
Definition: dofmanager.hh:1331
static ThisType & instance(const GridType &grid)
obtain a reference to the DofManager for a given grid
Definition: dofmanager.hh:1251
void addIndexSet(const IndexSetType &iset)
add index set to dof manager's list of index sets
Definition: dofmanager.hh:1296
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
Dune namespace.
Definition: alignedallocator.hh:13
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
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
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)