dune-fem 2.8.0
Loading...
Searching...
No Matches
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
16namespace 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;
34 typedef typename PetscDofBlock< PetscVectorType >::IndexType IndexType;
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
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