dune-fem  2.8-git
auxiliarydofs.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SPACE_COMMON_AUXILIARYDOFS_HH
2 #define DUNE_FEM_SPACE_COMMON_AUXILIARYDOFS_HH
3 
4 #include <cassert>
5 #include <cstddef>
6 
7 #include <algorithm>
8 #include <iostream>
9 #include <iterator>
10 #include <limits>
11 #include <map>
12 #include <memory>
13 #include <set>
14 
15 #include <dune/common/exceptions.hh>
16 #include <dune/common/genericiterator.hh>
17 #include <dune/common/ftraits.hh>
18 #include <dune/common/typetraits.hh>
19 
20 #include <dune/grid/common/gridenums.hh>
21 #include <dune/grid/common/datahandleif.hh>
22 
27 
28 namespace Dune
29 {
30 
31  namespace Fem
32  {
33 
45  template< class GridPart, class Mapper >
47  {
49 
50  class LinkBuilder;
51 
52  public:
54  typedef GridPart GridPartType;
55 
57  typedef Mapper MapperType;
58 
59  protected:
60  typedef Fem :: CommunicationIndexMap IndexMapType;
61 
64 
65  // type of communication indices
67 
68  public:
70  : public std::iterator< std::forward_iterator_tag, const int, int >
71  {
72  ConstIterator () = default;
73  ConstIterator ( const IndexMapType &auxiliarys, int index ) : auxiliarys_( &auxiliarys ), index_( index ) {}
74 
75  const int &operator* () const { return (*auxiliarys_)[ index_ ]; }
76  const int *operator-> () const { return &(*auxiliarys_)[ index_ ]; }
77 
78  const int &operator[] ( int n ) const noexcept { return (*auxiliarys_)[ index_ + n ]; }
79 
80  bool operator== ( const ConstIterator &other ) const { return (index_ == other.index_); }
81  bool operator!= ( const ConstIterator &other ) const { return (index_ != other.index_); }
82 
83  ConstIterator &operator++ () { ++index_; return *this; }
84  ConstIterator operator++ ( int ) { ConstIterator copy( *this ); ++(*this); return copy; }
85 
86  ThisType &operator-- () noexcept { --index_; return *this; }
87  ThisType operator-- ( int ) noexcept { ThisType copy( *this ); --(*this); return copy; }
88 
89  ThisType &operator+= ( int n ) noexcept { index_ += n; return *this; }
90  ThisType &operator-= ( int n ) noexcept { index_ -= n; return *this; }
91 
92  ThisType operator+ ( int n ) const noexcept { return ThisType( index_ + n ); }
93  ThisType operator- ( int n ) const noexcept { return ThisType( index_ - n ); }
94 
95  friend ThisType operator+ ( int n, const ThisType &i ) noexcept { return i + n; }
96 
97  int operator- ( const ThisType &other ) const noexcept { return (index_ - other.index_); }
98 
99  bool operator< ( const ThisType &other ) const noexcept { return (index_ < other.index_); }
100  bool operator<= ( const ThisType &other ) const noexcept { return (index_ <= other.index_); }
101  bool operator>= ( const ThisType &other ) const noexcept { return (index_ >= other.index_); }
102  bool operator> ( const ThisType &other ) const noexcept { return (index_ > other.index_); }
103 
104  private:
105  const IndexMapType *auxiliarys_ = nullptr;
106  int index_ = 0;
107  };
108 
109  AuxiliaryDofs ( const GridPartType &gridPart, const MapperType &mapper )
110  : gridPart_( gridPart ), mapper_( mapper )
111  {}
112 
113  AuxiliaryDofs ( const AuxiliaryDofs& ) = delete;
114 
116  int operator [] ( const int index ) const
117  {
118  return auxiliarys_[ index ];
119  }
120 
122  int size () const
123  {
124  return auxiliarys_.size();
125  }
126 
127  ConstIterator begin () const { return ConstIterator( auxiliarys_, 0 ); }
128  ConstIterator end () const { assert( size() > 0 ); return ConstIterator( auxiliarys_, size()-1 ); }
129 
131  bool contains ( int index ) const { return std::binary_search( begin(), end(), index ); }
132 
133  [[deprecated("Use contains instead")]]
134  bool isSlave ( int index ) const { return contains( index ); }
135 
136  void rebuild ()
137  {
138  std::set< int > auxiliarySet;
139  buildMaps( auxiliarySet );
140  auxiliarys_.clear();
141  auxiliarys_.set( auxiliarySet );
142  }
143 
144  const GridPartType &gridPart () const { return gridPart_; }
145 
146  protected:
147  void buildMaps ( std::set< int > &auxiliarySet )
148  {
149  // build linkage and index maps
150  for( int codim = 1; codim <= GridPartType::dimension; ++codim )
151  {
152  if( mapper_.contains( codim ) )
153  return buildCommunicatedMaps( auxiliarySet );
154  }
155  return buildDiscontinuousMaps( auxiliarySet );
156  }
157 
158  void buildDiscontinuousMaps ( std::set< int > &auxiliarySet )
159  {
160  // if DoFs are only attached to codimension 0, we do not have to communicate
161  const auto idxpitype = GridPartType::indexSetPartitionType;
162  for( auto it = gridPart().template begin< 0, idxpitype >(), end = gridPart().template end< 0, idxpitype >(); it != end; ++it )
163  {
164  const auto& entity = *it;
165  if( entity.partitionType() != Dune::InteriorEntity )
166  mapper_.mapEachEntityDof( entity, [ &auxiliarySet ] ( int, int value ) { auxiliarySet.insert( value ); } );
167  }
168  // insert overall size at the end
169  auxiliarySet.insert( mapper_.size() );
170  }
171 
172  void buildCommunicatedMaps ( std::set< int > &auxiliarySet )
173  {
174  // we have to skip communication when parallel program is executed only on one processor
175  // otherwise YaspGrid and Lagrange polorder=2 fails :(
176  if( gridPart().comm().size() > 1 )
177  {
178  try
179  {
180  LinkBuilder handle( auxiliarySet, gridPart(), mapper_ );
181  gridPart().communicate( handle, GridPartType::indexSetInterfaceType, ForwardCommunication );
182  }
183  catch( const Exception &e )
184  {
185  std::cerr << e << std::endl << "Exception thrown in: " << __FILE__ << " line:" << __LINE__ << std::endl;
186  std::abort();
187  }
188  }
189  // insert overall size at the end
190  auxiliarySet.insert( mapper_.size() );
191  }
192  };
193 
194 
195 
196  // AuxiliaryDofs::LinkBuilder
197  // ----------------------
198 
199  template< class GridPart, class Mapper >
200  class AuxiliaryDofs< GridPart, Mapper >::LinkBuilder
201  : public CommDataHandleIF< LinkBuilder, int >
202  {
203  public:
204  LinkBuilder( std::set< int > &auxiliarySet, const GridPartType &gridPart, const MapperType &mapper )
205  : myRank_( gridPart.comm().rank() ), mySize_( gridPart.comm().size() ),
206  auxiliarySet_( auxiliarySet ), mapper_( mapper )
207  {}
208 
209  bool contains ( int dim, int codim ) const { return mapper_.contains( codim ); }
210  bool fixedSize ( int dim, int codim ) const { return false; }
211 
213  template< class MessageBuffer, class Entity >
214  void gather ( MessageBuffer &buffer, const Entity &entity ) const
215  {
216  // for sending ranks write rank
217  if( sendRank( entity ) )
218  buffer.write( myRank_ );
219  }
220 
225  template< class MessageBuffer, class EntityType >
226  void scatter ( MessageBuffer &buffer, const EntityType &entity, std::size_t n )
227  {
228  if( n == 1 )
229  {
230  int rank;
231  buffer.read( rank );
232 
233  assert( (rank >= 0) && (rank < mySize_) );
234 
235  // if entity in not interiorBorder insert anyway
236  if ( rank < myRank_ || ! sendRank( entity ) )
237  mapper_.mapEachEntityDof( entity, [this]( const int , const auto& value ){auxiliarySet_.insert( value );} );
238  }
239  }
240 
242  template< class Entity >
243  std::size_t size ( const Entity &entity ) const
244  {
245  return (sendRank( entity )) ? 1 : 0;
246  }
247 
248  protected:
249  template <class Entity>
250  bool sendRank(const Entity& entity) const
251  {
252  const PartitionType ptype = entity.partitionType();
253  return (ptype == InteriorEntity) || (ptype == BorderEntity);
254  }
255 
256  private:
257  int myRank_, mySize_;
258  std::set< int > &auxiliarySet_;
259  const MapperType &mapper_;
260  };
261 
262 
263 
264  // PrimaryDofs
265  // -----------
266 
273  template< class AuxiliaryDofs >
274  struct PrimaryDofs;
275 
276  template< class GridPart, class Mapper >
277  struct PrimaryDofs< AuxiliaryDofs< GridPart, Mapper > >
278  {
280 
281  struct ConstIterator
282  : public std::iterator< std::forward_iterator_tag, int, std::ptrdiff_t, Envelope< int >, int >
283  {
284  ConstIterator () = default;
285 
286  ConstIterator ( int index, int auxiliary )
287  : index_( index ), auxiliary_( auxiliary )
288  {}
289 
290  ConstIterator ( const AuxiliaryDofsType &auxiliaryDofs, int index, int auxiliary )
291  : auxiliaryDofs_( &auxiliaryDofs ), index_( index ), auxiliary_( auxiliary )
292  {
293  skipAuxiliarys();
294  }
295 
296  int operator* () const { return index_; }
297  Envelope< int > operator-> () const { return Envelope< int >( index_ ); }
298 
299  bool operator== ( const ConstIterator &other ) const { return (index_ == other.index_); }
300  bool operator!= ( const ConstIterator &other ) const { return (index_ != other.index_); }
301 
302  ConstIterator &operator++ () { ++index_; skipAuxiliarys(); return *this; }
303  ConstIterator operator++ ( int ) { ConstIterator copy( *this ); ++(*this); return copy; }
304 
305  const AuxiliaryDofsType &auxiliaryDofs () const { assert( auxiliaryDofs_ ); return *auxiliaryDofs_; }
306 
307  bool contains( const int index ) const { return ! auxiliaryDofs().contains( index ); }
308 
309  private:
310  void skipAuxiliarys ()
311  {
312  const int aSize = auxiliaryDofs().size();
313  assert( auxiliary_ < aSize );
314  for( ; (index_ == auxiliaryDofs()[ auxiliary_ ]) && (++auxiliary_ != aSize); ++index_ )
315  continue;
316  }
317 
318  const AuxiliaryDofsType *auxiliaryDofs_ = nullptr;
319  int index_ = 0, auxiliary_ = 0;
320  };
321 
322  explicit PrimaryDofs ( const AuxiliaryDofsType &auxiliaryDofs )
323  : auxiliaryDofs_( auxiliaryDofs )
324  {}
325 
326  ConstIterator begin () const { return ConstIterator( auxiliaryDofs_, 0, 0 ); }
327  ConstIterator end () const { return ConstIterator( auxiliaryDofs_[ auxiliaryDofs_.size()-1 ], auxiliaryDofs_.size() ); }
328 
329  int size () const { return auxiliaryDofs_[ auxiliaryDofs_.size()-1 ] - (auxiliaryDofs_.size()-1); }
330 
331  private:
332  const AuxiliaryDofsType &auxiliaryDofs_;
333  };
334 
335 
336 
337  // primaryDofs
338  // -----------
339 
340  template< class AuxiliaryDofs >
341  inline static PrimaryDofs< AuxiliaryDofs > primaryDofs ( const AuxiliaryDofs &auxiliaryDofs )
342  {
343  return PrimaryDofs< AuxiliaryDofs >( auxiliaryDofs );
344  }
345 
346  template< class AuxiliaryDofs >
347  [[deprecated("Use primaryDofs instead!" )]]
348  inline static PrimaryDofs< AuxiliaryDofs > masterDofs ( const AuxiliaryDofs &auxiliaryDofs )
349  {
350  return PrimaryDofs< AuxiliaryDofs >( auxiliaryDofs );
351  }
352 
354 
355  } // end namespace Fem
356 
357 } // end namespace Dune
358 #endif // #ifndef DUNE_FEM_SPACE_COMMON_AUXILIARYDOFS_HH
GridPart GridPartType
type of grid part
Definition: auxiliarydofs.hh:50
const int * operator->() const
Definition: auxiliarydofs.hh:76
ConstIterator(const IndexMapType &auxiliarys, int index)
Definition: auxiliarydofs.hh:73
ThisType & operator--() noexcept
Definition: auxiliarydofs.hh:86
PrimaryDofs(const AuxiliaryDofsType &auxiliaryDofs)
Definition: auxiliarydofs.hh:322
ConstIterator end() const
Definition: auxiliarydofs.hh:128
ThisType operator+(int n) const noexcept
Definition: auxiliarydofs.hh:92
int operator[](const int index) const
return dof number of auxiliary for index
Definition: auxiliarydofs.hh:116
const AuxiliaryDofsType & auxiliaryDofs() const
Definition: auxiliarydofs.hh:305
static PrimaryDofs< AuxiliaryDofs > primaryDofs(const AuxiliaryDofs &auxiliaryDofs)
Definition: auxiliarydofs.hh:341
bool operator<(const ThisType &other) const noexcept
Definition: auxiliarydofs.hh:99
void rebuild()
Definition: auxiliarydofs.hh:136
bool operator!=(const ConstIterator &other) const
Definition: auxiliarydofs.hh:81
const GridPartType & gridPart() const
Definition: auxiliarydofs.hh:144
IndexMapType auxiliarys_
Definition: auxiliarydofs.hh:66
void buildMaps(std::set< int > &auxiliarySet)
Definition: auxiliarydofs.hh:147
ConstIterator begin() const
Definition: auxiliarydofs.hh:127
ThisType & operator+=(int n) noexcept
Definition: auxiliarydofs.hh:89
void buildCommunicatedMaps(std::set< int > &auxiliarySet)
Definition: auxiliarydofs.hh:172
void buildDiscontinuousMaps(std::set< int > &auxiliarySet)
Definition: auxiliarydofs.hh:158
bool contains(int index) const
return true if index is contained, meaning it is a auxiliary dof
Definition: auxiliarydofs.hh:131
bool contains(const int index) const
Definition: auxiliarydofs.hh:307
const int & operator*() const
Definition: auxiliarydofs.hh:75
std::size_t size(const Entity &entity) const
return local dof size to be communicated
Definition: auxiliarydofs.hh:243
Mapper MapperType
type of used mapper
Definition: auxiliarydofs.hh:57
ThisType & operator-=(int n) noexcept
Definition: auxiliarydofs.hh:90
const MapperType & mapper_
Definition: auxiliarydofs.hh:63
AuxiliaryDofs(const GridPartType &gridPart, const MapperType &mapper)
Definition: auxiliarydofs.hh:109
ConstIterator & operator++()
Definition: auxiliarydofs.hh:83
void scatter(MessageBuffer &buffer, const EntityType &entity, std::size_t n)
Definition: auxiliarydofs.hh:226
Fem ::CommunicationIndexMap IndexMapType
Definition: auxiliarydofs.hh:60
bool fixedSize(int dim, int codim) const
Definition: auxiliarydofs.hh:210
bool isSlave(int index) const
Definition: auxiliarydofs.hh:134
bool operator==(const ConstIterator &other) const
Definition: auxiliarydofs.hh:80
bool operator<=(const ThisType &other) const noexcept
Definition: auxiliarydofs.hh:100
LinkBuilder(std::set< int > &auxiliarySet, const GridPartType &gridPart, const MapperType &mapper)
Definition: auxiliarydofs.hh:204
bool contains(int dim, int codim) const
Definition: auxiliarydofs.hh:209
AuxiliaryDofs(const AuxiliaryDofs &)=delete
int size() const
Definition: auxiliarydofs.hh:329
ConstIterator end() const
Definition: auxiliarydofs.hh:327
bool sendRank(const Entity &entity) const
Definition: auxiliarydofs.hh:250
int size() const
return number of auxiliary dofs
Definition: auxiliarydofs.hh:122
ThisType operator-(int n) const noexcept
Definition: auxiliarydofs.hh:93
static PrimaryDofs< AuxiliaryDofs > masterDofs(const AuxiliaryDofs &auxiliaryDofs)
Definition: auxiliarydofs.hh:348
bool operator>(const ThisType &other) const noexcept
Definition: auxiliarydofs.hh:102
const GridPartType & gridPart_
Definition: auxiliarydofs.hh:62
ConstIterator(int index, int auxiliary)
Definition: auxiliarydofs.hh:286
AuxiliaryDofs< GridPart, Mapper > AuxiliaryDofsType
Definition: auxiliarydofs.hh:279
ConstIterator begin() const
Definition: auxiliarydofs.hh:326
ConstIterator(const AuxiliaryDofsType &auxiliaryDofs, int index, int auxiliary)
Definition: auxiliarydofs.hh:290
const int & operator[](int n) const noexcept
Definition: auxiliarydofs.hh:78
void gather(MessageBuffer &buffer, const Entity &entity) const
read buffer and apply operation
Definition: auxiliarydofs.hh:214
bool operator>=(const ThisType &other) const noexcept
Definition: auxiliarydofs.hh:101
Definition: bindguard.hh:11
Double operator*(const Double &a, const Double &b)
Definition: double.hh:506
bool operator==(const Double &a, const Double &b)
Definition: double.hh:600
bool operator!=(const Double &a, const Double &b)
Definition: double.hh:640
In parallel computations the dofs of a discrete function are made up by all primary dofs....
Definition: auxiliarydofs.hh:47
Definition: auxiliarydofs.hh:71
Definition: auxiliarydofs.hh:202
In parallel computations the dofs of a discrete function are made up by all primary dofs....
Definition: auxiliarydofs.hh:274
Definition: envelope.hh:11