dune-fem 2.8.0
Loading...
Searching...
No Matches
localmassmatrix.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_LOCALMASSMATRIX_HH
2#define DUNE_FEM_LOCALMASSMATRIX_HH
3
4//- dune-common includes
5#include <dune/common/dynmatrix.hh>
6#include <dune/common/fmatrix.hh>
7
8//- dune-geometry includes
9#include <dune/geometry/typeindex.hh>
10
11//- dune-fem includes
18#include <dune/fem/version.hh>
19
20namespace Dune
21{
22
23 namespace Fem
24 {
25
32 template< class DiscreteFunctionSpace, class VolumeQuadrature >
34 {
36
37 public:
39 typedef typename DiscreteFunctionSpaceType :: RangeFieldType ctype;
40 typedef typename DiscreteFunctionSpaceType :: RangeFieldType RangeFieldType;
41 typedef typename DiscreteFunctionSpaceType :: RangeType RangeType;
42
43 enum { dimRange = DiscreteFunctionSpaceType :: dimRange };
44 enum { localBlockSize = DiscreteFunctionSpaceType :: localBlockSize };
46
47 typedef Dune::FieldMatrix< ctype, dgNumDofs, dgNumDofs > DGMatrixType;
48 typedef Dune::FieldVector< ctype, dgNumDofs > DGVectorType;
49
50 typedef typename DiscreteFunctionSpaceType :: GridPartType GridPartType;
51
52 typedef typename DiscreteFunctionSpaceType :: IndexSetType IndexSetType;
53 typedef typename IndexSetType :: IndexType IndexType;
54
55 typedef typename DiscreteFunctionSpaceType :: BasisFunctionSetType BasisFunctionSetType;
56
57 typedef typename GridPartType :: GridType GridType;
58 typedef typename DiscreteFunctionSpaceType :: EntityType EntityType;
59 typedef typename EntityType :: Geometry Geometry;
60
61 typedef VolumeQuadrature VolumeQuadratureType;
62
64
66 enum { StructuredGrid = Dune::Capabilities::isCartesian< GridType >::v };
67
69 typedef typename GeometryInformationType :: DomainType DomainType;
70
71 // use dynamic matrix from dune-common
72 typedef Dune::DynamicMatrix< RangeFieldType > MatrixType;
73 typedef Dune::DynamicVector< RangeFieldType > VectorType;
74
75 protected:
76 std::shared_ptr< const DiscreteFunctionSpaceType > spc_;
78
80 const std::function<int(const int)> volumeQuadratureOrder_;
81 const bool affine_;
82
85
86 // use dynamic vector from dune-common
87
90
91 mutable std::vector< RangeType > phi_;
92 mutable std::vector< RangeType > phiMass_;
93
94 typedef std::pair< std::unique_ptr< MatrixType >, std::unique_ptr< VectorType > > MatrixPairType;
95 typedef std::map< const int, MatrixPairType > MassMatrixStorageType;
96 typedef std::vector< MassMatrixStorageType > LocalInverseMassMatrixStorageType;
97
99
100 // index of entity from index set, don't setup mass matrix for the same entity twice
102 mutable unsigned int lastTopologyId_ ;
103 // sequence number (obtained from DofManager via the space)
104 mutable int sequence_;
105
107 {
108 typedef Dune::FieldMatrix< ctype, dimRange, dimRange > MassFactorType;
109
110 // return false since we don;t have a mass term
111 bool hasMass () const { return false; }
112
113 void mass ( const EntityType &, const VolumeQuadratureType &, int, MassFactorType & ) const
114 {}
115 };
116
117
118 bool checkDiagonalMatrix( const MatrixType& matrix ) const
119 {
120 const int rows = matrix.rows();
121 for( int r=0; r<rows; ++r )
122 {
123 for( int c=0; c<r; ++c ) // the mass matrix is symmetric
124 {
125 // if we find one off diagonal non-zero return false
126 if( std::abs(matrix[r][c]) > 1e-12 )
127 return false;
128 }
129 }
130 return true;
131 }
132
133 template< class BasisFunctionSet >
135 getLocalInverseMassMatrix ( const EntityType &entity, const Geometry &geo,
136 const BasisFunctionSet &basisSet, int numBasisFct ) const
137 {
138 const GeometryType geomType = geo.type();
139 typedef typename MassMatrixStorageType::iterator iterator;
140 MassMatrixStorageType &massMap = localInverseMassMatrix_[ GlobalGeometryTypeIndex::index( geomType ) ];
141
142 auto it = massMap.find( numBasisFct );
143 if( it == massMap.end() )
144 {
145 std::pair< iterator, bool > insertPair = massMap.insert( std::make_pair( numBasisFct, MatrixPairType(nullptr,nullptr) ) );
146 it = insertPair.first;
147 insertPair.first->second.first.reset( new MatrixType( numBasisFct, numBasisFct, 0.0 ));
148 MatrixType& matrix = insertPair.first->second.first.operator *();
149 VolumeQuadratureType volQuad( entity, volumeQuadratureOrder( entity ) );
150 buildMatrixNoMassFactor( entity, geo, basisSet, volQuad, numBasisFct, matrix, false );
151 try {
152 matrix.invert();
153 }
154 catch ( Dune::FMatrixError &e )
155 {
156 std::cerr << "Matrix is singular:" << std::endl << matrix << std::endl;
157 std::terminate();
158 }
159 const bool diagonal = checkDiagonalMatrix( matrix );
160 // store diagonal if matrix is diagonal
161 if( diagonal )
162 {
163 insertPair.first->second.second.reset( new VectorType( matrix.rows() ) );
164 VectorType& diag = insertPair.first->second.second.operator *();
165 const int rows = matrix.rows();
166 for( int row=0; row<rows; ++row )
167 {
168 diag[ row ] = matrix[ row ][ row ];
169 }
170 }
171 }
172
173 return it->second;
174 }
175
176 template< class MassCaller, class BasisFunctionSet >
177 MatrixType &getLocalInverseMassMatrixDefault ( MassCaller &caller, const EntityType &entity,
178 const Geometry &geo, const BasisFunctionSet &basisSet ) const
179 {
180 const int numDofs = basisSet.size();
181 // if sequence changed or entity index changed: recompute mass matrix
182 if( entityHasChanged( entity ) || (numDofs != int( matrix_.rows())) )
183 {
184 // resize temporary memory if necessary
185 if( numDofs != int( matrix_.rows() ) )
186 matrix_.resize( numDofs, numDofs );
187
188 buildMatrix( caller, entity, geo, basisSet, numDofs, matrix_ );
189 matrix_.invert();
190 }
191
192 // make sure that rhs_ has the correct size
193 assert( numDofs == int( matrix_.rows() ) );
194 return matrix_;
195 }
196
197 // return number of max non blocked dofs
198 int maxNumDofs () const
199 {
200 return space().blockMapper().maxNumDofs() * localBlockSize;
201 }
202
205 {
206 return volumeQuadratureOrder_( space().order() );
207 }
208
209 public:
211 int volumeQuadratureOrder ( const EntityType &entity ) const
212 {
213 return volumeQuadratureOrder_( space().order( entity ) );
214 }
215
217 explicit LocalMassMatrixImplementation ( const DiscreteFunctionSpaceType &spc, int volQuadOrd )
218 : LocalMassMatrixImplementation( spc, [volQuadOrd](const int order) { return volQuadOrd; } )
219 {}
220
223 std::function<int(const int)> volQuadOrderFct = [](const int order) { return Capabilities::DefaultQuadrature< DiscreteFunctionSpaceType >::volumeOrder(order); } )
224 : spc_( referenceToSharedPtr( spc ) )
225 , indexSet_( space().indexSet() )
227 , volumeQuadratureOrder_ ( volQuadOrderFct )
228 , affine_ ( setup() )
229 , rhs_(), row_(), matrix_()
230 , phi_( maxNumDofs() )
231 , phiMass_( maxNumDofs() )
232 , localInverseMassMatrix_( GlobalGeometryTypeIndex :: size( GridType::dimension ) )
233 , lastEntityIndex_( -1 )
234 , lastTopologyId_( ~0u )
235 , sequence_( -1 )
236 {}
237
240 : spc_(other.spc_),
241 indexSet_( space().indexSet() ),
244 affine_( other.affine_ ),
245 rhs_( other.rhs_ ), row_( other.row_ ), matrix_( other.matrix_ ),
246 phi_( other.phi_ ),
247 phiMass_( other.phiMass_ ),
248 localInverseMassMatrix_( GlobalGeometryTypeIndex :: size( GridType::dimension ) ),
251 sequence_( other.sequence_ )
252 {}
253
254 public:
256 bool affine () const { return affine_; }
257
259 double getAffineMassFactor(const Geometry& geo) const
260 {
261 return geoInfo_.referenceVolume( geo.type() ) / geo.volume();
262 }
263
264 template< class BasisFunctionSet >
266 {
267 static const int quadPointSetId = SelectQuadraturePointSetId< VolumeQuadratureType >::value;
268 // if BasisFunctionSet does not have an static int member called pointSetId then this will be -1
269 static const int basePointSetId = detail::SelectPointSetId< BasisFunctionSet >::value;
270 // for Lagrange-type basis evaluated on interpolation points
271 // this is the Kronecker delta, so the mass matrix is diagonal even
272 // on non affine grids
273 if constexpr ( quadPointSetId == basePointSetId )
274 {
275 const unsigned int numShapeFunctions = bfs.size() / dimRange;
277 .isInterpolationQuadrature(numShapeFunctions);
278 }
279 return false;
280 }
281
284 template< class MassCaller, class BasisFunctionSet, class LocalFunction >
285 void applyInverse ( MassCaller &caller, const EntityType &entity, const BasisFunctionSet &basisFunctionSet, LocalFunction &lf ) const
286 {
287 Geometry geo = entity.geometry();
288 if( ( affine() || geo.affine() || checkInterpolationBFS(basisFunctionSet) )
289 && !caller.hasMass() )
290 applyInverseLocally( entity, geo, basisFunctionSet, lf );
291 else
292 applyInverseDefault( caller, entity, geo, basisFunctionSet, lf );
293 }
294
297 template< class MassCaller, class LocalFunction >
298 void applyInverse ( MassCaller &caller, const EntityType &entity, LocalFunction &lf ) const
299 {
300 applyInverse( caller, entity, lf.basisFunctionSet(), lf );
301 }
302
304 template< class LocalFunction >
305 void applyInverse ( const EntityType &entity, LocalFunction &lf ) const
306 {
307 NoMassDummyCaller caller;
308 applyInverse( caller, entity, lf.basisFunctionSet(), lf );
309 }
310 template< class BasisFunctionSet, class LocalFunction >
311 void applyInverse ( const EntityType &entity, const BasisFunctionSet &basisFunctionSet, LocalFunction &lf ) const
312 {
313 NoMassDummyCaller caller;
314 applyInverse( caller, entity, basisFunctionSet, lf );
315 }
316
317
319 template< class LocalFunction >
320 void applyInverse ( LocalFunction &lf ) const
321 {
322 applyInverse( lf.entity(), lf );
323 }
324
326 template< class LocalMatrix >
327 void rightMultiplyInverse ( LocalMatrix &localMatrix ) const
328 {
329 const EntityType &entity = localMatrix.rangeEntity();
330 Geometry geo = entity.geometry();
331 if( ( affine() || geo.affine() || checkInterpolationBFS(localMatrix.rangeBasisFunctionSet())) )
332 rightMultiplyInverseLocally( entity, geo, localMatrix );
333 else
334 rightMultiplyInverseDefault( entity, geo, localMatrix );
335 }
336
338 template< class LocalMatrix >
339 void leftMultiplyInverse ( LocalMatrix &localMatrix ) const
340 {
341 const EntityType &entity = localMatrix.domainEntity();
342 Geometry geo = entity.geometry();
343 if( ( affine() || geo.affine() || checkInterpolationBFS(localMatrix.rangeBasisFunctionSet())) )
344 leftMultiplyInverseLocally( entity, geo, localMatrix );
345 else
346 leftMultiplyInverseDefault( entity, geo, localMatrix );
347 }
348
349 const DiscreteFunctionSpaceType &space () const { return *spc_; }
350
352 // end of public methods
354
355 protected:
357 // applyInverse for DG spaces
359 template< class MassCaller, class BasisFunctionSet, class LocalFunction >
360 void applyInverseDgOrthoNormalBasis ( MassCaller &caller, const EntityType &entity,
361 const BasisFunctionSet &basisFunctionSet, LocalFunction &lf ) const
362 {
363 Geometry geo = entity.geometry();
364 assert( dgNumDofs == lf.size() );
365
366 // affine_ can be a static information
367 const bool isAffine = affine() || geo.affine();
368 // make sure that for affine grids the geometry info is also correct
369 assert( affine() ? geo.affine() : true );
370
371 // in case of affine mappings we only have to multiply with a factor
372 if( isAffine && !caller.hasMass() )
373 {
374 const double massVolInv = getAffineMassFactor( geo );
375
376 // apply inverse mass matrix
377 for( int l = 0; l < dgNumDofs; ++l )
378 lf[ l ] *= massVolInv;
379 }
380 else
381 {
382 // copy local function to right hand side
383 for( int l = 0; l < dgNumDofs; ++l )
384 dgRhs_[ l ] = lf[ l ];
385
386 buildMatrix( caller, entity, geo, basisFunctionSet, dgNumDofs, dgMatrix_ );
387 dgMatrix_.solve( dgX_, dgRhs_ );
388
389 // copy back to local function
390 for( int l = 0; l < dgNumDofs; ++l )
391 lf[ l ] = dgX_[ l ];
392 }
393 }
394
396 template< class LocalMatrix >
397 void rightMultiplyInverseDgOrthoNormalBasis ( LocalMatrix &localMatrix ) const
398 {
399 const EntityType &entity = localMatrix.rangeEntity();
400 Geometry geo = entity.geometry();
401 assert( dgNumDofs == localMatrix.columns() );
402
403 // in case of affine mappings we only have to multiply with a factor
404 if( affine() || geo.affine() )
405 {
406 localMatrix.scale( getAffineMassFactor( geo ) );
407 }
408 else
409 {
410 NoMassDummyCaller caller;
411 buildMatrix( caller, entity, geo, localMatrix.rangeBasisFunctionSet(), dgNumDofs, dgMatrix_ );
412 dgMatrix_.invert();
413
414 const int rows = localMatrix.rows();
415 for( int i = 0; i < rows; ++i )
416 {
417 for( int j = 0; j < dgNumDofs; ++j )
418 dgRhs_[ j ] = localMatrix.get( i, j );
419 dgMatrix_.mtv( dgRhs_, dgX_ );
420 for( int j = 0; j < dgNumDofs; ++j )
421 localMatrix.set( i, j, dgX_[ j ] );
422 }
423 }
424 }
425
427 template< class LocalMatrix >
428 void leftMultiplyInverseDgOrthoNormalBasis ( LocalMatrix &localMatrix ) const
429 {
430 const EntityType &entity = localMatrix.domainEntity();
431 Geometry geo = entity.geometry();
432 assert( dgNumDofs == localMatrix.columns() );
433
434 // in case of affine mappings we only have to multiply with a factor
435 if( affine() || geo.affine() )
436 {
437 localMatrix.scale( getAffineMassFactor( geo ) );
438 }
439 else
440 {
441 NoMassDummyCaller caller;
442 buildMatrix( caller, entity, geo, localMatrix.domainBasisFunctionSet(), dgNumDofs, dgMatrix_ );
443 dgMatrix_.invert();
444
445 const int rows = localMatrix.rows();
446 for( int i = 0; i < rows; ++i )
447 {
448 for( int j = 0; j < dgNumDofs; ++j )
449 dgRhs_[ j ] = localMatrix.get( i, j );
450 dgMatrix_.mv( dgRhs_, dgX_ );
451 for( int j = 0; j < dgNumDofs; ++j )
452 localMatrix.set( i, j, dgX_[ j ] );
453 }
454 }
455 }
456
458 bool entityHasChanged( const EntityType& entity ) const
459 {
460 // don't compute matrix new for the same entity
461 const int currentSequence = space().sequence();
462 const unsigned int topologyId = entity.type().id();
463 const IndexType entityIndex = indexSet_.index( entity ) ;
464
465 // check whether sequence has been updated
466 if( sequence_ != currentSequence ||
467 lastEntityIndex_ != entityIndex ||
468 lastTopologyId_ != topologyId )
469 {
470 // update identifiers
471 lastEntityIndex_ = entityIndex ;
472 sequence_ = currentSequence;
473 lastTopologyId_ = topologyId ;
474
475 return true ;
476 }
477 else
478 // the entity did not change
479 return false ;
480 }
481
483 // standard applyInverse method
487 template< class MassCaller, class BasisFunctionSet, class LocalFunction >
488 void applyInverseDefault ( MassCaller &caller, const EntityType &entity,
489 const Geometry &geo, const BasisFunctionSet &basisFunctionSet, LocalFunction &lf ) const
490 {
491 // get local inverted mass matrix
492 MatrixType &invMassMatrix
493 = getLocalInverseMassMatrixDefault ( caller, entity, geo, basisFunctionSet );
494
495 // copy local function to right hand side
496 const int numDofs = lf.size();
497 rhs_.resize( numDofs );
498 for( int l = 0; l < numDofs; ++l )
499 rhs_[ l ] = lf[ l ];
500
501 // apply inverse to right hand side and store in lf
502 multiply( numDofs, invMassMatrix, rhs_, lf );
503 }
504
506 template< class LocalMatrix >
507 void rightMultiplyInverseDefault ( const EntityType &entity, const Geometry &geo, LocalMatrix &localMatrix ) const
508 {
509 NoMassDummyCaller caller;
510 MatrixType &invMassMatrix
511 = getLocalInverseMassMatrixDefault ( caller, entity, geo, localMatrix.rangeBasisFunctionSet() );
512
513 const int cols = localMatrix.columns();
514 rhs_.resize( cols );
515 row_.resize( cols );
516
517 const int rows = localMatrix.rows();
518 for( int i = 0; i < rows; ++i )
519 {
520 // get i-th row from localMatrix
521 for( int j = 0; j < cols; ++j )
522 rhs_[ j ] = localMatrix.get( i, j );
523
524 // multiply with all columns of inverse mass matrix
525 invMassMatrix.mtv( rhs_, row_ );
526
527 // store as i-th row in localMatrix
528 for( int j = 0; j < cols; ++j )
529 localMatrix.set( i, j, row_[ j ] );
530 }
531 }
532
534 template< class LocalMatrix >
535 void leftMultiplyInverseDefault ( const EntityType &entity, const Geometry &geo, LocalMatrix &localMatrix ) const
536 {
537 NoMassDummyCaller caller;
538 MatrixType &invMassMatrix
539 = getLocalInverseMassMatrixDefault ( caller, entity, geo, localMatrix.rangeBasisFunctionSet() );
540
541 const int cols = localMatrix.columns();
542 rhs_.resize( cols );
543 row_.resize( cols );
544
545 const int rows = localMatrix.rows();
546 for( int i = 0; i < rows; ++i )
547 {
548 // get i-th column from localMatrix
549 for( int j = 0; j < cols; ++j )
550 rhs_[ j ] = localMatrix.get( j, i );
551
552 // multiply with all rows in inverse mass matrix
553 invMassMatrix.mv( rhs_, row_ );
554
555 // store as i-th column in localMatrix
556 for( int j = 0; j < cols; ++j )
557 localMatrix.set( j, i, row_[ j ] );
558 }
559 }
560
561
563 // local applyInverse method for affine geometries
566 template< class BasisFunctionSet, class LocalFunction >
567 void applyInverseLocally ( const EntityType &entity,
568 const Geometry &geo, const BasisFunctionSet &basisFunctionSet, LocalFunction &lf ) const
569 {
570 const int numDofs = lf.size();
571
572 // get local inverted mass matrix
573 MatrixPairType& matrixPair =
574 getLocalInverseMassMatrix( entity, geo, basisFunctionSet, numDofs );
575
576 // if diagonal exists then matrix is in diagonal form
577 if( matrixPair.second )
578 {
579 const VectorType& diagonal = *matrixPair.second;
580 assert( int(diagonal.size()) == numDofs );
581
582 VolumeQuadratureType volQuad( entity, volumeQuadratureOrder( entity ) );
583
584 const int nop = volQuad.nop();
585 assert(nop*dimRange == numDofs);
586 for( int l=0, qt = 0; qt < nop; ++qt )
587 {
588 const auto intel = geo.integrationElement( volQuad.point(qt) );
589 for (int r = 0; r < dimRange; ++r, ++l )
590 {
591 lf[ l ] *= diagonal[ l ] / intel;
592 }
593 }
594 }
595 else
596 {
597 const double massVolInv = getAffineMassFactor( geo );
598 // copy local function to right hand side
599 // and apply inverse mass volume fraction
600 rhs_.resize( numDofs );
601 for( int l = 0; l < numDofs; ++l )
602 rhs_[ l ] = lf[ l ] * massVolInv;
603
604 const MatrixType& invMassMatrix = *matrixPair.first;
605 // apply inverse local mass matrix and store in lf
606 multiply( numDofs, invMassMatrix, rhs_, lf );
607 }
608 }
609
610 template <class LocalMatrix>
611 const VectorType&
612 setupInverseDiagonal( const EntityType &entity, const Geometry &geo,
613 const VectorType& refElemDiagonal,
614 LocalMatrix &localMatrix ) const
615 {
616 const int cols = localMatrix.columns();
617
618 assert( int(refElemDiagonal.size()) == cols );
619
620 VolumeQuadratureType volQuad( entity, volumeQuadratureOrder( entity ) );
621
622 VectorType& elementDiagonal = rhs_;
623 elementDiagonal.resize( cols );
624
625 const int nop = volQuad.nop();
626 assert(nop*dimRange == cols);
627 for( int l = 0, qt = 0; qt < nop; ++qt )
628 {
629 const auto intel = geo.integrationElement( volQuad.point(qt) );
630 for (int r = 0; r < dimRange; ++r,++l )
631 {
632 elementDiagonal[ l ] = refElemDiagonal[ l ] / intel;
633 }
634 }
635 return elementDiagonal;
636 }
637
638 template< class LocalMatrix >
639 void rightMultiplyInverseLocally ( const EntityType &entity, const Geometry &geo, LocalMatrix &localMatrix ) const
640 {
641 const int cols = localMatrix.columns();
642 MatrixPairType& matrixPair =
643 getLocalInverseMassMatrix( entity, geo, localMatrix.rangeBasisFunctionSet(), cols );
644
645 // if diagonal exists then matrix is in diagonal form
646 // stored as inverse on the reference element
647 if( matrixPair.second )
648 {
649 const VectorType& elementDiagonal =
650 setupInverseDiagonal( entity, geo, *matrixPair.second, localMatrix );
651
652 row_.resize( cols );
653 const int rows = localMatrix.rows();
654 for( int i = 0; i < rows; ++i )
655 {
656 // get i-th row from localMatrix
657 // and multiply with diagonal of inverse mass matrix
658 for( int j = 0; j < cols; ++j )
659 row_[ j ] = elementDiagonal[ j ] * localMatrix.get( i, j );
660
661 // store as i-th row in localMatrix
662 for( int j = 0; j < cols; ++j )
663 localMatrix.set( i, j, row_[ j ] );
664 }
665 }
666 else
667 {
668 const MatrixType &invMassMatrix = *matrixPair.first;
669
670 const double massVolInv = getAffineMassFactor( geo );
671
672 rhs_.resize( cols );
673 row_.resize( cols );
674
675 const int rows = localMatrix.rows();
676 for( int i = 0; i < rows; ++i )
677 {
678 // get i-th row from localMatrix
679 // and multiply with diagonal of inverse mass matrix
680 for( int j = 0; j < cols; ++j )
681 rhs_[ j ] = localMatrix.get( i, j ) * massVolInv;
682
683 // multiply with all columns of inverse mass matrix
684 invMassMatrix.mtv( rhs_, row_ );
685
686 // store as i-th row of localMatrix
687 for( int j = 0; j < cols; ++j )
688 localMatrix.set( i, j, row_[ j ] );
689 }
690 }
691 }
692
694 template< class LocalMatrix >
695 void leftMultiplyInverseLocally ( const EntityType &entity, const Geometry &geo, LocalMatrix &localMatrix ) const
696 {
697 const int cols = localMatrix.columns();
698 MatrixPairType& matrixPair =
699 getLocalInverseMassMatrix( entity, geo, localMatrix.rangeBasisFunctionSet(), cols );
700
701 // if diagonal exists then matrix is in diagonal form
702 // stored as inverse on the reference element
703 if( matrixPair.second )
704 {
705 const VectorType& elementDiagonal =
706 setupInverseDiagonal( entity, geo, *matrixPair.second, localMatrix );
707
708 row_.resize( cols );
709 const int rows = localMatrix.rows();
710 for( int i = 0; i < rows; ++i )
711 {
712 // get i-th column from localMatrix
713 // and multiply with diagonal of inverse mass matrix
714 for( int j = 0; j < cols; ++j )
715 row_[ j ] = elementDiagonal[ j ] * localMatrix.get( j, i );
716
717 // store as i-th column of localMatrix
718 for( int j = 0; j < cols; ++j )
719 localMatrix.set( j, i, row_[ j ] );
720 }
721 }
722 else
723 {
724 const MatrixType &invMassMatrix = *matrixPair.first;
725
726 const double massVolInv = getAffineMassFactor( geo );
727
728 rhs_.resize( cols );
729 row_.resize( cols );
730
731 const int rows = localMatrix.rows();
732 for( int i = 0; i < rows; ++i )
733 {
734 // get i-th column from localMatrix
735 for( int j = 0; j < cols; ++j )
736 rhs_[ j ] = localMatrix.get( j, i ) * massVolInv;
737
738 // apply to all rows of inverse mass matrix
739 invMassMatrix.mv( rhs_, row_ );
740
741 // store as i-th column of localMatrix
742 for( int j = 0; j < cols; ++j )
743 localMatrix.set( j, i, row_[ j ] );
744 }
745 }
746 }
747
748
750 bool setup () const
751 {
752 // for structured grids this is always true
753 if( StructuredGrid )
754 return true;
755
756 // get types for codim 0
757 const std::vector<Dune::GeometryType>& geomTypes = geoInfo_.geomTypes(0);
758
759 // for simplices we also have affine mappings
760 if( (geomTypes.size() == 1) && geomTypes[0].isSimplex() )
761 {
762 return true;
763 }
764
765 // otherwise use geometry affinity
766 return false ;
767 }
768
770 template< class MassCaller, class Matrix >
771 void buildMatrix ( MassCaller &caller, const EntityType &entity,
772 const Geometry &geo, const BasisFunctionSetType &set,
773 std::size_t numDofs, Matrix &matrix ) const
774 {
775 assert( numDofs == set.size() );
776
777 // clear matrix
778 matrix = 0;
779
780 // create quadrature
781 VolumeQuadratureType volQuad( entity, volumeQuadratureOrder( entity ) );
782
783 if( caller.hasMass() )
784 buildMatrixWithMassFactor( caller, entity, geo, set, volQuad, numDofs, matrix );
785 else
786 buildMatrixNoMassFactor( entity, geo, set, volQuad, numDofs, matrix );
787 }
788
790 template <class Matrix>
792 const EntityType& en,
793 const Geometry& geo,
794 const BasisFunctionSetType& set,
795 const VolumeQuadratureType& volQuad,
796 const int numDofs,
797 Matrix& matrix,
798 const bool applyIntegrationElement = true ) const
799 {
800 const int volNop = volQuad.nop();
801 for(int qp=0; qp<volNop; ++qp)
802 {
803 // calculate integration weight
804 const double intel = ( applyIntegrationElement ) ?
805 ( volQuad.weight(qp) * geo.integrationElement(volQuad.point(qp)) ) : volQuad.weight(qp) ;
806
807 // eval basis functions
808 set.evaluateAll(volQuad[qp], phi_);
809
810 for(int m=0; m<numDofs; ++m)
811 {
812 const RangeType& phi_m = phi_[m];
813 const ctype val = intel * (phi_m * phi_m);
814 matrix[m][m] += val;
815
816 for(int k=m+1; k<numDofs; ++k)
817 {
818 const ctype val = intel * (phi_m * phi_[k]);
819 matrix[m][k] += val;
820 matrix[k][m] += val;
821 }
822 }
823 }
824 }
825
827 template <class MassCallerType, class Matrix>
829 MassCallerType& caller,
830 const EntityType& en,
831 const Geometry& geo,
832 const BasisFunctionSetType& set,
833 const VolumeQuadratureType& volQuad,
834 const int numDofs,
835 Matrix& matrix) const
836 {
837 typedef typename MassCallerType :: MassFactorType MassFactorType;
838 MassFactorType mass;
839
840 const int volNop = volQuad.nop();
841 for(int qp=0; qp<volNop; ++qp)
842 {
843 // calculate integration weight
844 const double intel = volQuad.weight(qp)
845 * geo.integrationElement(volQuad.point(qp));
846
847 // eval basis functions
848 set.evaluateAll( volQuad[qp], phi_);
849
850 // call mass factor
851 caller.mass( en, volQuad, qp, mass);
852
853 // apply mass matrix to all basis functions
854 for(int m=0; m<numDofs; ++m)
855 {
856 mass.mv( phi_[m], phiMass_[m] );
857 }
858
859 // add values to matrix
860 for(int m=0; m<numDofs; ++m)
861 {
862 for(int k=0; k<numDofs; ++k)
863 {
864 matrix[m][k] += intel * (phiMass_[m] * phi_[k]);
865 }
866 }
867 }
868 }
869
870 // implement matvec with matrix (mv of densematrix is too stupid)
871 template <class Matrix, class Rhs, class X>
872 void multiply( const int size,
873 const Matrix& matrix,
874 const Rhs& rhs,
875 X& x ) const
876 {
877 assert( (int) matrix.rows() == size );
878 assert( (int) matrix.cols() == size );
879 assert( (int) rhs.size() == size );
880
881 for( int row = 0; row < size; ++ row )
882 {
883 RangeFieldType sum = 0;
884 // get matrix row
885 typedef typename Matrix :: const_row_reference MatRow;
886 MatRow matRow = matrix[ row ];
887
888 // multiply row with right hand side
889 for( int col = 0; col < size; ++ col )
890 {
891 sum += matRow[ col ] * rhs[ col ];
892 }
893
894 // set to result to result vector
895 x[ row ] = sum;
896 }
897 }
898 };
899
900
901
902 // LocalMassMatrix
903 // ---------------
904
906 template< class DiscreteFunctionSpace, class VolumeQuadrature >
908 : public LocalMassMatrixImplementation< DiscreteFunctionSpace, VolumeQuadrature >
909 {
911
912 public:
913 // copy base class constructors
914 using BaseType :: LocalMassMatrixImplementation;
915 };
916
917
918
920 //
921 // DG LocalMassMatrix Implementation
922 //
924
926 template< class DiscreteFunctionSpace, class VolumeQuadrature >
928 : public LocalMassMatrixImplementation< DiscreteFunctionSpace, VolumeQuadrature >
929 {
931
932 public:
933 typedef typename BaseType :: EntityType EntityType;
934
935 // copy base class constructors
936 using BaseType :: LocalMassMatrixImplementation;
937
940 template <class MassCallerType, class BasisFunction, class LocalFunctionType>
941 void applyInverse(MassCallerType& caller,
942 const EntityType& en,
943 const BasisFunction &basisFunction,
944 LocalFunctionType& lf) const
945 {
946 BaseType :: applyInverseDgOrthoNormalBasis( caller, en, basisFunction, lf );
947 }
948 template <class MassCallerType, class LocalFunctionType>
949 void applyInverse(MassCallerType& caller,
950 const EntityType& en,
951 LocalFunctionType& lf) const
952 {
953 BaseType :: applyInverseDgOrthoNormalBasis( caller, en, lf.basisFunctionSet(), lf );
954 }
955
957 template <class LocalFunctionType>
958 void applyInverse(const EntityType& en,
959 LocalFunctionType& lf) const
960 {
961 typename BaseType :: NoMassDummyCaller caller;
962 applyInverse(caller, en, lf );
963 }
964
965 template <class BasisFunction, class LocalFunctionType>
966 void applyInverse(const EntityType& en,
967 const BasisFunction &basisFunction,
968 LocalFunctionType& lf) const
969 {
970 typename BaseType :: NoMassDummyCaller caller;
971 applyInverse(caller, en, basisFunction, lf );
972 }
973
975 template< class LocalFunction >
976 void applyInverse ( LocalFunction &lf ) const
977 {
978 applyInverse( lf.entity(), lf.basisFunctionSet(), lf );
979 }
980
982 template< class LocalMatrix >
983 void rightMultiplyInverse ( LocalMatrix &localMatrix ) const
984 {
986 }
987
989 template< class LocalMatrix >
990 void leftMultiplyInverse ( LocalMatrix &localMatrix ) const
991 {
993 }
994 };
995
997
998 } // namespace Fem
999
1000} // namespace Dune
1001
1002#endif // #ifndef DUNE_FEM_LOCALMASSMATRIX_HH
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
Definition: bindguard.hh:11
detail::SelectPointSetId< Quadrature, -Dune::QuadratureType::size > SelectQuadraturePointSetId
Definition: quadrature.hh:541
static std::shared_ptr< T > referenceToSharedPtr(T &t)
Definition: memory.hh:19
interface for local functions
Definition: localfunction.hh:77
const EntityType & entity() const
obtain the entity, this local function lives on
Definition: localfunction.hh:302
const BasisFunctionSetType & basisFunctionSet() const
obtain the basis function set for this local function
Definition: localfunction.hh:296
SizeType size() const
obtain the number of local DoFs
Definition: localfunction.hh:360
Helper class to check affinity of the grid's geometries.
Definition: checkgeomaffinity.hh:22
Local Mass Matrix inversion implementation, select the correct method in your implementation.
Definition: localmassmatrix.hh:34
DiscreteFunctionSpaceType::EntityType EntityType
Definition: localmassmatrix.hh:58
void applyInverse(const EntityType &entity, const BasisFunctionSet &basisFunctionSet, LocalFunction &lf) const
Definition: localmassmatrix.hh:311
void applyInverse(MassCaller &caller, const EntityType &entity, LocalFunction &lf) const
Definition: localmassmatrix.hh:298
std::vector< RangeType > phi_
Definition: localmassmatrix.hh:91
MatrixType & getLocalInverseMassMatrixDefault(MassCaller &caller, const EntityType &entity, const Geometry &geo, const BasisFunctionSet &basisSet) const
Definition: localmassmatrix.hh:177
AllGeomTypes< typename GridPartType ::IndexSetType, GridType > GeometryInformationType
Definition: localmassmatrix.hh:68
DiscreteFunctionSpaceType::IndexSetType IndexSetType
Definition: localmassmatrix.hh:52
bool checkInterpolationBFS(const BasisFunctionSet &bfs) const
Definition: localmassmatrix.hh:265
int volumeQuadratureOrder(const EntityType &entity) const
return appropriate quadrature order, default is 2 * order(entity)
Definition: localmassmatrix.hh:211
void applyInverseDefault(MassCaller &caller, const EntityType &entity, const Geometry &geo, const BasisFunctionSet &basisFunctionSet, LocalFunction &lf) const
Definition: localmassmatrix.hh:488
void applyInverseLocally(const EntityType &entity, const Geometry &geo, const BasisFunctionSet &basisFunctionSet, LocalFunction &lf) const
apply local mass matrix to local function lf
Definition: localmassmatrix.hh:567
const bool affine_
Definition: localmassmatrix.hh:81
const DiscreteFunctionSpaceType & space() const
Definition: localmassmatrix.hh:349
bool setup() const
setup and return affinity
Definition: localmassmatrix.hh:750
DiscreteFunctionSpaceType::RangeFieldType ctype
Definition: localmassmatrix.hh:39
GridPartType::GridType GridType
Definition: localmassmatrix.hh:57
Dune::FieldVector< ctype, dgNumDofs > DGVectorType
Definition: localmassmatrix.hh:48
const IndexSetType & indexSet_
Definition: localmassmatrix.hh:77
void applyInverse(const EntityType &entity, LocalFunction &lf) const
apply local dg mass matrix to local function lf without mass factor
Definition: localmassmatrix.hh:305
@ dimRange
Definition: localmassmatrix.hh:43
DGMatrixType dgMatrix_
Definition: localmassmatrix.hh:83
LocalMassMatrixImplementation(const DiscreteFunctionSpaceType &spc, int volQuadOrd)
constructor taking space and volume quadrature order
Definition: localmassmatrix.hh:217
Fem::GeometryAffinityCheck< VolumeQuadratureType > GeometryAffinityCheckType
Definition: localmassmatrix.hh:63
double getAffineMassFactor(const Geometry &geo) const
return mass factor for diagonal mass matrix
Definition: localmassmatrix.hh:259
@ StructuredGrid
Definition: localmassmatrix.hh:66
int maxNumDofs() const
Definition: localmassmatrix.hh:198
IndexSetType::IndexType IndexType
Definition: localmassmatrix.hh:53
VolumeQuadrature VolumeQuadratureType
Definition: localmassmatrix.hh:61
std::vector< MassMatrixStorageType > LocalInverseMassMatrixStorageType
Definition: localmassmatrix.hh:96
std::vector< RangeType > phiMass_
Definition: localmassmatrix.hh:92
LocalMassMatrixImplementation(const ThisType &other)
copy constructor
Definition: localmassmatrix.hh:239
const VectorType & setupInverseDiagonal(const EntityType &entity, const Geometry &geo, const VectorType &refElemDiagonal, LocalMatrix &localMatrix) const
Definition: localmassmatrix.hh:612
void applyInverse(LocalFunction &lf) const
apply local dg mass matrix to local function lf without mass factor
Definition: localmassmatrix.hh:320
Dune::DynamicVector< RangeFieldType > VectorType
Definition: localmassmatrix.hh:73
bool checkDiagonalMatrix(const MatrixType &matrix) const
Definition: localmassmatrix.hh:118
DGVectorType dgRhs_
Definition: localmassmatrix.hh:84
void applyInverse(MassCaller &caller, const EntityType &entity, const BasisFunctionSet &basisFunctionSet, LocalFunction &lf) const
Definition: localmassmatrix.hh:285
void leftMultiplyInverseDefault(const EntityType &entity, const Geometry &geo, LocalMatrix &localMatrix) const
compute M^-1 * localMatrix
Definition: localmassmatrix.hh:535
@ localBlockSize
Definition: localmassmatrix.hh:44
std::pair< std::unique_ptr< MatrixType >, std::unique_ptr< VectorType > > MatrixPairType
Definition: localmassmatrix.hh:94
MatrixPairType & getLocalInverseMassMatrix(const EntityType &entity, const Geometry &geo, const BasisFunctionSet &basisSet, int numBasisFct) const
Definition: localmassmatrix.hh:135
void multiply(const int size, const Matrix &matrix, const Rhs &rhs, X &x) const
Definition: localmassmatrix.hh:872
bool entityHasChanged(const EntityType &entity) const
returns true if the entity has been changed
Definition: localmassmatrix.hh:458
DiscreteFunctionSpaceType::RangeType RangeType
Definition: localmassmatrix.hh:41
DGVectorType dgX_
Definition: localmassmatrix.hh:84
LocalInverseMassMatrixStorageType localInverseMassMatrix_
Definition: localmassmatrix.hh:98
void buildMatrixWithMassFactor(MassCallerType &caller, const EntityType &en, const Geometry &geo, const BasisFunctionSetType &set, const VolumeQuadratureType &volQuad, const int numDofs, Matrix &matrix) const
build local mass matrix with mass factor
Definition: localmassmatrix.hh:828
DiscreteFunctionSpaceType::BasisFunctionSetType BasisFunctionSetType
Definition: localmassmatrix.hh:55
DiscreteFunctionSpace DiscreteFunctionSpaceType
Definition: localmassmatrix.hh:38
DiscreteFunctionSpaceType::RangeFieldType RangeFieldType
Definition: localmassmatrix.hh:40
MatrixType matrix_
Definition: localmassmatrix.hh:89
void rightMultiplyInverseDefault(const EntityType &entity, const Geometry &geo, LocalMatrix &localMatrix) const
compute localMatrix * M^-1
Definition: localmassmatrix.hh:507
void leftMultiplyInverseDgOrthoNormalBasis(LocalMatrix &localMatrix) const
compute M^-1 * localMatrix
Definition: localmassmatrix.hh:428
bool affine() const
returns true if geometry mapping is affine
Definition: localmassmatrix.hh:256
VectorType rhs_
Definition: localmassmatrix.hh:88
IndexType lastEntityIndex_
Definition: localmassmatrix.hh:101
@ dgNumDofs
Definition: localmassmatrix.hh:45
Dune::FieldMatrix< ctype, dgNumDofs, dgNumDofs > DGMatrixType
Definition: localmassmatrix.hh:47
void rightMultiplyInverseLocally(const EntityType &entity, const Geometry &geo, LocalMatrix &localMatrix) const
Definition: localmassmatrix.hh:639
GeometryInformationType::DomainType DomainType
Definition: localmassmatrix.hh:69
void applyInverseDgOrthoNormalBasis(MassCaller &caller, const EntityType &entity, const BasisFunctionSet &basisFunctionSet, LocalFunction &lf) const
Definition: localmassmatrix.hh:360
LocalMassMatrixImplementation(const DiscreteFunctionSpaceType &spc, std::function< int(const int)> volQuadOrderFct=[](const int order) { return Capabilities::DefaultQuadrature< DiscreteFunctionSpaceType >::volumeOrder(order);})
constructor taking space and volume quadrature order
Definition: localmassmatrix.hh:222
void rightMultiplyInverse(LocalMatrix &localMatrix) const
compute localMatrix * M^-1
Definition: localmassmatrix.hh:327
const std::function< int(const int)> volumeQuadratureOrder_
Definition: localmassmatrix.hh:80
unsigned int lastTopologyId_
Definition: localmassmatrix.hh:102
void leftMultiplyInverse(LocalMatrix &localMatrix) const
compute M^-1 * localMatrix
Definition: localmassmatrix.hh:339
void buildMatrix(MassCaller &caller, const EntityType &entity, const Geometry &geo, const BasisFunctionSetType &set, std::size_t numDofs, Matrix &matrix) const
build local mass matrix
Definition: localmassmatrix.hh:771
GeometryInformationType geoInfo_
Definition: localmassmatrix.hh:79
VectorType row_
Definition: localmassmatrix.hh:88
int maxVolumeQuadratureOrder() const
return appropriate quadrature order, default is 2 * order()
Definition: localmassmatrix.hh:204
void buildMatrixNoMassFactor(const EntityType &en, const Geometry &geo, const BasisFunctionSetType &set, const VolumeQuadratureType &volQuad, const int numDofs, Matrix &matrix, const bool applyIntegrationElement=true) const
build local mass matrix with mass factor
Definition: localmassmatrix.hh:791
std::shared_ptr< const DiscreteFunctionSpaceType > spc_
Definition: localmassmatrix.hh:76
EntityType::Geometry Geometry
Definition: localmassmatrix.hh:59
void leftMultiplyInverseLocally(const EntityType &entity, const Geometry &geo, LocalMatrix &localMatrix) const
compute M^-1 * localMatrix
Definition: localmassmatrix.hh:695
void rightMultiplyInverseDgOrthoNormalBasis(LocalMatrix &localMatrix) const
compute localMatrix * M^-1
Definition: localmassmatrix.hh:397
Dune::DynamicMatrix< RangeFieldType > MatrixType
Definition: localmassmatrix.hh:72
std::map< const int, MatrixPairType > MassMatrixStorageType
Definition: localmassmatrix.hh:95
int sequence_
Definition: localmassmatrix.hh:104
DiscreteFunctionSpaceType::GridPartType GridPartType
Definition: localmassmatrix.hh:50
bool hasMass() const
Definition: localmassmatrix.hh:111
Dune::FieldMatrix< ctype, dimRange, dimRange > MassFactorType
Definition: localmassmatrix.hh:108
void mass(const EntityType &, const VolumeQuadratureType &, int, MassFactorType &) const
Definition: localmassmatrix.hh:113
Local Mass Matrix for arbitrary spaces.
Definition: localmassmatrix.hh:909
DG Local Mass Matrix for arbitrary spaces.
Definition: localmassmatrix.hh:929
BaseType::EntityType EntityType
Definition: localmassmatrix.hh:933
void applyInverse(LocalFunction &lf) const
apply local dg mass matrix to local function lf without mass factor
Definition: localmassmatrix.hh:976
void applyInverse(const EntityType &en, LocalFunctionType &lf) const
apply local dg mass matrix to local function lf without mass factor
Definition: localmassmatrix.hh:958
void applyInverse(MassCallerType &caller, const EntityType &en, LocalFunctionType &lf) const
Definition: localmassmatrix.hh:949
void applyInverse(MassCallerType &caller, const EntityType &en, const BasisFunction &basisFunction, LocalFunctionType &lf) const
Definition: localmassmatrix.hh:941
void applyInverse(const EntityType &en, const BasisFunction &basisFunction, LocalFunctionType &lf) const
Definition: localmassmatrix.hh:966
void rightMultiplyInverse(LocalMatrix &localMatrix) const
compute localMatrix * M^-1
Definition: localmassmatrix.hh:983
void leftMultiplyInverse(LocalMatrix &localMatrix) const
compute M^-1 * localMatrix
Definition: localmassmatrix.hh:990
Interface class for basis function sets.
Definition: basisfunctionset/basisfunctionset.hh:31
const EntityType & entity() const
return entity
std::size_t size() const
return size of basis function set
ctype referenceVolume(const GeometryType &type) const
return volume of reference element for geometry of type type
Definition: allgeomtypes.hh:72
default implementation uses method geomTypes of given index set. Used in DiscreteFunctionSpaces.
Definition: allgeomtypes.hh:99
const std ::vector< GeometryType > & geomTypes(unsigned int codim) const
returns vector with geometry tpyes this index set has indices for
Definition: allgeomtypes.hh:171
static int volumeOrder(const int k)
return quadrature order for volume quadratures for given polynomial order k
Definition: space/common/capabilities.hh:138
discrete function space