DUNE-FEM (unstable)

densematrix.hh
1#ifndef DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
2#define DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
3
4#include <memory>
5#include <iostream>
6
7#include <dune/common/dynmatrixev.hh>
8
9#include <dune/fem/function/adaptivefunction.hh>
10#include <dune/fem/misc/functor.hh>
11#include <dune/fem/operator/common/localmatrix.hh>
12#include <dune/fem/operator/common/localmatrixwrapper.hh>
13#include <dune/fem/solver/krylovinverseoperators.hh>
14#include <dune/fem/storage/objectstack.hh>
15
16namespace Dune
17{
18
19 namespace Fem
20 {
21
22 // DenseRowMatrix
23 // --------------
24
25 template< class F >
26 class DenseRowMatrix
27 {
28 typedef DenseRowMatrix< F > ThisType;
29
30 public:
31 typedef F Field;
32
33 template< class RF >
34 class Row;
35
36 DenseRowMatrix () : rows_( 0 ), cols_( 0 ) {}
37
38 DenseRowMatrix ( unsigned int rows, unsigned int cols )
39 : rows_( 0 ), cols_( 0 )
40 {
41 reserve( rows, cols );
42 }
43
44 unsigned int rows () const { return rows_; }
45 unsigned int cols () const { return cols_; }
46
47 const Field &operator() ( unsigned int row, unsigned int col ) const
48 {
49 assert( (row < rows()) && (col < cols()) );
50 return fields_[ row*cols() + col ];
51 }
52
53 Field &operator() ( unsigned int row, unsigned int col )
54 {
55 assert( (row < rows()) && (col < cols()) );
56 return fields_[ row*cols() + col ];
57 }
58
59 void add ( unsigned int row, unsigned int col, const Field &value )
60 {
61 assert( (row < rows()) && (col < cols()) );
62 fields_[ row*cols() + col ] += value;
63 }
64
65 Row< const Field > operator[] ( unsigned int row ) const
66 {
67 assert( row < rows() );
68 return Row< const Field >( cols(), fields_.get() + row*cols() );
69 }
70
71 Row< Field > operator[] ( unsigned int row )
72 {
73 assert( row < rows() );
74 return Row< Field >( cols(), fields_.get() + row*cols() );
75 }
76
77 void clear () { std::fill( fields_.get(), fields_.get() + (rows()*cols()), Field( 0 ) ); }
78
79 void multiply ( const Field *x, Field *y ) const
80 {
81 for( unsigned int row = 0; row < rows(); ++row )
82 {
83 const Field *fields = fields_.get() + row*cols();
84 y[ row ] = Field( 0 );
85 for( unsigned int col = 0; col < cols(); ++col )
86 y[ row ] += fields[ col ] * x[ col ];
87 }
88 }
89
99 std::vector< std::complex< double > > eigenValues ()
100 {
101 if( rows() != cols() )
102 DUNE_THROW( InvalidStateException, "Requiring square matrix for eigenvalue computation" );
103
104 const long int N = rows();
105 const char jobvl = 'n';
106 const char jobvr = 'n';
107
108 // working memory
109 std::unique_ptr< double[] > work = std::make_unique< double[] >( 5*N );
110
111 // return value information
112 long int info = 0;
113 long int lwork = 3*N;
114
115 // call LAPACK routine (see fmatrixev_ext.cc)
116 DynamicMatrixHelp::eigenValuesNonsymLapackCall( &jobvl, &jobvr, &N, fields_, &N, work.get(), work.get()+N, 0, &N, 0, &N, work.get()+2*N, &lwork, &info );
117
118 if( info != 0 )
119 DUNE_THROW( MathError, "DenseRowMatrix: Eigenvalue computation failed" );
120
121 std::vector< std::complex< double > > eigenValues( N );
122 std::transform( work.get(), work.get()+N, work.get()+N, eigenValues.begin(), [] ( double r, double i ) { return std::complex< double >( r, i ); } );
123 return eigenValues;
124 }
125
126 void reserve ( unsigned int rows, unsigned int cols )
127 {
128 if( (rows != rows_) || (cols != cols_) )
129 {
130 fields_.reset( new Field[ rows*cols ] );
131 rows_ = rows;
132 cols_ = cols;
133 }
134 // Martin: Is clearing required, here?
135 clear();
136 }
137
138 void print( std::ostream& s=std::cout ) const
139 {
140 s.precision( 6 );
141 for( unsigned int row = 0; row < rows(); ++row )
142 {
143 const Field *fields = fields_ + row*cols();
144 for( unsigned int col = 0; col < cols(); ++col )
145 s << fields[ col ] << " ";
146
147 s << std::endl;
148 }
149 }
150
151 private:
152 unsigned int rows_, cols_;
153 std::unique_ptr< Field[] > fields_;
154 };
155
156
157
158 // DenseRowMatrix::Row
159 // -------------------
160
161 template< class F >
162 template< class RF >
163 class DenseRowMatrix< F >::Row
164 {
165 typedef Row< RF > ThisType;
166
167 template< class > friend class Row;
168
169 public:
170 Row ( unsigned int cols, RF *fields )
171 : cols_( cols ),
172 fields_( fields )
173 {}
174
175 Row ( const Row< F > &row )
176 : cols_( row.cols_ ),
177 fields_( row.fields_ )
178 {}
179
180 const RF &operator[] ( const unsigned int col ) const
181 {
182 assert( col < size() );
183 return fields_[ col ];
184 }
185
186 RF &operator[] ( const unsigned int col )
187 {
188 assert( col < size() );
189 return fields_[ col ];
190 }
191
192 void clear ()
193 {
194 Field *const end = fields_ + size();
195 for( Field *it = fields_; it != end; ++it )
196 *it = Field( 0 );
197 }
198
199 unsigned int size () const
200 {
201 return cols_;
202 }
203
204 private:
205 unsigned int cols_;
206 RF *fields_;
207 };
208
209
210
211 // DenseRowMatrixObject
212 // --------------------
213
214 template< class DomainSpace, class RangeSpace >
215 class DenseRowMatrixObject
216 {
217 typedef DenseRowMatrixObject< DomainSpace, RangeSpace > ThisType;
218
219 public:
220 typedef DomainSpace DomainSpaceType;
221 typedef RangeSpace RangeSpaceType;
222
223 typedef typename RangeSpaceType::RangeFieldType Field;
224
225 typedef typename DomainSpaceType::BlockMapperType DomainBlockMapperType;
226 typedef NonBlockMapper< DomainBlockMapperType, DomainSpaceType::localBlockSize > DomainMapperType;
227 typedef typename RangeSpaceType::BlockMapperType RangeBlockMapperType;
228 typedef NonBlockMapper< RangeBlockMapperType, RangeSpaceType::localBlockSize > RangeMapperType;
229
230 typedef typename DomainSpaceType::EntityType DomainEntityType;
231 typedef typename RangeSpaceType::EntityType RangeEntityType;
232 typedef typename DomainSpace::GridType::template Codim< 0 >::Entity ColEntityType;
233 typedef typename RangeSpace::GridType::template Codim< 0 >::Entity RowEntityType;
234
235 typedef DenseRowMatrix< Field > MatrixType;
236
237 private:
238 class LocalMatrixTraits;
239 class LocalMatrix;
240 class LocalMatrixFactory;
241
242 typedef Fem::ObjectStack< LocalMatrixFactory > LocalMatrixStack;
243
244 public:
245 typedef LocalMatrixWrapper< LocalMatrixStack > LocalMatrixType;
246
247 DenseRowMatrixObject ( const DomainSpaceType &domainSpace,
248 const RangeSpaceType &rangeSpace )
249 : domainSpace_( domainSpace ),
250 rangeSpace_( rangeSpace ),
251 domainMapper_( domainSpace.blockMapper() ),
252 rangeMapper_( rangeSpace.blockMapper() ),
253 domainSequence_( -1 ),
254 rangeSequence_( -1 ),
255 localMatrixFactory_( *this ),
256 localMatrixStack_( localMatrixFactory_ )
257 {}
258
259 MatrixType &matrix ()
260 {
261 return matrix_;
262 }
263
264 LocalMatrixType localMatrix ( const RowEntityType &rowEntity,
265 const ColEntityType &colEntity )
266 {
267 return LocalMatrixType( localMatrixStack_, rowEntity, colEntity );
268 }
269
270
271 LocalMatrixType localMatrix () const
272 {
273 return LocalMatrixType( localMatrixStack_ );
274 }
275
276 void clear ()
277 {
278 matrix_.clear();
279 }
280
281 template< class Stencil >
282 void reserve ( const Stencil &stencil, bool verbose = false )
283 {
284 if( (domainSequence_ != domainSpace().sequence()) || (rangeSequence_ != rangeSpace().sequence()) )
285 {
286 matrix_.reserve( rangeSpace().size(), domainSpace().size() );
287 domainSequence_ = domainSpace().sequence();
288 rangeSequence_ = rangeSpace().sequence();
289 }
290 }
291
292 template< class DomainFunction, class RangeFunction >
293 void apply ( const DomainFunction &u, RangeFunction &w ) const
294 {
295 matrix_.multiply( u.leakPointer(), w.leakPointer() );
296 rangeSpace().communicate( w );
297 }
298
299 Field ddotOEM ( const Field *v, const Field *w ) const
300 {
301 typedef AdaptiveDiscreteFunction< RangeSpaceType > RangeFunction;
302 RangeFunction vFunction( "v (ddotOEM)", rangeSpace(), v );
303 RangeFunction wFunction( "w (ddotOEM)", rangeSpace(), w );
304 return vFunction.scalarProductDofs( wFunction );
305 }
306
307 void multOEM ( const Field *u, Field *w ) const
308 {
309 matrix_.multiply( u, w );
310
311 typedef AdaptiveDiscreteFunction< RangeSpaceType > RangeFunction;
312 RangeFunction wFunction( "w (multOEM)", rangeSpace(), w );
313 rangeSpace().communicate( wFunction );
314 }
315
316 template< class DiscreteFunctionType >
317 void extractDiagonal( DiscreteFunctionType &diag ) const
318 {
319 assert( matrix_.rows() == matrix_.cols() );
320 const auto dofEnd = diag.dend();
321 unsigned int row = 0;
322 for( auto dofIt = diag.dbegin(); dofIt != dofEnd; ++dofIt, ++row )
323 {
324 assert( row < matrix_.rows() );
325 (*dofIt) = matrix_( row, row );
326 }
327 }
328
329 const DomainSpace &domainSpace () const { return domainSpace_; }
330 const RangeSpace &rangeSpace () const { return rangeSpace_; }
331
332 private:
333 const DomainSpaceType &domainSpace_;
334 const RangeSpaceType &rangeSpace_;
335 DomainMapperType domainMapper_;
336 RangeMapperType rangeMapper_;
337
338 int domainSequence_;
339 int rangeSequence_;
340
341 MatrixType matrix_;
342
343 LocalMatrixFactory localMatrixFactory_;
344 mutable LocalMatrixStack localMatrixStack_;
345 };
346
347
348
349 // DenseRowMatrixObject::LocalMatrixTraits
350 // ---------------------------------------
351
352 template< class DomainSpace, class RangeSpace >
353 class DenseRowMatrixObject< DomainSpace, RangeSpace >::LocalMatrixTraits
354 {
355 typedef DenseRowMatrixObject< DomainSpace, RangeSpace > MatrixObject;
356
357 public:
358 typedef typename MatrixObject::DomainSpaceType DomainSpaceType;
359 typedef typename MatrixObject::RangeSpaceType RangeSpaceType;
360
361 typedef typename MatrixObject::Field RangeFieldType;
362 typedef RangeFieldType LittleBlockType;
363
364 typedef typename MatrixObject::LocalMatrix LocalMatrixType;
365 };
366
367
368
369 // DenseRowMatrixObject::LocalMatrix
370 // ---------------------------------
371
372 template< class DomainSpace, class RangeSpace >
373 class DenseRowMatrixObject< DomainSpace, RangeSpace >::LocalMatrix
374 : public LocalMatrixDefault< LocalMatrixTraits >
375 {
376 typedef DenseRowMatrixObject< DomainSpace, RangeSpace > MatrixObject;
377
378 typedef LocalMatrix ThisType;
379 typedef LocalMatrixDefault< LocalMatrixTraits > BaseType;
380
381 public:
382 typedef LocalMatrixTraits Traits;
383
384 typedef typename MatrixObject::MatrixType MatrixType;
385
386 typedef typename Traits::RangeFieldType RangeFieldType;
387 typedef typename Traits::LittleBlockType LittleBlockType;
388 typedef RangeFieldType DofType;
389
390 LocalMatrix ( MatrixType &matrix, const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const DomainMapperType &domainMapper, const RangeMapperType &rangeMapper )
391 : BaseType( domainSpace, rangeSpace ),
392 matrix_( matrix ),
393 domainMapper_( domainMapper ),
394 rangeMapper_( rangeMapper )
395 {}
396
397 LocalMatrix ( const ThisType & ) = delete;
398 ThisType &operator= ( const ThisType & ) = delete;
399
400 void init ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
401 {
402 BaseType::init( domainEntity, rangeEntity );
403
404 rowIndices_.resize( rangeMapper_.numDofs( rangeEntity ) );
405 rangeMapper_.mapEach( rangeEntity, Fem::AssignFunctor< std::vector< unsigned int > >( rowIndices_ ) );
406 colIndices_.resize( domainMapper_.numDofs( domainEntity ) );
407 domainMapper_.mapEach( domainEntity, Fem::AssignFunctor< std::vector< unsigned int > >( colIndices_ ) );
408 }
409
410 int rows () const { return rowIndices_.size(); }
411 int cols () const { return colIndices_.size(); }
412
413 void add ( const int row, const int col, const DofType &value )
414 {
415 assert( (row >= 0) && (row < rows()) );
416 assert( (col >= 0) && (col < cols()) );
417 matrix_( rowIndices_[ row ], colIndices_[ col ] ) += value;
418 }
419
420 const DofType &get ( const int row, const int col ) const
421 {
422 assert( (row >= 0) && (row < rows()) );
423 assert( (col >= 0) && (col < cols()) );
424 return matrix_( rowIndices_[ row ], colIndices_[ col ] );
425 }
426
427 void set ( const int row, const int col, const DofType &value )
428 {
429 assert( (row >= 0) && (row < rows()) );
430 assert( (col >= 0) && (col < cols()) );
431 matrix_( rowIndices_[ row ], colIndices_[ col ] ) = value;
432 }
433
434 void clearRow ( const int row )
435 {
436 assert( (row >= 0) && (row < rows()) );
437 const unsigned int rowIndex = rowIndices_[ row ];
438 matrix_[ rowIndex ].clear();
439 }
440
441 void unitRow ( const int row )
442 {
443 clearRow( row );
444 set( row, row, DofType( 1 ) );
445 }
446
447 void clear ()
448 {
449 typedef std::vector< unsigned int >::const_iterator Iterator;
450 const Iterator rowEnd = rowIndices_.end();
451 for( Iterator rowIt = rowIndices_.begin(); rowIt != rowEnd; ++rowIt )
452 {
453 const Iterator colEnd = colIndices_.end();
454 for( Iterator colIt = colIndices_.begin(); colIt != colEnd; ++colIt )
455 matrix_( *rowIt, *colIt ) = DofType( 0 );
456 }
457 }
458
459 void scale ( const DofType &value )
460 {
461 typedef std::vector< unsigned int >::const_iterator Iterator;
462 const Iterator rowEnd = rowIndices_.end();
463 for( Iterator rowIt = rowIndices_.begin(); rowIt != rowEnd; ++rowIt )
464 {
465 const Iterator colEnd = colIndices_.end();
466 for( Iterator colIt = colIndices_.begin(); colIt != colEnd; ++colIt )
467 matrix_( *rowIt, *colIt ) *= value;
468 }
469 }
470
471 protected:
472 using BaseType::domainSpace_;
473 using BaseType::rangeSpace_;
474
475 private:
476 MatrixType &matrix_;
477 const DomainMapperType &domainMapper_;
478 const RangeMapperType &rangeMapper_;
479 std::vector< unsigned int > rowIndices_;
480 std::vector< unsigned int > colIndices_;
481 };
482
483
484
485 // DenseRowMatrixObject::LocalMatrixFactory
486 // ----------------------------------------
487
488 template< class DomainSpace, class RangeSpace >
489 class DenseRowMatrixObject< DomainSpace, RangeSpace >::LocalMatrixFactory
490 {
491 typedef DenseRowMatrixObject< DomainSpace, RangeSpace > MatrixObject;
492
493 public:
494 typedef LocalMatrix ObjectType;
495
496 LocalMatrixFactory ( MatrixObject &matrixObject )
497 : matrixObject_( &matrixObject )
498 {}
499
500 ObjectType *newObject () const
501 {
502 return new ObjectType( matrixObject_->matrix_, matrixObject_->domainSpace_, matrixObject_->rangeSpace_, matrixObject_->domainMapper_, matrixObject_->rangeMapper_ );
503 }
504
505 private:
506 MatrixObject *matrixObject_;
507 };
508
509 } // namespace Fem
510
511} // namespace Dune
512
513#endif // #ifndef DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
Dune::Fem::DiscreteFunctionDefault::dbegin
ConstDofIteratorType dbegin() const
Obtain the constant iterator pointing to the first dof.
Definition: discretefunction.hh:761
Dune::Fem::DiscreteFunctionDefault::dend
ConstDofIteratorType dend() const
Obtain the constant iterator pointing to the last dof.
Definition: discretefunction.hh:773
Dune::Fem::TupleDiscreteFunction
forward declaration
Definition: discretefunction.hh:51
dynmatrixev.hh
utility functions to compute eigenvalues for dense matrices.
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::size
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Dune::get
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 27, 22:29, 2024)