1#ifndef DUNE_FEM_PETSCDOFBLOCK_HH
2#define DUNE_FEM_PETSCDOFBLOCK_HH
21 template <
class PVector >
28 template<
class PVector >
32 typedef PVector PetscVectorType;
33 typedef typename PetscDofBlock< PetscVectorType >::DofProxy ThisType;
34 typedef typename PetscDofBlock< PetscVectorType >::IndexType
IndexType;
36 static const unsigned int blockSize = PetscVectorType::blockSize;
39 PetscDofProxy () =
default;
42 PetscDofProxy ( PetscScalar s ) {}
45 PetscDofProxy ( PetscVectorType &petscVector, IndexType blockIndex, PetscInt indexInBlock )
46 : petscVector_( &petscVector ),
47 blockIndex_( blockIndex ),
48 indexInBlock_( indexInBlock )
52 void assign (
const ThisType &other )
54 petscVector_ = other.petscVector_;
55 blockIndex_ = other.blockIndex_;
56 indexInBlock_ = other.indexInBlock_;
59 const ThisType& operator= ( ThisType &other )
61 setValue( other.getValue() );
65 const ThisType& operator= ( PetscScalar val )
71 const ThisType& operator*= (
const ThisType& other )
73 PetscScalar value = getValue() * other.getValue();
78 const ThisType& operator*= (
const PetscScalar scalar )
80 PetscScalar value = getValue() * scalar ;
85 const ThisType& operator+= (
const ThisType &other )
87 setValue( other.getValue(), ADD_VALUES );
91 const ThisType& operator+= (
const PetscScalar scalar )
93 setValue( scalar, ADD_VALUES );
97 const ThisType& operator-= (
const ThisType &other )
99 setValue( -other.getValue(), ADD_VALUES );
103 const ThisType& operator-= (
const PetscScalar scalar )
105 setValue( -scalar, ADD_VALUES );
110 operator PetscScalar ()
const
112 return valid() ? getValue() : PetscScalar( 0 );
115 bool valid()
const {
return bool( petscVector_ ); }
118 PetscScalar getValue ()
const
121 PetscInt index = blockSize * petscVector().mappers().ghostIndex( blockIndex_ ) + indexInBlock_;
123 ::Dune::Petsc::VecGetValues( *petscVector().getGhostedVector(), 1, &index, &ret );
127 void setValue (
const PetscScalar &val, InsertMode mode = INSERT_VALUES )
129 PetscInt index = blockSize * petscVector().mappers().ghostIndex( blockIndex_ ) + indexInBlock_;
130 ::Dune::Petsc::VecSetValue( *( petscVector().getGhostedVector() ), index, val, mode );
131 petscVector().hasBeenModified();
134 PetscVectorType& petscVector ()
136 assert( petscVector_ );
137 return *petscVector_;
140 const PetscVectorType& petscVector ()
const
142 assert( petscVector_ );
143 return *petscVector_;
147 PetscVectorType *petscVector_ =
nullptr;
149 PetscInt indexInBlock_ = 0;
152 template<
class Scalar,
class PVector >
153 void axpy (
const Scalar &a,
const Scalar &x, PetscDofProxy< PVector > proxy )
165 template<
class PVector >
168 typedef PetscDofBlock< PVector > ThisType;
171 typedef PVector PetscVectorType;
174 static const unsigned int blockSize = PetscVectorType::blockSize;
176 typedef PetscDofProxy< PVector > DofProxy;
180 typedef std::pair< PetscVectorType&, IndexType > UnaryConstructorParamType;
183 PetscDofBlock ( PetscVectorType &petscVector, IndexType blockIndex )
184 : petscVector_( petscVector ),
185 blockIndex_( blockIndex )
189 PetscDofBlock (
const PetscDofBlock& other )
190 : petscVector_( other.petscVector_ ),
191 blockIndex_( other.blockIndex_ )
195 explicit PetscDofBlock ( UnaryConstructorParamType arg )
196 : petscVector_( arg.first ),
197 blockIndex_( arg.second )
200 const PetscDofBlock& operator*= (
const PetscScalar value )
202 for(
unsigned int i=0; i<blockSize; ++i )
203 (*
this)[ i ] *= value ;
207 PetscDofBlock () =
delete;
209 IndexType size()
const {
return blockSize; }
211 DofProxy operator [] (
unsigned int index )
213 assert( index < blockSize );
214 return DofProxy( petscVector_, blockIndex_, index );
217 const DofProxy operator [] (
unsigned int index )
const
219 assert( index < blockSize );
220 return DofProxy( petscVector_, blockIndex_, index );
224 PetscVectorType &petscVector_;
230 template<
class Traits,
class PVector >
231 inline OutStreamInterface< Traits >&
233 const PetscDofProxy< PVector >& value )
236 out << PetscScalar( value );
242 template<
class Traits,
class PVector >
243 inline InStreamInterface< Traits >&
245 PetscDofProxy< PVector > value )
265 template<
class PVector >
266 class PetscDofBlock< PVector >::DofIterator
267 :
public std::iterator< std::input_iterator_tag, PetscScalar >
269 typedef typename PetscDofBlock< PVector >::DofIterator ThisType;
270 typedef PetscDofBlock< PVector > DofBlockType;
272 typedef std::iterator< std::input_iterator_tag, PetscScalar > BaseType;
275 typedef std::shared_ptr< DofBlockType > DofBlockSharedPointer;
278 typedef PVector PetscVectorType;
279 typedef typename DofBlockType::DofProxy value_type;
282 DofIterator ( PetscVectorType &petscVector,
unsigned int blockIndex, PetscInt indexInBlock = 0 )
283 : petscVector_( petscVector ),
284 blockIndex_( blockIndex ),
285 indexInBlock_( indexInBlock )
288 assert(
static_cast< std::size_t
>( blockIndex ) <= petscVector_.mappers().size() );
291 if(
static_cast< std::size_t
>( blockIndex ) < petscVector_.mappers().size() )
297 bool operator== (
const ThisType &other )
const
299 return (blockIndex_ == other.blockIndex_) && (indexInBlock_ == other.indexInBlock_);
302 bool operator!= (
const ThisType &other )
const {
return !this->
operator==( other ); };
304 value_type operator* () {
return block()[ indexInBlock_]; }
305 const value_type operator* ()
const {
return block()[ indexInBlock_ ]; }
308 ThisType& operator++ () { increment();
return *
this; }
311 ThisType& operator-- () { decrement();
return *
this; }
313 DofIterator () =
delete;
316 PetscVectorType& petscVector () {
return petscVector_; }
321 if(
static_cast< std::size_t
>( indexInBlock_ ) >= DofBlockType::blockSize )
325 if(
static_cast< std::size_t
>( blockIndex_ ) < petscVector().mappers().size() )
332 assert( blockIndex_ > 0 );
333 if( indexInBlock_ == 0 )
335 indexInBlock_ = DofBlockType::blockSize - 1;
346 void resetBlockPtr () { blockPtr_.reset(
new DofBlockType( *petscVector().block( blockIndex_ ) ) ); }
348 DofBlockType& block ()
const {
return *blockPtr_.get(); }
350 PetscVectorType &petscVector_;
351 unsigned int blockIndex_;
352 PetscInt indexInBlock_;
353 DofBlockSharedPointer blockPtr_;
bool operator==(const DiscreteFunctionInterface< ImplX > &x, const DiscreteFunctionInterface< ImplY > &y)
Definition: common/discretefunction.hh:1053
Definition: bindguard.hh:11
OutStreamInterface< StreamTraits > & operator<<(OutStreamInterface< StreamTraits > &out, const DiscreteFunctionInterface< Impl > &df)
write a discrete function into an output stream
Definition: discretefunction_inline.hh:396
void axpy(const T &a, const T &x, T &y)
Definition: space/basisfunctionset/functor.hh:38
InStreamInterface< StreamTraits > & operator>>(InStreamInterface< StreamTraits > &in, DiscreteFunctionInterface< Impl > &df)
read a discrete function from an input stream
Definition: discretefunction_inline.hh:416
typename Impl::Index< Range >::Type IndexType
Definition: hybrid.hh:69