dune-fem 2.8.0
Loading...
Searching...
No Matches
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 <unordered_map>
10
11#include <dune/grid/common/gridenums.hh>
12#include <dune/grid/common/rangegenerators.hh>
15
16namespace 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 typename std::vector< std::size_t > :: size_type IndexType;
55
56 static const bool indexIsSimple = std::is_convertible< RangeGlobalKeyType, IndexType >::value;
57 //static const bool indexIsPOD = false ;//std::is_convertible< RangeGlobalKeyType, IndexType >::value;
58 typedef typename std::conditional< indexIsSimple,
59 std::unordered_map< RangeGlobalKeyType, LocalStencilType >,
60 std::map< RangeGlobalKeyType, LocalStencilType > > :: type GlobalStencilType;
61
62 public:
69 Stencil(const DomainSpace &dSpace, const RangeSpace &rSpace)
70 : domainSpace_(dSpace), rangeSpace_(rSpace)
71 , domainBlockMapper_( dSpace.blockMapper() )
72 , rangeBlockMapper_( rSpace.blockMapper() )
73 {
74 }
75
76 const DomainSpace &domainSpace() const
77 {
78 return domainSpace_;
79 }
80 const RangeSpace &rangeSpace() const
81 {
82 return rangeSpace_;
83 }
84
92 void fill ( const DomainEntityType &dEntity, const RangeEntityType &rEntity,
93 bool fillGhost=true ) const
94 {
95 if( (dEntity.partitionType() == GhostEntity) && !fillGhost )
96 return;
97
98 rangeBlockMapper_.mapEach( rEntity, [ this, &dEntity ] ( int localRow, auto globalRow ) {
99 domainBlockMapper_.mapEach( dEntity, RowFillFunctor( globalStencil_[ globalRow ] ) );
100 } );
101 }
102
109 {
110 return globalStencil()[ key ];
111 }
112
116 {
117 return globalStencil_;
118 }
119
123 {
124 int maxNZ = 0;
125 for( const auto& entry : globalStencil_ )
126 {
127 maxNZ = std::max( maxNZ, static_cast<int>( entry.second.size() ) );
128 }
129 return maxNZ;
130 }
131
132 int rows () const { return rangeBlockMapper_.size(); }
133 int cols () const { return domainBlockMapper_.size(); }
134
136 void update()
137 {
138 globalStencil_.clear();
139 // compute stencil based on overloaded implementation
140 setupStencil();
141 }
142
143 [[deprecated("Use Stencil::update instead()")]]
144 void setup() {
145 update();
146 }
147
148 protected:
150 virtual void setupStencil() const = 0;
151
152 private:
153 struct RowFillFunctor
154 {
155 explicit RowFillFunctor ( LocalStencilType &localStencil )
156 : localStencil_( localStencil )
157 {}
158
159 void operator() ( const std::size_t, const DomainGlobalKeyType &domainGlobal ) const
160 {
161 localStencil_.insert( domainGlobal );
162 }
163
164 private:
165 LocalStencilType &localStencil_;
166 };
167
168 protected:
169 const DomainSpace &domainSpace_;
170 const RangeSpace &rangeSpace_;
171 const DomainBlockMapper &domainBlockMapper_;
172 const RangeBlockMapper &rangeBlockMapper_;
173
175 };
176
184 template <class DomainSpace, class RangeSpace>
186 {
188 public:
195
196 SimpleStencil(int maxNZ)
197 : maxNZ_(maxNZ)
198 {}
200 {
201 maxNZ_ = 1;
202 }
204 {
205 return maxNZ_;
206 }
208 {
209 DUNE_THROW( Dune::NotImplemented, "SimpleStencil: exact stencil information is unavailable." );
210 return localStencil_;
211 }
213 {
214 DUNE_THROW( Dune::NotImplemented, "SimpleStencil: global stencil is unavailable." );
215 return stencil_;
216 }
217 protected:
221 };
222
223
231 template <class DomainSpace, class RangeSpace, class LocalStencil>
233 {
236 public:
241
242 static_assert( Std::is_pod< DomainGlobalKeyType > :: value, "StencilWrapper only works for POD DomainGlobalKeyType");
243
244 typedef LocalStencil LocalStencilType;
245 typedef std::vector< LocalStencilType > GlobalStencilType;
246
247 struct Iterator
248 {
249 typedef std::pair< DomainGlobalKeyType, const LocalStencilType& > ItemType;
250
253
254
255 Iterator( const DomainGlobalKeyType& init, const GlobalStencilType& stencil )
256 : index_( init ), stencil_( stencil ) {}
257
259 {
260 ++ index_ ;
261 return *this;
262 }
263
265 {
266 assert( index_ < stencil_.size() );
267 return ItemType( index_, stencil_[ index_ ] );
268 }
269
270 bool operator == ( const Iterator& other ) const
271 {
272 return index_ == other.index_;
273 }
274
275 bool operator != ( const Iterator& other ) const
276 {
277 return !this->operator==( other );
278 }
279 };
280
282 : stencil_( stencil )
283 , maxNZ_( computeMaxNZ() )
284 {
285 }
286
288 {
289 return maxNZ_;
290 }
291
293 {
294 return stencil_[ key ];
295 }
296
297 const ThisType& globalStencil() const
298 {
299 return *this;
300 }
301
309 void fill ( const DomainEntityType &dEntity, const RangeEntityType &rEntity,
310 bool fillGhost=true )
311 {
312 }
313
314 Iterator begin() const { return Iterator(0, stencil_); }
315 Iterator end() const { return Iterator(stencil_.size(), stencil_); }
317 {
318 assert( key < stencil_.size() );
319 return Iterator( key, stencil_);
320 }
321
322 protected:
323 int computeMaxNZ() const
324 {
325 int maxNZ = 0;
326 for( const auto& row : stencil_ )
327 {
328 maxNZ = std::max( maxNZ, int(row.size()) );
329 }
330 return maxNZ;
331 }
332
335 };
336
337
346 template <class DomainSpace, class RangeSpace, class Partition = Dune::Partitions::InteriorBorder>
347 struct DiagonalStencil : public Stencil<DomainSpace,RangeSpace>
348 {
350 public:
351 typedef Partition PartitionType;
358
359 DiagonalStencil(const DomainSpace &dSpace, const RangeSpace &rSpace)
360 : BaseType( dSpace, rSpace )
361 {
362 setupStencil();
363 }
364
365 protected:
366 virtual void setupStencil () const override
367 {
368 const DomainSpace &dSpace = BaseType::domainSpace();
369 for (const auto& entity : elements( dSpace.gridPart(), PartitionType{} ) )
370 BaseType::fill(entity,entity);
371 }
372 };
373
383 template <class DomainSpace, class RangeSpace, class Partition = Dune::Partitions::InteriorBorder>
384 struct DiagonalAndNeighborStencil : public Stencil<DomainSpace,RangeSpace>
385 {
387 public:
388 typedef Partition PartitionType;
395
396 DiagonalAndNeighborStencil(const DomainSpace &dSpace, const RangeSpace &rSpace,
397 bool onlyNonContinuousNeighbors = false)
398 : BaseType( dSpace, rSpace ),
399 onlyNonContinuousNeighbors_(onlyNonContinuousNeighbors)
400 {
401 setupStencil();
402 }
403
404 protected:
405 virtual void setupStencil() const override
406 {
407 const DomainSpace &dSpace = BaseType::domainSpace();
408 const RangeSpace &rSpace = BaseType::rangeSpace();
409 for (const auto & entity: elements( dSpace.gridPart(), PartitionType{} ) )
410 {
411 BaseType::fill(entity,entity);
412 for (const auto & intersection: intersections(dSpace.gridPart(), entity) )
413 {
415 && rSpace.continuous(intersection) && dSpace.continuous(intersection) )
416 continue;
417 if( intersection.neighbor() )
418 {
419 auto neighbor = intersection.outside();
420 BaseType::fill(neighbor,entity);
421 }
422 }
423 }
424 }
425
427 };
428
429 } // namespace Fem
430
431} // namespace Dune
432
433#endif // #if defined DUNE_FEM_STENCIL_HH
434
double max(const Dune::Fem::Double &v, const double p)
Definition: double.hh:965
Definition: bindguard.hh:11
Definition: utility.hh:162
default implementation for a general operator stencil
Definition: stencil.hh:35
const DomainSpace & domainSpace_
Definition: stencil.hh:169
static const bool indexIsSimple
Definition: stencil.hh:56
const RangeBlockMapper & rangeBlockMapper_
Definition: stencil.hh:172
GlobalStencilType globalStencil_
Definition: stencil.hh:174
int rows() const
Definition: stencil.hh:132
RangeIteratorType::Entity RangeEntityType
Definition: stencil.hh:46
void update()
clear previously computed entries such that a re-compute happens when used again
Definition: stencil.hh:136
DomainBlockMapper::GlobalKeyType DomainGlobalKeyType
Definition: stencil.hh:47
const RangeSpace & rangeSpace_
Definition: stencil.hh:170
void setup()
Definition: stencil.hh:144
const RangeSpace & rangeSpace() const
Definition: stencil.hh:80
virtual void setupStencil() const =0
method to setup stencil depending on entity set defined in derived class
RangeBlockMapper::GlobalKeyType RangeGlobalKeyType
Definition: stencil.hh:48
const DomainSpace & domainSpace() const
Definition: stencil.hh:76
const GlobalStencilType & globalStencil() const
Return the full stencil.
Definition: stencil.hh:115
DomainIteratorType::Entity DomainEntityType
Definition: stencil.hh:45
std::set< DomainGlobalKeyType > LocalStencilType
type for storing the stencil of one row
Definition: stencil.hh:51
void fill(const DomainEntityType &dEntity, const RangeEntityType &rEntity, bool fillGhost=true) const
Create stencil entries for (dEntity,rEntity) pair.
Definition: stencil.hh:92
std::vector< std::size_t >::size_type IndexType
type of std::vector for indexing
Definition: stencil.hh:54
std::conditional< indexIsSimple, std::unordered_map< RangeGlobalKeyType, LocalStencilType >, std::map< RangeGlobalKeyType, LocalStencilType > >::type GlobalStencilType
Definition: stencil.hh:60
int maxNonZerosEstimate() const
Return an upper bound for the maximum number of non-zero entries in all rows.
Definition: stencil.hh:122
Stencil(const DomainSpace &dSpace, const RangeSpace &rSpace)
Constructor.
Definition: stencil.hh:69
const LocalStencilType & localStencil(const RangeGlobalKeyType &key) const
Return stencil for a given row of the matrix.
Definition: stencil.hh:108
const DomainBlockMapper & domainBlockMapper_
Definition: stencil.hh:171
int cols() const
Definition: stencil.hh:133
a watered down stencil providing only the upper bound for the non-zero entries per row.
Definition: stencil.hh:186
StencilType::RangeGlobalKeyType RangeGlobalKeyType
Definition: stencil.hh:192
int maxNZ_
Definition: stencil.hh:218
StencilType::RangeEntityType RangeEntityType
Definition: stencil.hh:190
LocalStencilType localStencil_
Definition: stencil.hh:220
const LocalStencilType & localStencil(const DomainGlobalKeyType &key) const
Definition: stencil.hh:207
GlobalStencilType stencil_
Definition: stencil.hh:219
SimpleStencil()
Definition: stencil.hh:199
int maxNonZerosEstimate() const
Definition: stencil.hh:203
SimpleStencil(int maxNZ)
Definition: stencil.hh:196
const GlobalStencilType & globalStencil() const
Definition: stencil.hh:212
StencilType::DomainEntityType DomainEntityType
Definition: stencil.hh:189
StencilType::GlobalStencilType GlobalStencilType
Definition: stencil.hh:194
StencilType::LocalStencilType LocalStencilType
Definition: stencil.hh:193
StencilType::DomainGlobalKeyType DomainGlobalKeyType
Definition: stencil.hh:191
a simple wrapper class for sparsity patterns provide as vector< set< size_t > >
Definition: stencil.hh:233
int computeMaxNZ() const
Definition: stencil.hh:323
void fill(const DomainEntityType &dEntity, const RangeEntityType &rEntity, bool fillGhost=true)
Create stencil entries for (dEntity,rEntity) pair.
Definition: stencil.hh:309
int maxNZ_
Definition: stencil.hh:334
StencilType::DomainEntityType DomainEntityType
Definition: stencil.hh:237
std::vector< LocalStencilType > GlobalStencilType
Definition: stencil.hh:245
const LocalStencilType & localStencil(const DomainGlobalKeyType &key) const
Definition: stencil.hh:292
const ThisType & globalStencil() const
Definition: stencil.hh:297
StencilType::RangeGlobalKeyType RangeGlobalKeyType
Definition: stencil.hh:240
StencilType::DomainGlobalKeyType DomainGlobalKeyType
Definition: stencil.hh:239
LocalStencil LocalStencilType
Definition: stencil.hh:244
Iterator find(const DomainGlobalKeyType &key) const
Definition: stencil.hh:316
StencilWrapper(const GlobalStencilType &stencil)
Definition: stencil.hh:281
StencilType::RangeEntityType RangeEntityType
Definition: stencil.hh:238
Iterator begin() const
Definition: stencil.hh:314
int maxNonZerosEstimate() const
Definition: stencil.hh:287
const GlobalStencilType & stencil_
Definition: stencil.hh:333
Iterator end() const
Definition: stencil.hh:315
Definition: stencil.hh:248
bool operator==(const Iterator &other) const
Definition: stencil.hh:270
const GlobalStencilType & stencil_
Definition: stencil.hh:252
bool operator!=(const Iterator &other) const
Definition: stencil.hh:275
Iterator(const DomainGlobalKeyType &init, const GlobalStencilType &stencil)
Definition: stencil.hh:255
ItemType operator*() const
Definition: stencil.hh:264
std::pair< DomainGlobalKeyType, const LocalStencilType & > ItemType
Definition: stencil.hh:249
DomainGlobalKeyType index_
Definition: stencil.hh:251
Iterator & operator++()
Definition: stencil.hh:258
Stencil contaning the entries (en,en) for all entities in the space. Defailt for an operator over a L...
Definition: stencil.hh:348
BaseType::DomainGlobalKeyType DomainGlobalKeyType
Definition: stencil.hh:354
BaseType::GlobalStencilType GlobalStencilType
Definition: stencil.hh:357
BaseType::RangeEntityType RangeEntityType
Definition: stencil.hh:353
BaseType::DomainEntityType DomainEntityType
Definition: stencil.hh:352
DiagonalStencil(const DomainSpace &dSpace, const RangeSpace &rSpace)
Definition: stencil.hh:359
virtual void setupStencil() const override
method to setup stencil depending on entity set defined in derived class
Definition: stencil.hh:366
Partition PartitionType
Definition: stencil.hh:351
BaseType::LocalStencilType LocalStencilType
Definition: stencil.hh:356
BaseType::RangeGlobalKeyType RangeGlobalKeyType
Definition: stencil.hh:355
Stencil< DomainSpace, RangeSpace > BaseType
Definition: stencil.hh:349
Stencil contaning the entries (en,en) and (en,nb) for all entities en in the space and neighbors nb o...
Definition: stencil.hh:385
virtual void setupStencil() const override
method to setup stencil depending on entity set defined in derived class
Definition: stencil.hh:405
bool onlyNonContinuousNeighbors_
Definition: stencil.hh:426
BaseType::RangeEntityType RangeEntityType
Definition: stencil.hh:390
BaseType::DomainGlobalKeyType DomainGlobalKeyType
Definition: stencil.hh:391
DiagonalAndNeighborStencil(const DomainSpace &dSpace, const RangeSpace &rSpace, bool onlyNonContinuousNeighbors=false)
Definition: stencil.hh:396
Stencil< DomainSpace, RangeSpace > BaseType
Definition: stencil.hh:386
BaseType::LocalStencilType LocalStencilType
Definition: stencil.hh:393
BaseType::GlobalStencilType GlobalStencilType
Definition: stencil.hh:394
Partition PartitionType
Definition: stencil.hh:388
BaseType::DomainEntityType DomainEntityType
Definition: stencil.hh:389
BaseType::RangeGlobalKeyType RangeGlobalKeyType
Definition: stencil.hh:392