dune-fem 2.8.0
Loading...
Searching...
No Matches
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
35
39
40namespace 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 const auto spEnd = sparsityPattern.end();
142
143 // insert map of indices into matrix
144 auto endcreate = this->createend();
145 for(auto create = this->createbegin(); create != endcreate; ++create)
146 {
147 const auto row = sparsityPattern.find( create.index() );
148 // if a row is empty then do nothing
149 if( row == spEnd ) continue;
150
151 // insert all indices for this row
152 for ( const auto& col : row->second )
153 create.insert( col );
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 proposeImplicit = true )
886 {
887 // if grid sequence number changed, rebuild matrix
888 if(sequence_ != domainSpace().sequence())
889 {
890 removeObj();
891
892 // do not use implicit build mode when multi threading is enabled
893 const bool implicit = proposeImplicit && MPIManager::numThreads() == 1;
894 if( implicit )
895 {
896 auto nnz = stencil.maxNonZerosEstimate();
897 if( nnz == 0 )
898 {
899 Stencil tmpStencil( stencil );
900 tmpStencil.fill( *(domainSpace_.begin()), *(rangeSpace_.begin()) );
901 nnz = tmpStencil.maxNonZerosEstimate();
902 }
903 matrix_.reset( new MatrixType( rowMapper_.size(), colMapper_.size(), nnz, overflowFraction_ ) );
904 }
905 else
906 {
907 matrix_.reset( new MatrixType( rowMapper_.size(), colMapper_.size() ) );
908 matrix().createEntries( stencil.globalStencil() );
909 }
910
911 sequence_ = domainSpace().sequence();
912 }
913 }
914
916 template <class HangingNodesType>
917 void changeHangingNodes(const HangingNodesType& hangingNodes)
918 {
919 // create new matrix
920 MatrixType* newMatrix = new MatrixType(rowMapper_.size(), colMapper_.size());
921
922 // setup with hanging rows
923 newMatrix->setup( *matrix_ , hangingNodes );
924
925 // remove old matrix
926 removeObj();
927
928 // store new matrix
929 matrix_.reset( newMatrix );
930 }
931
934 void extractDiagonal( ColumnDiscreteFunctionType& diag ) const
935 {
936 // extract diagonal entries
937 matrix().extractDiagonal( diag.blockVector() );
938 }
939
941 bool rightPrecondition() const
942 {
943 return false;
944 }
945
947 void apply(const ColumnDiscreteFunctionType& arg, RowDiscreteFunctionType& dest) const
948 {
949 matrixAdapter().apply( arg.blockVector(), dest.blockVector() );
950 }
951
953 void resort()
954 {}
955
958 void createPreconditionMatrix()
959 {}
960
962 void print(std::ostream & s) const
963 {
964 matrix().print(s);
965 }
966
967 const DomainSpaceType& domainSpace() const
968 {
969 return domainSpace_;
970 }
971 const RangeSpaceType& rangeSpace() const
972 {
973 return rangeSpace_;
974 }
975
976 const RowMapperType& rowMapper() const
977 {
978 return rowMapper_;
979 }
980 const ColMapperType& colMapper() const
981 {
982 return colMapper_;
983 }
984
986 ObjectType* newObject() const
987 {
988 return new ObjectType(*this, domainSpace(), rangeSpace());
989 }
990
995 [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
996 LocalMatrixType localMatrix( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity ) const
997 {
998 return LocalMatrixType( localMatrixStack_, domainEntity, rangeEntity );
999 }
1000
1005 [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
1006 LocalMatrixType localMatrix() const
1007 {
1008 return LocalMatrixType( localMatrixStack_ );
1009 }
1010
1011 LocalColumnObjectType localColumn( const DomainEntityType &domainEntity ) const
1012 {
1013 return LocalColumnObjectType ( *this, domainEntity );
1014 }
1015
1016 protected:
1017 template< class LocalBlock, class Operation >
1018 void applyToBlock ( const size_t row, const size_t col,
1019 const LocalBlock &localBlock,
1020 Operation& operation )
1021 {
1022 LittleBlockType& block = ( implicitModeActive() ) ? matrix().entry( row, col ) : matrix()[ row ][ col ];
1023 for( int i = 0; i < littleRows; ++i )
1024 for( int j = 0; j < littleCols; ++j )
1025 operation( block[ i ][ j ], localBlock[ i ][ j ] );
1026 }
1027
1028 template< class LocalBlock >
1029 void setBlock ( const size_t row, const size_t col,
1030 const LocalBlock &localBlock )
1031 {
1032 typedef typename DomainSpaceType :: RangeFieldType Field;
1033 auto copy = [] ( Field& a, const typename LocalBlock::field_type& b ) { a = b; };
1034 applyToBlock( row, col, localBlock, copy );
1035 }
1036
1037 template< class LocalBlock >
1038 void addBlock ( const size_t row, const size_t col,
1039 const LocalBlock &localBlock )
1040 {
1041 typedef typename DomainSpaceType :: RangeFieldType Field;
1042 auto add = [] ( Field& a, const typename LocalBlock::field_type& b ) { a += b; };
1043 applyToBlock( row, col, localBlock, add );
1044 }
1045
1046 public:
1047 template< class LocalMatrix, class Operation >
1048 void applyToLocalMatrix ( const DomainEntityType &domainEntity,
1049 const RangeEntityType &rangeEntity,
1050 const LocalMatrix &localMat,
1051 Operation& operation )
1052 {
1053 typedef typename MatrixType::size_type Index;
1054 if( implicitModeActive() )
1055 {
1056 auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
1057 {
1058 return matrix().entry( index.first, index.second );
1059 };
1060 auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local, const std::pair< Index, Index > &index )
1061 {
1062 LittleBlockType& block = blockAccess( index );
1063 for( int i = 0; i < littleRows; ++i )
1064 for( int j = 0; j < littleCols; ++j )
1065 operation( block[ i ][ j ], localMat.get( local.first * littleRows + i, local.second *littleCols + j ) );
1066 };
1067 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1068 }
1069 else
1070 {
1071 auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
1072 {
1073 return matrix()[ index.first][ index.second ];
1074 };
1075 auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local, const std::pair< Index, Index > &index )
1076 {
1077 LittleBlockType& block = blockAccess( index );
1078 for( int i = 0; i < littleRows; ++i )
1079 for( int j = 0; j < littleCols; ++j )
1080 operation( block[ i ][ j ], localMat.get( local.first * littleRows + i, local.second *littleCols + j ) );
1081 };
1082 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1083 }
1084 }
1085
1086 template< class LocalMatrix >
1087 void setLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
1088 {
1089 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1090 auto operation = [] ( RangeFieldType& a, const RangeFieldType& b )
1091 {
1092 a = b;
1093 };
1094 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1095 }
1096
1097 template< class LocalMatrix >
1098 void addLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
1099 {
1100 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1101 auto operation = [] ( RangeFieldType& a, const RangeFieldType& b )
1102 {
1103 a += b;
1104 };
1105 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1106 }
1107
1108 template< class LocalMatrix, class Scalar >
1109 void addScaledLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat, const Scalar &s )
1110 {
1111 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1112 auto operation = [ &s ] ( RangeFieldType& a, const RangeFieldType& b )
1113 {
1114 a += s * b;
1115 };
1116 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1117 }
1118
1119 template< class LocalMatrix >
1120 void getLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMat ) const
1121 {
1122 // make sure that matrix is in compressed state if build mode is implicit
1123 finalizeAssembly();
1124
1125 typedef typename MatrixType::size_type Index;
1126 auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< Index, Index > &global )
1127 {
1128 for( std::size_t i = 0; i < littleRows; ++i )
1129 for( std::size_t j = 0; j < littleCols; ++j )
1130 localMat.set( local.first * littleRows + i, local.second *littleCols + j,
1131 matrix()[ global.first ][ global.second ][i][j] );
1132 };
1133
1134 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1135 }
1136
1137 protected:
1138 void preConErrorMsg(int preCon) const
1139 {
1140 exit(1);
1141 }
1142
1143 void removeObj ()
1144 {
1145 Dest_.reset( nullptr );
1146 Arg_.reset( nullptr );
1147 matrixAdap_.reset( nullptr );
1148 }
1149
1150 void createBlockVectors () const
1151 {
1152 if( !Arg_ )
1153 Arg_.reset( new RowBlockVectorType( rowMapper_.size() ) );
1154 if( !Dest_ )
1155 Dest_.reset( new ColumnBlockVectorType( colMapper_.size() ) );
1156
1157 createMatrixAdapter ();
1158 }
1159 };
1160
1161 } // namespace Fem
1162
1163} // namespace Dune
1164
1165#endif // #if HAVE_DUNE_ISTL
1166
1167#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