dune-fem  2.8-git
istlmatrix.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_ISTLMATRIXWRAPPER_HH
2 #define DUNE_FEM_ISTLMATRIXWRAPPER_HH
3 
4 #if HAVE_DUNE_ISTL
5 
6 //- system includes
7 #include <vector>
8 #include <iostream>
9 #include <fstream>
10 #include <algorithm>
11 #include <set>
12 #include <map>
13 #include <string>
14 
15 //- Dune common includes
16 #include <dune/common/exceptions.hh>
17 #include <dune/common/fmatrix.hh>
18 
19 //- Dune istl includes
20 #include <dune/istl/bvector.hh>
21 #include <dune/istl/bcrsmatrix.hh>
22 #include <dune/istl/preconditioners.hh>
23 
24 //- Dune fem includes
32 #include <dune/fem/io/parameter.hh>
35 
39 
40 namespace Dune
41 {
42 
43  namespace Fem
44  {
45 
46 
48  template< class MatrixObject >
49  class ISTLLocalMatrix;
50 
51  template <class DomainSpaceImp, class RangeSpaceImp,
52  class DomainBlock = Dune::FieldVector< typename DomainSpaceImp :: RangeFieldType, DomainSpaceImp:: localBlockSize >,
53  class RangeBlock = Dune::FieldVector< typename RangeSpaceImp :: RangeFieldType, RangeSpaceImp :: localBlockSize > >
54  class ISTLMatrixObject;
55 
57  // --ISTLMatrixHandle
59  template <class LittleBlockType, class RowDiscreteFunctionImp, class ColDiscreteFunctionImp = RowDiscreteFunctionImp>
60  class ImprovedBCRSMatrix : public BCRSMatrix<LittleBlockType>
61  {
62  friend struct MatrixDimension<ImprovedBCRSMatrix>;
63  public:
64  typedef RowDiscreteFunctionImp RowDiscreteFunctionType;
65  typedef ColDiscreteFunctionImp ColDiscreteFunctionType;
66 
67  typedef BCRSMatrix<LittleBlockType> BaseType;
69  typedef BaseType MatrixBaseType;
70  typedef typename BaseType :: RowIterator RowIteratorType ;
71  typedef typename BaseType :: ColIterator ColIteratorType ;
72 
73  typedef ImprovedBCRSMatrix< LittleBlockType, RowDiscreteFunctionImp, ColDiscreteFunctionImp > ThisType;
74 
75  typedef typename BaseType :: size_type size_type;
76 
77  //===== type definitions and constants
78 
80  typedef typename BaseType::field_type field_type;
81 
83  typedef typename BaseType::block_type block_type;
84 
86  typedef typename BaseType:: allocator_type allocator_type;
87 
89  typedef typename BaseType :: row_type row_type;
90 
92  typedef typename BaseType :: ColIterator ColIterator;
93 
95  typedef typename BaseType :: ConstColIterator ConstColIterator;
96 
98  typedef typename BaseType :: RowIterator RowIterator;
99 
101  typedef typename BaseType :: ConstRowIterator ConstRowIterator;
102 
104  typedef typename ColDiscreteFunctionType :: DiscreteFunctionSpaceType RangeSpaceType;
105 
107  typedef typename RowDiscreteFunctionType :: DofStorageType RowBlockVectorType;
108 
110  typedef typename ColDiscreteFunctionType :: DofStorageType ColBlockVectorType;
111 
113  typedef typename RangeSpaceType :: GridType :: Traits :: CollectiveCommunication CollectiveCommunictionType ;
114 
115  typedef typename BaseType :: BuildMode BuildMode ;
116 
117  public:
119  ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nnz, double overflowFraction) :
120  BaseType (rows, cols, nnz, overflowFraction, BaseType::implicit)
121  {}
122 
124  ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nz = 0 ) :
125  BaseType (rows, cols, BaseType :: row_wise)
126  {}
127 
129  ImprovedBCRSMatrix( ) :
130  BaseType ()
131  {}
132 
134  ImprovedBCRSMatrix(const ImprovedBCRSMatrix& org) :
135  BaseType(org)
136  {}
137 
138  template < class SparsityPattern >
139  void createEntries(const SparsityPattern& sparsityPattern)
140  {
141  // not insert map of indices into matrix
142  auto endcreate = this->createend();
143  const auto endsp = sparsityPattern.end();
144  for(auto create = this->createbegin(); create != endcreate; ++create)
145  {
146  const auto row = sparsityPattern.find( create.index() );
147  if ( row == endsp )
148  continue;
149  const auto& localIndices = ( *row ).second;
150  const auto end = localIndices.end();
151  // insert all indices for this row
152  for (auto it = localIndices.begin(); it != end; ++it)
153  create.insert( *it );
154  }
155  }
156 
158  void clear()
159  {
160  for (auto& row : *this)
161  for (auto& entry : row)
162  entry = 0;
163  }
164 
166  void unitRow( const size_t row )
167  {
168  block_type idBlock( 0 );
169  for (int i = 0; i < idBlock.rows; ++i)
170  idBlock[i][i] = 1.0;
171 
172  auto& matRow = (*this)[ row ];
173  auto colIt = matRow.begin();
174  const auto& colEndIt = matRow.end();
175  for (; colIt != colEndIt; ++colIt)
176  {
177  if( colIt.index() == row )
178  *colIt = idBlock;
179  else
180  *colIt = 0.0;
181  }
182  }
183 
185  template <class HangingNodesType>
186  void setup(ThisType& oldMatrix, const HangingNodesType& hangingNodes)
187  {
188  // necessary because element traversal not necessarily is in ascending order
189  typedef std::set< std::pair<int, block_type> > LocalEntryType;
190  typedef std::map< int , LocalEntryType > EntriesType;
191  EntriesType entries;
192 
193  // map of indices
194  std::map< int , std::set<int> > indices;
195  // not insert map of indices into matrix
196  auto rowend = oldMatrix.end();
197  for(auto it = oldMatrix.begin(); it != rowend; ++it)
198  {
199  const auto row = it.index();
200  auto& localIndices = indices[ row ];
201 
202  if( hangingNodes.isHangingNode( row ) )
203  {
204  // insert columns into other columns
205  const auto& cols = hangingNodes.associatedDofs( row );
206  const auto colSize = cols.size();
207  for(auto i=0; i<colSize; ++i)
208  {
209  assert( ! hangingNodes.isHangingNode( cols[i].first ) );
210 
211  // get local indices of col
212  auto& localColIndices = indices[ cols[i].first ];
213  auto& localEntry = entries[ cols[i].first ];
214 
215  // copy from old matrix
216  auto endj = (*it).end();
217  for (auto j= (*it).begin(); j!=endj; ++j)
218  {
219  localColIndices.insert( j.index () );
220  localEntry.insert( std::make_pair( j.index(), (cols[i].second * (*j)) ));
221  }
222  }
223 
224  // insert diagonal and hanging columns
225  localIndices.insert( row );
226  for(auto i=0; i<colSize; ++i)
227  localIndices.insert( cols[i].first );
228  }
229  else
230  {
231  // copy from old matrix
232  auto endj = (*it).end();
233  for (auto j= (*it).begin(); j!=endj; ++j)
234  localIndices.insert( j.index () );
235  }
236  }
237 
238  // create matrix from entry map
239  createEntries( indices );
240 
241  // not insert map of indices into matrix
242  auto rowit = oldMatrix.begin();
243 
244  auto endcreate = this->end();
245  for(auto create = this->begin(); create != endcreate; ++create, ++rowit )
246  {
247  assert( rowit != oldMatrix.end() );
248 
249  const auto row = create.index();
250  if( hangingNodes.isHangingNode( row ) )
251  {
252  const auto& cols = hangingNodes.associatedDofs( row );
253 
254  std::map< const int , block_type > colMap;
255  // only working for block size 1 ath the moment
256  assert( block_type :: rows == 1 );
257  // insert columns into map
258  const auto colSize = cols.size();
259  for( auto i=0; i<colSize; ++i)
260  colMap[ cols[i].first ] = -cols[i].second;
261  // insert diagonal into map
262  colMap[ row ] = 1;
263 
264  auto endj = (*create).end();
265  for (auto j= (*create).begin(); j!=endj; ++j)
266  {
267  assert( colMap.find( j.index() ) != colMap.end() );
268  (*j) = colMap[ j.index() ];
269  }
270  }
271  // if entries are equal, just copy
272  else if ( entries.find( row ) == entries.end() )
273  {
274  auto colit = (*rowit).begin();
275  auto endj = (*create).end();
276  for (auto j= (*create).begin(); j!=endj; ++j, ++colit )
277  {
278  assert( colit != (*rowit).end() );
279  (*j) = (*colit);
280  }
281  }
282  else
283  {
284  std::map< int , block_type > oldCols;
285 
286  {
287  auto colend = (*rowit).end();
288  for(auto colit = (*rowit).begin(); colit != colend; ++colit)
289  oldCols[ colit.index() ] = 0;
290  }
291 
292  auto entry = entries.find( row );
293  assert( entry != entries.end ());
294 
295  {
296  auto endcol = (*entry).second.end();
297  for( auto co = (*entry).second.begin(); co != endcol; ++co)
298  oldCols[ (*co).first ] = 0;
299  }
300 
301  {
302  auto colend = (*rowit).end();
303  for(auto colit = (*rowit).begin(); colit != colend; ++colit)
304  oldCols[ colit.index() ] += (*colit);
305  }
306 
307  {
308  auto endcol = (*entry).second.end();
309  for( auto co = (*entry).second.begin(); co != endcol; ++co)
310  oldCols[ (*co).first ] += (*co).second;
311  }
312 
313  auto endj = (*create).end();
314  for (auto j= (*create).begin(); j!=endj; ++j )
315  {
316  auto colEntry = oldCols.find( j.index() );
317  if( colEntry != oldCols.end() )
318  (*j) = (*colEntry).second;
319  else
320  abort();
321  }
322  }
323  }
324  }
325 
327  void extractDiagonal( ColBlockVectorType& diag ) const
328  {
329  const auto endi = this->end();
330  for (auto i = this->begin(); i!=endi; ++i)
331  {
332  // get diagonal entry of matrix
333  const auto row = i.index();
334  auto entry = (*i).find( row );
335  const LittleBlockType& block = (*entry);
336  enum { blockSize = LittleBlockType :: rows };
337  for( auto l=0; l<blockSize; ++l )
338  diag[ row ][ l ] = block[ l ][ l ];
339  }
340  }
341 
344  field_type operator()(const std::size_t row, const std::size_t col) const
345  {
346  const std::size_t blockRow(row/(LittleBlockType :: rows));
347  const std::size_t localRowIdx(row%(LittleBlockType :: rows));
348  const std::size_t blockCol(col/(LittleBlockType :: cols));
349  const std::size_t localColIdx(col%(LittleBlockType :: cols));
350 
351  const auto& matrixRow(this->operator[](blockRow));
352  auto entry = matrixRow.find( blockCol );
353  const LittleBlockType& block = (*entry);
354  return block[localRowIdx][localColIdx];
355  }
356 
359  void set(const std::size_t row, const std::size_t col, field_type value)
360  {
361  const std::size_t blockRow(row/(LittleBlockType :: rows));
362  const std::size_t localRowIdx(row%(LittleBlockType :: rows));
363  const std::size_t blockCol(col/(LittleBlockType :: cols));
364  const std::size_t localColIdx(col%(LittleBlockType :: cols));
365 
366  auto& matrixRow(this->operator[](blockRow));
367  auto entry = matrixRow.find( blockCol );
368  LittleBlockType& block = (*entry);
369  block[localRowIdx][localColIdx] = value;
370  }
371 
373  void print(std::ostream& s=std::cout, unsigned int offset=0) const
374  {
375  s.precision( 6 );
376  const auto endi=this->end();
377  for (auto i=this->begin(); i!=endi; ++i)
378  {
379  const auto endj = (*i).end();
380  for (auto j=(*i).begin(); j!=endj; ++j)
381  if( (*j).infinity_norm() > 1.e-15)
382  s << i.index()+offset << " " << j.index()+offset << " " << *j << std::endl;
383  }
384  }
385  };
386 
387 
388 
389  // ISTLLocalMatrixTraits
390  // ---------------------
391 
392  template< class MatrixObject >
393  struct ISTLLocalMatrixTraits
394  {
395  typedef typename MatrixObject::DomainSpaceType DomainSpaceType;
396  typedef typename MatrixObject::RangeSpaceType RangeSpaceType;
397  typedef typename DomainSpaceType::RangeFieldType RangeFieldType;
398 
399  typedef ISTLLocalMatrix< MatrixObject > LocalMatrixType;
400  typedef typename MatrixObject::MatrixType::block_type LittleBlockType;
401  };
402 
403 
404  // ISTLLocalMatrix
405  // ---------------
406 
407  template< class MatrixObject >
408  class ISTLLocalMatrix
409  : public LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > >
410  {
411  public:
413  typedef LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > > BaseType;
414 
416  typedef MatrixObject MatrixObjectType;
418  typedef typename MatrixObjectType::MatrixType MatrixType;
420  typedef typename MatrixType::block_type LittleBlockType;
421  // type of row and column indices
422  typedef typename MatrixType :: size_type size_type;
423 
424  typedef typename MatrixObjectType::DomainSpaceType DomainSpaceType;
425  typedef typename MatrixObjectType::RangeSpaceType RangeSpaceType;
426 
427  typedef typename MatrixObjectType::DomainEntityType DomainEntityType;
428  typedef typename MatrixObjectType::RangeEntityType RangeEntityType;
429 
431  typedef typename DomainSpaceType::RangeFieldType DofType;
432  typedef typename MatrixType::row_type RowType;
433 
435  typedef typename DomainSpaceType::BlockMapperType ColMapperType;
437  typedef typename RangeSpaceType::BlockMapperType RowMapperType;
438 
439  static const int littleCols = MatrixObjectType::littleCols;
440  static const int littleRows = MatrixObjectType::littleRows;
441 
442  typedef typename MatrixType::size_type Index;
443 
447  [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
448  ISTLLocalMatrix ( const MatrixObjectType& mObj, const DomainSpaceType& domainSpace, const RangeSpaceType& rangeSpace )
449  : BaseType( domainSpace, rangeSpace ),
450  rowMapper_( rangeSpace.blockMapper() ),
451  colMapper_( domainSpace.blockMapper() ),
452  numRows_( rowMapper_.maxNumDofs() ),
453  numCols_( colMapper_.maxNumDofs() ),
454  matrixObj_( mObj )
455  {}
456 
460  [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
461  ISTLLocalMatrix ( const ISTLLocalMatrix& org )
462  : BaseType( org ),
463  rowMapper_(org.rowMapper_),
464  colMapper_(org.colMapper_),
465  numRows_( org.numRows_ ),
466  numCols_( org.numCols_ ),
467  matrices_(org.matrices_),
468  matrixObj_(org.matrixObj_)
469  {}
470 
472  void init ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
473  {
474  // initialize base functions sets
475  BaseType :: init ( domainEntity, rangeEntity );
476 
477  numRows_ = rowMapper_.numDofs( rangeEntity );
478  numCols_ = colMapper_.numDofs( domainEntity );
479 
480  // resize matrix pointer storage for new row/col numbers
481  matrices_.resize( numRows_ );
482  for( auto& row : matrices_ )
483  {
484  row.resize( numCols_, nullptr );
485  }
486 
487  if( matrixObj_.implicitModeActive() )
488  {
489  auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
490  {
491  return matrixObj_.matrix().entry( index.first, index.second );
492  };
493  initBlocks( blockAccess, domainEntity, rangeEntity );
494  }
495  else
496  {
497  auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
498  {
499  return matrixObj_.matrix()[ index.first ][ index.second ];
500  };
501  initBlocks( blockAccess, domainEntity, rangeEntity );
502  }
503  }
504 
505  template <class BlockAccess>
506  void initBlocks( BlockAccess& blockAccess, const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
507  {
508  auto functor = [ this, &blockAccess ] ( std::pair< int, int > local, const std::pair< Index, Index > &index )
509  {
510  matrices_[ local.first ][ local.second ] = &blockAccess( index );
511  };
512 
513  rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
514  }
515 
516  private:
517  // check whether given (row,col) pair is valid
518  void check(int localRow, int localCol) const
519  {
520 #ifndef NDEBUG
521  const std::size_t row = (int) localRow / littleRows;
522  const std::size_t col = (int) localCol / littleCols;
523  const int lRow = localRow%littleRows;
524  const int lCol = localCol%littleCols;
525  assert( row < matrices_.size() ) ;
526  assert( col < matrices_[row].size() );
527  assert( lRow < littleRows );
528  assert( lCol < littleCols );
529 #endif
530  }
531 
532  DofType& getValue(const int localRow, const int localCol)
533  {
534  const int row = (int) localRow / littleRows;
535  const int col = (int) localCol / littleCols;
536  const int lRow = localRow%littleRows;
537  const int lCol = localCol%littleCols;
538  return (*matrices_[row][col])[lRow][lCol];
539  }
540 
541  public:
542  const DofType get(const int localRow, const int localCol) const
543  {
544  const int row = (int) localRow / littleRows;
545  const int col = (int) localCol / littleCols;
546  const int lRow = localRow%littleRows;
547  const int lCol = localCol%littleCols;
548  return (*matrices_[row][col])[lRow][lCol];
549  }
550 
551  void scale (const DofType& scalar)
552  {
553  const size_type nrows = matrices_.size();
554  for(size_type i=0; i<nrows; ++i)
555  {
556  const size_type ncols = matrices_[i].size();
557  for(size_type j=0; j<ncols; ++j)
558  {
559  (*matrices_[i][j]) *= scalar;
560  }
561  }
562  }
563 
564  void add(const int localRow, const int localCol , const DofType value)
565  {
566 #ifndef NDEBUG
567  check(localRow,localCol);
568 #endif
569  getValue(localRow,localCol) += value;
570  }
571 
572  void set(const int localRow, const int localCol , const DofType value)
573  {
574 #ifndef NDEBUG
575  check(localRow,localCol);
576 #endif
577  getValue(localRow,localCol) = value;
578  }
579 
581  void unitRow(const int localRow)
582  {
583  const int row = (int) localRow / littleRows;
584  const int lRow = localRow%littleRows;
585 
586  // clear row
587  doClearRow( row, lRow );
588 
589  // set diagonal entry to 1
590  (*matrices_[row][row])[lRow][lRow] = 1;
591  }
592 
594  void clear ()
595  {
596  const size_type nrows = matrices_.size();
597  for(size_type i=0; i<nrows; ++i)
598  {
599  const size_type ncols = matrices_[i].size();
600  for(size_type j=0; j<ncols; ++j)
601  {
602  (*matrices_[i][j]) = (DofType) 0;
603  }
604  }
605  }
606 
608  void clearRow ( const int localRow )
609  {
610  const int row = (int) localRow / littleRows;
611  const int lRow = localRow%littleRows;
612 
613  // clear the row
614  doClearRow( row, lRow );
615  }
616 
618  void resort ()
619  {}
620 
621  protected:
623  void doClearRow ( const int row, const int lRow )
624  {
625  // get number of columns
626  const auto col = this->columns();
627  for(auto localCol=0; localCol<col; ++localCol)
628  {
629  const int col = (int) localCol / littleCols;
630  const int lCol = localCol%littleCols;
631  (*matrices_[row][col])[lRow][lCol] = 0;
632  }
633  }
634 
635  private:
636  // special mapper omitting block size
637  const RowMapperType& rowMapper_;
638  const ColMapperType& colMapper_;
639 
640  // number of local matrices
641  int numRows_;
642  int numCols_;
643 
644  // dynamic matrix with pointers to block matrices
645  std::vector< std::vector< LittleBlockType* > > matrices_;
646 
647  // matrix to build
648  const MatrixObjectType& matrixObj_;
649  };
650 
651 
653  template <class DomainSpaceImp, class RangeSpaceImp, class DomainBlock, class RangeBlock >
654  class ISTLMatrixObject
655  {
656  public:
658  typedef DomainSpaceImp DomainSpaceType;
660  typedef RangeSpaceImp RangeSpaceType;
661 
663  typedef ISTLMatrixObject< DomainSpaceImp, RangeSpaceImp, DomainBlock, RangeBlock > ThisType;
664  typedef ThisType PreconditionMatrixType;
665 
666  typedef typename DomainSpaceType::GridType GridType;
667 
668  typedef typename RangeSpaceType :: EntityType RangeEntityType ;
669  typedef typename DomainSpaceType :: EntityType DomainEntityType ;
670 
671  enum { littleCols = DomainSpaceType :: localBlockSize };
672  enum { littleRows = RangeSpaceType :: localBlockSize };
673 
674  typedef FieldMatrix<typename DomainSpaceType :: RangeFieldType, littleRows, littleCols> LittleBlockType;
675  typedef LittleBlockType block_type;
676  typedef LittleBlockType MatrixBlockType;
677 
678  typedef ISTLBlockVectorDiscreteFunction< RangeSpaceType, RangeBlock > RowDiscreteFunctionType;
679  typedef ISTLBlockVectorDiscreteFunction< DomainSpaceType, DomainBlock > ColumnDiscreteFunctionType;
680 
681  typedef typename RowDiscreteFunctionType :: DofStorageType RowBlockVectorType;
682  typedef typename ColumnDiscreteFunctionType :: DofStorageType ColumnBlockVectorType;
683 
684  typedef typename RangeSpaceType :: BlockMapperType RowMapperType;
685  typedef typename DomainSpaceType :: BlockMapperType ColMapperType;
686 
688  typedef ImprovedBCRSMatrix< LittleBlockType , ColumnDiscreteFunctionType , RowDiscreteFunctionType > MatrixType;
689  typedef ISTLParallelMatrixAdapterInterface< MatrixType > MatrixAdapterType;
690 
692  typedef ISTLLocalMatrix<ThisType> ObjectType;
693  friend class ISTLLocalMatrix<ThisType>;
694 
695  typedef ThisType LocalMatrixFactoryType;
696  typedef ObjectStack< LocalMatrixFactoryType > LocalMatrixStackType;
698  typedef LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType;
699  typedef ColumnObject< ThisType > LocalColumnObjectType;
700 
701  protected:
702  const DomainSpaceType& domainSpace_;
703  const RangeSpaceType& rangeSpace_;
704 
705  // special row mapper
706  RowMapperType& rowMapper_;
707  // special col mapper
708  ColMapperType& colMapper_;
709 
710  int size_;
711  int sequence_;
712 
713  mutable std::unique_ptr< MatrixType > matrix_;
714 
715  mutable LocalMatrixStackType localMatrixStack_;
716 
717  mutable std::unique_ptr< MatrixAdapterType > matrixAdap_;
718  mutable std::unique_ptr< ColumnBlockVectorType > Arg_;
719  mutable std::unique_ptr< RowBlockVectorType > Dest_;
720  // overflow fraction for implicit build mode
721  const double overflowFraction_;
722  ISTLSolverParameter param_;
723  public:
724  ISTLMatrixObject(const ISTLMatrixObject&) = delete;
725 
729  ISTLMatrixObject ( const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const ISTLSolverParameter& param = ISTLSolverParameter() ) :
730  domainSpace_(domainSpace)
731  , rangeSpace_(rangeSpace)
732  , rowMapper_( rangeSpace.blockMapper() )
733  , colMapper_( domainSpace.blockMapper() )
734  , size_(-1)
735  , sequence_(-1)
736  , localMatrixStack_( *this )
737  , overflowFraction_( param.overflowFraction() )
738  , param_( param )
739  {}
740 
741  protected:
744  MatrixType& matrix() const
745  {
746  assert( matrix_ );
747  return *matrix_;
748  }
749 
750  public:
751  void printTexInfo(std::ostream& out) const
752  {
753  out << "ISTL MatrixObj: ";
754  out << "\\\\ \n";
755  }
756 
758  std::string preconditionName() const
759  {
760  return "";
761  }
762 
763  void createMatrixAdapter ( const ISTLSolverParameter& param ) const
764  {
765  if( !matrixAdap_ )
766  {
767  matrixAdap_ = ISTLMatrixAdapterFactory< ThisType > :: matrixAdapter( *this, param );
768  }
769  assert( matrixAdap_ );
770  }
771 
773  MatrixAdapterType& matrixAdapter( const ISTLSolverParameter& parameter ) const
774  {
775  const ISTLSolverParameter* param = dynamic_cast< const ISTLSolverParameter* > (&parameter);
776  if( param )
777  createMatrixAdapter( *param );
778  else
779  {
780  ISTLSolverParameter newparam( parameter );
781  createMatrixAdapter( newparam );
782  }
783 
784  finalizeAssembly();
785  return *matrixAdap_;
786  }
787 
789  MatrixAdapterType& matrixAdapter() const
790  {
791  return matrixAdapter( param_ );
792  }
793 
794  public:
795  MatrixType& exportMatrix() const
796  {
797  finalizeAssembly();
798  return matrix();
799  }
800 
801  bool implicitModeActive() const
802  {
803  // implicit build mode is only active when the
804  // build mode of the matrix is implicit and the
805  // matrix is currently being build
806 
807  if( matrix().buildMode() == MatrixType::implicit && matrix().buildStage() == MatrixType::building )
808  return true;
809  else
810  return false;
811  }
812 
813  void finalizeAssembly() const
814  {
815  const_cast< ThisType& > (*this).compress();
816  }
817 
818  // compress matrix if not already done before and only in implicit build mode
819  void compress( )
820  {
821  if( implicitModeActive() )
822  {
823  matrix().compress();
824  }
825  }
826 
828  void clear()
829  {
830  matrix().clear();
831  // clean matrix adapter and other helper classes
832  removeObj();
833  }
834 
835  void unitRow( const size_t row )
836  {
837  matrix().unitRow( row );
838  }
839 
840  template <class Vector>
841  void setUnitRows( const Vector &rows )
842  {
843  const auto &auxiliaryDofs = domainSpace().auxiliaryDofs();
844 
845  for (auto r : rows)
846  {
847  const std::size_t blockRow( r/(LittleBlockType :: rows) );
848  const std::size_t localRowIdx( r%(LittleBlockType :: rows) );
849  auto& row = matrix()[blockRow];
850  const auto endcol = row.end();
851 #ifndef NDEBUG
852  bool set = false;
853 #endif
854  for (auto col=row.begin(); col!=endcol; ++col)
855  {
856  for (auto& entry : (*col)[localRowIdx])
857  entry = 0;
858  if (col.index() == blockRow)
859  {
860  (*col)[localRowIdx][localRowIdx] = auxiliaryDofs.contains( r )? 0.0 : 1.0;
861 #ifndef NDEBUG
862  set = true;
863 #endif
864  }
865  }
866  assert(set);
867  }
868  }
869 
871  void reserve()
872  {
873  reserve( Stencil<DomainSpaceType,RangeSpaceType>(domainSpace_, rangeSpace_), true );
874  }
875 
877  template <class Set>
878  void reserve (const std::vector< Set >& sparsityPattern )
879  {
880  reserve( StencilWrapper< DomainSpaceType,RangeSpaceType, Set >( sparsityPattern ), false );
881  }
882 
884  template <class Stencil>
885  void reserve(const Stencil &stencil, const bool implicit = true )
886  {
887  // if grid sequence number changed, rebuild matrix
888  if(sequence_ != domainSpace().sequence())
889  {
890  removeObj();
891 
892  if( implicit )
893  {
894  auto nnz = stencil.maxNonZerosEstimate();
895  if( nnz == 0 )
896  {
897  Stencil tmpStencil( stencil );
898  tmpStencil.fill( *(domainSpace_.begin()), *(rangeSpace_.begin()) );
899  nnz = tmpStencil.maxNonZerosEstimate();
900  }
901  matrix_.reset( new MatrixType( rowMapper_.size(), colMapper_.size(), nnz, overflowFraction_ ) );
902  }
903  else
904  {
905  matrix_.reset( new MatrixType( rowMapper_.size(), colMapper_.size() ) );
906  matrix().createEntries( stencil.globalStencil() );
907  }
908 
909  sequence_ = domainSpace().sequence();
910  }
911  }
912 
914  template <class HangingNodesType>
915  void changeHangingNodes(const HangingNodesType& hangingNodes)
916  {
917  // create new matrix
918  MatrixType* newMatrix = new MatrixType(rowMapper_.size(), colMapper_.size());
919 
920  // setup with hanging rows
921  newMatrix->setup( *matrix_ , hangingNodes );
922 
923  // remove old matrix
924  removeObj();
925 
926  // store new matrix
927  matrix_.reset( newMatrix );
928  }
929 
932  void extractDiagonal( ColumnDiscreteFunctionType& diag ) const
933  {
934  // extract diagonal entries
935  matrix().extractDiagonal( diag.blockVector() );
936  }
937 
939  bool rightPrecondition() const
940  {
941  return false;
942  }
943 
945  void apply(const ColumnDiscreteFunctionType& arg, RowDiscreteFunctionType& dest) const
946  {
947  matrixAdapter().apply( arg.blockVector(), dest.blockVector() );
948  }
949 
951  void resort()
952  {}
953 
956  void createPreconditionMatrix()
957  {}
958 
960  void print(std::ostream & s) const
961  {
962  matrix().print(s);
963  }
964 
965  const DomainSpaceType& domainSpace() const
966  {
967  return domainSpace_;
968  }
969  const RangeSpaceType& rangeSpace() const
970  {
971  return rangeSpace_;
972  }
973 
974  const RowMapperType& rowMapper() const
975  {
976  return rowMapper_;
977  }
978  const ColMapperType& colMapper() const
979  {
980  return colMapper_;
981  }
982 
984  ObjectType* newObject() const
985  {
986  return new ObjectType(*this, domainSpace(), rangeSpace());
987  }
988 
993  [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
994  LocalMatrixType localMatrix( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity ) const
995  {
996  return LocalMatrixType( localMatrixStack_, domainEntity, rangeEntity );
997  }
998 
1003  [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
1004  LocalMatrixType localMatrix() const
1005  {
1006  return LocalMatrixType( localMatrixStack_ );
1007  }
1008 
1009  LocalColumnObjectType localColumn( const DomainEntityType &domainEntity ) const
1010  {
1011  return LocalColumnObjectType ( *this, domainEntity );
1012  }
1013 
1014  protected:
1015  template< class LocalBlock, class Operation >
1016  void applyToBlock ( const size_t row, const size_t col,
1017  const LocalBlock &localBlock,
1018  Operation& operation )
1019  {
1020  LittleBlockType& block = ( implicitModeActive() ) ? matrix().entry( row, col ) : matrix()[ row ][ col ];
1021  for( int i = 0; i < littleRows; ++i )
1022  for( int j = 0; j < littleCols; ++j )
1023  operation( block[ i ][ j ], localBlock[ i ][ j ] );
1024  }
1025 
1026  template< class LocalBlock >
1027  void setBlock ( const size_t row, const size_t col,
1028  const LocalBlock &localBlock )
1029  {
1030  typedef typename DomainSpaceType :: RangeFieldType Field;
1031  auto copy = [] ( Field& a, const typename LocalBlock::field_type& b ) { a = b; };
1032  applyToBlock( row, col, localBlock, copy );
1033  }
1034 
1035  template< class LocalBlock >
1036  void addBlock ( const size_t row, const size_t col,
1037  const LocalBlock &localBlock )
1038  {
1039  typedef typename DomainSpaceType :: RangeFieldType Field;
1040  auto add = [] ( Field& a, const typename LocalBlock::field_type& b ) { a += b; };
1041  applyToBlock( row, col, localBlock, add );
1042  }
1043 
1044  public:
1045  template< class LocalMatrix, class Operation >
1046  void applyToLocalMatrix ( const DomainEntityType &domainEntity,
1047  const RangeEntityType &rangeEntity,
1048  const LocalMatrix &localMat,
1049  Operation& operation )
1050  {
1051  typedef typename MatrixType::size_type Index;
1052  if( implicitModeActive() )
1053  {
1054  auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
1055  {
1056  return matrix().entry( index.first, index.second );
1057  };
1058  auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local, const std::pair< Index, Index > &index )
1059  {
1060  LittleBlockType& block = blockAccess( index );
1061  for( int i = 0; i < littleRows; ++i )
1062  for( int j = 0; j < littleCols; ++j )
1063  operation( block[ i ][ j ], localMat.get( local.first * littleRows + i, local.second *littleCols + j ) );
1064  };
1065  rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1066  }
1067  else
1068  {
1069  auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
1070  {
1071  return matrix()[ index.first][ index.second ];
1072  };
1073  auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local, const std::pair< Index, Index > &index )
1074  {
1075  LittleBlockType& block = blockAccess( index );
1076  for( int i = 0; i < littleRows; ++i )
1077  for( int j = 0; j < littleCols; ++j )
1078  operation( block[ i ][ j ], localMat.get( local.first * littleRows + i, local.second *littleCols + j ) );
1079  };
1080  rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1081  }
1082  }
1083 
1084  template< class LocalMatrix >
1085  void setLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
1086  {
1087  typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1088  auto operation = [] ( RangeFieldType& a, const RangeFieldType& b )
1089  {
1090  a = b;
1091  };
1092  applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1093  }
1094 
1095  template< class LocalMatrix >
1096  void addLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
1097  {
1098  typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1099  auto operation = [] ( RangeFieldType& a, const RangeFieldType& b )
1100  {
1101  a += b;
1102  };
1103  applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1104  }
1105 
1106  template< class LocalMatrix, class Scalar >
1107  void addScaledLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat, const Scalar &s )
1108  {
1109  typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1110  auto operation = [ &s ] ( RangeFieldType& a, const RangeFieldType& b )
1111  {
1112  a += s * b;
1113  };
1114  applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1115  }
1116 
1117  template< class LocalMatrix >
1118  void getLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMat ) const
1119  {
1120  // make sure that matrix is in compressed state if build mode is implicit
1121  finalizeAssembly();
1122 
1123  typedef typename MatrixType::size_type Index;
1124  auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< Index, Index > &global )
1125  {
1126  for( std::size_t i = 0; i < littleRows; ++i )
1127  for( std::size_t j = 0; j < littleCols; ++j )
1128  localMat.set( local.first * littleRows + i, local.second *littleCols + j,
1129  matrix()[ global.first ][ global.second ][i][j] );
1130  };
1131 
1132  rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1133  }
1134 
1135  protected:
1136  void preConErrorMsg(int preCon) const
1137  {
1138  exit(1);
1139  }
1140 
1141  void removeObj ()
1142  {
1143  Dest_.reset( nullptr );
1144  Arg_.reset( nullptr );
1145  matrixAdap_.reset( nullptr );
1146  }
1147 
1148  void createBlockVectors () const
1149  {
1150  if( !Arg_ )
1151  Arg_.reset( new RowBlockVectorType( rowMapper_.size() ) );
1152  if( !Dest_ )
1153  Dest_.reset( new ColumnBlockVectorType( colMapper_.size() ) );
1154 
1155  createMatrixAdapter ();
1156  }
1157  };
1158 
1159  } // namespace Fem
1160 
1161 } // namespace Dune
1162 
1163 #endif // #if HAVE_DUNE_ISTL
1164 
1165 #endif // #ifndef DUNE_FEM_ISTLMATRIXWRAPPER_HH
Definition: bindguard.hh:11
std::tuple_element< i, Tuple >::type & get(Dune::TypeIndexedTuple< Tuple, Types > &tuple)
Definition: typeindexedtuple.hh:122
PairFunctor< Mapper, Entity, Functor > makePairFunctor(const Mapper &mapper, const Entity &entity, Functor functor)
Definition: operator/matrix/functor.hh:65