diff --git a/dune/fem-dg/misc/Makefile.am b/dune/fem-dg/misc/Makefile.am
index 5417c1d5e48169fc732d65c20915e99b6785dbf2..969ed3789ad33df58d641e8d822cf5dade01bad2 100644
--- a/dune/fem-dg/misc/Makefile.am
+++ b/dune/fem-dg/misc/Makefile.am
@@ -1,4 +1,4 @@
miscdir = $(includedir)/dune/fem-dg/misc
-misc_HEADERS = diagnostics.hh cons2prim.hh streams.hh
+misc_HEADERS = diagnostics.hh cons2prim.hh crs.hh streams.hh
include $(top_srcdir)/am/global-rules
diff --git a/dune/fem-dg/misc/crs.hh b/dune/fem-dg/misc/crs.hh
new file mode 100644
index 0000000000000000000000000000000000000000..44ce49d4a8a7363a3c4ea9bbef3b4cdd08bf41bc
--- /dev/null
+++ b/dune/fem-dg/misc/crs.hh
@@ -0,0 +1,148 @@
+#ifndef DUNE_CRSMATRIX_HH
+#define DUNE_CRSMATRIX_HH
+
+#include <cmath>
+#include <cstddef>
+#include <iostream>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/densematrix.hh>
+#include <dune/common/precision.hh>
+#include <dune/common/static_assert.hh>
+#include <dune/common/std/constexpr.hh>
+
+namespace Dune
+{
+
+ namespace Fem {
+
+
+ template<class K, int ROWS, int COLS>
+ class BlockCRSMatrix
+ {
+ std::vector< double > data_;
+ std::vector< int > rows_;
+ std::vector< int > cols_;
+ public:
+
+ //! export size
+ enum {
+ //! The number of rows.
+ rows = ROWS,
+ //! The number of columns.
+ cols = COLS
+ };
+
+ typedef K Field;
+ typedef K value_type ;
+
+ typedef std::vector< value_type > row_type ;
+
+ //===== constructors
+ /** \brief Default constructor
+ */
+ BlockCRSMatrix ()
+ : data_(), rows_(rows + 1, int(0) ), cols_() {}
+
+ BlockCRSMatrix ( const BlockCRSMatrix& other )
+ : data_( other.data_ ), rows_( other.rows_ ), cols_( other.cols_ ) {}
+
+ template <class Matrix>
+ void set( const Matrix& matrix )
+ {
+ data_.clear();
+ cols_.clear();
+
+ assert( matrix.size() == rows );
+
+ int count = 0;
+ for( int row=0; row<rows; ++row )
+ for( int col=0; col<cols; ++col )
+ if( std::abs( matrix[ row ][ col ] ) > 0.0 )
+ ++ count;
+
+ data_.resize( count, Field( 0 ) );
+ cols_.resize( count, int( 0 ) );
+
+ count = 0 ;
+ for( int row=0; row<rows; ++row )
+ {
+ rows_[ row ] = count ;
+ for( int col=0; col<cols; ++col )
+ {
+ const Field val = matrix[ row ][ col ];
+ if( std::abs( val ) > 0.0 )
+ {
+ assert( count < int( data_.size() ) );
+ cols_[ count ] = col ;
+ data_[ count ] = val ;
+ ++ count ;
+ }
+ }
+ }
+ rows_[ rows ] = count ;
+ }
+
+ const int N () const { return rows; }
+ const int M () const { return cols; }
+
+ template <class X, class Y>
+ void mv ( const X& x, Y& y ) const
+ {
+ for( int row=0; row<rows; ++row )
+ {
+ Field sum = 0;
+ for( int c = rows_[ row ]; c < rows_[ row+1 ]; ++ c )
+ {
+ sum += data_[ c ] * x[ cols_[ c ] ];
+ }
+ y[ row ] = sum;
+ }
+ }
+
+ template <class X, class Y>
+ void mvb( const int blockSize, const X& x, Y& y ) const
+ {
+ // clear result
+ y.clear();
+
+ for( int row=0; row<rows; ++row )
+ {
+ for( int c = rows_[ row ]; c < rows_[ row+1 ]; ++ c )
+ {
+ Field matrix = data_[ c ];
+ for( int r=0, ir = row * blockSize, jr = cols_[ c ] * blockSize ;
+ r<blockSize; ++ r, ++ ir, ++ jr )
+ {
+ y[ ir ] += matrix * x[ jr ] ;
+ }
+ }
+ }
+ }
+
+ //! Multiplies M from the right to this matrix
+ template < class M >
+ BlockCRSMatrix& rightmultiply (const M& )
+ {
+ abort();
+ return *this;
+ }
+
+ template < class M >
+ BlockCRSMatrix& leftmultiply (const M& )
+ {
+ abort();
+ return *this;
+ }
+
+ void invert()
+ {
+ abort();
+ }
+ };
+
+ } // end namespace Fem
+
+} // end namespace
+#endif
diff --git a/dune/fem-dg/pass/Makefile.am b/dune/fem-dg/pass/Makefile.am
index 51e4c605c86dac6189d207385c91aedb58d9d61e..ab5dde7ce7712fd62753fa0aad17d7791e0f18a5 100644
--- a/dune/fem-dg/pass/Makefile.am
+++ b/dune/fem-dg/pass/Makefile.am
@@ -1,3 +1,4 @@
passdir = $(includedir)/dune/fem-dg/pass
-pass_HEADERS = assembleddiffusionpass.hh dgmodelcaller.hh dgpass.hh ellipticmodelcaller.hh threadhandle.hh threadpass.hh
+pass_HEADERS = assembleddiffusionpass.hh dgmodelcaller.hh dgmasspass.hh \
+ dgpass.hh ellipticmodelcaller.hh threadhandle.hh threadpass.hh
include $(top_srcdir)/am/global-rules
diff --git a/dune/fem-dg/pass/dgmasspass.hh b/dune/fem-dg/pass/dgmasspass.hh
index 41b078985514289f01fc467378c989aadb61e16f..97ef309ccdc4058563d91421c42d794753c1f18f 100644
--- a/dune/fem-dg/pass/dgmasspass.hh
+++ b/dune/fem-dg/pass/dgmasspass.hh
@@ -12,6 +12,7 @@
#include <dune/fem/pass/common/local.hh>
#include <dune/fem/quadrature/cachingquadrature.hh>
#include <dune/fem/pass/dginversemass.hh>
+#include <dune/fem-dg/misc/crs.hh>
namespace Dune
{
@@ -69,7 +70,8 @@ namespace Dune
private:
typedef typename ScalarDiscreteFunctionSpaceType :: RangeFieldType ctype;
- typedef Dune::FieldMatrix< ctype, numScalarBasis, numScalarBasis > LocalMatrixType ;
+ typedef Dune::FieldMatrix< ctype, numScalarBasis, numScalarBasis > LocalMatrixType ;
+ typedef Dune::Fem::BlockCRSMatrix< ctype, numScalarBasis, numScalarBasis > SparseLocalMatrixType ;
typedef Fem::CachingQuadrature< GridPartType, 0 > VolumeQuadratureType;
@@ -88,9 +90,9 @@ namespace Dune
struct VectorEntry
{
- LocalMatrixType matrix_ ;
+ SparseLocalMatrixType matrix_ ;
GlobalCoordinate center_ ;
- VectorEntry () : matrix_( 0 ), center_( std::numeric_limits< double > :: max() )
+ VectorEntry () : matrix_(), center_( std::numeric_limits< double > :: max() )
{}
};
@@ -223,13 +225,15 @@ namespace Dune
protected:
template < class ConstLocalFunction, class LocalFunctionType >
void multiplyBlock( const int dimRange,
- const LocalMatrixType& matrix,
+ const SparseLocalMatrixType& matrix,
const ConstLocalFunction& arg,
LocalFunctionType& dest ) const
{
// clear destination
- dest.clear();
+ // dest.clear();
+ matrix.mvb( dimRange, arg, dest );
+ /*
static const int rows = LocalMatrixType::rows;
static const int cols = LocalMatrixType::cols;
@@ -247,6 +251,7 @@ namespace Dune
}
}
}
+ */
}
void setup() const
@@ -306,10 +311,10 @@ namespace Dune
const EntityType& entity,
const Geometry& geo,
const BasisFunctionSetType& set,
- LocalMatrixType& matrix ) const
+ SparseLocalMatrixType& storedMatrix ) const
{
// clear matrix
- matrix = 0;
+ LocalMatrixType matrix( 0 );
// vector holding basis function evaluation
Dune::array< ScalarRangeType, numScalarBasis > phi ;
@@ -356,6 +361,9 @@ namespace Dune
// if we are looking for the inverse, then invert matrix
if( inverse )
matrix.invert();
+
+ // store matrix
+ storedMatrix.set( matrix );
}
protected: