1#ifndef DUNE_FEM_PETSCLINEAROPERATOR_HH
2#define DUNE_FEM_PETSCLINEAROPERATOR_HH
13#include <dune/common/dynmatrix.hh>
43 struct PetscSolverParameter :
public LocalParameter< SolverParameter, PetscSolverParameter >
45 typedef LocalParameter< SolverParameter, PetscSolverParameter > BaseType;
48 using BaseType :: parameter ;
49 using BaseType :: keyPrefix ;
52 : BaseType( parameter )
55 PetscSolverParameter(
const SolverParameter& sp )
56 : PetscSolverParameter( sp.keyPrefix(), sp.parameter() )
60 : BaseType( keyPrefix, parameter )
63 bool isPetscSolverParameter()
const {
return true; }
65 static const int boomeramg = 0;
66 static const int parasails = 1;
67 static const int pilut = 2;
69 int hypreMethod()
const
71 const std::string hyprePCNames[] = {
"boomer-amg",
"parasails",
"pilu-t" };
73 if (parameter().exists(
"petsc.preconditioning.method"))
75 hypreType = parameter().getEnum(
"petsc.preconditioning.hypre.method", hyprePCNames, 0 );
76 std::cout <<
"WARNING: using deprecated parameter 'petsc.preconditioning.hypre.method' use "
77 << keyPrefix() <<
"preconditioning.hypre.method instead\n";
80 hypreType = parameter().getEnum( keyPrefix()+
"hypre.method", hyprePCNames, 0 );
84 int superluMethod()
const
86 const std::string factorizationNames[] = {
"petsc",
"superlu",
"mumps" };
88 if (parameter().exists(
"petsc.preconditioning.lu.method"))
90 factorType = parameter().getEnum(
"petsc.preconditioning.lu.method", factorizationNames, 0 );
91 std::cout <<
"WARNING: using deprecated parameter 'petsc.preconditioning.lu.method' use "
92 << keyPrefix() <<
"preconditioning.lu.method instead\n";
95 factorType = parameter().getEnum( keyPrefix()+
"preconditioning.lu.method", factorizationNames, 0 );
99 bool viennaCL ()
const {
100 return parameter().getValue<
bool > ( keyPrefix() +
"petsc.viennacl", false );
102 bool blockedMode ()
const {
103 return parameter().getValue<
bool > ( keyPrefix() +
"petsc.blockedmode", true );
112 template<
typename DomainFunction,
typename RangeFunction >
113 class PetscLinearOperator
114 :
public Fem::AssembledOperator< DomainFunction, RangeFunction >
116 typedef PetscLinearOperator< DomainFunction, RangeFunction > ThisType;
118 typedef Mat MatrixType;
119 typedef DomainFunction DomainFunctionType;
120 typedef RangeFunction RangeFunctionType;
121 typedef typename DomainFunctionType::RangeFieldType DomainFieldType;
122 typedef typename RangeFunctionType::RangeFieldType RangeFieldType;
123 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType;
124 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType;
126 typedef PetscDiscreteFunction< DomainSpaceType > PetscDomainFunctionType;
127 typedef PetscDiscreteFunction< RangeSpaceType > PetscRangeFunctionType;
129 typedef typename DomainSpaceType::GridPartType::template Codim< 0 >::EntityType DomainEntityType;
130 typedef typename RangeSpaceType::GridPartType::template Codim< 0 >::EntityType RangeEntityType;
132 static const unsigned int domainLocalBlockSize = DomainSpaceType::localBlockSize;
133 static const unsigned int rangeLocalBlockSize = RangeSpaceType::localBlockSize;
135 static constexpr bool squareBlocked = domainLocalBlockSize == rangeLocalBlockSize ;
136 static constexpr bool blockedMatrix = domainLocalBlockSize > 1 && squareBlocked;
140 typedef FlatFieldMatrix< RangeFieldType, domainLocalBlockSize, rangeLocalBlockSize > MatrixBlockType;
141 typedef MatrixBlockType block_type;
144 enum Status {statAssembled=0,statAdd=1,statInsert=2,statGet=3,statNothing=4};
146 typedef PetscMappers< DomainSpaceType > DomainMappersType;
147 typedef PetscMappers< RangeSpaceType > RangeMappersType;
149 typedef AuxiliaryDofs< typename DomainSpaceType::GridPartType, typename DomainMappersType::GhostMapperType > DomainAuxiliaryDofsType;
150 typedef AuxiliaryDofs< typename RangeSpaceType::GridPartType, typename RangeMappersType::GhostMapperType > RangeAuxiliaryDofsType;
156 struct LocalMatrixTraits
158 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainSpaceType;
159 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeSpaceType;
160 typedef LocalMatrix LocalMatrixType;
161 typedef typename RangeSpaceType::RangeFieldType RangeFieldType;
164 typedef RangeFieldType LittleBlockType;
168 typedef LocalMatrix ObjectType;
169 typedef ThisType LocalMatrixFactoryType;
173 typedef LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType;
174 typedef ColumnObject< ThisType > LocalColumnObjectType;
176 PetscLinearOperator (
const DomainSpaceType &domainSpace,
const RangeSpaceType &rangeSpace,
177 const PetscSolverParameter& param = PetscSolverParameter() )
178 : domainMappers_( domainSpace ),
179 rangeMappers_( rangeSpace ),
180 localMatrixStack_( *this ),
181 status_(statNothing),
182 viennaCL_( param.viennaCL() ),
183 blockedMode_( blockedMatrix && (!viennaCL_) && param.blockedMode() )
186 PetscLinearOperator (
const std::string &,
const DomainSpaceType &domainSpace,
const RangeSpaceType &rangeSpace,
187 const PetscSolverParameter& param = PetscSolverParameter() )
188 : PetscLinearOperator( domainSpace, rangeSpace, param )
192 ~PetscLinearOperator ()
199 ::Dune::Petsc::MatAssemblyBegin( petscMatrix_, MAT_FLUSH_ASSEMBLY );
200 ::Dune::Petsc::MatAssemblyEnd ( petscMatrix_, MAT_FLUSH_ASSEMBLY );
201 setStatus( statAssembled );
206 if( ! finalizedAlready() )
208 ::Dune::Petsc::MatAssemblyBegin( petscMatrix_, MAT_FINAL_ASSEMBLY );
209 ::Dune::Petsc::MatAssemblyEnd ( petscMatrix_, MAT_FINAL_ASSEMBLY );
211 if( ! unitRows_.empty() )
213 ::Dune::Petsc::MatZeroRows( petscMatrix_, unitRows_.size(), unitRows_.data(), 1. );
221 bool finalizedAlready()
const
223 PetscBool assembled = PETSC_FALSE ;
224 ::Dune::Petsc::MatAssembled( petscMatrix_, &assembled );
225 return assembled == PETSC_TRUE;
228 void finalizeAssembly ()
const
230 const_cast< ThisType&
> (*this).finalize();
234 const DomainSpaceType &domainSpace ()
const {
return domainMappers_.space(); }
235 const RangeSpaceType &rangeSpace ()
const {
return rangeMappers_.space(); }
240 template <
class DF,
class RF>
241 void apply (
const DF &arg, RF &dest )
const
244 petscArg_.reset(
new PetscDomainFunctionType(
"PetscOp-arg", domainSpace() ) );
246 petscDest_.reset(
new PetscRangeFunctionType(
"PetscOp-arg", rangeSpace() ) );
248 petscArg_->assign( arg );
249 apply( *petscArg_, *petscDest_ );
250 dest.assign( *petscDest_ );
254 void apply (
const PetscDomainFunctionType &arg, PetscRangeFunctionType &dest )
const
258 ::Dune::Petsc::MatMult( petscMatrix_, *arg.petscVec() , *dest.petscVec() );
261 void operator() (
const DomainFunctionType &arg, RangeFunctionType &dest )
const
268 reserve( SimpleStencil<DomainSpaceType,RangeSpaceType>(0),
true );
272 void reserve (
const std::vector< Set >& sparsityPattern )
274 reserve( StencilWrapper< DomainSpaceType,RangeSpaceType, Set >( sparsityPattern ),
true );
278 template <
class StencilType>
279 void reserve (
const StencilType &stencil,
const bool isSimpleStencil =
false )
281 domainMappers_.update();
282 rangeMappers_.update();
284 if(sequence_ != domainSpace().sequence())
296 ::Dune::Petsc::MatCreate( domainSpace().gridPart().comm(), &petscMatrix_ );
299 const PetscInt localCols = domainMappers_.ghostMapper().interiorSize() * domainLocalBlockSize;
300 const PetscInt localRows = rangeMappers_.ghostMapper().interiorSize() * rangeLocalBlockSize;
302 const PetscInt globalCols = domainMappers_.parallelMapper().size() * domainLocalBlockSize;
303 const PetscInt globalRows = rangeMappers_.parallelMapper().size() * rangeLocalBlockSize;
305 ::Dune::Petsc::MatSetSizes( petscMatrix_, localRows, localCols, globalRows, globalCols );
310 ::Dune::Petsc::MatSetType( petscMatrix_, MATAIJVIENNACL );
312 else if( blockedMode_ )
314 bs = domainLocalBlockSize ;
315 ::Dune::Petsc::MatSetType( petscMatrix_, MATBAIJ );
317 ::Dune::Petsc::MatSetBlockSize( petscMatrix_, bs );
321 ::Dune::Petsc::MatSetType( petscMatrix_, MATAIJ );
339 if( isSimpleStencil || std::is_same< StencilType,SimpleStencil<DomainSpaceType,RangeSpaceType> >::value )
341 ::Dune::Petsc::MatSetUp( petscMatrix_, bs, domainLocalBlockSize * stencil.maxNonZerosEstimate() );
345 DomainAuxiliaryDofsType domainAuxiliaryDofs( domainMappers_.ghostMapper() );
346 RangeAuxiliaryDofsType rangeAuxiliaryDofs( rangeMappers_.ghostMapper() );
348 std::vector< PetscInt > d_nnz( localRows / bs, 0 );
349 std::vector< PetscInt > o_nnz( localRows / bs, 0 );
350 for(
const auto entry : stencil.globalStencil() )
352 const int petscIndex = rangeMappers_.ghostIndex( entry.first );
353 if( rangeAuxiliaryDofs.contains( petscIndex ) )
356 for (
unsigned int rb = 0; rb<rangeLocalBlockSize/bs; ++rb)
358 const size_t row = petscIndex*rangeLocalBlockSize/bs + rb;
361 assert( row <
size_t(localRows/bs) );
362 d_nnz[ row ] = o_nnz[ row ] = 0;
363 for(
const auto local : entry.second )
365 if( !domainAuxiliaryDofs.contains( domainMappers_.ghostIndex( local ) ) )
366 d_nnz[ row ] += domainLocalBlockSize/bs;
368 o_nnz[ row ] += domainLocalBlockSize/bs;
372 ::Dune::Petsc::MatSetUp( petscMatrix_, bs, &d_nnz[0], &o_nnz[0] );
374 sequence_ = domainSpace().sequence();
378 status_ = statAssembled ;
383 ::Dune::Petsc::MatZeroEntries( petscMatrix_ );
387 template <
class Vector>
388 void setUnitRows(
const Vector &rows )
390 std::vector< PetscInt > r( rows.size() );
391 for( std::size_t i =0 ; i< rows.size(); ++i )
393 const PetscInt block = rangeMappers_.parallelIndex( rows[ i ] / rangeLocalBlockSize );
394 r[ i ] = block * rangeLocalBlockSize + (rows[ i ] % rangeLocalBlockSize);
397 if( finalizedAlready() )
399 ::Dune::Petsc::MatZeroRows( petscMatrix_, r.size(), r.data(), 1. );
403 unitRows_.reserve( unitRows_.size() + r.size() );
404 for(
const auto& row : r )
405 unitRows_.push_back( row );
410 ObjectType* newObject()
const
412 return new ObjectType( *
this, domainSpace(), rangeSpace() );
419 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
420 LocalMatrixType localMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
const
422 return LocalMatrixType(localMatrixStack_, domainEntity, rangeEntity);
425 LocalColumnObjectType localColumn(
const DomainEntityType &colEntity )
const
427 return LocalColumnObjectType ( *
this, colEntity );
431 void unitRow(
const PetscInt localRow,
const PetscScalar diag = 1.0 )
433 std::array< PetscInt, domainLocalBlockSize > rows;
434 const PetscInt row = rangeMappers_.parallelIndex( localRow );
435 for(
unsigned int i=0, r = row * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r )
440 if( finalizedAlready() )
443 ::Dune::Petsc::MatZeroRows( petscMatrix_, domainLocalBlockSize, rows.data(), diag );
448 assert(
std::abs( diag - 1. ) < 1e-12 );
449 unitRows_.reserve( unitRows_.size() + domainLocalBlockSize );
450 for(
const auto& r : rows )
452 unitRows_.push_back( r );
457 bool blockedMode()
const {
return blockedMode_; }
460 template<
class PetscOp >
461 void applyToBlock (
const PetscInt localRow,
const PetscInt localCol,
const MatrixBlockType& block, PetscOp op )
464 const PetscInt localCols = domainMappers_.ghostMapper().interiorSize() * domainLocalBlockSize;
465 const PetscInt localRows = rangeMappers_.ghostMapper().interiorSize() * rangeLocalBlockSize;
466 assert( localRow < localRows );
467 assert( localCol < localCols );
473 const PetscInt row = rangeMappers_.parallelIndex( localRow );
474 const PetscInt col = rangeMappers_.parallelIndex( localCol );
475 ::Dune::Petsc::MatSetValuesBlocked( petscMatrix_, 1, &row, 1, &col, block.data(), op );
480 const PetscInt row = rangeMappers_.parallelIndex( localRow );
481 const PetscInt col = rangeMappers_.parallelIndex( localCol );
482 std::array< PetscInt, domainLocalBlockSize > rows;
483 std::array< PetscInt, domainLocalBlockSize > cols;
484 for(
unsigned int i=0, r = row * domainLocalBlockSize, c = col * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r, ++c )
491 ::Dune::Petsc::MatSetValues( petscMatrix_, domainLocalBlockSize, rows.data(), domainLocalBlockSize, cols.data(), block.data(), op );
493 setStatus( statAssembled );
496 template<
class LocalBlock,
class PetscOp >
497 void applyToBlock (
const size_t row,
const size_t col,
const LocalBlock& block, PetscOp op )
499 assert( block.rows() == rangeLocalBlockSize );
500 assert( block.cols() == domainLocalBlockSize );
503 MatrixBlockType matBlock( block );
504 applyToBlock( row, col, matBlock, op );
507 template<
class LocalBlock >
508 void setBlock (
const size_t row,
const size_t col,
const LocalBlock& block )
511 assert( status_==statAssembled || status_==statInsert );
516 setStatus( statInsert );
517 applyToBlock(
static_cast< PetscInt
> (row),
static_cast< PetscInt
> (col), block, INSERT_VALUES );
520 template<
class LocalBlock >
521 void addBlock (
const size_t row,
const size_t col,
const LocalBlock& block )
524 assert( status_==statAssembled || status_==statInsert );
529 setStatus( statAdd );
530 applyToBlock(
static_cast< PetscInt
> (row),
static_cast< PetscInt
> (col), block, ADD_VALUES );
534 typedef TemporaryLocalMatrix< DomainSpaceType, RangeSpaceType > TemporaryLocalMatrixType;
537 template <
class Operation>
538 const PetscScalar* getValues(
const unsigned int rSize,
const unsigned int cSize,
539 const TemporaryLocalMatrixType &localMat,
const Operation&,
540 const std::integral_constant< bool, false> nonscaled )
542 return localMat.data();
546 template <
class LM,
class Operation>
547 const PetscScalar* getValues(
const unsigned int rSize,
const unsigned int cSize,
548 const Assembly::Impl::LocalMatrixGetter< LM >& localMat,
const Operation&,
549 const std::integral_constant< bool, false> nonscaled )
551 return localMat.localMatrix().data();
555 template <
class LocalMatrix,
class Operation,
bool T>
556 const PetscScalar* getValues(
const unsigned int rSize,
const unsigned int cSize,
557 const LocalMatrix &localMat,
const Operation& operation,
558 const std::integral_constant< bool, T> scaled )
560 std::vector< PetscScalar >& v = v_;
561 v.resize( rSize * cSize );
562 for(
unsigned int i = 0, ic = 0 ; i< rSize; ++i )
564 for(
unsigned int j =0; j< cSize; ++j, ++ic )
566 v[ ic ] = operation( localMat.get( i, j ) );
572 template<
class LocalMatrix,
class Operation,
class PetscOp,
bool T >
573 void applyLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
574 const LocalMatrix &localMat,
const Operation& operation,
576 const std::integral_constant<bool, T> scaled )
578 std::vector< PetscInt >& r = r_;
579 std::vector< PetscInt >& c = c_;
583 setupIndicesBlocked( rangeMappers_, rangeEntity, r );
584 setupIndicesBlocked( domainMappers_, domainEntity, c );
587 const unsigned int rSize = r.size() * domainLocalBlockSize ;
588 const unsigned int cSize = c.size() * domainLocalBlockSize ;
590 const PetscScalar* values = getValues( rSize, cSize, localMat, operation, scaled );
591 ::Dune::Petsc::MatSetValuesBlocked( petscMatrix_, r.size(), r.data(), c.size(), c.data(), values, petscOp );
595 setupIndices( rangeMappers_, rangeEntity, r );
596 setupIndices( domainMappers_, domainEntity, c );
598 const unsigned int rSize = r.size();
599 const unsigned int cSize = c.size();
601 const PetscScalar* values = getValues( rSize, cSize, localMat, operation, scaled );
602 ::Dune::Petsc::MatSetValues( petscMatrix_, rSize, r.data(), cSize, c.data(), values, petscOp );
604 setStatus( statAssembled );
608 template<
class LocalMatrix >
609 void addLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
611 assert( status_==statAssembled || status_==statAdd );
612 setStatus( statAdd );
614 auto operation = [] (
const PetscScalar& value ) -> PetscScalar {
return value; };
616 applyLocalMatrix( domainEntity, rangeEntity, localMat, operation, ADD_VALUES, std::integral_constant< bool, false>() );
619 template<
class LocalMatrix,
class Scalar >
620 void addScaledLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat,
const Scalar &s )
622 assert( status_==statAssembled || status_==statAdd );
623 setStatus( statAdd );
625 auto operation = [ &s ] (
const PetscScalar& value ) -> PetscScalar {
return s * value; };
627 applyLocalMatrix( domainEntity, rangeEntity, localMat, operation, ADD_VALUES, std::integral_constant< bool, true>() );
630 template<
class LocalMatrix >
631 void setLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
633 assert( status_==statAssembled || status_==statInsert );
634 setStatus( statInsert );
636 auto operation = [] (
const PetscScalar& value ) -> PetscScalar {
return value; };
638 applyLocalMatrix( domainEntity, rangeEntity, localMat, operation, INSERT_VALUES, std::integral_constant< bool, false>() );
641 template<
class LocalMatrix >
642 void getLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity, LocalMatrix &localMat )
const
647 assert( status_==statAssembled || status_==statGet );
648 setStatus( statGet );
650 std::vector< PetscInt >& r = r_;
651 std::vector< PetscInt >& c = c_;
652 std::vector< PetscScalar >& v = v_;
654 setupIndices( rangeMappers_, rangeEntity, r );
655 setupIndices( domainMappers_, domainEntity, c );
657 v.resize( r.size() * c.size() );
658 ::Dune::Petsc::MatGetValues( petscMatrix_, r.size(), r.data(), c.size(), c.data(), v.data() );
660 for( std::size_t i =0 ; i< r.size(); ++i )
661 for( std::size_t j =0; j< c.size(); ++j )
662 localMat.set( i, j, v[ i * c.size() +j ] );
664 setStatus( statAssembled );
667 void scaleLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const RangeFieldType &s )
672 assert( status_==statAssembled || status_==statGet );
673 setStatus( statGet );
675 std::vector< PetscInt >& r = r_;
676 std::vector< PetscInt >& c = c_;
677 std::vector< PetscScalar >& v = v_;
679 setupIndices( rangeMappers_, rangeEntity, r );
680 setupIndices( domainMappers_, domainEntity, c );
682 const unsigned int rSize = r.size();
683 const unsigned int cSize = c.size();
684 const std::size_t size = rSize * cSize;
688 ::Dune::Petsc::MatGetValues( petscMatrix_, r.size(), r.data(), c.size(), c.data(), v.data() );
691 for( std::size_t i=0; i<size; ++i )
695 ::Dune::Petsc::MatSetValues( petscMatrix_, rSize, r.data(), cSize, c.data(), v.data(), INSERT_VALUES );
696 setStatus( statAssembled );
702 ::Dune::Petsc::MatView( petscMatrix_, PETSC_VIEWER_STDOUT_WORLD );
706 void print( std::ostream& s )
const
708 if( &s == &std::cout || &s == &std::cerr )
715 Mat& exportMatrix ()
const
723 PetscLinearOperator ();
728 if( status_ != statNothing )
730 ::Dune::Petsc::MatDestroy( &petscMatrix_ );
731 status_ = statNothing ;
736 void setStatus (
const Status &newstatus)
const
744 template<
class DFS,
class Entity >
745 static void setupIndices (
const PetscMappers< DFS > &mappers,
const Entity &entity, std::vector< PetscInt > &indices )
747 NonBlockMapper< const typename PetscMappers< DFS >::ParallelMapperType, DFS::localBlockSize > nonBlockMapper( mappers.parallelMapper() );
748 nonBlockMapper.map( entity, indices );
751 template<
class DFS,
class Entity >
752 static void setupIndicesBlocked (
const PetscMappers< DFS > &mappers,
const Entity &entity, std::vector< PetscInt > &indices )
754 mappers.parallelMapper().map( entity, indices );
760 DomainMappersType domainMappers_;
761 RangeMappersType rangeMappers_;
764 mutable Mat petscMatrix_;
766 mutable LocalMatrixStackType localMatrixStack_;
767 mutable Status status_;
769 const bool viennaCL_;
770 const bool blockedMode_;
772 mutable std::unique_ptr< PetscDomainFunctionType > petscArg_;
773 mutable std::unique_ptr< PetscRangeFunctionType > petscDest_;
775 mutable std::vector< PetscScalar > v_;
776 mutable std::vector< PetscInt > r_;
777 mutable std::vector< PetscInt > c_;
779 mutable std::vector< PetscInt > unitRows_;
787 template<
typename DomainFunction,
typename RangeFunction >
788 class PetscLinearOperator< DomainFunction, RangeFunction >::LocalMatrix
789 :
public LocalMatrixDefault< typename PetscLinearOperator< DomainFunction, RangeFunction >::LocalMatrixTraits >
791 typedef LocalMatrix ThisType;
792 typedef LocalMatrixDefault< typename PetscLinearOperator< DomainFunction, RangeFunction >::LocalMatrixTraits > BaseType;
794 typedef PetscLinearOperator< DomainFunction, RangeFunction > PetscLinearOperatorType;
798 typedef PetscInt DofIndexType;
799 typedef std::vector< DofIndexType > IndexVectorType;
800 typedef typename DomainFunction::DiscreteFunctionSpaceType DomainSpaceType;
801 typedef typename RangeFunction::DiscreteFunctionSpaceType RangeSpaceType;
802 typedef typename DomainSpaceType::BasisFunctionSetType DomainBasisFunctionSetType;
803 typedef typename RangeSpaceType::BasisFunctionSetType RangeBasisFunctionSetType;
805 enum { rangeBlockSize = RangeSpaceType::localBlockSize };
806 enum { domainBlockSize = DomainSpaceType ::localBlockSize };
811 template<
typename PetscMapping >
812 struct PetscAssignFunctor
814 explicit PetscAssignFunctor (
const PetscMapping &petscMapping, IndexVectorType &indices )
815 : petscMapping_( petscMapping ),
819 template<
typename T >
820 void operator() (
const std::size_t localIndex, T globalIndex ) { indices_[ localIndex ] = petscMapping_.globalMapping( globalIndex ); }
823 const PetscMapping &petscMapping_;
824 IndexVectorType &indices_;
828 LocalMatrix (
const PetscLinearOperatorType &petscLinOp,
829 const DomainSpaceType &domainSpace,
830 const RangeSpaceType &rangeSpace )
831 : BaseType( domainSpace, rangeSpace ),
832 petscLinearOperator_( petscLinOp )
839 void init (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
842 BaseType :: init( domainEntity, rangeEntity );
852 setupIndices( petscLinearOperator_.rangeMappers_, rangeEntity, rowIndices_ );
853 setupIndices( petscLinearOperator_.domainMappers_, domainEntity, colIndices_ );
855 status_ = statAssembled;
856 petscLinearOperator_.setStatus(status_);
859 inline void add (
const int localRow,
const int localCol,
const RangeFieldType &value )
861 assert( status_==statAssembled || status_==statAdd );
863 petscLinearOperator_.setStatus(status_);
864 ::Dune::Petsc::MatSetValue( petscMatrix(), globalRowIndex( localRow ), globalColIndex( localCol ) , value, ADD_VALUES );
867 inline void set(
const int localRow,
const int localCol,
const RangeFieldType &value )
869 assert( status_==statAssembled || status_==statInsert );
870 status_ = statInsert;
871 petscLinearOperator_.setStatus(status_);
872 ::Dune::Petsc::MatSetValue( petscMatrix(), globalRowIndex( localRow ), globalColIndex( localCol ) , value, INSERT_VALUES );
876 void clearRow (
const int localRow )
878 assert( status_==statAssembled || status_==statInsert );
879 status_ = statInsert;
880 petscLinearOperator_.setStatus(status_);
881 const int col = this->columns();
882 const int globalRowIdx = globalRowIndex( localRow );
883 for(
int localCol=0; localCol<col; ++localCol)
885 ::Dune::Petsc::MatSetValue( petscMatrix(), globalRowIdx, globalColIndex( localCol ), 0.0, INSERT_VALUES );
889 inline const RangeFieldType
get (
const int localRow,
const int localCol )
const
891 assert( status_==statAssembled || status_==statGet );
893 petscLinearOperator_.setStatus(status_);
895 const int r[] = {globalRowIndex( localRow )};
896 const int c[] = {globalColIndex( localCol )};
897 ::Dune::Petsc::MatGetValues( petscMatrix(), 1, r, 1, c, v );
901 inline void scale (
const RangeFieldType &factor )
903 const_cast< PetscLinearOperatorType&
> (petscLinearOperator_).scaleLocalMatrix( this->domainEntity(), this->rangeEntity(), factor );
909 Mat& petscMatrix () {
return petscLinearOperator_.petscMatrix_; }
910 const Mat& petscMatrix ()
const {
return petscLinearOperator_.petscMatrix_; }
913 int rows()
const {
return rowIndices_.size(); }
914 int columns()
const {
return colIndices_.size(); }
917 DofIndexType globalRowIndex(
const int localRow )
const
919 assert( localRow <
static_cast< int >( rowIndices_.size() ) );
920 return rowIndices_[ localRow ];
923 DofIndexType globalColIndex(
const int localCol )
const
925 assert( localCol <
static_cast< int >( colIndices_.size() ) );
926 return colIndices_[ localCol ];
932 const PetscLinearOperatorType &petscLinearOperator_;
933 IndexVectorType rowIndices_;
934 IndexVectorType colIndices_;
935 mutable Status status_;
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
Definition: bindguard.hh:11
std::tuple_element< i, Tuple >::type & get(Dune::TypeIndexedTuple< Tuple, Types > &tuple)
Definition: typeindexedtuple.hh:122
BasicParameterReader< std::function< const std::string *(const std::string &, const std::string *) > > ParameterReader
Definition: reader.hh:315
static constexpr T max(T a)
Definition: utility.hh:77
static ParameterContainer & container()
Definition: io/parameter.hh:193