diff --git a/dumux/implicit/staggered/newtoncontroller.hh b/dumux/implicit/staggered/newtoncontroller.hh index edb72ac64eac1940448e25f4fcde986658c358d1..1e4cf1dd4a867d24d564b3fdcde479b0a21f9ae1 100644 --- a/dumux/implicit/staggered/newtoncontroller.hh +++ b/dumux/implicit/staggered/newtoncontroller.hh @@ -31,6 +31,7 @@ #include #include +#include #include "newtonconvergencewriter.hh" namespace Dumux { @@ -53,6 +54,7 @@ namespace Properties * which allows the newton method to abort quicker if the solution is * way out of bounds. */ + template class StaggeredNewtonController : public NewtonController { @@ -101,116 +103,34 @@ public: if (this->numSteps_ == 0) this->initialResidual_ = b.two_norm(); - // copy the matrix and the vector to types the IterativeSolverBackend can handle - using MatrixBlock = typename Dune::FieldMatrix; - using SparseMatrix = typename Dune::BCRSMatrix; - - // get the new matrix sizes - std::size_t numRows = numEqCellCenter*A[cellCenterIdx][cellCenterIdx].N() + numEqFace*A[faceIdx][cellCenterIdx].N(); - std::size_t numCols = numEqCellCenter*A[cellCenterIdx][cellCenterIdx].M() + numEqFace*A[cellCenterIdx][faceIdx].M(); - // check matrix sizes assert(A[cellCenterIdx][cellCenterIdx].N() == A[cellCenterIdx][faceIdx].N()); assert(A[faceIdx][cellCenterIdx].N() == A[faceIdx][faceIdx].N()); - assert(numRows == numCols); // create the bcrs matrix the IterativeSolver backend can handle - auto M = SparseMatrix(numRows, numCols, SparseMatrix::random); - - // set the rowsizes - // A11 and A12 - for (auto row = A[cellCenterIdx][cellCenterIdx].begin(); row != A[cellCenterIdx][cellCenterIdx].end(); ++row) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - M.setrowsize(numEqCellCenter*row.index() + i, row->size()*numEqCellCenter); - for (auto row = A[cellCenterIdx][faceIdx].begin(); row != A[cellCenterIdx][faceIdx].end(); ++row) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - M.setrowsize(numEqCellCenter*row.index() + i, M.getrowsize(numEqCellCenter*row.index() + i) + row->size()*numEqFace); - // A21 and A22 - for (auto row = A[faceIdx][cellCenterIdx].begin(); row != A[faceIdx][cellCenterIdx].end(); ++row) - for (std::size_t i = 0; i < numEqFace; ++i) - M.setrowsize(numEqFace*row.index() + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter, row->size()*numEqCellCenter); - for (auto row = A[faceIdx][faceIdx].begin(); row != A[faceIdx][faceIdx].end(); ++row) - for (std::size_t i = 0; i < numEqFace; ++i) - M.setrowsize(numEqFace*row.index() + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter, M.getrowsize(numEqFace*row.index() + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter) + row->size()*numEqFace); - M.endrowsizes(); - - // set the indices - for (auto row = A[cellCenterIdx][cellCenterIdx].begin(); row != A[cellCenterIdx][cellCenterIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - M.addindex(row.index()*numEqCellCenter + i, col.index()*numEqCellCenter + j); - - for (auto row = A[cellCenterIdx][faceIdx].begin(); row != A[cellCenterIdx][faceIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - M.addindex(row.index()*numEqCellCenter + i, col.index()*numEqFace + j + A[cellCenterIdx][cellCenterIdx].M()*numEqCellCenter); - - for (auto row = A[faceIdx][cellCenterIdx].begin(); row != A[faceIdx][cellCenterIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqFace; ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - M.addindex(row.index()*numEqFace + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter, col.index()*numEqCellCenter + j); - - for (auto row = A[faceIdx][faceIdx].begin(); row != A[faceIdx][faceIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqFace; ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - M.addindex(row.index()*numEqFace + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter, col.index()*numEqFace + j + A[cellCenterIdx][cellCenterIdx].M()*numEqCellCenter); - M.endindices(); - - // copy values - for (auto row = A[cellCenterIdx][cellCenterIdx].begin(); row != A[cellCenterIdx][cellCenterIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - M[row.index()*numEqCellCenter + i][col.index()*numEqCellCenter + j] = A[cellCenterIdx][cellCenterIdx][row.index()][col.index()][i][j]; - - for (auto row = A[cellCenterIdx][faceIdx].begin(); row != A[cellCenterIdx][faceIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqCellCenter; ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - M[row.index()*numEqCellCenter + i][col.index()*numEqFace + j + A[cellCenterIdx][cellCenterIdx].M()*numEqCellCenter] = A[cellCenterIdx][faceIdx][row.index()][col.index()][i][j]; - - for (auto row = A[faceIdx][cellCenterIdx].begin(); row != A[faceIdx][cellCenterIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqFace; ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - M[row.index()*numEqFace + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter][col.index()*numEqCellCenter + j] = A[faceIdx][cellCenterIdx][row.index()][col.index()][i][j]; - - for (auto row = A[faceIdx][faceIdx].begin(); row != A[faceIdx][faceIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqFace; ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - M[row.index()*numEqFace + i + A[cellCenterIdx][cellCenterIdx].N()*numEqCellCenter][col.index()*numEqFace + j + A[cellCenterIdx][cellCenterIdx].M()*numEqCellCenter] = A[faceIdx][faceIdx][row.index()][col.index()][i][j]; + const auto M = MatrixConverter::multiTypeToBCRSMatrix(A); + + // get the new matrix sizes + const std::size_t numRows = M.N(); + assert(numRows == M.M()); // create the vector the IterativeSolver backend can handle + const auto bTmp = VectorConverter::multiTypeToBlockVector(b); + assert(bTmp.size() == numRows); + + // create a blockvector to which the linear solver writes the solution using VectorBlock = typename Dune::FieldVector; using BlockVector = typename Dune::BlockVector; - - BlockVector y, bTmp; + BlockVector y; y.resize(numRows); - bTmp.resize(numCols); - for (std::size_t i = 0; i < b[cellCenterIdx].N(); ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - bTmp[i*numEqCellCenter + j] = b[cellCenterIdx][i][j]; - for (std::size_t i = 0; i < b[faceIdx].N(); ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - bTmp[i*numEqFace + j + b[cellCenterIdx].N()*numEqCellCenter] = b[faceIdx][i][j]; // printmatrix(std::cout, M, "", ""); // solve - bool converged = this->linearSolver_.solve(M, y, bTmp); + const bool converged = this->linearSolver_.solve(M, y, bTmp); // copy back the result y into x - for (std::size_t i = 0; i < x[cellCenterIdx].N(); ++i) - for (std::size_t j = 0; j < numEqCellCenter; ++j) - x[cellCenterIdx][i][j] = y[i*numEqCellCenter + j]; - for (std::size_t i = 0; i < x[faceIdx].N(); ++i) - for (std::size_t j = 0; j < numEqFace; ++j) - x[faceIdx][i][j] = y[i*numEqFace + j + x[cellCenterIdx].N()*numEqCellCenter]; + VectorConverter::retrieveValues(x, y); if (!converged) DUNE_THROW(NumericalProblem, "Linear solver did not converge"); diff --git a/dumux/linear/matrixconverter.hh b/dumux/linear/matrixconverter.hh new file mode 100644 index 0000000000000000000000000000000000000000..0615754cd0968cd04a4b77018dec2f8e918e15bb --- /dev/null +++ b/dumux/linear/matrixconverter.hh @@ -0,0 +1,266 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \brief A helper classe that converts a Dune::MultiTypeBlockMatrix into a plain Dune::BCRSMatrix + * + */ +#ifndef DUMUX_MATRIX_CONVERTER +#define DUMUX_MATRIX_CONVERTER + +#include +#include +#include +#include +#include +#include + +namespace Dumux { + +/*! + * A helper classe that converts a Dune::MultiTypeBlockMatrix into a plain Dune::BCRSMatrix + * TODO: allow block sizes for BCRSMatrix other than 1x1 ? + * + */ +template +class MatrixConverter +{ + using MatrixBlock = typename Dune::FieldMatrix; + using BCRSMatrix = typename Dune::BCRSMatrix; + +public: + + /*! + * \brief Converts the matrix to a type the IterativeSolverBackend can handle + * + * \param M The converted matrix + * \param A The original multitype blockmatrix + */ + static auto multiTypeToBCRSMatrix(const MultiTypeBlockMatrix &A) + { + // get the size for the converted matrix + const auto numRows = getNumRows_(A); + + // create an empty BCRS matrix with 1x1 blocks + auto M = BCRSMatrix(numRows, numRows, BCRSMatrix::random); + + // set the occupation pattern and copy the values + setOccupationPattern_(M, A); + copyValues_(M, A); + + return M; + } + +private: + + /*! + * \brief Sets the occupation pattern and indices for the converted matrix + * + * \param M The converted matrix + * \param A The original multitype blockmatrix + */ + static void setOccupationPattern_(BCRSMatrix& M, const MultiTypeBlockMatrix& A) + { + // prepare the occupation pattern + const auto numRows = M.N(); + Dune::MatrixIndexSet occupationPattern; + occupationPattern.resize(numRows, numRows); + + // lambda function to fill the occupation pattern + auto addIndices = [&occupationPattern](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol) + { + using BlockType = typename std::decay_t::block_type; + const auto blockSizeI = BlockType::rows; + const auto blockSizeJ = BlockType::cols; + for(auto row = subMatrix.begin(); row != subMatrix.end(); ++row) + for(auto col = row->begin(); col != row->end(); ++col) + for(std::size_t i = 0; i < blockSizeI; ++i) + for(std::size_t j = 0; j < blockSizeJ; ++j) + occupationPattern.add(startRow + row.index()*blockSizeI + i, startCol + col.index()*blockSizeJ + j); + }; + + // fill the pattern + std::size_t rowIndex = 0; + Dune::Hybrid::forEach(A, [&addIndices, &rowIndex, &occupationPattern, numRows](const auto& rowOfMultiTypeMatrix) + { + std::size_t colIndex = 0; + Dune::Hybrid::forEach(rowOfMultiTypeMatrix, [&addIndices, &occupationPattern, &colIndex, &rowIndex, numRows](const auto& subMatrix) + { + addIndices(subMatrix, rowIndex, colIndex); + + using SubBlockType = typename std::decay_t::block_type; + + colIndex += SubBlockType::cols * subMatrix.M(); + + // if we have arrived at the right side of the matrix, increase the row index + if(colIndex == numRows) + rowIndex += SubBlockType::rows * subMatrix.N(); + }); + }); + + occupationPattern.exportIdx(M); + } + + /*! + * \brief Sets the occupation pattern (i.e. indices) for the converted matrix + * + * \param M The converted matrix + * \param A The original subMatrix of the multitype blockmatrix + */ + static void copyValues_(BCRSMatrix& M, const MultiTypeBlockMatrix& A) + { + // get number of rows + const auto numRows = M.N(); + + // lambda function to copy the values + auto copyValues = [&M](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol) + { + using BlockType = typename std::decay_t::block_type; + const auto blockSizeI = BlockType::rows; + const auto blockSizeJ = BlockType::cols; + for (auto row = subMatrix.begin(); row != subMatrix.end(); ++row) + for (auto col = row->begin(); col != row->end(); ++col) + for (std::size_t i = 0; i < blockSizeI; ++i) + for (std::size_t j = 0; j < blockSizeJ; ++j) + M[startRow + row.index()*blockSizeI + i][startCol + col.index()*blockSizeJ + j] = subMatrix[row.index()][col.index()][i][j]; + + }; + + std::size_t rowIndex = 0; + Dune::Hybrid::forEach(A, [©Values, &rowIndex, numRows](const auto& rowOfMultiTypeMatrix) + { + std::size_t colIndex = 0; + Dune::Hybrid::forEach(rowOfMultiTypeMatrix, [©Values, &colIndex, &rowIndex, numRows](const auto& subMatrix) + { + copyValues(subMatrix, rowIndex, colIndex); + + using SubBlockType = typename std::decay_t::block_type; + + colIndex += SubBlockType::cols * subMatrix.M(); + + // if we have arrived at the right side of the matrix, increase the row index + if(colIndex == numRows) + rowIndex += SubBlockType::rows * subMatrix.N(); + }); + }); + } + + /*! + * \brief Calculates the total number of rows (== number of cols) for the converted matrix, assuming a block size of 1x1 + * + * \param A The original multitype blockmatrix + */ + static std::size_t getNumRows_(const MultiTypeBlockMatrix& A) + { + // iterate over the first row of the multitype blockmatrix + std::size_t numRows = 0; + Dune::Hybrid::forEach(Dune::Hybrid::elementAt(A, Dune::Indices::_0), [&numRows](const auto& subMatrixInFirstRow) + { + // the number of cols of the individual submatrice's block equals the respective number of equations. + const auto numEq = std::decay_t::block_type::cols; + numRows += numEq * subMatrixInFirstRow.M(); + }); + + return numRows; + } + +}; + +/*! + * A helper classe that converts a Dune::MultiTypeBlockVector into a plain Dune::BlockVector and transfers back values + */ +template +class VectorConverter +{ + using VectorBlock = typename Dune::FieldVector; + using BlockVector = typename Dune::BlockVector; + +public: + + /*! + * \brief Converts a Dune::MultiTypeBlockVector to a plain 1x1 Dune::BlockVector + * + * \param A The original multitype blockvector + */ + static auto multiTypeToBlockVector(const MultiTypeBlockVector& b) + { + const auto size = getSize_(b); + + BlockVector bTmp; + bTmp.resize(size); + + std::size_t startIndex = 0; + Dune::Hybrid::forEach(b, [&bTmp, &startIndex](const auto& subVector) + { + const auto numEq = std::decay_t::block_type::size(); + + for(std::size_t i = 0; i < subVector.size(); ++i) + for(std::size_t j = 0; j < numEq; ++j) + bTmp[startIndex + i*numEq + j] = subVector[i][j]; + + startIndex += numEq*subVector.size(); + }); + + return bTmp; + } + + /*! + * \brief Copys the entries of a Dune::BlockVector to a Dune::MultiTypeBlockVector + * + * \param x The multitype blockvector where the values are copied to + * \param y The regular blockvector where the values are copied from + */ + static void retrieveValues(MultiTypeBlockVector& x, const BlockVector& y) + { + std::size_t startIndex = 0; + Dune::Hybrid::forEach(x, [&y, &startIndex](auto& subVector) + { + const auto numEq = std::decay_t::block_type::size(); + + for(std::size_t i = 0; i < subVector.size(); ++i) + for(std::size_t j = 0; j < numEq; ++j) + subVector[i][j] = y[startIndex + i*numEq + j]; + + startIndex += numEq*subVector.size(); + }); + } + +private: + + /*! + * \brief Returns the size of the expanded multitype block vector + * + * \param b The multitype blockvector + */ + static std::size_t getSize_(const MultiTypeBlockVector& b) + { + std::size_t size = 0; + Dune::Hybrid::forEach(b, [&size](const auto& subVector) + { + // the size of the individual vector blocks equals the respective number of equations. + const auto numEq = std::decay_t::block_type::size(); + size += numEq * subVector.size(); + }); + + return size; + } +}; +} + +#endif diff --git a/dumux/mixeddimension/newtoncontroller.hh b/dumux/mixeddimension/newtoncontroller.hh index e938fa248f5c44b9b9c3ad7bf84c0507cf383cd8..a95a7bd1f33ef4f268da8ac455e54dffd98d5c46 100644 --- a/dumux/mixeddimension/newtoncontroller.hh +++ b/dumux/mixeddimension/newtoncontroller.hh @@ -31,6 +31,7 @@ #include #include #include +#include namespace Dumux { @@ -346,110 +347,36 @@ public: using MatrixBlock = typename Dune::FieldMatrix; using SparseMatrix = typename Dune::BCRSMatrix; - // get the new matrix sizes - std::size_t numRows = numEqBulk*A[bulkIdx][bulkIdx].N() + numEqLowDim*A[lowDimIdx][bulkIdx].N(); - std::size_t numCols = numEqBulk*A[bulkIdx][bulkIdx].M() + numEqLowDim*A[bulkIdx][lowDimIdx].M(); - // check matrix sizes assert(A[bulkIdx][bulkIdx].N() == A[bulkIdx][lowDimIdx].N()); assert(A[lowDimIdx][bulkIdx].N() == A[lowDimIdx][lowDimIdx].N()); - assert(numRows == numCols); // create the bcrs matrix the IterativeSolver backend can handle - auto M = SparseMatrix(numRows, numCols, SparseMatrix::random); - - // set the rowsizes - // A11 and A12 - for (auto row = A[bulkIdx][bulkIdx].begin(); row != A[bulkIdx][bulkIdx].end(); ++row) - for (std::size_t i = 0; i < numEqBulk; ++i) - M.setrowsize(numEqBulk*row.index() + i, row->size()*numEqBulk); - for (auto row = A[bulkIdx][lowDimIdx].begin(); row != A[bulkIdx][lowDimIdx].end(); ++row) - for (std::size_t i = 0; i < numEqBulk; ++i) - M.setrowsize(numEqBulk*row.index() + i, M.getrowsize(numEqBulk*row.index() + i) + row->size()*numEqLowDim); - // A21 and A22 - for (auto row = A[lowDimIdx][bulkIdx].begin(); row != A[lowDimIdx][bulkIdx].end(); ++row) - for (std::size_t i = 0; i < numEqLowDim; ++i) - M.setrowsize(numEqLowDim*row.index() + i + A[bulkIdx][bulkIdx].N()*numEqBulk, row->size()*numEqBulk); - for (auto row = A[lowDimIdx][lowDimIdx].begin(); row != A[lowDimIdx][lowDimIdx].end(); ++row) - for (std::size_t i = 0; i < numEqLowDim; ++i) - M.setrowsize(numEqLowDim*row.index() + i + A[bulkIdx][bulkIdx].N()*numEqBulk, M.getrowsize(numEqLowDim*row.index() + i + A[bulkIdx][bulkIdx].N()*numEqBulk) + row->size()*numEqLowDim); - M.endrowsizes(); - - // set the indices - for (auto row = A[bulkIdx][bulkIdx].begin(); row != A[bulkIdx][bulkIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqBulk; ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - M.addindex(row.index()*numEqBulk + i, col.index()*numEqBulk + j); - - for (auto row = A[bulkIdx][lowDimIdx].begin(); row != A[bulkIdx][lowDimIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqBulk; ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - M.addindex(row.index()*numEqBulk + i, col.index()*numEqLowDim + j + A[bulkIdx][bulkIdx].M()*numEqBulk); - - for (auto row = A[lowDimIdx][bulkIdx].begin(); row != A[lowDimIdx][bulkIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqLowDim; ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - M.addindex(row.index()*numEqLowDim + i + A[bulkIdx][bulkIdx].N()*numEqBulk, col.index()*numEqBulk + j); - - for (auto row = A[lowDimIdx][lowDimIdx].begin(); row != A[lowDimIdx][lowDimIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqLowDim; ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - M.addindex(row.index()*numEqLowDim + i + A[bulkIdx][bulkIdx].N()*numEqBulk, col.index()*numEqLowDim + j + A[bulkIdx][bulkIdx].M()*numEqBulk); - M.endindices(); - - // copy values - for (auto row = A[bulkIdx][bulkIdx].begin(); row != A[bulkIdx][bulkIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqBulk; ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - M[row.index()*numEqBulk + i][col.index()*numEqBulk + j] = A[bulkIdx][bulkIdx][row.index()][col.index()][i][j]; - - for (auto row = A[bulkIdx][lowDimIdx].begin(); row != A[bulkIdx][lowDimIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqBulk; ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - M[row.index()*numEqBulk + i][col.index()*numEqLowDim + j + A[bulkIdx][bulkIdx].M()*numEqBulk] = A[bulkIdx][lowDimIdx][row.index()][col.index()][i][j]; - - for (auto row = A[lowDimIdx][bulkIdx].begin(); row != A[lowDimIdx][bulkIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqLowDim; ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - M[row.index()*numEqLowDim + i + A[bulkIdx][bulkIdx].N()*numEqBulk][col.index()*numEqBulk + j] = A[lowDimIdx][bulkIdx][row.index()][col.index()][i][j]; - - for (auto row = A[lowDimIdx][lowDimIdx].begin(); row != A[lowDimIdx][lowDimIdx].end(); ++row) - for (auto col = row->begin(); col != row->end(); ++col) - for (std::size_t i = 0; i < numEqLowDim; ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - M[row.index()*numEqLowDim + i + A[bulkIdx][bulkIdx].N()*numEqBulk][col.index()*numEqLowDim + j + A[bulkIdx][bulkIdx].M()*numEqBulk] = A[lowDimIdx][lowDimIdx][row.index()][col.index()][i][j]; + const auto M = MatrixConverter::multiTypeToBCRSMatrix(A); + + // get the new matrix sizes + const std::size_t numRows = M.N(); + assert(numRows == M.M()); + + // create the vector the IterativeSolver backend can handle + const auto bTmp = VectorConverter::multiTypeToBlockVector(b); + assert(bTmp.size() == numRows); // create the vector the IterativeSolver backend can handle using VectorBlock = typename Dune::FieldVector; using BlockVector = typename Dune::BlockVector; - BlockVector y, bTmp; + // create a blockvector to which the linear solver writes the solution + using VectorBlock = typename Dune::FieldVector; + using BlockVector = typename Dune::BlockVector; + BlockVector y; y.resize(numRows); - bTmp.resize(numCols); - for (std::size_t i = 0; i < b[bulkIdx].N(); ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - bTmp[i*numEqBulk + j] = b[bulkIdx][i][j]; - for (std::size_t i = 0; i < b[lowDimIdx].N(); ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - bTmp[i*numEqLowDim + j + b[bulkIdx].N()*numEqBulk] = b[lowDimIdx][i][j]; // solve bool converged = linearSolver_.solve(M, y, bTmp); // copy back the result y into x - for (std::size_t i = 0; i < x[bulkIdx].N(); ++i) - for (std::size_t j = 0; j < numEqBulk; ++j) - x[bulkIdx][i][j] = y[i*numEqBulk + j]; - for (std::size_t i = 0; i < x[lowDimIdx].N(); ++i) - for (std::size_t j = 0; j < numEqLowDim; ++j) - x[lowDimIdx][i][j] = y[i*numEqLowDim + j + x[bulkIdx].N()*numEqBulk]; + VectorConverter::retrieveValues(x, y); if (!converged) DUNE_THROW(NumericalProblem, "Linear solver did not converge");