Add a method that combines a vector of matrices to one big blockmatrix

dune/fufem/blockmatrixview.hh

@@ -43,6 +43,7 @@ public:
         return nCols_;
     }
 
+    static void setupBlockMatrix(const std::vector<const MatrixType*>& submat, MatrixType& mat);
 protected:
     std::vector<Dune::array<int,2> > offsets_;
 
@@ -50,4 +51,45 @@ protected:
 
 };
 
+template <class MatrixType>
+void BlockMatrixView<MatrixType>::setupBlockMatrix(const std::vector<const MatrixType*>& submat, MatrixType& mat) {
+
+    BlockMatrixView<MatrixType> blockView(submat);
+    Dune::MatrixIndexSet indexSet(blockView.nRows(),blockView.nCols());
+
+    for (size_t i=0; i<submat.size(); i++) {
+
+        const auto row = submat[i]->begin();
+        const auto rowEnd = submat[i]->end();
+
+        for (; row != rowEnd; row++) {
+
+            const auto col = row.begin();
+            const auto colEnd = row.end();
+
+            for (; col != colEnd; col++)
+                indexSet.add(blockView.row(i,row.index()),
+                             blockView.col(i,col.index()));
+        }
+    }
+
+    indexSet.exportIdx(mat);
+    mat = 0;
+
+    for (size_t i=0; i<submat.size(); i++) {
+
+        const auto row = submat[i]->begin();
+        const auto rowEnd = submat[i]->end();
+
+        for (; row != rowEnd; row++) {
+
+            const auto col = row.begin();
+            const auto colEnd = row.end();
+
+            for (; col != colEnd; col++)
+                mat[blockView.row(i,row.index())][blockView.col(i,col.index())] = *col;
+        }
+    }
+}
+
 #endif