dune-fem  2.6-git
blockmapper.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_HPDG_SPACE_DISCONTINUOUSGALERKIN_BLOCKMAPPER_HH
2 #define DUNE_FEM_HPDG_SPACE_DISCONTINUOUSGALERKIN_BLOCKMAPPER_HH
3 
4 #include <cassert>
5 #include <cstddef>
6 
7 #include <algorithm>
8 #include <functional>
9 #include <limits>
10 #include <type_traits>
11 #include <utility>
12 #include <vector>
13 
14 #include <dune/common/exceptions.hh>
15 
16 #include <dune/geometry/dimension.hh>
17 
18 #include <dune/grid/common/gridenums.hh>
19 #include <dune/grid/utility/persistentcontainer.hh>
20 
25 
27 
28 #include "datahandle.hh"
29 #include "localdofstorage.hh"
30 
31 namespace Dune
32 {
33 
34  namespace Fem
35  {
36 
37  namespace hpDG
38  {
39 
40  // Internal forward declaration
41  // ----------------------------
42 
43  template< class GridPart, class LocalKeys >
44  struct DiscontinuousGalerkinBlockMapper;
45 
46 
47 
48 #ifndef DOXYGEN
49 
50  // DiscontinuousGalerkinBlockMapperTraits
51  // --------------------------------------
52 
53  template< class GridPart, class LocalKeys >
54  struct DiscontinuousGalerkinBlockMapperTraits
55  {
56  using DofMapperType = DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys >;
57  using ElementType = typename GridPart::template Codim< 0 >::EntityType;
58  using SizeType = std::size_t;
59  };
60 
61 #endif // #ifndef DOXYGEN
62 
63 
64 
65  // DiscontinuousGalerkinBlockMapper
66  // --------------------------------
67 
99  template< class GridPart, class LocalKeys >
101  : public AdaptiveDofMapper< DiscontinuousGalerkinBlockMapperTraits< GridPart, LocalKeys > >
102  {
105 
106  public:
108  using SizeType = typename BaseType::SizeType;
111 
113  using GridPartType = GridPart;
116 
118  using LocalKeysType = LocalKeys;
120  using KeyType = typename LocalKeysType::KeyType;
121 
122  private:
124  friend class DataHandle< DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys > >;
125 
126  struct Reserve;
127  struct Resize;
128 
130 
131  using GridType = typename GridPartType::GridType;
132  using GridElementType = typename GridType::template Codim< 0 >::Entity;
133 
134  public:
137 
142  template< class Function >
144  const LocalKeysType &localKeys,
145  const KeyType &value, Function function )
146  : gridPart_( gridPart ),
147  localKeys_( localKeys ),
148  key_( value ),
149  keys_( gridPart_.get().grid(), 0, std::make_pair( key_, key_ ) ),
150  size_( 0u ),
151  dofs_( gridPart_.get().grid(), 0 ),
152  dofManager_( DofManagerType::instance( gridPart.grid() ) )
153  {
154  auto first = gridPart.template begin< 0, All_Partition >();
155  auto last = gridPart.template end< 0, All_Partition >();
156  for( ; first != last; ++first )
157  {
158  const ElementType &element = *first;
159  const KeyType key = function( element );
160 
161  const GridElementType &gridElement = gridEntity( element );
162  keys_[ gridElement ] = std::make_pair( key, key );
163 
164  auto &dofs = dofs_[ gridElement ];
165  resize( dofs, localKeys.blocks( gridElement.type(), key ) );
166  }
167 
168  dofManager_.get().addIndexSet( *this );
169  }
170 
172  const LocalKeysType &localKeys,
173  const KeyType &value )
174  : DiscontinuousGalerkinBlockMapper( gridPart, localKeys, value, [&value]( const ElementType & ){ return value; } )
175  {}
176 
185 
188 
190  ThisType &operator= ( const ThisType & ) = delete;
191 
193  ThisType &operator= ( ThisType && ) = default;
194 
197 #ifndef DOXYGEN
198 
200  {
201  dofManager_.get().removeIndexSet( *this );
202  }
203 
204 #endif // #ifndef DOXYGEN
205 
211  SizeType size () const { return size_; }
212 
214  int maxNumDofs () const { return localKeys().maxBlocks(); }
215 
217  SizeType numDofs ( const ElementType &element ) const
218  {
219  return numDofs( element, Codim< ElementType::codimension >() );
220  }
221 
223  template< class Entity >
224  SizeType numEntityDofs ( const Entity &entity ) const
225  {
226  return numDofs( entity, Codim< Entity::codimension >() );
227  }
228 
230  static constexpr bool contains ( const int codim ) { return (codim == 0); }
231 
233  static constexpr bool fixedDataSize ( int codim ) { return (codim != 0); }
234 
236  template< class Function >
237  void mapEach ( const ElementType &element, Function function ) const
238  {
239  mapEach( element, function, Codim< ElementType::codimension >() );
240  }
241 
243  template< class Entity, class Function >
244  void mapEachEntityDof ( const Entity &entity, Function function ) const
245  {
246  mapEach( entity, function, Codim< Entity::codimension >() );
247  }
248 
250  SizeType numberOfHoles ( const int block ) const
251  {
252  assert( block == 0 );
253  return indices_.size();
254  }
255 
257  GlobalKeyType oldIndex ( const int hole, const int block ) const
258  {
259  assert( block == 0 );
260  GlobalKeyType dof = indices_[ hole ].first;
262  return std::move( dof );
263  }
264 
266  GlobalKeyType newIndex ( const int hole, const int block ) const
267  {
268  assert( block == 0 );
269  return indices_[ hole ].second;
270  }
271 
273  static constexpr bool consecutive () { return true; }
274 
276  static constexpr SizeType oldOffSet ( const int block ) { return 0u; }
277 
279  static constexpr SizeType offSet ( const int block ) { return 0u; }
280 
282  static constexpr SizeType numBlocks () { return 1u; }
283 
291  template< class Element >
292  typename std::enable_if< (std::is_same< Element, ElementType >::value || std::is_same< Element, GridElementType >::value), const KeyType & >::type
293  key ( const Element &element ) const
294  {
295  return keys_[ gridEntity( element ) ].first;
296  }
297 
299  void mark ( const KeyType &key, const ElementType &element )
300  {
301  assert( gridPart().indexSet().contains( element ) );
302  keys_[ gridEntity( element ) ].second = key;
303  }
304 
306  KeyType getMark ( const ElementType &element ) const
307  {
308  assert( gridPart().indexSet().contains( element ) );
309  return keys_[ gridEntity( element ) ].second;
310  }
311 
313  template< class Function >
314  bool adapt ( Function function );
315 
317  bool adapt ()
318  {
319  auto function = []( const ElementType &,
320  const KeyType &, const KeyType &,
321  const std::vector< std::size_t > &,
322  const std::vector< std::size_t > & )
323  {};
324  return adapt( function );
325  }
326 
334  const DofManagerType &dofManager () const { return dofManager_.get(); }
335 
337  void insertEntity ( const GridElementType &gridElement )
338  {
339  resize( dofs_ );
340  resize( keys_, std::make_pair( key_, key_ ) );
341 
342  resize( dofs_[ gridElement ], localKeys().blocks( gridElement.type(), key( gridElement ) ) );
343  }
344 
346  void removeEntity ( const GridElementType &gridElement )
347  {
348  auto &dofs = dofs_[ gridElement ];
349 
350  size_ = dofs.reserve( 0u, Reserve( size_ ) );
351  Resize function( holes_ );
352  for( auto dof : dofs )
353  function( dof );
354  }
355 
357  void insertNewEntity ( const GridElementType &gridElement )
358  {
359  if( gridElement.isNew() )
360  {
361  resize( dofs_[ gridElement ], localKeys().blocks( gridElement.type(), key( gridElement ) ) );
362  assert( gridElement.hasFather() );
363  const GridElementType &father = gridElement.father();
364  insertNewEntity( father );
365  removeEntity( father );
366  }
367  else
368  {
369  auto &dofs = dofs_[ gridElement ];
370  resize( dofs, localKeys().blocks( gridElement.type(), key( gridElement ) ) );
371  for( auto dof : dofs )
372  {
373  auto iterator = std::lower_bound( holes_.begin(), holes_.end(), dof );
374  if( iterator != holes_.end() && *iterator == dof )
375  holes_.erase( iterator );
376  }
377  }
378  }
379 
380  void resize ()
381  {
382  resize( dofs_ );
383  resize( keys_, std::make_pair( key_, key_ ) );
384 
385  auto first = gridPart().template begin< 0, All_Partition >();
386  auto last = gridPart().template end< 0, All_Partition >();
387  for( ; first != last; ++first )
388  {
389  const ElementType &element = *first;
390  insertNewEntity( gridEntity( element ) );
391  }
392 #if 0
393  auto &dofs = dofs_[ gridEntity( element ) ];
394  resize( dofs, localKeys().blocks( element.type(), key( element ) ) );
395 #endif
396  }
397 
399  bool compress ();
400 
402  template< class Traits >
404  {
405  DUNE_THROW( NotImplemented, "Method write() not implemented yet" );
406  }
407 
409  template< class Traits >
411  {
412  DUNE_THROW( NotImplemented, "Method read() not implemented yet" );
413  }
414 
416  void backup () const
417  {
418  DUNE_THROW( NotImplemented, "Method backup() not implemented" );
419  }
420 
422  void restore ()
423  {
424  DUNE_THROW( NotImplemented, "Method restore() not implemented" );
425  }
426 
427  private:
428  bool compressed () const { return holes_.empty(); }
429 
430  template< class Element, class Function >
431  void mapEach ( const Element &element, Function function, Codim< 0 > ) const
432  {
433  const auto &dofs = dofs_[ gridEntity( element ) ];
434  std::size_t size = dofs.size();
435  for( std::size_t i = 0; i < size; ++i )
436  function( i, dofs[ i ] );
437  }
438 
439  template< class Entity, class Function, int codim >
440  void mapEach ( const Entity &entity, Function function, Codim< codim > ) const
441  {}
442 
443  template< class Element >
444  SizeType numDofs ( const Element &element, Codim< 0 > ) const
445  {
446  return dofs_[ gridEntity( element ) ].size();
447  }
448 
449  template< class Entity, int codim >
450  SizeType numDofs ( const Entity &entity, Codim< codim > ) const
451  {
452  return 0u;
453  }
454 
455  void resize ( LocalDofStorageType &dofs, SizeType n )
456  {
457  size_ = dofs.reserve( n, Reserve( size_ ) );
458  dofs.resize( Resize( holes_ ) );
459  }
460 
461  template< class T >
462  static void resize ( PersistentContainer< GridType, T > &container, const T &value = T() )
463  {
464  container.resize( value );
465  container.shrinkToFit();
466  }
467 
468  const GridPartType &gridPart () const { return gridPart_.get(); }
469 
470  const LocalKeysType &localKeys () const { return localKeys_.get(); }
471 
472  std::reference_wrapper< const GridPartType > gridPart_;
473  std::reference_wrapper< const LocalKeysType > localKeys_;
474 
475  const KeyType key_;
476  PersistentContainer< GridType, std::pair< KeyType, KeyType > > keys_;
477 
478  SizeType size_;
479  PersistentContainer< GridType, LocalDofStorageType > dofs_;
480 
481  std::vector< GlobalKeyType > holes_;
482  std::vector< std::pair< GlobalKeyType, GlobalKeyType > > indices_;
483 
484  std::reference_wrapper< DofManagerType > dofManager_;
485  };
486 
487 
488 
489  // DiscontinuousGalerkinBlockMapper::Reserve
490  // -----------------------------------------
491 
492  template< class GridPart, class LocalKeys >
493  struct DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys >::Reserve
494  {
495  explicit Reserve ( SizeType &size ) : size_( size ) {}
496 
497  GlobalKeyType operator() () { return size_++; }
498 
499  operator SizeType () const { return size_; }
500 
501  private:
502  SizeType &size_;
503  };
504 
505 
506 
507  // DiscontinuousGalerkinBlockMapper::Resize
508  // ----------------------------------------
509 
510  template< class GridPart, class LocalKeys >
511  struct DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys >::Resize
512  {
513  explicit Resize ( std::vector< GlobalKeyType > &holes )
514  : holes_( holes )
515  {
516  assert( std::is_sorted( holes_.begin(), holes_.end() ) );
517  }
518 
520  {
521  assert( std::is_sorted( holes_.begin(), holes_.end() ) );
522  }
523 
524  void operator() ( const GlobalKeyType &dof )
525  {
526  auto iterator = std::lower_bound( holes_.begin(), holes_.end(), dof );
527  if( iterator == holes_.end() || *iterator != dof )
528  holes_.insert( iterator, dof );
529  }
530 
531  private:
532  std::vector< GlobalKeyType > &holes_;
533  };
534 
535 
536 
537  // DiscontinuousGalerkinBlockMapper::adapt
538  // ---------------------------------------
539 
540  template< class GridPart, class LocalKeys >
541  template< class Function >
543  {
544  DataHandleType dataHandle( *this );
545  gridPart().communicate( dataHandle, InteriorBorder_All_Interface, ForwardCommunication );
546 
547  std::vector< std::size_t > origin, destination;
548  origin.reserve( maxNumDofs() );
549  destination.reserve( maxNumDofs() );
550 
551  std::size_t count = 0;
552 
553  auto first = gridPart().template begin< 0, All_Partition >();
554  auto last = gridPart().template end< 0, All_Partition >();
555  for( ; first != last; ++first )
556  {
557  const ElementType &element = *first;
558  const GridElementType &gridElement = gridEntity( element );
559  auto &keys = keys_[ gridElement ];
560  if( keys.first == keys.second )
561  continue;
562 
563  // remember old state
564  const KeyType prior = keys.first;
565  auto &dofs = dofs_[ gridElement ];
566  origin.resize( dofs.size() );
567  std::copy( dofs.begin(), dofs.end(), origin.begin() );
568 
569  // reset to new state
570  keys.first = keys.second;
571  resize( dofs, localKeys().blocks( gridElement.type(), keys.first ) );
572 
573  // call function
574  destination.resize( dofs.size() );
575  std::copy( dofs.begin(), dofs.end(), destination.begin() );
576  function( element, prior, keys.first, origin, destination );
577 
578  ++count;
579  }
580  return (count > 0u);
581  }
582 
583 
584 
585  // DiscontinuousGalerkinBlockMapper::compress
586  // ------------------------------------------
587 
588  template< class GridPart, class LocalKeys >
590  {
591  // resize persistent containers
592  resize( dofs_ );
593  resize( keys_, std::make_pair( key_, key_ ) );
594 
595  for( auto &dofs : dofs_ )
596  dofs.resize( Resize( holes_ ) );
597 
598  // check if compress is necessary
599  if( holes_.empty() )
600  {
601  indices_.clear();
602  return false;
603  }
604 
605  // get number of dofs after compress
606  assert( size_ >= holes_.size() );
607  size_ -= holes_.size();
608 
609  // remove trailing indices
610  auto iterator = std::lower_bound( holes_.begin(), holes_.end(), size_ );
611  holes_.erase( iterator, holes_.end() );
612 
613  // initialize vector of old and new indices
614  auto assign = []( GlobalKeyType newIndex ) -> std::pair< GlobalKeyType, GlobalKeyType >
615  {
617  return std::make_pair( oldIndex, newIndex );
618  };
619  indices_.resize( holes_.size() );
620  std::transform( holes_.begin(), holes_.end(), indices_.begin(), assign );
621 
622  // update local dof storages
623  std::size_t hole = 0u;
624  for( auto &dofs : dofs_ )
625  {
626  for( auto &dof : dofs )
627  {
628  if( dof >= size_ )
629  {
630  auto &indices = indices_.at( hole++ );
631  indices.first = dof;
632  dof = indices.second;
633  }
634  }
635  }
636 
637  assert( hole == holes_.size() );
638  holes_.clear();
639 
640  // replace all but markers currently used by default value
641  PersistentContainer< GridType, std::pair< KeyType, KeyType > > tmp( gridPart().grid(), 0, std::make_pair( key_, key_ ) );
642  auto first = gridPart().template begin< 0, All_Partition >();
643  auto last = gridPart().template end< 0, All_Partition >();
644  for( ; first != last; ++first )
645  {
646  const ElementType &element = *first;
647  const GridElementType &gridElement = gridEntity( element );
648  tmp[ gridElement ] = keys_[ gridElement ];
649  }
650  keys_.swap( tmp );
651 
652  return true;
653  }
654 
655  } // namespace hpDG
656 
657  namespace Capabilities
658  {
659  // isConsecutiveIndexSet
660  // ---------------------
661 
662  template< class GridPart, class LocalKeys >
663  struct isConsecutiveIndexSet< hpDG::DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys > >
664  {
665  static const bool v = true;
666  };
667 
668  // isAdaptiveDofMapper
669  // -------------------
670 
671  template< class GridPart, class LocalKeys >
672  struct isAdaptiveDofMapper< hpDG::DiscontinuousGalerkinBlockMapper< GridPart, LocalKeys > >
673  {
674  static const bool v = true;
675  };
676 
677  } // namespace Capabilities
678 
679  } // namespace Fem
680 
681 } // namespace Dune
682 
683 #endif // #ifndef DUNE_FEM_HPDG_SPACE_DISCONTINUOUSGALERKIN_BLOCKMAPPER_HH
Definition: bindguard.hh:11
std::tuple_element< i, Tuple >::type & get(Dune::TypeIndexedTuple< Tuple, Types > &tuple)
Definition: typeindexedtuple.hh:122
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
const GridEntityAccess< Entity >::GridEntityType & gridEntity(const Entity &entity)
Definition: gridpart.hh:412
static constexpr T max(T a)
Definition: utility.hh:77
Abstract class representing a function.
Definition: common/function.hh:50
specialize with true if index set implements the interface for consecutive index sets
Definition: common/indexset.hh:42
static const bool v
Definition: common/indexset.hh:49
abstract interface for an output stream
Definition: streams.hh:46
abstract interface for an input stream
Definition: streams.hh:179
Definition: dofmanager.hh:825
An -adaptive Dune::Fem::DofMapper.
Definition: blockmapper.hh:102
void backup() const
this mapper has no I/O capabilities
Definition: blockmapper.hh:416
static constexpr SizeType oldOffSet(const int block)
return 0 (this mapper has no offset)
Definition: blockmapper.hh:276
DiscontinuousGalerkinBlockMapper(const GridPartType &gridPart, const LocalKeysType &localKeys, const KeyType &value)
Definition: blockmapper.hh:171
SizeType size() const
return number of dofs
Definition: blockmapper.hh:211
GlobalKeyType newIndex(const int hole, const int block) const
return new index of given hole during compression
Definition: blockmapper.hh:266
typename LocalKeysType::KeyType KeyType
key type
Definition: blockmapper.hh:120
static constexpr bool consecutive()
return true (this mapper yields a consecutive DOF numbering)
Definition: blockmapper.hh:273
GridPart GridPartType
grid part type
Definition: blockmapper.hh:113
void write(OutStreamInterface< Traits > &)
this mapper has no I/O capabilities
Definition: blockmapper.hh:403
void mapEachEntityDof(const Entity &entity, Function function) const
map local dof to global key
Definition: blockmapper.hh:244
DiscontinuousGalerkinBlockMapper(const ThisType &)=delete
copy constructor
void mapEach(const ElementType &element, Function function) const
map local dof to global key
Definition: blockmapper.hh:237
LocalKeys LocalKeysType
basis function sets type
Definition: blockmapper.hh:118
ThisType & operator=(const ThisType &)=delete
assignment operator
SizeType numEntityDofs(const Entity &entity) const
return number of dofs for given element
Definition: blockmapper.hh:224
typename BaseType::ElementType ElementType
element type
Definition: blockmapper.hh:115
void insertNewEntity(const GridElementType &gridElement)
add DOFs for new element
Definition: blockmapper.hh:357
bool adapt()
please doc me
Definition: blockmapper.hh:317
void insertEntity(const GridElementType &gridElement)
add DOFs for element
Definition: blockmapper.hh:337
static constexpr bool fixedDataSize(int codim)
return true if number of dofs is fixed for given codimension
Definition: blockmapper.hh:233
const DofManagerType & dofManager() const
return DOF manager
Definition: blockmapper.hh:334
DiscontinuousGalerkinBlockMapper(const GridPartType &gridPart, const LocalKeysType &localKeys, const KeyType &value, Function function)
Definition: blockmapper.hh:143
KeyType getMark(const ElementType &element) const
get key to be assigned to an entity after next call to adapt()
Definition: blockmapper.hh:306
bool compress()
compress DOF mapping
Definition: blockmapper.hh:589
DiscontinuousGalerkinBlockMapper(ThisType &&)=default
move constructor
typename BaseType::GlobalKeyType GlobalKeyType
global key type
Definition: blockmapper.hh:110
void read(InStreamInterface< Traits > &)
this mapper has no I/O capabilities
Definition: blockmapper.hh:410
SizeType numDofs(const ElementType &element) const
return number of dofs for given element
Definition: blockmapper.hh:217
int maxNumDofs() const
return upper bound for number of dofs
Definition: blockmapper.hh:214
SizeType numberOfHoles(const int block) const
return number of holes during compression
Definition: blockmapper.hh:250
static constexpr bool contains(const int codim)
return true if dofs are associated to codimension
Definition: blockmapper.hh:230
void removeEntity(const GridElementType &gridElement)
mark DOFs for removal
Definition: blockmapper.hh:346
std::enable_if<(std::is_same< Element, ElementType >::value||std::is_same< Element, GridElementType >::value), const KeyType & >::type key(const Element &element) const
get key currently assigned to an entity
Definition: blockmapper.hh:293
GlobalKeyType oldIndex(const int hole, const int block) const
return old index of given hole during compression
Definition: blockmapper.hh:257
static constexpr SizeType offSet(const int block)
return 0 (this mapper has no offset)
Definition: blockmapper.hh:279
void resize()
Definition: blockmapper.hh:380
void restore()
this mapper has no I/O capabilities
Definition: blockmapper.hh:422
void mark(const KeyType &key, const ElementType &element)
set key to be assigned to an entity after next call to adapt()
Definition: blockmapper.hh:299
static constexpr SizeType numBlocks()
return 1 (this mapper has one block)
Definition: blockmapper.hh:282
typename BaseType::SizeType SizeType
size type
Definition: blockmapper.hh:108
Reserve(SizeType &size)
Definition: blockmapper.hh:495
Resize(std::vector< GlobalKeyType > &holes)
Definition: blockmapper.hh:513
Definition: space/hpdg/datahandle.hh:32
Definition: localdofstorage.hh:25
Definition: space/mapper/capabilities.hh:22
static const bool v
Definition: space/mapper/capabilities.hh:23
Traits::ElementType ElementType
type of codimension 0 entities
Definition: mapper/dofmapper.hh:54
Extended interface for adaptive DoF mappers.
Definition: mapper/dofmapper.hh:219
SizeType GlobalKeyType
at the moment this should be similar to SizeType
Definition: mapper/dofmapper.hh:230
BaseType::SizeType SizeType
type of size integer
Definition: mapper/dofmapper.hh:227