Add a blockmatrix view that computes global indices from block-indices

This class also provides a static method to combine submatrices into
one block matrix by copying the entries.

dune/matrix-vector/CMakeLists.txt

@@ -5,6 +5,7 @@ install(FILES
   addtodiagonal.hh
   axpy.hh
   axy.hh
+  blockmatrixview.hh
   componentwisematrixmap.hh
   crossproduct.hh
   genericvectortools.hh

dune/matrix-vector/blockmatrixview.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=4 sw=2 et sts=2:
+#ifndef DUNE_MATRIX_VECTOR_BLOCK_MATRIX_VIEW_HH
+#define DUNE_MATRIX_VECTOR_BLOCK_MATRIX_VIEW_HH
+
+#include <array>
+#include <vector>
+
+#include <dune/istl/matrixindexset.hh>
+
+namespace Dune {
+namespace MatrixVector {
+
+/** \brief Class that provides global indices for a block matrix given local indices of the blocks. Furthermore, it provides a method to construct a block matrix from given sub matrices.
+ *
+ *  This class provides a method to combine matrices \f A_1,\ldots,A_n\f into one single block matrix
+ *  \f B = \begin{pmatrix} A_1 &  0  & \ldots & 0 \\
+ *                          0 &  A_2 & \ldots & 0 \\
+ *                                \ldots          \\
+ *                          0 & \ldots & 0 & A_n \end{pmatrix} \f,
+ *  by copying the corresponding entries.
+ *  A BlockMatrixView object provides methods that return global indices for \f B\f given a local index of \f A_i\f and the block index \f i\f.
+ *
+ *  \tparam MatrixType The matrix type of the matrix blocks.
+ */
+template <class MatrixType>
+class BlockMatrixView
+{
+public:
+  //! Construct block view from vector of matrices corr. to the blocks
+  BlockMatrixView(const std::vector<const MatrixType*>& submat) {
+
+    offsets_.resize(submat.size());
+    offsets_[0] = {0, 0};
+
+    for (size_t i=0; i<submat.size()-1; i++)
+      offsets_[i+1] = {offsets_[i][0] + submat[i]->N(),
+                        offsets_[i][1] + submat[i]->M()};
+
+    nRows_ = offsets_.back()[0] + submat.back()->N();
+    nCols_ = offsets_.back()[1] + submat.back()->M();
+  }
+
+  //! Get global row index for given block and a local row index
+  size_t row(size_t block, size_t index) const {
+    return offsets_[block][0] + index;
+  }
+
+  //! Get global column index given a block and a local column index
+  size_t col(size_t block, size_t index) const {
+    return offsets_[block][1] + index;
+  }
+
+  //! Return the total number of rows
+  size_t nRows() const {
+    return nRows_;
+  }
+
+  //! Return the total number of columns
+  size_t nCols() const {
+    return nCols_;
+  }
+
+  /** \brief Combine submatrices in a big matrix. */
+  static void setupBlockMatrix(const std::vector<MatrixType>& submat, MatrixType& mat) {
+    std::vector<const MatrixType*> dummy(submat.size());
+    for (size_t i=0; i<dummy.size(); i++)
+      dummy[i] = &submat[i];
+    setupBlockMatrix(dummy,mat);
+  }
+
+  /** \brief Combine submatrices in a big matrix. */
+  static void setupBlockMatrix(const std::vector<const MatrixType*>& submat, MatrixType& mat)
+  {
+    BlockMatrixView<MatrixType> blockView(submat);
+    MatrixIndexSet indexSet(blockView.nRows(), blockView.nCols());
+
+    for (size_t i=0; i<submat.size(); i++)
+      indexSet.import(*submat[i], blockView.row(i,0), blockView.col(i,0));
+
+    indexSet.exportIdx(mat);
+
+    for (size_t i=0; i<submat.size(); i++) {
+
+      auto row = submat[i]->begin();
+      auto rowEnd = submat[i]->end();
+
+      for (; row != rowEnd; row++) {
+
+        auto col = row->begin();
+        auto colEnd = row->end();
+
+        for (; col != colEnd; col++)
+          mat[blockView.row(i,row.index())][blockView.col(i,col.index())] = *col;
+      }
+    }
+  }
+
+protected:
+  //! Offsets for the different blocks
+  std::vector<std::array<size_t, 2> > offsets_;
+  //! Total numbers of rows and columns
+  size_t nRows_, nCols_;
+};
+
+} // end namespace MatrixVector
+} // end namespace Dune
+
+#endif