1#ifndef DUNE_FEM_ISTLMATRIXWRAPPER_HH
2#define DUNE_FEM_ISTLMATRIXWRAPPER_HH
16#include <dune/common/exceptions.hh>
17#include <dune/common/fmatrix.hh>
20#include <dune/istl/bvector.hh>
21#include <dune/istl/bcrsmatrix.hh>
22#include <dune/istl/preconditioners.hh>
48 template<
class MatrixObject >
49 class ISTLLocalMatrix;
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;
59 template <
class LittleBlockType,
class RowDiscreteFunctionImp,
class ColDiscreteFunctionImp = RowDiscreteFunctionImp>
60 class ImprovedBCRSMatrix :
public BCRSMatrix<LittleBlockType>
62 friend struct MatrixDimension<ImprovedBCRSMatrix>;
64 typedef RowDiscreteFunctionImp RowDiscreteFunctionType;
65 typedef ColDiscreteFunctionImp ColDiscreteFunctionType;
67 typedef BCRSMatrix<LittleBlockType> BaseType;
69 typedef BaseType MatrixBaseType;
70 typedef typename BaseType :: RowIterator RowIteratorType ;
71 typedef typename BaseType :: ColIterator ColIteratorType ;
73 typedef ImprovedBCRSMatrix< LittleBlockType, RowDiscreteFunctionImp, ColDiscreteFunctionImp > ThisType;
75 typedef typename BaseType :: size_type size_type;
80 typedef typename BaseType::field_type field_type;
83 typedef typename BaseType::block_type block_type;
86 typedef typename BaseType:: allocator_type allocator_type;
89 typedef typename BaseType :: row_type row_type;
92 typedef typename BaseType :: ColIterator ColIterator;
95 typedef typename BaseType :: ConstColIterator ConstColIterator;
98 typedef typename BaseType :: RowIterator RowIterator;
101 typedef typename BaseType :: ConstRowIterator ConstRowIterator;
104 typedef typename ColDiscreteFunctionType :: DiscreteFunctionSpaceType RangeSpaceType;
107 typedef typename RowDiscreteFunctionType :: DofStorageType RowBlockVectorType;
110 typedef typename ColDiscreteFunctionType :: DofStorageType ColBlockVectorType;
113 typedef typename RangeSpaceType :: GridType :: Traits :: CollectiveCommunication CollectiveCommunictionType ;
115 typedef typename BaseType :: BuildMode BuildMode ;
119 ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nnz,
double overflowFraction) :
120 BaseType (rows, cols, nnz, overflowFraction, BaseType::implicit)
124 ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nz = 0 ) :
125 BaseType (rows, cols, BaseType :: row_wise)
129 ImprovedBCRSMatrix( ) :
134 ImprovedBCRSMatrix(
const ImprovedBCRSMatrix& org) :
138 template <
class SparsityPattern >
139 void createEntries(
const SparsityPattern& sparsityPattern)
141 const auto spEnd = sparsityPattern.end();
144 auto endcreate = this->createend();
145 for(
auto create = this->createbegin(); create != endcreate; ++create)
147 const auto row = sparsityPattern.find( create.index() );
149 if( row == spEnd )
continue;
152 for (
const auto& col : row->second )
153 create.insert( col );
160 for (
auto& row : *
this)
161 for (
auto& entry : row)
166 void unitRow(
const size_t row )
168 block_type idBlock( 0 );
169 for (
int i = 0; i < idBlock.rows; ++i)
172 auto& matRow = (*this)[ row ];
173 auto colIt = matRow.begin();
174 const auto& colEndIt = matRow.end();
175 for (; colIt != colEndIt; ++colIt)
177 if( colIt.index() == row )
185 template <
class HangingNodesType>
186 void setup(ThisType& oldMatrix,
const HangingNodesType& hangingNodes)
189 typedef std::set< std::pair<int, block_type> > LocalEntryType;
190 typedef std::map< int , LocalEntryType > EntriesType;
194 std::map< int , std::set<int> > indices;
196 auto rowend = oldMatrix.end();
197 for(
auto it = oldMatrix.begin(); it != rowend; ++it)
199 const auto row = it.index();
200 auto& localIndices = indices[ row ];
202 if( hangingNodes.isHangingNode( row ) )
205 const auto& cols = hangingNodes.associatedDofs( row );
206 const auto colSize = cols.size();
207 for(
auto i=0; i<colSize; ++i)
209 assert( ! hangingNodes.isHangingNode( cols[i].first ) );
212 auto& localColIndices = indices[ cols[i].first ];
213 auto& localEntry = entries[ cols[i].first ];
216 auto endj = (*it).end();
217 for (
auto j= (*it).begin(); j!=endj; ++j)
219 localColIndices.insert( j.index () );
220 localEntry.insert( std::make_pair( j.index(), (cols[i].second * (*j)) ));
225 localIndices.insert( row );
226 for(
auto i=0; i<colSize; ++i)
227 localIndices.insert( cols[i].first );
232 auto endj = (*it).end();
233 for (
auto j= (*it).begin(); j!=endj; ++j)
234 localIndices.insert( j.index () );
239 createEntries( indices );
242 auto rowit = oldMatrix.begin();
244 auto endcreate = this->end();
245 for(
auto create = this->begin(); create != endcreate; ++create, ++rowit )
247 assert( rowit != oldMatrix.end() );
249 const auto row = create.index();
250 if( hangingNodes.isHangingNode( row ) )
252 const auto& cols = hangingNodes.associatedDofs( row );
254 std::map< const int , block_type > colMap;
256 assert( block_type :: rows == 1 );
258 const auto colSize = cols.size();
259 for(
auto i=0; i<colSize; ++i)
260 colMap[ cols[i].first ] = -cols[i].second;
264 auto endj = (*create).end();
265 for (
auto j= (*create).begin(); j!=endj; ++j)
267 assert( colMap.find( j.index() ) != colMap.end() );
268 (*j) = colMap[ j.index() ];
272 else if ( entries.find( row ) == entries.end() )
274 auto colit = (*rowit).begin();
275 auto endj = (*create).end();
276 for (
auto j= (*create).begin(); j!=endj; ++j, ++colit )
278 assert( colit != (*rowit).end() );
284 std::map< int , block_type > oldCols;
287 auto colend = (*rowit).end();
288 for(
auto colit = (*rowit).begin(); colit != colend; ++colit)
289 oldCols[ colit.index() ] = 0;
292 auto entry = entries.find( row );
293 assert( entry != entries.end ());
296 auto endcol = (*entry).second.end();
297 for(
auto co = (*entry).second.begin(); co != endcol; ++co)
298 oldCols[ (*co).first ] = 0;
302 auto colend = (*rowit).end();
303 for(
auto colit = (*rowit).begin(); colit != colend; ++colit)
304 oldCols[ colit.index() ] += (*colit);
308 auto endcol = (*entry).second.end();
309 for(
auto co = (*entry).second.begin(); co != endcol; ++co)
310 oldCols[ (*co).first ] += (*co).second;
313 auto endj = (*create).end();
314 for (
auto j= (*create).begin(); j!=endj; ++j )
316 auto colEntry = oldCols.find( j.index() );
317 if( colEntry != oldCols.end() )
318 (*j) = (*colEntry).second;
327 void extractDiagonal( ColBlockVectorType& diag )
const
329 const auto endi = this->end();
330 for (
auto i = this->begin(); i!=endi; ++i)
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 ];
344 field_type operator()(
const std::size_t row,
const std::size_t col)
const
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));
351 const auto& matrixRow(this->
operator[](blockRow));
352 auto entry = matrixRow.find( blockCol );
353 const LittleBlockType& block = (*entry);
354 return block[localRowIdx][localColIdx];
359 void set(
const std::size_t row,
const std::size_t col, field_type value)
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));
366 auto& matrixRow(this->
operator[](blockRow));
367 auto entry = matrixRow.find( blockCol );
368 LittleBlockType& block = (*entry);
369 block[localRowIdx][localColIdx] = value;
373 void print(std::ostream& s=std::cout,
unsigned int offset=0)
const
376 const auto endi=this->end();
377 for (
auto i=this->begin(); i!=endi; ++i)
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;
392 template<
class MatrixObject >
393 struct ISTLLocalMatrixTraits
395 typedef typename MatrixObject::DomainSpaceType DomainSpaceType;
396 typedef typename MatrixObject::RangeSpaceType RangeSpaceType;
397 typedef typename DomainSpaceType::RangeFieldType RangeFieldType;
399 typedef ISTLLocalMatrix< MatrixObject > LocalMatrixType;
400 typedef typename MatrixObject::MatrixType::block_type LittleBlockType;
407 template<
class MatrixObject >
408 class ISTLLocalMatrix
409 :
public LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > >
413 typedef LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > > BaseType;
416 typedef MatrixObject MatrixObjectType;
418 typedef typename MatrixObjectType::MatrixType MatrixType;
420 typedef typename MatrixType::block_type LittleBlockType;
422 typedef typename MatrixType :: size_type size_type;
424 typedef typename MatrixObjectType::DomainSpaceType DomainSpaceType;
425 typedef typename MatrixObjectType::RangeSpaceType RangeSpaceType;
427 typedef typename MatrixObjectType::DomainEntityType DomainEntityType;
428 typedef typename MatrixObjectType::RangeEntityType RangeEntityType;
431 typedef typename DomainSpaceType::RangeFieldType DofType;
432 typedef typename MatrixType::row_type RowType;
435 typedef typename DomainSpaceType::BlockMapperType ColMapperType;
437 typedef typename RangeSpaceType::BlockMapperType RowMapperType;
439 static const int littleCols = MatrixObjectType::littleCols;
440 static const int littleRows = MatrixObjectType::littleRows;
442 typedef typename MatrixType::size_type Index;
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() ),
460 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
461 ISTLLocalMatrix (
const ISTLLocalMatrix& 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_)
472 void init (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
475 BaseType :: init ( domainEntity, rangeEntity );
477 numRows_ = rowMapper_.numDofs( rangeEntity );
478 numCols_ = colMapper_.numDofs( domainEntity );
481 matrices_.resize( numRows_ );
482 for(
auto& row : matrices_ )
484 row.resize( numCols_,
nullptr );
487 if( matrixObj_.implicitModeActive() )
489 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
491 return matrixObj_.matrix().entry( index.first, index.second );
493 initBlocks( blockAccess, domainEntity, rangeEntity );
497 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
499 return matrixObj_.matrix()[ index.first ][ index.second ];
501 initBlocks( blockAccess, domainEntity, rangeEntity );
505 template <
class BlockAccess>
506 void initBlocks( BlockAccess& blockAccess,
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
508 auto functor = [
this, &blockAccess ] ( std::pair< int, int > local,
const std::pair< Index, Index > &index )
510 matrices_[ local.first ][ local.second ] = &blockAccess( index );
513 rowMapper_.mapEach( rangeEntity,
makePairFunctor( colMapper_, domainEntity, functor ) );
518 void check(
int localRow,
int localCol)
const
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 );
532 DofType& getValue(
const int localRow,
const int localCol)
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];
542 const DofType
get(
const int localRow,
const int localCol)
const
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];
551 void scale (
const DofType& scalar)
553 const size_type nrows = matrices_.size();
554 for(size_type i=0; i<nrows; ++i)
556 const size_type ncols = matrices_[i].size();
557 for(size_type j=0; j<ncols; ++j)
559 (*matrices_[i][j]) *= scalar;
564 void add(
const int localRow,
const int localCol ,
const DofType value)
567 check(localRow,localCol);
569 getValue(localRow,localCol) += value;
572 void set(
const int localRow,
const int localCol ,
const DofType value)
575 check(localRow,localCol);
577 getValue(localRow,localCol) = value;
581 void unitRow(
const int localRow)
583 const int row = (int) localRow / littleRows;
584 const int lRow = localRow%littleRows;
587 doClearRow( row, lRow );
590 (*matrices_[row][row])[lRow][lRow] = 1;
596 const size_type nrows = matrices_.size();
597 for(size_type i=0; i<nrows; ++i)
599 const size_type ncols = matrices_[i].size();
600 for(size_type j=0; j<ncols; ++j)
602 (*matrices_[i][j]) = (DofType) 0;
608 void clearRow (
const int localRow )
610 const int row = (int) localRow / littleRows;
611 const int lRow = localRow%littleRows;
614 doClearRow( row, lRow );
623 void doClearRow (
const int row,
const int lRow )
626 const auto col = this->columns();
627 for(
auto localCol=0; localCol<col; ++localCol)
629 const int col = (int) localCol / littleCols;
630 const int lCol = localCol%littleCols;
631 (*matrices_[row][col])[lRow][lCol] = 0;
637 const RowMapperType& rowMapper_;
638 const ColMapperType& colMapper_;
645 std::vector< std::vector< LittleBlockType* > > matrices_;
648 const MatrixObjectType& matrixObj_;
653 template <
class DomainSpaceImp,
class RangeSpaceImp,
class DomainBlock,
class RangeBlock >
654 class ISTLMatrixObject
658 typedef DomainSpaceImp DomainSpaceType;
660 typedef RangeSpaceImp RangeSpaceType;
663 typedef ISTLMatrixObject< DomainSpaceImp, RangeSpaceImp, DomainBlock, RangeBlock > ThisType;
664 typedef ThisType PreconditionMatrixType;
666 typedef typename DomainSpaceType::GridType GridType;
668 typedef typename RangeSpaceType :: EntityType RangeEntityType ;
669 typedef typename DomainSpaceType :: EntityType DomainEntityType ;
671 enum { littleCols = DomainSpaceType :: localBlockSize };
672 enum { littleRows = RangeSpaceType :: localBlockSize };
674 typedef FieldMatrix<typename DomainSpaceType :: RangeFieldType, littleRows, littleCols> LittleBlockType;
675 typedef LittleBlockType block_type;
676 typedef LittleBlockType MatrixBlockType;
678 typedef ISTLBlockVectorDiscreteFunction< RangeSpaceType, RangeBlock > RowDiscreteFunctionType;
679 typedef ISTLBlockVectorDiscreteFunction< DomainSpaceType, DomainBlock > ColumnDiscreteFunctionType;
681 typedef typename RowDiscreteFunctionType :: DofStorageType RowBlockVectorType;
682 typedef typename ColumnDiscreteFunctionType :: DofStorageType ColumnBlockVectorType;
684 typedef typename RangeSpaceType :: BlockMapperType RowMapperType;
685 typedef typename DomainSpaceType :: BlockMapperType ColMapperType;
688 typedef ImprovedBCRSMatrix< LittleBlockType , ColumnDiscreteFunctionType , RowDiscreteFunctionType > MatrixType;
689 typedef ISTLParallelMatrixAdapterInterface< MatrixType > MatrixAdapterType;
692 typedef ISTLLocalMatrix<ThisType> ObjectType;
693 friend class ISTLLocalMatrix<ThisType>;
695 typedef ThisType LocalMatrixFactoryType;
698 typedef LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType;
699 typedef ColumnObject< ThisType > LocalColumnObjectType;
702 const DomainSpaceType& domainSpace_;
703 const RangeSpaceType& rangeSpace_;
706 RowMapperType& rowMapper_;
708 ColMapperType& colMapper_;
713 mutable std::unique_ptr< MatrixType > matrix_;
715 mutable LocalMatrixStackType localMatrixStack_;
717 mutable std::unique_ptr< MatrixAdapterType > matrixAdap_;
718 mutable std::unique_ptr< ColumnBlockVectorType > Arg_;
719 mutable std::unique_ptr< RowBlockVectorType > Dest_;
721 const double overflowFraction_;
722 ISTLSolverParameter param_;
724 ISTLMatrixObject(
const ISTLMatrixObject&) =
delete;
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() )
736 , localMatrixStack_( *this )
737 , overflowFraction_( param.overflowFraction() )
744 MatrixType& matrix()
const
751 void printTexInfo(std::ostream& out)
const
753 out <<
"ISTL MatrixObj: ";
758 std::string preconditionName()
const
763 void createMatrixAdapter (
const ISTLSolverParameter& param )
const
767 matrixAdap_ = ISTLMatrixAdapterFactory< ThisType > :: matrixAdapter( *
this, param );
769 assert( matrixAdap_ );
773 MatrixAdapterType& matrixAdapter(
const ISTLSolverParameter& parameter )
const
775 const ISTLSolverParameter* param =
dynamic_cast< const ISTLSolverParameter*
> (¶meter);
777 createMatrixAdapter( *param );
780 ISTLSolverParameter newparam( parameter );
781 createMatrixAdapter( newparam );
789 MatrixAdapterType& matrixAdapter()
const
791 return matrixAdapter( param_ );
795 MatrixType& exportMatrix()
const
801 bool implicitModeActive()
const
807 if( matrix().buildMode() == MatrixType::implicit && matrix().buildStage() == MatrixType::building )
813 void finalizeAssembly()
const
815 const_cast< ThisType&
> (*this).compress();
821 if( implicitModeActive() )
835 void unitRow(
const size_t row )
837 matrix().unitRow( row );
840 template <
class Vector>
841 void setUnitRows(
const Vector &rows )
843 const auto &auxiliaryDofs = domainSpace().auxiliaryDofs();
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();
854 for (
auto col=row.begin(); col!=endcol; ++col)
856 for (
auto& entry : (*col)[localRowIdx])
858 if (col.index() == blockRow)
860 (*col)[localRowIdx][localRowIdx] = auxiliaryDofs.contains( r )? 0.0 : 1.0;
873 reserve( Stencil<DomainSpaceType,RangeSpaceType>(domainSpace_, rangeSpace_),
true );
878 void reserve (
const std::vector< Set >& sparsityPattern )
880 reserve( StencilWrapper< DomainSpaceType,RangeSpaceType, Set >( sparsityPattern ),
false );
884 template <
class Stencil>
885 void reserve(
const Stencil &stencil,
const bool proposeImplicit =
true )
888 if(sequence_ != domainSpace().sequence())
893 const bool implicit = proposeImplicit && MPIManager::numThreads() == 1;
896 auto nnz = stencil.maxNonZerosEstimate();
899 Stencil tmpStencil( stencil );
900 tmpStencil.fill( *(domainSpace_.begin()), *(rangeSpace_.begin()) );
901 nnz = tmpStencil.maxNonZerosEstimate();
903 matrix_.reset(
new MatrixType( rowMapper_.size(), colMapper_.size(), nnz, overflowFraction_ ) );
907 matrix_.reset(
new MatrixType( rowMapper_.size(), colMapper_.size() ) );
908 matrix().createEntries( stencil.globalStencil() );
911 sequence_ = domainSpace().sequence();
916 template <
class HangingNodesType>
917 void changeHangingNodes(
const HangingNodesType& hangingNodes)
920 MatrixType* newMatrix =
new MatrixType(rowMapper_.size(), colMapper_.size());
923 newMatrix->setup( *matrix_ , hangingNodes );
929 matrix_.reset( newMatrix );
934 void extractDiagonal( ColumnDiscreteFunctionType& diag )
const
937 matrix().extractDiagonal( diag.blockVector() );
941 bool rightPrecondition()
const
947 void apply(
const ColumnDiscreteFunctionType& arg, RowDiscreteFunctionType& dest)
const
949 matrixAdapter().apply( arg.blockVector(), dest.blockVector() );
958 void createPreconditionMatrix()
962 void print(std::ostream & s)
const
967 const DomainSpaceType& domainSpace()
const
971 const RangeSpaceType& rangeSpace()
const
976 const RowMapperType& rowMapper()
const
980 const ColMapperType& colMapper()
const
986 ObjectType* newObject()
const
988 return new ObjectType(*
this, domainSpace(), rangeSpace());
995 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
996 LocalMatrixType localMatrix(
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
const
998 return LocalMatrixType( localMatrixStack_, domainEntity, rangeEntity );
1005 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
1006 LocalMatrixType localMatrix()
const
1008 return LocalMatrixType( localMatrixStack_ );
1011 LocalColumnObjectType localColumn(
const DomainEntityType &domainEntity )
const
1013 return LocalColumnObjectType ( *
this, domainEntity );
1017 template<
class LocalBlock,
class Operation >
1018 void applyToBlock (
const size_t row,
const size_t col,
1019 const LocalBlock &localBlock,
1020 Operation& operation )
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 ] );
1028 template<
class LocalBlock >
1029 void setBlock (
const size_t row,
const size_t col,
1030 const LocalBlock &localBlock )
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 );
1037 template<
class LocalBlock >
1038 void addBlock (
const size_t row,
const size_t col,
1039 const LocalBlock &localBlock )
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 );
1047 template<
class LocalMatrix,
class Operation >
1048 void applyToLocalMatrix (
const DomainEntityType &domainEntity,
1049 const RangeEntityType &rangeEntity,
1050 const LocalMatrix &localMat,
1051 Operation& operation )
1053 typedef typename MatrixType::size_type Index;
1054 if( implicitModeActive() )
1056 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
1058 return matrix().entry( index.first, index.second );
1060 auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local,
const std::pair< Index, Index > &index )
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 ) );
1067 rowMapper_.mapEach( rangeEntity,
makePairFunctor( colMapper_, domainEntity, functor ) );
1071 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
1073 return matrix()[ index.first][ index.second ];
1075 auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local,
const std::pair< Index, Index > &index )
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 ) );
1082 rowMapper_.mapEach( rangeEntity,
makePairFunctor( colMapper_, domainEntity, functor ) );
1086 template<
class LocalMatrix >
1087 void setLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
1089 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1090 auto operation = [] ( RangeFieldType& a,
const RangeFieldType& b )
1094 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1097 template<
class LocalMatrix >
1098 void addLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
1100 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1101 auto operation = [] ( RangeFieldType& a,
const RangeFieldType& b )
1105 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1108 template<
class LocalMatrix,
class Scalar >
1109 void addScaledLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat,
const Scalar &s )
1111 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1112 auto operation = [ &s ] ( RangeFieldType& a,
const RangeFieldType& b )
1116 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1119 template<
class LocalMatrix >
1120 void getLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity, LocalMatrix &localMat )
const
1125 typedef typename MatrixType::size_type Index;
1126 auto functor = [ &localMat, this ] ( std::pair< int, int > local,
const std::pair< Index, Index > &global )
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] );
1134 rowMapper_.mapEach( rangeEntity,
makePairFunctor( colMapper_, domainEntity, functor ) );
1138 void preConErrorMsg(
int preCon)
const
1145 Dest_.reset(
nullptr );
1146 Arg_.reset(
nullptr );
1147 matrixAdap_.reset(
nullptr );
1150 void createBlockVectors ()
const
1153 Arg_.reset(
new RowBlockVectorType( rowMapper_.size() ) );
1155 Dest_.reset(
new ColumnBlockVectorType( colMapper_.size() ) );
1157 createMatrixAdapter ();
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