Syncronize with dune-fufem

[[Imported from SVN: r6115]]

dune/solvers/common/staticmatrixtools.hh

@@ -129,15 +129,47 @@ class StaticMatrix
                 {
                     for (size_t k=0; k<B.N(); ++k)
                     {
-                        typename MatrixB::row_type::ConstIterator lIt=B[k].begin();
-                        typename MatrixB::row_type::ConstIterator lEnd=B[k].end();
-                        for(; lIt!=lEnd; ++lIt)
-                            A[i][jIt.index()] += T1[k][i] * B[k][lIt.index()] * T2[lIt.index()][jIt.index()];
+                        if (T1[k][i]!=0)
+                        {
+                            typename MatrixB::row_type::ConstIterator lIt=B[k].begin();
+                            typename MatrixB::row_type::ConstIterator lEnd=B[k].end();
+                            for(; lIt!=lEnd; ++lIt)
+                                A[i][jIt.index()] += T1[k][i] * B[k][lIt.index()] * T2[lIt.index()][jIt.index()];
+                        }
                     }
                 }
             }
         }
 
+        template<class K1, class K2, class K3, int n, int m>
+        static void addTransformedMatrix(
+            typename Dune::FieldMatrix<K1, m, m>& A,
+            const typename Dune::FieldMatrix<K2, n, m>& T1,
+            const typename Dune::FieldMatrix<K3, n, n>& B,
+            const typename Dune::FieldMatrix<K2, n, m>& T2)
+        {
+            typename Dune::FieldMatrix<K1, m, n> T2transposedB;
+            T2transposedB = 0;
+            for (size_t i=0; i<A.N(); ++i)
+            {
+                for (size_t k=0; k<B.N(); ++k)
+                {
+                    if (T1[k][i]!=0)
+                        for (size_t l=0; l<B.M(); ++l)
+                            T2transposedB[i][l] += T1[k][i] * B[k][l];
+                }
+            }
+            for (size_t k=0; k<T2.N(); ++k)
+            {
+                for (size_t l=0; l<T2.M(); ++l)
+                {
+                    if (T2[k][l]!=0)
+                        for (size_t i=0; i<A.N(); ++i)
+                            A[i][l] += T2transposedB[i][k] * T2[k][l];
+                }
+            }
+        }
+
         template<class MatrixA, class MatrixB>
         static void addTransformedMatrix(MatrixA& X, const double& A, const MatrixB& Y, const double& B)
         {