dune-fem 2.8.0
Loading...
Searching...
No Matches
blockdiagonal.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_OPERATOR_LINEAR_BLOCKDIAGONAL_HH
2#define DUNE_FEM_OPERATOR_LINEAR_BLOCKDIAGONAL_HH
3
4// system includes
5#include <cassert>
6#include <string>
7#include <vector>
8
9// local includes
10#include <dune/common/fmatrix.hh>
17
18namespace Dune
19{
20
21 namespace Fem
22 {
23
25 template< class DiscreteFunctionSpace,
26 class LocalBlock = Dune::FieldMatrix< typename DiscreteFunctionSpace ::
27 RangeFieldType, DiscreteFunctionSpace::localBlockSize, DiscreteFunctionSpace::localBlockSize > >
29 : public Fem::AssembledOperator< AdaptiveDiscreteFunction< DiscreteFunctionSpace > >
30 {
33
34 public:
37
40
41 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType;
42 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType;
43
44 typedef typename DomainSpaceType::EntityType DomainEntityType;
45 typedef typename RangeSpaceType::EntityType RangeEntityType;
46
47 static const int localBlockSize = DomainSpaceType::localBlockSize;
48
49 typedef LocalBlock LocalBlockType;
50
51 // types needed for CommunicationManager to fake DiscreteFunction interface
54 typedef typename LocalBlockType::row_type DofType ;
56
57 template< class Operation >
59 {
60 typedef typename DiscreteFunctionSpaceType
63 };
64
65 private:
66
68 class LocalMatrix;
69 struct LocalMatrixFactory;
70
71 public:
74
76
77
78 BlockDiagonalLinearOperator ( const std::string &name,
81 : name_( name ),
83 localMatrixFactory_( *this ),
85 {
86 if( &domainSpace != &rangeSpace )
87 DUNE_THROW( InvalidStateException, "BlockDiagonalLinearOperator must be created with identical spaces." );
88 }
89
91 {
92 multiply( u, w );
93 }
94
95 template < class DomainSpace, class RangeSpace >
98 {
99 multiply( u, w );
100 }
101
102 template < class DomainSpace, class RangeSpace >
105 {
106 const auto uit = u.leakPointer();
107 auto wit = w.leakPointer();
108 for( auto& entry : diagonal_ )
109 {
110 entry.mv( uit, wit );
111 uit += entry.M();
112 wit += entry.N();
113 }
114 assert( uit == u.leakPointer() + u.size() );
115 assert( wit == w.leakPointer() + w.size() );
116 }
117
118 void clear ()
119 {
120 for( auto& entry : diagonal_ )
121 entry = RangeFieldType( 0 );
122 }
123
124 template< class Functor >
125 void forEach ( const Functor &functor )
126 {
127 for( auto& entry : diagonal_ )
128 functor( entry );
129 }
130
131 void invert ()
132 {
133 for( auto& entry : diagonal_ )
134 entry.invert();
135 }
136
137 void rightmultiply( const ThisType& other )
138 {
139 assert( other.diagonal_.size() == diagonal_.size() );
140 auto it = other.diagonal_.begin();
141 for( auto& entry : diagonal_ )
142 {
143 entry.rightmultiply( *it );
144 ++it;
145 }
146 }
147
148 void leftmultiply( const ThisType& other )
149 {
150 assert( other.diagonal_.size() == diagonal_.size() );
151 auto it = other.diagonal_.begin();
152 for( auto& entry : diagonal_ )
153 {
154 entry.leftmultiply( *it );
155 ++it;
156 }
157 }
158
160 DofBlockPtrType block( const std::size_t block )
161 {
162 assert( block < diagonal_.size() );
163 return &diagonal_[ block ];
164 }
165
167 ConstDofBlockPtrType block( const std::size_t block ) const
168 {
169 assert( block < diagonal_.size() );
170 return &diagonal_[ block ];
171 }
172
178 {
179 domainSpace().communicate( *this );
180 }
181
185 template< class Operation >
186 typename CommDataHandle< Operation > :: Type
187 dataHandle( const Operation &operation )
188 {
189 return space().createDataHandle( *this, operation );
190 }
191
192 template< class Stencil >
193 void reserve ( const Stencil &stencil, bool verbose = false )
194 {
195 // note: to use DynamicMatrix as LocalBlockType, also resize local blocks
196 diagonal_.resize( domainSpace().blockMapper().size() );
197 }
198
200 {
201 return LocalColumnObjectType( *this, domainEntity );
202 }
203
204 LocalMatrixType localMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity ) const;
205
207 {
208 return space_;
209 }
211 {
212 return space_;
213 }
214
216 const DomainSpaceType &space () const
217 {
218 return space_;
219 }
220
221 const std::string &name () const
222 {
223 return name_;
224 }
225
226 protected:
227 std::string name_;
229 std::vector< LocalBlockType > diagonal_;
232 };
233
234
235
236 // BlockDiagonalLinearOperator::LocalMatrixTraits
237 // ----------------------------------------------
238
239 template< class DiscreteFunctionSpace, class LocalBlock >
241 {
243
244 public:
246
248
251
253 };
254
255
256
257 // BlockDiagonalLinearOperator::LocalMatrix
258 // ----------------------------------------
259
260 template< class DiscreteFunctionSpace, class LocalBlock >
262 : public LocalMatrixInterface< LocalMatrixTraits >
263 {
264 typedef LocalMatrix ThisType;
266
267 public:
269
271
274
277
278 private:
279 typedef DomainBasisFunctionSetType BasisFunctionSetType;
281
282 struct SetLocalBlockFunctor
283 {
284 SetLocalBlockFunctor ( std::vector< LocalBlockType > &diagonal,
285 LocalBlockType *&localBlock )
286 : diagonal_( diagonal ),
287 localBlock_( localBlock )
288 {}
289
290 void operator() ( int localDoF, std::size_t globalDoF )
291 {
292 assert( localDoF == 0 );
293 assert( globalDoF < diagonal_.size() );
294 localBlock_ = &diagonal_[ globalDoF ];
295 }
296
297 private:
298 std::vector< LocalBlockType > &diagonal_;
299 LocalBlockType *&localBlock_;
300 };
301
302 public:
303 explicit LocalMatrix ( OperatorType &op )
304 : op_( &op ),
305 localBlock_( nullptr )
306 {}
307
308 void init ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
309 {
310 basisFunctionSet_ = domainSpace().basisFunctionSet( domainEntity );
311 SetLocalBlockFunctor f( op_->diagonal_, localBlock_ );
312 domainSpace().blockMapper().mapEach( domainEntity, f );
313 if( &domainEntity != &rangeEntity )
314 {
315 static LocalBlockType dummyBlock( 0 );
316
317 LocalBlockType *otherBlock = 0;
318 SetLocalBlockFunctor f( op_->diagonal_, otherBlock );
319 rangeSpace().blockMapper().mapEach( rangeEntity, f );
320 // check whether the blocks match, otherwise off-diagonal
321 // for off-diagonal we simply use a dummy local matrix
322 if( otherBlock != localBlock_ )
323 localBlock_ = &dummyBlock ;
324 }
325 }
326
327 void clear ()
328 {
329 localBlock() = RangeFieldType( 0 );
330 }
331 void scale ( const RangeFieldType &a )
332 {
333 localBlock() *= a;
334 }
335
336 RangeFieldType get ( int i, int j ) const
337 {
338 return localBlock()[ i ][ j ];
339 }
340 void add ( int i, int j, const RangeFieldType &value )
341 {
342 localBlock()[ i ][ j ] += value;
343 }
344 void set ( int i, int j, const RangeFieldType &value )
345 {
346 localBlock()[ i ][ j ] = value;
347 }
348
349 void clearRow ( int i )
350 {
351 localBlock()[ i ] = RangeFieldType( 0 );
352 }
353
354 void clearCol ( int j )
355 {
356 for( int i = 0; i < rows(); ++i )
357 localBlock()[ i ][ j ] = RangeFieldType( 0 );
358 }
359
360 template< class DomainLocalFunction, class RangeLocalFunction >
361 void multiplyAdd ( const DomainLocalFunction &x, RangeLocalFunction &y ) const
362 {
363 localBlock().umv( x, y );
364 }
365
366 void finalize ()
367 {}
368 void resort ()
369 {}
370
371 int rows () const
372 {
373 return localBlock().N();
374 }
375 int columns () const
376 {
377 return localBlock().M();
378 }
379
381 {
382 return op_->domainSpace();
383 }
385 { return op_->rangeSpace();
386 }
387
389 {
390 return basisFunctionSet_;
391 }
393 {
394 return basisFunctionSet_;
395 }
396
398 {
399 return domainBasisFunctionSet().entity();
400 }
402 {
403 return rangeBasisFunctionSet().entity();
404 }
405
406 private:
407 const LocalBlockType &localBlock () const
408 {
409 assert( localBlock_ );
410 return *localBlock_;
411 }
412 LocalBlockType &localBlock ()
413 {
414 assert( localBlock_ );
415 return *localBlock_;
416 }
417
418 OperatorType *op_;
419 BasisFunctionSetType basisFunctionSet_;
420 LocalBlockType *localBlock_;
421 };
422
423
424
425 // BlockDiagonalLinearOperator::LocalMatrixFactory
426 // -----------------------------------------------
427
428 template< class DiscreteFunctionSpace, class LocalBlock >
430 {
433
435 : op_( &op )
436 {}
437
439 {
440 return new ObjectType( *op_ );
441 }
442
443 private:
444 OperatorType *op_;
445 };
446
447
448
449 // Implementation of BlockDiagonalLinearOperator
450 // ---------------------------------------------
451
452 template< class DiscreteFunctionSpace, class LocalBlock >
455 ::localMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity ) const
456 {
457 return LocalMatrixType( localMatrixStack_, domainEntity, rangeEntity );
458 }
459
460 } // namespace Fem
461
462} // namespace Dune
463
464#endif // #ifndef DUNE_FEM_OPERATOR_LINEAR_BLOCKDIAGONAL_HH
Definition: bindguard.hh:11
Definition: adaptivefunction/adaptivefunction.hh:48
DofType * leakPointer()
Definition: adaptivefunction/adaptivefunction.hh:130
SizeType size() const
Return the number of blocks in the block vector.
Definition: common/discretefunction.hh:755
Interface for local matrix classes.
Definition: localmatrix.hh:32
Traits::DomainSpaceType DomainSpaceType
type of domain discrete function space
Definition: localmatrix.hh:52
Traits::RangeSpaceType RangeSpaceType
type of range discrete function space
Definition: localmatrix.hh:55
RangeSpaceType::EntityType RangeEntityType
Definition: localmatrix.hh:66
RangeSpaceType::BasisFunctionSetType RangeBasisFunctionSetType
type of base function sets within range function space
Definition: localmatrix.hh:63
DomainSpaceType::BasisFunctionSetType DomainBasisFunctionSetType
type of base function sets within domain function space
Definition: localmatrix.hh:59
Traits::RangeFieldType RangeFieldType
type of range field
Definition: localmatrix.hh:49
DomainSpaceType::EntityType DomainEntityType
Definition: localmatrix.hh:65
Definition: localmatrixwrapper.hh:48
DomainFunction DomainFunctionType
type of discrete function in the operator's domain
Definition: operator.hh:36
RangeFunction::RangeFieldType RangeFieldType
field type of the operator's range
Definition: operator.hh:43
DomainFunction::RangeFieldType DomainFieldType
field type of the operator's domain
Definition: operator.hh:41
RangeFunction RangeFunctionType
type of discrete function in the operator's range
Definition: operator.hh:38
abstract matrix operator
Definition: operator.hh:124
default implementation for a general operator stencil
Definition: stencil.hh:35
BlockDiagonalLinearOperator.
Definition: blockdiagonal.hh:30
LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType
Definition: blockdiagonal.hh:73
CommDataHandle< Operation >::Type dataHandle(const Operation &operation)
return reference to data handle object (needed to make this class work with CommunicationManager)
Definition: blockdiagonal.hh:187
ConstDofBlockPtrType block(const std::size_t block) const
return block matrix for given block number (== entity number)
Definition: blockdiagonal.hh:167
BaseType::RangeFieldType RangeFieldType
Definition: blockdiagonal.hh:39
DomainSpaceType DiscreteFunctionSpaceType
Definition: blockdiagonal.hh:55
void forEach(const Functor &functor)
Definition: blockdiagonal.hh:125
void leftmultiply(const ThisType &other)
Definition: blockdiagonal.hh:148
const RangeSpaceType & space_
Definition: blockdiagonal.hh:228
LocalColumnObjectType localColumn(const DomainEntityType &domainEntity) const
Definition: blockdiagonal.hh:199
static const int localBlockSize
Definition: blockdiagonal.hh:47
void operator()(const DomainFunctionType &u, RangeFunctionType &w) const
application operator
Definition: blockdiagonal.hh:90
ColumnObject< ThisType > LocalColumnObjectType
Definition: blockdiagonal.hh:75
void communicate()
copy matrices to ghost cells to make this class work in parallel
Definition: blockdiagonal.hh:177
BlockDiagonalLinearOperator(const std::string &name, const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace)
Definition: blockdiagonal.hh:78
const DomainSpaceType & space() const
return reference to space (needed to make this class work with CommunicationManager)
Definition: blockdiagonal.hh:216
void invert()
Definition: blockdiagonal.hh:131
RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType
Definition: blockdiagonal.hh:42
void multiply(const AdaptiveDiscreteFunction< DomainSpace > &u, AdaptiveDiscreteFunction< RangeSpace > &w) const
Definition: blockdiagonal.hh:103
void clear()
Definition: blockdiagonal.hh:118
BaseType::DomainFunctionType DomainFunctionType
Definition: blockdiagonal.hh:35
std::vector< LocalBlockType > diagonal_
Definition: blockdiagonal.hh:229
const RangeSpaceType & rangeSpace() const
Definition: blockdiagonal.hh:210
DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType
Definition: blockdiagonal.hh:41
LocalBlockType::row_type DofType
Definition: blockdiagonal.hh:54
DomainSpaceType::EntityType DomainEntityType
Definition: blockdiagonal.hh:44
const std::string & name() const
Definition: blockdiagonal.hh:221
LocalMatrixType localMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity) const
Definition: blockdiagonal.hh:455
void rightmultiply(const ThisType &other)
Definition: blockdiagonal.hh:137
BaseType::RangeFunctionType RangeFunctionType
Definition: blockdiagonal.hh:36
ObjectStack< LocalMatrixFactory > LocalMatrixStackType
Definition: blockdiagonal.hh:72
const DomainSpaceType & domainSpace() const
Definition: blockdiagonal.hh:206
BaseType::DomainFieldType DomainFieldType
Definition: blockdiagonal.hh:38
LocalMatrixFactory localMatrixFactory_
Definition: blockdiagonal.hh:230
std::string name_
Definition: blockdiagonal.hh:227
void reserve(const Stencil &stencil, bool verbose=false)
Definition: blockdiagonal.hh:193
DofBlockPtrType block(const std::size_t block)
return block matrix for given block number (== entity number)
Definition: blockdiagonal.hh:160
const LocalBlockType * ConstDofBlockPtrType
Definition: blockdiagonal.hh:53
RangeSpaceType::EntityType RangeEntityType
Definition: blockdiagonal.hh:45
LocalMatrixStackType localMatrixStack_
Definition: blockdiagonal.hh:231
LocalBlockType * DofBlockPtrType
Definition: blockdiagonal.hh:52
LocalBlock LocalBlockType
Definition: blockdiagonal.hh:49
DiscreteFunctionSpaceType::template CommDataHandle< ThisType, Operation >::Type Type
Definition: blockdiagonal.hh:62
OperatorType::RangeFieldType RangeFieldType
Definition: blockdiagonal.hh:247
OperatorType::DomainSpaceType DomainSpaceType
Definition: blockdiagonal.hh:249
RangeFieldType LittleBlockType
Definition: blockdiagonal.hh:252
OperatorType::LocalMatrix LocalMatrixType
Definition: blockdiagonal.hh:245
OperatorType::RangeSpaceType RangeSpaceType
Definition: blockdiagonal.hh:250
void finalize()
Definition: blockdiagonal.hh:366
RangeFieldType get(int i, int j) const
Definition: blockdiagonal.hh:336
int columns() const
Definition: blockdiagonal.hh:375
void add(int i, int j, const RangeFieldType &value)
Definition: blockdiagonal.hh:340
const DomainBasisFunctionSetType & domainBasisFunctionSet() const
Definition: blockdiagonal.hh:388
BlockDiagonalLinearOperator< DiscreteFunctionSpace, LocalBlock > OperatorType
Definition: blockdiagonal.hh:268
void scale(const RangeFieldType &a)
Definition: blockdiagonal.hh:331
BaseType::DomainEntityType DomainEntityType
Definition: blockdiagonal.hh:275
void resort()
Definition: blockdiagonal.hh:368
void init(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity)
Definition: blockdiagonal.hh:308
const DomainSpaceType & domainSpace() const
Definition: blockdiagonal.hh:380
void set(int i, int j, const RangeFieldType &value)
Definition: blockdiagonal.hh:344
void clearRow(int i)
Definition: blockdiagonal.hh:349
BaseType::RangeFieldType RangeFieldType
Definition: blockdiagonal.hh:270
void clearCol(int j)
Definition: blockdiagonal.hh:354
BaseType::DomainBasisFunctionSetType DomainBasisFunctionSetType
Definition: blockdiagonal.hh:272
const RangeBasisFunctionSetType & rangeBasisFunctionSet() const
Definition: blockdiagonal.hh:392
void multiplyAdd(const DomainLocalFunction &x, RangeLocalFunction &y) const
Definition: blockdiagonal.hh:361
void clear()
Definition: blockdiagonal.hh:327
int rows() const
Definition: blockdiagonal.hh:371
BaseType::RangeBasisFunctionSetType RangeBasisFunctionSetType
Definition: blockdiagonal.hh:273
const RangeEntityType & rangeEntity() const
Definition: blockdiagonal.hh:401
BaseType::RangeEntityType RangeEntityType
Definition: blockdiagonal.hh:276
const DomainEntityType & domainEntity() const
Definition: blockdiagonal.hh:397
LocalMatrix(OperatorType &op)
Definition: blockdiagonal.hh:303
const RangeSpaceType & rangeSpace() const
Definition: blockdiagonal.hh:384
LocalMatrix ObjectType
Definition: blockdiagonal.hh:432
ObjectType * newObject() const
Definition: blockdiagonal.hh:438
BlockDiagonalLinearOperator< DiscreteFunctionSpace, LocalBlock > OperatorType
Definition: blockdiagonal.hh:431
LocalMatrixFactory(OperatorType &op)
Definition: blockdiagonal.hh:434
Definition: columnobject.hh:12
discrete function space