Skip to content
Snippets Groups Projects
Commit 20e7cd9a authored by graeser's avatar graeser Committed by graeser
Browse files

Sync with dune-fufem

[[Imported from SVN: r6058]]
parent f542de4a
No related branches found
No related tags found
No related merge requests found
...@@ -24,6 +24,18 @@ class StaticMatrix ...@@ -24,6 +24,18 @@ class StaticMatrix
typedef Matrix Type; typedef Matrix Type;
}; };
template<typename FieldType, int n>
struct Promote<Dune::FieldMatrix<FieldType,n,n>, Dune::DiagonalMatrix<FieldType,n> >
{
typedef Dune::FieldMatrix<FieldType,n,n> Type;
};
template<typename FieldType, int n>
struct Promote<Dune::DiagonalMatrix<FieldType,n>, Dune::FieldMatrix<FieldType,n,n> >
{
typedef Dune::FieldMatrix<FieldType,n,n> Type;
};
template<typename FieldType, int n> template<typename FieldType, int n>
struct Promote<Dune::DiagonalMatrix<FieldType,n>, Dune::ScaledIdentityMatrix<FieldType,n> > struct Promote<Dune::DiagonalMatrix<FieldType,n>, Dune::ScaledIdentityMatrix<FieldType,n> >
{ {
...@@ -47,7 +59,14 @@ class StaticMatrix ...@@ -47,7 +59,14 @@ class StaticMatrix
} }
template <typename FieldType, int n> template <typename FieldType, int n>
static void addToDiagonal(Dune::ScaledIdentityMatrix<FieldType,n>& x, const FieldType a) // static void addToDiagonal(Dune::ScaledIdentityMatrix<FieldType,n>& x, const FieldType a)
// the commented line should NOT be used as it may lead to ambiguities in which case the general method above will be used.
// an example is the line taken from massassembler.hh (line 83 at the time):
// StaticMatrix::addToDiagonal(localMatrix[i][i], values[i] * zi);
// where values[i] is of type FielVector<FieldType,1> and zi is a double. This call then does not exactly fit this specialization (without the implicit cast of FV<FT,1>)
// and hence some wild template voodoo decides which of the templates is to be taken - in this case with a gcc 4.4.5 that was the one above leading to a wrong factor n
// in the diagonal value
static void addToDiagonal(Dune::ScaledIdentityMatrix<FieldType,n>& x, const typename Dune::ScaledIdentityMatrix<FieldType,n>::field_type a)
{ {
x.scalar() += a; x.scalar() += a;
} }
...@@ -58,7 +77,7 @@ class StaticMatrix ...@@ -58,7 +77,7 @@ class StaticMatrix
template<class MatrixA, class MatrixB> template<class MatrixA, class MatrixB>
static void axpy(MatrixA& A, const typename MatrixA::field_type a, const MatrixB& B) static void axpy(MatrixA& A, const typename MatrixA::field_type a, const MatrixB& B)
{ {
for(int i=0; i<B.N(); ++i) for(size_t i=0; i<B.N(); ++i)
{ {
typename MatrixB::row_type::ConstIterator it=B[i].begin(); typename MatrixB::row_type::ConstIterator it=B[i].begin();
typename MatrixB::row_type::ConstIterator end=B[i].end(); typename MatrixB::row_type::ConstIterator end=B[i].end();
...@@ -76,13 +95,13 @@ class StaticMatrix ...@@ -76,13 +95,13 @@ class StaticMatrix
// } // }
template <typename FieldType, int n, class MatrixB> template <typename FieldType, int n, class MatrixB>
static void axpy(Dune::ScaledIdentityMatrix<FieldType,n>& x, const FieldType a, const MatrixB& y) static void axpy(Dune::ScaledIdentityMatrix<FieldType,n>& x, const typename Dune::ScaledIdentityMatrix<FieldType,n>::field_type a, const MatrixB& y)
{ {
x.scalar() += a * y[0][0]; x.scalar() += a * y[0][0];
} }
template <typename FieldType, int n, class MatrixB> template <typename FieldType, int n, class MatrixB>
static void axpy(Dune::DiagonalMatrix<FieldType,n>& x, const FieldType a, const MatrixB& y) static void axpy(Dune::DiagonalMatrix<FieldType,n>& x, const typename Dune::DiagonalMatrix<FieldType,n>::field_type a, const MatrixB& y)
{ {
for(int i=0; i<n; ++i) for(int i=0; i<n; ++i)
x.diagonal()[i] += a * y[i][i]; x.diagonal()[i] += a * y[i][i];
...@@ -95,17 +114,20 @@ class StaticMatrix ...@@ -95,17 +114,20 @@ class StaticMatrix
} }
// add transformed matrix X = A^T*Y*B ****************************************************** // add transformed matrix X = A^T*Y*B ******************************************************
template<class MatrixA, class MatrixB, class TransformationMatrix> template<class MatrixA, class MatrixB, class TransformationMatrix>
static void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2) static void addTransformedMatrix(MatrixA& A, const TransformationMatrix& T1, const MatrixB& B, const TransformationMatrix& T2)
{ {
for (int i=0; i<A.N(); ++i) for (size_t i=0; i<A.N(); ++i)
{ {
typename MatrixA::row_type::ConstIterator jIt=A[i].begin(); typename MatrixA::row_type::ConstIterator jIt=A[i].begin();
typename MatrixA::row_type::ConstIterator jEnd=A[i].end(); typename MatrixA::row_type::ConstIterator jEnd=A[i].end();
for(; jIt!=jEnd; ++jIt) for(; jIt!=jEnd; ++jIt)
{ {
for (int k=0; k<B.N(); ++k) for (size_t k=0; k<B.N(); ++k)
{ {
typename MatrixB::row_type::ConstIterator lIt=B[k].begin(); typename MatrixB::row_type::ConstIterator lIt=B[k].begin();
typename MatrixB::row_type::ConstIterator lEnd=B[k].end(); typename MatrixB::row_type::ConstIterator lEnd=B[k].end();
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment