dune-fem  2.6-git
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 
18 namespace 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 
67  class LocalMatrixTraits;
68  class LocalMatrix;
69  struct LocalMatrixFactory;
70 
71  public:
74 
76 
77 
78  BlockDiagonalLinearOperator ( const std::string &name,
80  const RangeSpaceType &rangeSpace )
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 
177  void communicate ()
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 
199  LocalColumnObjectType localColumn ( const DomainEntityType &domainEntity ) const
200  {
201  return LocalColumnObjectType( *this, domainEntity );
202  }
203 
204  LocalMatrixType localMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity ) const;
205 
206  const DomainSpaceType &domainSpace () const
207  {
208  return space_;
209  }
210  const RangeSpaceType &rangeSpace () const
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 
273  typedef typename BaseType::RangeBasisFunctionSetType RangeBasisFunctionSetType;
274 
275  typedef typename BaseType::DomainEntityType DomainEntityType;
276  typedef typename BaseType::RangeEntityType RangeEntityType;
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 
380  const DomainSpaceType &domainSpace () const
381  {
382  return op_->domainSpace();
383  }
384  const RangeSpaceType &rangeSpace () const
385  { return op_->rangeSpace();
386  }
387 
389  {
390  return basisFunctionSet_;
391  }
393  {
394  return basisFunctionSet_;
395  }
396 
398  {
399  return domainBasisFunctionSet().entity();
400  }
401  const RangeEntityType &rangeEntity () const
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: function/adaptivefunction/adaptivefunction.hh:45
DofType * leakPointer()
Definition: function/adaptivefunction/adaptivefunction.hh:118
SizeType size() const
Return the number of blocks in the block vector.
Definition: common/discretefunction.hh:704
Definition: bartonnackmaninterface.hh:17
Interface for local matrix classes.
Definition: localmatrix.hh:32
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
RangeSpaceType ::BasisFunctionSetType RangeBasisFunctionSetType
type of base function sets within range function space
Definition: localmatrix.hh:63
Traits ::RangeSpaceType RangeSpaceType
type of range discrete function space
Definition: localmatrix.hh:55
Traits ::DomainSpaceType DomainSpaceType
type of domain discrete function space
Definition: localmatrix.hh:52
RangeSpaceType::EntityType RangeEntityType
Definition: localmatrix.hh:66
DomainSpaceType::EntityType DomainEntityType
Definition: localmatrix.hh:65
Definition: localmatrixwrapper.hh:48
AdaptiveDiscreteFunction< DiscreteFunctionSpace > DomainFunctionType
type of discrete function in the operator's domain
Definition: operator.hh:28
RangeFunction::RangeFieldType RangeFieldType
field type of the operator's range
Definition: operator.hh:35
DomainFunction::RangeFieldType DomainFieldType
field type of the operator's domain
Definition: operator.hh:33
AdaptiveDiscreteFunction< DiscreteFunctionSpace > RangeFunctionType
type of discrete function in the operator's range
Definition: operator.hh:30
abstract matrix operator
Definition: operator.hh:110
default implementation for a general operator stencil
Definition: stencil.hh:32
BlockDiagonalLinearOperator.
Definition: blockdiagonal.hh:30
LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType
Definition: blockdiagonal.hh:73
const DomainSpaceType & space() const
return reference to space (needed to make this class work with CommunicationManager)
Definition: blockdiagonal.hh:216
const RangeSpaceType & rangeSpace() const
Definition: blockdiagonal.hh:210
ConstDofBlockPtrType block(const std::size_t block) const
return block matrix for given block number (== entity number)
Definition: blockdiagonal.hh:167
const std::string & name() const
Definition: blockdiagonal.hh:221
BaseType::RangeFieldType RangeFieldType
Definition: blockdiagonal.hh:39
DomainSpaceType DiscreteFunctionSpaceType
Definition: blockdiagonal.hh:55
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
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
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
const DomainSpaceType & domainSpace() const
Definition: blockdiagonal.hh:206
BlockDiagonalLinearOperator(const std::string &name, const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace)
Definition: blockdiagonal.hh:78
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
DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType
Definition: blockdiagonal.hh:41
LocalBlockType::row_type DofType
Definition: blockdiagonal.hh:54
DomainSpaceType::EntityType DomainEntityType
Definition: blockdiagonal.hh:44
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:69
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
const DomainSpaceType & domainSpace() const
Definition: blockdiagonal.hh:380
RangeFieldType get(int i, int j) const
Definition: blockdiagonal.hh:336
const DomainBasisFunctionSetType & domainBasisFunctionSet() const
Definition: blockdiagonal.hh:388
int columns() const
Definition: blockdiagonal.hh:375
void add(int i, int j, const RangeFieldType &value)
Definition: blockdiagonal.hh:340
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 DomainEntityType & domainEntity() const
Definition: blockdiagonal.hh:397
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
void multiplyAdd(const DomainLocalFunction &x, RangeLocalFunction &y) const
Definition: blockdiagonal.hh:361
const RangeEntityType & rangeEntity() const
Definition: blockdiagonal.hh:401
void clear()
Definition: blockdiagonal.hh:327
int rows() const
Definition: blockdiagonal.hh:371
BaseType::RangeBasisFunctionSetType RangeBasisFunctionSetType
Definition: blockdiagonal.hh:273
const RangeSpaceType & rangeSpace() const
Definition: blockdiagonal.hh:384
BaseType::RangeEntityType RangeEntityType
Definition: blockdiagonal.hh:276
const RangeBasisFunctionSetType & rangeBasisFunctionSet() const
Definition: blockdiagonal.hh:392
LocalMatrix(OperatorType &op)
Definition: blockdiagonal.hh:303
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