dune-fem  2.8-git
petscdofblock.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_PETSCDOFBLOCK_HH
2 #define DUNE_FEM_PETSCDOFBLOCK_HH
3 
4 #include <algorithm>
5 #include <memory>
6 
8 
9 #if HAVE_PETSC
10 
13 
15 
16 namespace Dune
17 {
18 
19  namespace Fem
20  {
21  template < class PVector >
22  class PetscDofBlock;
23 
24  /* ========================================
25  * class PetscDofBlock::DofProxy
26  * =======================================
27  */
28  template< class PVector >
29  class PetscDofProxy
30  {
31  public:
32  typedef PVector PetscVectorType;
33  typedef typename PetscDofBlock< PetscVectorType >::DofProxy ThisType;
35 
36  static const unsigned int blockSize = PetscVectorType::blockSize;
37 
38  // This is needed to put DofProxies in STL (or STL-like) containers...
39  PetscDofProxy () = default;
40 
41  // ???
42  PetscDofProxy ( PetscScalar s ) {}
43 
44  // this is called by a friend
45  PetscDofProxy ( PetscVectorType &petscVector, IndexType blockIndex, PetscInt indexInBlock )
46  : petscVector_( &petscVector ),
47  blockIndex_( blockIndex ),
48  indexInBlock_( indexInBlock )
49  {}
50 
51  // Does not change the value of the DoF which this proxy, but it assigns this proxy instance...
52  void assign ( const ThisType &other )
53  {
54  petscVector_ = other.petscVector_;
55  blockIndex_ = other.blockIndex_;
56  indexInBlock_ = other.indexInBlock_;
57  }
58 
59  const ThisType& operator= ( ThisType &other )
60  {
61  setValue( other.getValue() );
62  return *this;
63  }
64 
65  const ThisType& operator= ( PetscScalar val )
66  {
67  setValue( val );
68  return *this;
69  }
70 
71  const ThisType& operator*= ( const ThisType& other )
72  {
73  PetscScalar value = getValue() * other.getValue();
74  setValue( value );
75  return *this;
76  }
77 
78  const ThisType& operator*= ( const PetscScalar scalar )
79  {
80  PetscScalar value = getValue() * scalar ;
81  setValue( value );
82  return *this;
83  }
84 
85  const ThisType& operator+= ( const ThisType &other )
86  {
87  setValue( other.getValue(), ADD_VALUES );
88  return *this;
89  }
90 
91  const ThisType& operator+= ( const PetscScalar scalar )
92  {
93  setValue( scalar, ADD_VALUES );
94  return *this;
95  }
96 
97  const ThisType& operator-= ( const ThisType &other )
98  {
99  setValue( -other.getValue(), ADD_VALUES );
100  return *this;
101  }
102 
103  const ThisType& operator-= ( const PetscScalar scalar )
104  {
105  setValue( -scalar, ADD_VALUES );
106  return *this;
107  }
108 
109  // conversion operators
110  operator PetscScalar () const
111  {
112  return valid() ? getValue() : PetscScalar( 0 );
113  }
114 
115  bool valid() const { return bool( petscVector_ ); }
116 
117  protected:
118  PetscScalar getValue () const
119  {
120  assert( valid() );
121  PetscInt index = blockSize * petscVector().mappers().ghostIndex( blockIndex_ ) + indexInBlock_;
122  PetscScalar ret;
123  ::Dune::Petsc::VecGetValues( *petscVector().getGhostedVector(), 1, &index, &ret );
124  return ret;
125  }
126 
127  void setValue ( const PetscScalar &val, InsertMode mode = INSERT_VALUES )
128  {
129  PetscInt index = blockSize * petscVector().mappers().ghostIndex( blockIndex_ ) + indexInBlock_;
130  ::Dune::Petsc::VecSetValue( *( petscVector().getGhostedVector() ), index, val, mode );
131  petscVector().hasBeenModified();
132  }
133 
134  PetscVectorType& petscVector ()
135  {
136  assert( petscVector_ );
137  return *petscVector_;
138  }
139 
140  const PetscVectorType& petscVector () const
141  {
142  assert( petscVector_ );
143  return *petscVector_;
144  }
145 
146  // data fields
147  PetscVectorType *petscVector_ = nullptr;
148  IndexType blockIndex_ = 0;
149  PetscInt indexInBlock_ = 0;
150  };
151 
152  template< class Scalar, class PVector >
153  void axpy ( const Scalar &a, const Scalar &x, PetscDofProxy< PVector > proxy )
154  {
155  proxy += a*x;
156  }
157 
158 
159 
160 
161  /* ========================================
162  * class PetscDofBlock
163  * =======================================
164  */
165  template< class PVector >
166  class PetscDofBlock
167  {
168  typedef PetscDofBlock< PVector > ThisType;
169 
170  public:
171  typedef PVector PetscVectorType;
172  typedef PetscInt IndexType;
173 
174  static const unsigned int blockSize = PetscVectorType::blockSize;
175 
176  typedef PetscDofProxy< PVector > DofProxy;
177  class DofIterator;
178 
179  // Needed so that we can put this class in a Dune::Envelope
180  typedef std::pair< PetscVectorType&, IndexType > UnaryConstructorParamType;
181 
182  // this is the ctor to be used
183  PetscDofBlock ( PetscVectorType &petscVector, IndexType blockIndex )
184  : petscVector_( petscVector ),
185  blockIndex_( blockIndex )
186  {}
187 
188  // this is the ctor to be used
189  PetscDofBlock ( const PetscDofBlock& other )
190  : petscVector_( other.petscVector_ ),
191  blockIndex_( other.blockIndex_ )
192  {}
193 
194  // ..or this one, which is semantically equivalent to the above ctor
195  explicit PetscDofBlock ( UnaryConstructorParamType arg )
196  : petscVector_( arg.first ),
197  blockIndex_( arg.second )
198  {}
199 
200  const PetscDofBlock& operator*= ( const PetscScalar value )
201  {
202  for( unsigned int i=0; i<blockSize; ++i )
203  (*this)[ i ] *= value ;
204  return *this;
205  }
206 
207  PetscDofBlock () = delete;
208 
209  IndexType size() const { return blockSize; }
210 
211  DofProxy operator [] ( unsigned int index )
212  {
213  assert( index < blockSize );
214  return DofProxy( petscVector_, blockIndex_, index );
215  }
216 
217  const DofProxy operator [] ( unsigned int index ) const
218  {
219  assert( index < blockSize );
220  return DofProxy( petscVector_, blockIndex_, index );
221  }
222 
223  private:
224  PetscVectorType &petscVector_;
225  IndexType blockIndex_;
226  };
227 
229  //template< class Traits, class PVector >
230  template< class Traits, class PVector >
231  inline OutStreamInterface< Traits >&
232  operator<< ( OutStreamInterface< Traits > &out,
233  const PetscDofProxy< PVector >& value )
234  {
235  // write to stream
236  out << PetscScalar( value );
237  return out;
238  }
239 
241  //template< class Traits, class PVector >
242  template< class Traits, class PVector >
243  inline InStreamInterface< Traits >&
244  operator>> ( InStreamInterface< Traits > &in,
245  PetscDofProxy< PVector > value )
246  {
247  PetscScalar val;
248  // get value from stream
249  in >> val;
250  // write value to discrete function
251  value = val;
252  return in;
253  }
254 
255 
256  /*
257  * This is almost a bidirectional iterator but does not completely satisfy the required
258  * interface (see the C++2003 standard, 24.1.4) [no default ctor, no operator->].
259  */
260 
261  /* ========================================
262  * class PetscDofBlock::DofIterator
263  * =======================================
264  */
265  template< class PVector >
266  class PetscDofBlock< PVector >::DofIterator
267  : public std::iterator< std::input_iterator_tag, PetscScalar >
268  {
269  typedef typename PetscDofBlock< PVector >::DofIterator ThisType;
270  typedef PetscDofBlock< PVector > DofBlockType;
271 
272  typedef std::iterator< std::input_iterator_tag, PetscScalar > BaseType;
273 
274  // TODO: get rid of this! we don't like shared pointers. Own a real instance instead!
275  typedef std::shared_ptr< DofBlockType > DofBlockSharedPointer;
276 
277  public:
278  typedef PVector PetscVectorType;
279  typedef typename DofBlockType::DofProxy value_type; // (this overrides the 2nd template parameter of BaseType...)
280 
281  // standard ctor
282  DofIterator ( PetscVectorType &petscVector, unsigned int blockIndex, PetscInt indexInBlock = 0 )
283  : petscVector_( petscVector ),
284  blockIndex_( blockIndex ),
285  indexInBlock_( indexInBlock )
286  {
287  // blockIndex == size denotes the end iterator
288  assert( static_cast< std::size_t >( blockIndex ) <= petscVector_.mappers().size() );
289 
290  // Is this not the end iterator?
291  if( static_cast< std::size_t >( blockIndex ) < petscVector_.mappers().size() )
292  {
293  resetBlockPtr();
294  }
295  }
296 
297  bool operator== ( const ThisType &other ) const
298  {
299  return (blockIndex_ == other.blockIndex_) && (indexInBlock_ == other.indexInBlock_);
300  }
301 
302  bool operator!= ( const ThisType &other ) const { return !this->operator==( other ); };
303 
304  value_type operator* () { return block()[ indexInBlock_]; }
305  const value_type operator* () const { return block()[ indexInBlock_ ]; }
306 
307  // prefix increment
308  ThisType& operator++ () { increment(); return *this; }
309 
310  // prefix decrement
311  ThisType& operator-- () { decrement(); return *this; }
312 
313  DofIterator () = delete;
314 
315  private:
316  PetscVectorType& petscVector () { return petscVector_; }
317 
318  void increment ()
319  {
320  ++indexInBlock_;
321  if( static_cast< std::size_t >( indexInBlock_ ) >= DofBlockType::blockSize )
322  {
323  ++blockIndex_;
324  indexInBlock_ = 0;
325  if( static_cast< std::size_t >( blockIndex_ ) < petscVector().mappers().size() )
326  resetBlockPtr();
327  }
328  }
329 
330  void decrement ()
331  {
332  assert( blockIndex_ > 0 );
333  if( indexInBlock_ == 0 )
334  {
335  indexInBlock_ = DofBlockType::blockSize - 1;
336  --blockIndex_;
337  resetBlockPtr();
338  }
339  else
340  {
341  --blockIndex_;
342  }
343  }
344 
345  // TODO: iterator should own the block - _no_ new...
346  void resetBlockPtr () { blockPtr_.reset( new DofBlockType( *petscVector().block( blockIndex_ ) ) ); }
347 
348  DofBlockType& block () const { return *blockPtr_.get(); }
349 
350  PetscVectorType &petscVector_;
351  unsigned int blockIndex_;
352  PetscInt indexInBlock_;
353  DofBlockSharedPointer blockPtr_;
354  };
355 
356  } // namespace Fem
357 
358 } // namespace Dune
359 
360 #endif // #if HAVE_PETSC
361 
362 #endif // DUNE_FEM_PETSCDOFBLOCK_HH
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:397
Double operator*(const Double &a, const Double &b)
Definition: double.hh:506
bool operator==(const Double &a, const Double &b)
Definition: double.hh:600
bool operator!=(const Double &a, const Double &b)
Definition: double.hh:640
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:417
typename Impl::Index< Range >::Type IndexType
Definition: hybrid.hh:69