dune-fem  2.8-git
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 
20 namespace 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;
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 
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 
88  mutable VectorType rhs_, row_;
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() )
226  , geoInfo_( 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() ),
242  geoInfo_( 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;
276  return VolumeQuadratureType( bfs.entity(), volumeQuadratureOrder( bfs.entity() ) )
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
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:
934 
935  // copy base class constructors
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
Definition: bindguard.hh:11
detail::SelectPointSetId< Quadrature, -Dune::QuadratureType::size > SelectQuadraturePointSetId
Definition: quadrature.hh:541
Double abs(const Double &a)
Definition: double.hh:871
static std::shared_ptr< T > referenceToSharedPtr(T &t)
Definition: memory.hh:19
typename Impl::Index< Range >::Type IndexType
Definition: hybrid.hh:69
static constexpr std::decay_t< T > sum(T a)
Definition: utility.hh:33
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
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
AllGeomTypes< typename GridPartType ::IndexSetType, GridType > GeometryInformationType
Definition: localmassmatrix.hh:68
GeometryInformationType ::DomainType DomainType
Definition: localmassmatrix.hh:69
DiscreteFunctionSpaceType ::RangeFieldType RangeFieldType
Definition: localmassmatrix.hh:40
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
const DiscreteFunctionSpaceType & space() const
Definition: localmassmatrix.hh:349
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
DiscreteFunctionSpaceType ::RangeType RangeType
Definition: localmassmatrix.hh:41
const bool affine_
Definition: localmassmatrix.hh:81
bool setup() const
setup and return affinity
Definition: localmassmatrix.hh:750
GridPartType ::GridType GridType
Definition: localmassmatrix.hh:57
Dune::FieldVector< ctype, dgNumDofs > DGVectorType
Definition: localmassmatrix.hh:48
MatrixType & getLocalInverseMassMatrixDefault(MassCaller &caller, const EntityType &entity, const Geometry &geo, const BasisFunctionSet &basisSet) const
Definition: localmassmatrix.hh:177
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
VolumeQuadrature VolumeQuadratureType
Definition: localmassmatrix.hh:61
std::vector< MassMatrixStorageType > LocalInverseMassMatrixStorageType
Definition: localmassmatrix.hh:96
DiscreteFunctionSpaceType ::IndexSetType IndexSetType
Definition: localmassmatrix.hh:52
MatrixPairType & getLocalInverseMassMatrix(const EntityType &entity, const Geometry &geo, const BasisFunctionSet &basisSet, int numBasisFct) const
Definition: localmassmatrix.hh:135
std::vector< RangeType > phiMass_
Definition: localmassmatrix.hh:92
LocalMassMatrixImplementation(const ThisType &other)
copy constructor
Definition: localmassmatrix.hh:239
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
IndexSetType ::IndexType IndexType
Definition: localmassmatrix.hh:53
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
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
DiscreteFunctionSpace DiscreteFunctionSpaceType
Definition: localmassmatrix.hh:38
MatrixType matrix_
Definition: localmassmatrix.hh:89
DiscreteFunctionSpaceType ::RangeFieldType ctype
Definition: localmassmatrix.hh:39
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
const VectorType & setupInverseDiagonal(const EntityType &entity, const Geometry &geo, const VectorType &refElemDiagonal, LocalMatrix &localMatrix) const
Definition: localmassmatrix.hh:612
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
DiscreteFunctionSpaceType ::BasisFunctionSetType BasisFunctionSetType
Definition: localmassmatrix.hh:55
void rightMultiplyInverseLocally(const EntityType &entity, const Geometry &geo, LocalMatrix &localMatrix) const
Definition: localmassmatrix.hh:639
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
EntityType ::Geometry Geometry
Definition: localmassmatrix.hh:59
DiscreteFunctionSpaceType ::GridPartType GridPartType
Definition: localmassmatrix.hh:50
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
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
DiscreteFunctionSpaceType ::EntityType EntityType
Definition: localmassmatrix.hh:58
int sequence_
Definition: localmassmatrix.hh:104
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
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
BaseType ::EntityType EntityType
Definition: localmassmatrix.hh:933
Interface class for basis function sets.
Definition: basisfunctionset/basisfunctionset.hh:31
std::size_t size() const
return size of basis function set
const EntityType & entity() const
return entity
ctype referenceVolume(const GeometryType &type) const
return volume of reference element for geometry of type type
Definition: allgeomtypes.hh:72
FieldVector< ctype, dim > DomainType
type of domain vector
Definition: allgeomtypes.hh:51
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