dune-fem  2.8-git
stencil.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_STENCIL_HH
2 #define DUNE_FEM_STENCIL_HH
3 
4 #include <algorithm>
5 #include <map>
6 #include <numeric>
7 #include <set>
8 #include <type_traits>
9 #include <vector>
10 
11 #include <dune/grid/common/gridenums.hh>
12 #include <dune/grid/common/rangegenerators.hh>
13 #include <dune/fem/misc/functor.hh>
15 
16 namespace Dune
17 {
18  namespace Fem
19  {
33  template <class DomainSpace, class RangeSpace>
34  class Stencil
35  {
36  // Domain = Row
37  typedef typename DomainSpace::IteratorType DomainIteratorType;
38  typedef typename DomainSpace::BlockMapperType DomainBlockMapper;
39 
40  // Range = Column
41  typedef typename RangeSpace::IteratorType RangeIteratorType;
42  typedef typename RangeSpace::BlockMapperType RangeBlockMapper;
43 
44  public:
45  typedef typename DomainIteratorType::Entity DomainEntityType;
46  typedef typename RangeIteratorType::Entity RangeEntityType;
47  typedef typename DomainBlockMapper::GlobalKeyType DomainGlobalKeyType;
48  typedef typename RangeBlockMapper::GlobalKeyType RangeGlobalKeyType;
49 
51  typedef std::set<DomainGlobalKeyType> LocalStencilType;
52 
54  typedef std::map<RangeGlobalKeyType,LocalStencilType> GlobalStencilType;
55 
56  public:
63  Stencil(const DomainSpace &dSpace, const RangeSpace &rSpace)
64  : domainBlockMapper_( dSpace.blockMapper() )
65  , rangeBlockMapper_( rSpace.blockMapper() )
66  {
67  }
68 
76  void fill ( const DomainEntityType &dEntity, const RangeEntityType &rEntity,
77  bool fillGhost=true )
78  {
79  if( (dEntity.partitionType() == GhostEntity) && !fillGhost )
80  return;
81 
82  rangeBlockMapper_.mapEach( rEntity, [ this, &dEntity ] ( int localRow, auto globalRow ) {
83  domainBlockMapper_.mapEach( dEntity, RowFillFunctor( globalStencil_[ globalRow ] ) );
84  } );
85  }
86 
93  {
94  return globalStencil_[key];
95  }
99  {
100  return globalStencil_;
101  }
105  {
106  return std::accumulate( globalStencil().begin(), globalStencil().end(), 0,
107  []( int ret, const auto& entry ){ return std::max( ret, static_cast<int>( entry.second.size() ) ); } );
108  }
109 
110  int rows () const { return rangeBlockMapper_.size(); }
111  int cols () const { return domainBlockMapper_.size(); }
112 
113  private:
114  struct RowFillFunctor
115  {
116  explicit RowFillFunctor ( LocalStencilType &localStencil )
117  : localStencil_( localStencil )
118  {}
119 
120  void operator() ( const std::size_t, const DomainGlobalKeyType &domainGlobal ) const
121  {
122  localStencil_.insert( domainGlobal );
123  }
124 
125  private:
126  LocalStencilType &localStencil_;
127  };
128 
129  protected:
130  const DomainBlockMapper &domainBlockMapper_;
131  const RangeBlockMapper &rangeBlockMapper_;
133  };
134 
142  template <class DomainSpace, class RangeSpace>
144  {
146  public:
153 
154  SimpleStencil(int maxNZ)
155  : maxNZ_(maxNZ)
156  {}
158  {
159  maxNZ_ = 1;
160  }
162  {
163  return maxNZ_;
164  }
166  {
167  DUNE_THROW( Dune::NotImplemented, "SimpleStencil: exact stencil information is unavailable." );
168  return localStencil_;
169  }
171  {
172  DUNE_THROW( Dune::NotImplemented, "SimpleStencil: global stencil is unavailable." );
173  return stencil_;
174  }
175  protected:
176  int maxNZ_;
179  };
180 
181 
189  template <class DomainSpace, class RangeSpace, class LocalStencil>
191  {
194  public:
199 
200  static_assert( Std::is_pod< DomainGlobalKeyType > :: value, "StencilWrapper only works for POD DomainGlobalKeyType");
201 
202  typedef LocalStencil LocalStencilType;
203  typedef std::vector< LocalStencilType > GlobalStencilType;
204 
205  struct Iterator
206  {
207  typedef std::pair< DomainGlobalKeyType, const LocalStencilType& > ItemType;
208 
211 
212 
213  Iterator( const DomainGlobalKeyType& init, const GlobalStencilType& stencil )
214  : index_( init ), stencil_( stencil ) {}
215 
217  {
218  ++ index_ ;
219  return *this;
220  }
221 
223  {
224  assert( index_ < stencil_.size() );
225  return ItemType( index_, stencil_[ index_ ] );
226  }
227 
228  bool operator == ( const Iterator& other ) const
229  {
230  return index_ == other.index_;
231  }
232 
233  bool operator != ( const Iterator& other ) const
234  {
235  return !this->operator==( other );
236  }
237  };
238 
240  : stencil_( stencil )
241  , maxNZ_( computeNNZ() )
242  {
243  }
244 
246  {
247  return maxNZ_;
248  }
249 
251  {
252  return stencil_[ key ];
253  }
254 
255  const ThisType& globalStencil() const
256  {
257  return *this;
258  }
259 
267  void fill ( const DomainEntityType &dEntity, const RangeEntityType &rEntity,
268  bool fillGhost=true )
269  {
270  }
271 
272  Iterator begin() const { return Iterator(0, stencil_); }
273  Iterator end() const { return Iterator(stencil_.size(), stencil_); }
274  Iterator find( const DomainGlobalKeyType &key) const
275  {
276  assert( key < stencil_.size() );
277  return Iterator( key, stencil_);
278  }
279 
280  protected:
281  int computeNNZ() const
282  {
283  int maxNNZ = 0;
284  for( const auto& row : stencil_ )
285  {
286  maxNNZ = std::max( maxNNZ, int(row.size()) );
287  }
288  return maxNNZ;
289  }
290 
292  int maxNZ_;
293  };
294 
295 
304  template <class DomainSpace, class RangeSpace, class Partition = Dune::Partitions::InteriorBorder>
305  struct DiagonalStencil : public Stencil<DomainSpace,RangeSpace>
306  {
308  public:
309  typedef Partition PartitionType;
316 
317  DiagonalStencil(const DomainSpace &dSpace, const RangeSpace &rSpace)
318  : BaseType( dSpace, rSpace )
319  {
320  for (const auto& entity : elements( dSpace.gridPart(), PartitionType{} ) )
321  BaseType::fill(entity,entity);
322  }
323  };
324 
334  template <class DomainSpace, class RangeSpace, class Partition = Dune::Partitions::InteriorBorder>
335  struct DiagonalAndNeighborStencil : public Stencil<DomainSpace,RangeSpace>
336  {
338  public:
339  typedef Partition PartitionType;
346 
347  DiagonalAndNeighborStencil(const DomainSpace &dSpace, const RangeSpace &rSpace,
348  bool onlyNonContinuousNeighbors = false)
349  : BaseType( dSpace, rSpace )
350  {
351  for (const auto & entity: elements( dSpace.gridPart(), PartitionType{} ) )
352  {
353  BaseType::fill(entity,entity);
354  for (const auto & intersection: intersections(dSpace.gridPart(), entity) )
355  {
356  if ( onlyNonContinuousNeighbors
357  && rSpace.continuous(intersection) && dSpace.continuous(intersection) )
358  continue;
359  if( intersection.neighbor() )
360  {
361  auto neighbor = intersection.outside();
362  BaseType::fill(neighbor,entity);
363  }
364  }
365  }
366  }
367  };
368 
369  } // namespace Fem
370 
371 } // namespace Dune
372 
373 #endif // #if defined DUNE_FEM_STENCIL_HH
374 
Definition: bindguard.hh:11
static constexpr T max(T a)
Definition: utility.hh:77
Definition: utility.hh:162
default implementation for a general operator stencil
Definition: stencil.hh:35
const RangeBlockMapper & rangeBlockMapper_
Definition: stencil.hh:131
GlobalStencilType globalStencil_
Definition: stencil.hh:132
int rows() const
Definition: stencil.hh:110
const LocalStencilType & localStencil(const RangeGlobalKeyType &key) const
Return stencil for a given row of the matrix.
Definition: stencil.hh:92
RangeIteratorType::Entity RangeEntityType
Definition: stencil.hh:46
DomainBlockMapper::GlobalKeyType DomainGlobalKeyType
Definition: stencil.hh:47
const GlobalStencilType & globalStencil() const
Return the full stencil.
Definition: stencil.hh:98
std::map< RangeGlobalKeyType, LocalStencilType > GlobalStencilType
type for storing the full stencil
Definition: stencil.hh:54
RangeBlockMapper::GlobalKeyType RangeGlobalKeyType
Definition: stencil.hh:48
DomainIteratorType::Entity DomainEntityType
Definition: stencil.hh:45
void fill(const DomainEntityType &dEntity, const RangeEntityType &rEntity, bool fillGhost=true)
Create stencil entries for (dEntity,rEntity) pair.
Definition: stencil.hh:76
int maxNonZerosEstimate() const
Return an upper bound for the maximum number of non-zero entries in all row.
Definition: stencil.hh:104
Stencil(const DomainSpace &dSpace, const RangeSpace &rSpace)
Constructor.
Definition: stencil.hh:63
const DomainBlockMapper & domainBlockMapper_
Definition: stencil.hh:130
int cols() const
Definition: stencil.hh:111
std::set< DomainGlobalKeyType > LocalStencilType
type for storing the stencil of one row
Definition: stencil.hh:51
a watered down stencil providing only the upper bound for the non-zero entries per row.
Definition: stencil.hh:144
const LocalStencilType & localStencil(const DomainGlobalKeyType &key) const
Definition: stencil.hh:165
StencilType::RangeGlobalKeyType RangeGlobalKeyType
Definition: stencil.hh:150
int maxNZ_
Definition: stencil.hh:176
StencilType::RangeEntityType RangeEntityType
Definition: stencil.hh:148
LocalStencilType localStencil_
Definition: stencil.hh:178
GlobalStencilType stencil_
Definition: stencil.hh:177
SimpleStencil()
Definition: stencil.hh:157
int maxNonZerosEstimate() const
Definition: stencil.hh:161
SimpleStencil(int maxNZ)
Definition: stencil.hh:154
const GlobalStencilType & globalStencil() const
Definition: stencil.hh:170
StencilType::DomainEntityType DomainEntityType
Definition: stencil.hh:147
StencilType::GlobalStencilType GlobalStencilType
Definition: stencil.hh:152
StencilType::LocalStencilType LocalStencilType
Definition: stencil.hh:151
StencilType::DomainGlobalKeyType DomainGlobalKeyType
Definition: stencil.hh:149
a simple wrapper class for sparsity patterns provide as vector< set< size_t > >
Definition: stencil.hh:191
void fill(const DomainEntityType &dEntity, const RangeEntityType &rEntity, bool fillGhost=true)
Create stencil entries for (dEntity,rEntity) pair.
Definition: stencil.hh:267
int maxNZ_
Definition: stencil.hh:292
StencilType::DomainEntityType DomainEntityType
Definition: stencil.hh:195
std::vector< LocalStencilType > GlobalStencilType
Definition: stencil.hh:203
StencilType::RangeGlobalKeyType RangeGlobalKeyType
Definition: stencil.hh:198
StencilType::DomainGlobalKeyType DomainGlobalKeyType
Definition: stencil.hh:197
int computeNNZ() const
Definition: stencil.hh:281
LocalStencil LocalStencilType
Definition: stencil.hh:200
Iterator find(const DomainGlobalKeyType &key) const
Definition: stencil.hh:274
StencilWrapper(const GlobalStencilType &stencil)
Definition: stencil.hh:239
StencilType::RangeEntityType RangeEntityType
Definition: stencil.hh:196
Iterator begin() const
Definition: stencil.hh:272
const LocalStencilType & localStencil(const DomainGlobalKeyType &key) const
Definition: stencil.hh:250
int maxNonZerosEstimate() const
Definition: stencil.hh:245
const GlobalStencilType & stencil_
Definition: stencil.hh:291
const ThisType & globalStencil() const
Definition: stencil.hh:255
Iterator end() const
Definition: stencil.hh:273
Definition: stencil.hh:206
bool operator==(const Iterator &other) const
Definition: stencil.hh:228
const GlobalStencilType & stencil_
Definition: stencil.hh:210
bool operator!=(const Iterator &other) const
Definition: stencil.hh:233
Iterator(const DomainGlobalKeyType &init, const GlobalStencilType &stencil)
Definition: stencil.hh:213
ItemType operator*() const
Definition: stencil.hh:222
std::pair< DomainGlobalKeyType, const LocalStencilType & > ItemType
Definition: stencil.hh:207
DomainGlobalKeyType index_
Definition: stencil.hh:209
Iterator & operator++()
Definition: stencil.hh:216
Stencil contaning the entries (en,en) for all entities in the space. Defailt for an operator over a L...
Definition: stencil.hh:306
BaseType::DomainGlobalKeyType DomainGlobalKeyType
Definition: stencil.hh:312
BaseType::GlobalStencilType GlobalStencilType
Definition: stencil.hh:315
BaseType::RangeEntityType RangeEntityType
Definition: stencil.hh:311
BaseType::DomainEntityType DomainEntityType
Definition: stencil.hh:310
DiagonalStencil(const DomainSpace &dSpace, const RangeSpace &rSpace)
Definition: stencil.hh:317
Partition PartitionType
Definition: stencil.hh:309
BaseType::LocalStencilType LocalStencilType
Definition: stencil.hh:314
BaseType::RangeGlobalKeyType RangeGlobalKeyType
Definition: stencil.hh:313
Stencil< DomainSpace, RangeSpace > BaseType
Definition: stencil.hh:307
Stencil contaning the entries (en,en) and (en,nb) for all entities en in the space and neighbors nb o...
Definition: stencil.hh:336
BaseType::RangeEntityType RangeEntityType
Definition: stencil.hh:341
BaseType::DomainGlobalKeyType DomainGlobalKeyType
Definition: stencil.hh:342
DiagonalAndNeighborStencil(const DomainSpace &dSpace, const RangeSpace &rSpace, bool onlyNonContinuousNeighbors=false)
Definition: stencil.hh:347
Stencil< DomainSpace, RangeSpace > BaseType
Definition: stencil.hh:337
BaseType::LocalStencilType LocalStencilType
Definition: stencil.hh:344
BaseType::GlobalStencilType GlobalStencilType
Definition: stencil.hh:345
Partition PartitionType
Definition: stencil.hh:339
BaseType::DomainEntityType DomainEntityType
Definition: stencil.hh:340
BaseType::RangeGlobalKeyType RangeGlobalKeyType
Definition: stencil.hh:343