3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
matrixconverter.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_MATRIX_CONVERTER
25#define DUMUX_MATRIX_CONVERTER
26
27#include <cmath>
28#include <utility>
29#include <dune/common/indices.hh>
30#include <dune/common/hybridutilities.hh>
31#include <dune/istl/bvector.hh>
32#include <dune/istl/bcrsmatrix.hh>
33#include <dune/istl/matrixindexset.hh>
34
36
37namespace Dumux {
38
45template <class MultiTypeBlockMatrix, class Scalar=double>
47{
48 using MatrixBlock = typename Dune::FieldMatrix<Scalar, 1, 1>;
49 using BCRSMatrix = typename Dune::BCRSMatrix<MatrixBlock>;
50
51public:
52
58 static auto multiTypeToBCRSMatrix(const MultiTypeBlockMatrix &A)
59 {
60 // get the size for the converted matrix
61 const auto numRows = getNumRows_(A);
62
63 // create an empty BCRS matrix with 1x1 blocks
64 auto M = BCRSMatrix(numRows, numRows, BCRSMatrix::random);
65
66 // set the occupation pattern and copy the values
67 setOccupationPattern_(M, A);
68 copyValues_(M, A);
69
70 return M;
71 }
72
73private:
74
81 static void setOccupationPattern_(BCRSMatrix& M, const MultiTypeBlockMatrix& A)
82 {
83 // prepare the occupation pattern
84 const auto numRows = M.N();
85 Dune::MatrixIndexSet occupationPattern;
86 occupationPattern.resize(numRows, numRows);
87
88 // lambda function to fill the occupation pattern
89 auto addIndices = [&](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol)
90 {
91 using std::abs;
92 static const Scalar eps = getParam<Scalar>("MatrixConverter.DeletePatternEntriesBelowAbsThreshold", -1.0);
93
94 using BlockType = typename std::decay_t<decltype(subMatrix)>::block_type;
95 const auto blockSizeI = BlockType::rows;
96 const auto blockSizeJ = BlockType::cols;
97 for(auto row = subMatrix.begin(); row != subMatrix.end(); ++row)
98 for(auto col = row->begin(); col != row->end(); ++col)
99 for(std::size_t i = 0; i < blockSizeI; ++i)
100 for(std::size_t j = 0; j < blockSizeJ; ++j)
101 if(abs(subMatrix[row.index()][col.index()][i][j]) > eps)
102 occupationPattern.add(startRow + row.index()*blockSizeI + i, startCol + col.index()*blockSizeJ + j);
103 };
104
105 // fill the pattern
106 using namespace Dune::Hybrid;
107 std::size_t rowIndex = 0;
108 forEach(std::make_index_sequence<MultiTypeBlockMatrix::N()>(), [&A, &addIndices, &rowIndex, numRows](const auto i)
109 {
110 std::size_t colIndex = 0;
111 forEach(A[i], [&](const auto& subMatrix)
112 {
113 addIndices(subMatrix, rowIndex, colIndex);
114
115 using SubBlockType = typename std::decay_t<decltype(subMatrix)>::block_type;
116
117 colIndex += SubBlockType::cols * subMatrix.M();
118
119 // if we have arrived at the right side of the matrix, increase the row index
120 if(colIndex == numRows)
121 rowIndex += SubBlockType::rows * subMatrix.N();
122 });
123 });
124
125 occupationPattern.exportIdx(M);
126 }
127
134 static void copyValues_(BCRSMatrix& M, const MultiTypeBlockMatrix& A)
135 {
136 // get number of rows
137 const auto numRows = M.N();
138
139 // lambda function to copy the values
140 auto copyValues = [&](const auto& subMatrix, const std::size_t startRow, const std::size_t startCol)
141 {
142 using std::abs;
143 static const Scalar eps = getParam<Scalar>("MatrixConverter.DeletePatternEntriesBelowAbsThreshold", -1.0);
144
145 using BlockType = typename std::decay_t<decltype(subMatrix)>::block_type;
146 const auto blockSizeI = BlockType::rows;
147 const auto blockSizeJ = BlockType::cols;
148 for (auto row = subMatrix.begin(); row != subMatrix.end(); ++row)
149 for (auto col = row->begin(); col != row->end(); ++col)
150 for (std::size_t i = 0; i < blockSizeI; ++i)
151 for (std::size_t j = 0; j < blockSizeJ; ++j)
152 if(abs(subMatrix[row.index()][col.index()][i][j]) > eps)
153 M[startRow + row.index()*blockSizeI + i][startCol + col.index()*blockSizeJ + j] = subMatrix[row.index()][col.index()][i][j];
154
155 };
156
157 using namespace Dune::Hybrid;
158 std::size_t rowIndex = 0;
159 forEach(std::make_index_sequence<MultiTypeBlockMatrix::N()>(), [&A, &copyValues, &rowIndex, numRows](const auto i)
160 {
161 std::size_t colIndex = 0;
162 forEach(A[i], [&](const auto& subMatrix)
163 {
164 copyValues(subMatrix, rowIndex, colIndex);
165
166 using SubBlockType = typename std::decay_t<decltype(subMatrix)>::block_type;
167
168 colIndex += SubBlockType::cols * subMatrix.M();
169
170 // if we have arrived at the right side of the matrix, increase the row index
171 if(colIndex == numRows)
172 rowIndex += SubBlockType::rows * subMatrix.N();
173 });
174 });
175 }
176
182 static std::size_t getNumRows_(const MultiTypeBlockMatrix& A)
183 {
184 // iterate over the first row of the multitype blockmatrix
185 std::size_t numRows = 0;
186 Dune::Hybrid::forEach(Dune::Hybrid::elementAt(A, Dune::Indices::_0), [&numRows](const auto& subMatrixInFirstRow)
187 {
188 // the number of cols of the individual submatrice's block equals the respective number of equations.
189 const auto numEq = std::decay_t<decltype(subMatrixInFirstRow)>::block_type::cols;
190 numRows += numEq * subMatrixInFirstRow.M();
191 });
192
193 return numRows;
194 }
195
196};
197
202template<class MultiTypeBlockVector, class Scalar=double>
204{
205 using VectorBlock = typename Dune::FieldVector<Scalar, 1>;
206 using BlockVector = typename Dune::BlockVector<VectorBlock>;
207
208public:
209
215 static auto multiTypeToBlockVector(const MultiTypeBlockVector& b)
216 {
217 BlockVector bTmp;
218 bTmp.resize(b.dim());
219
220 std::size_t startIndex = 0;
221 Dune::Hybrid::forEach(b, [&](const auto& subVector)
222 {
223 const auto numEq = std::decay_t<decltype(subVector)>::block_type::size();
224
225 for(std::size_t i = 0; i < subVector.size(); ++i)
226 for(std::size_t j = 0; j < numEq; ++j)
227 bTmp[startIndex + i*numEq + j] = subVector[i][j];
228
229 startIndex += numEq*subVector.size();
230 });
231
232 return bTmp;
233 }
234
241 static void retrieveValues(MultiTypeBlockVector& x, const BlockVector& y)
242 {
243 std::size_t startIndex = 0;
244 Dune::Hybrid::forEach(x, [&](auto& subVector)
245 {
246 const auto numEq = std::decay_t<decltype(subVector)>::block_type::size();
247
248 for(std::size_t i = 0; i < subVector.size(); ++i)
249 for(std::size_t j = 0; j < numEq; ++j)
250 subVector[i][j] = y[startIndex + i*numEq + j];
251
252 startIndex += numEq*subVector.size();
253 });
254 }
255};
256
257} // end namespace Dumux
258
259#endif
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
Definition: adapt.hh:29
Definition: matrix.hh:32
A helper classe that converts a Dune::MultiTypeBlockMatrix into a plain Dune::BCRSMatrix TODO: allow ...
Definition: matrixconverter.hh:47
static auto multiTypeToBCRSMatrix(const MultiTypeBlockMatrix &A)
Converts the matrix to a type the IterativeSolverBackend can handle.
Definition: matrixconverter.hh:58
A helper classe that converts a Dune::MultiTypeBlockVector into a plain Dune::BlockVector and transfe...
Definition: matrixconverter.hh:204
static void retrieveValues(MultiTypeBlockVector &x, const BlockVector &y)
Copys the entries of a Dune::BlockVector to a Dune::MultiTypeBlockVector.
Definition: matrixconverter.hh:241
static auto multiTypeToBlockVector(const MultiTypeBlockVector &b)
Converts a Dune::MultiTypeBlockVector to a plain 1x1 Dune::BlockVector.
Definition: matrixconverter.hh:215