matrixblock.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 This file is part of the Open Porous Media project (OPM).
5
6 OPM 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 2 of the License, or
9 (at your option) any later version.
10
11 OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
23#ifndef EWOMS_MATRIX_BLOCK_HH
24#define EWOMS_MATRIX_BLOCK_HH
25
26#include <dune/common/dynmatrix.hh>
27#include <dune/common/fmatrix.hh>
28#include <dune/common/typetraits.hh>
29
30#include <dune/istl/superlu.hh>
31#include <dune/istl/umfpack.hh>
32#include <dune/istl/matrixutils.hh>
33
34#include <opm/common/Exceptions.hpp>
35
36#include <limits>
37
38namespace Opm {
39namespace detail {
40
41template <typename K, int m, int n>
42static inline void invertMatrix(Dune::FieldMatrix<K,m,n>& matrix)
43{
44 matrix.invert();
45}
46
47template <typename K>
48static inline void invertMatrix(Dune::FieldMatrix<K,1,1>& matrix)
49{
50 Dune::FieldMatrix<K,1,1> tmp(matrix);
51 Dune::FMatrixHelp::invertMatrix(tmp,matrix);
52}
53
54template <typename K>
55static inline void invertMatrix(Dune::FieldMatrix<K,2,2>& matrix)
56{
57 Dune::FieldMatrix<K,2,2> tmp(matrix);
58 Dune::FMatrixHelp::invertMatrix(tmp,matrix);
59}
60
61template <typename K>
62static inline void invertMatrix(Dune::FieldMatrix<K,3,3>& matrix)
63{
64 Dune::FieldMatrix<K,3,3> tmp(matrix);
65 Dune::FMatrixHelp::invertMatrix(tmp,matrix);
66}
67
69template <template<class K> class Matrix, typename K>
70static inline K invertMatrix4(const Matrix<K>& matrix, Matrix<K>& inverse)
71{
72 inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] -
73 matrix[1][1] * matrix[2][3] * matrix[3][2] -
74 matrix[2][1] * matrix[1][2] * matrix[3][3] +
75 matrix[2][1] * matrix[1][3] * matrix[3][2] +
76 matrix[3][1] * matrix[1][2] * matrix[2][3] -
77 matrix[3][1] * matrix[1][3] * matrix[2][2];
78
79 inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] +
80 matrix[1][0] * matrix[2][3] * matrix[3][2] +
81 matrix[2][0] * matrix[1][2] * matrix[3][3] -
82 matrix[2][0] * matrix[1][3] * matrix[3][2] -
83 matrix[3][0] * matrix[1][2] * matrix[2][3] +
84 matrix[3][0] * matrix[1][3] * matrix[2][2];
85
86 inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] -
87 matrix[1][0] * matrix[2][3] * matrix[3][1] -
88 matrix[2][0] * matrix[1][1] * matrix[3][3] +
89 matrix[2][0] * matrix[1][3] * matrix[3][1] +
90 matrix[3][0] * matrix[1][1] * matrix[2][3] -
91 matrix[3][0] * matrix[1][3] * matrix[2][1];
92
93 inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] +
94 matrix[1][0] * matrix[2][2] * matrix[3][1] +
95 matrix[2][0] * matrix[1][1] * matrix[3][2] -
96 matrix[2][0] * matrix[1][2] * matrix[3][1] -
97 matrix[3][0] * matrix[1][1] * matrix[2][2] +
98 matrix[3][0] * matrix[1][2] * matrix[2][1];
99
100 inverse[0][1]= -matrix[0][1] * matrix[2][2] * matrix[3][3] +
101 matrix[0][1] * matrix[2][3] * matrix[3][2] +
102 matrix[2][1] * matrix[0][2] * matrix[3][3] -
103 matrix[2][1] * matrix[0][3] * matrix[3][2] -
104 matrix[3][1] * matrix[0][2] * matrix[2][3] +
105 matrix[3][1] * matrix[0][3] * matrix[2][2];
106
107 inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] -
108 matrix[0][0] * matrix[2][3] * matrix[3][2] -
109 matrix[2][0] * matrix[0][2] * matrix[3][3] +
110 matrix[2][0] * matrix[0][3] * matrix[3][2] +
111 matrix[3][0] * matrix[0][2] * matrix[2][3] -
112 matrix[3][0] * matrix[0][3] * matrix[2][2];
113
114 inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] +
115 matrix[0][0] * matrix[2][3] * matrix[3][1] +
116 matrix[2][0] * matrix[0][1] * matrix[3][3] -
117 matrix[2][0] * matrix[0][3] * matrix[3][1] -
118 matrix[3][0] * matrix[0][1] * matrix[2][3] +
119 matrix[3][0] * matrix[0][3] * matrix[2][1];
120
121 inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] -
122 matrix[0][0] * matrix[2][2] * matrix[3][1] -
123 matrix[2][0] * matrix[0][1] * matrix[3][2] +
124 matrix[2][0] * matrix[0][2] * matrix[3][1] +
125 matrix[3][0] * matrix[0][1] * matrix[2][2] -
126 matrix[3][0] * matrix[0][2] * matrix[2][1];
127
128 inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] -
129 matrix[0][1] * matrix[1][3] * matrix[3][2] -
130 matrix[1][1] * matrix[0][2] * matrix[3][3] +
131 matrix[1][1] * matrix[0][3] * matrix[3][2] +
132 matrix[3][1] * matrix[0][2] * matrix[1][3] -
133 matrix[3][1] * matrix[0][3] * matrix[1][2];
134
135 inverse[1][2] = -matrix[0][0] * matrix[1][2] * matrix[3][3] +
136 matrix[0][0] * matrix[1][3] * matrix[3][2] +
137 matrix[1][0] * matrix[0][2] * matrix[3][3] -
138 matrix[1][0] * matrix[0][3] * matrix[3][2] -
139 matrix[3][0] * matrix[0][2] * matrix[1][3] +
140 matrix[3][0] * matrix[0][3] * matrix[1][2];
141
142 inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] -
143 matrix[0][0] * matrix[1][3] * matrix[3][1] -
144 matrix[1][0] * matrix[0][1] * matrix[3][3] +
145 matrix[1][0] * matrix[0][3] * matrix[3][1] +
146 matrix[3][0] * matrix[0][1] * matrix[1][3] -
147 matrix[3][0] * matrix[0][3] * matrix[1][1];
148
149 inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] +
150 matrix[0][0] * matrix[1][2] * matrix[3][1] +
151 matrix[1][0] * matrix[0][1] * matrix[3][2] -
152 matrix[1][0] * matrix[0][2] * matrix[3][1] -
153 matrix[3][0] * matrix[0][1] * matrix[1][2] +
154 matrix[3][0] * matrix[0][2] * matrix[1][1];
155
156 inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] +
157 matrix[0][1] * matrix[1][3] * matrix[2][2] +
158 matrix[1][1] * matrix[0][2] * matrix[2][3] -
159 matrix[1][1] * matrix[0][3] * matrix[2][2] -
160 matrix[2][1] * matrix[0][2] * matrix[1][3] +
161 matrix[2][1] * matrix[0][3] * matrix[1][2];
162
163 inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] -
164 matrix[0][0] * matrix[1][3] * matrix[2][2] -
165 matrix[1][0] * matrix[0][2] * matrix[2][3] +
166 matrix[1][0] * matrix[0][3] * matrix[2][2] +
167 matrix[2][0] * matrix[0][2] * matrix[1][3] -
168 matrix[2][0] * matrix[0][3] * matrix[1][2];
169
170 inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] +
171 matrix[0][0] * matrix[1][3] * matrix[2][1] +
172 matrix[1][0] * matrix[0][1] * matrix[2][3] -
173 matrix[1][0] * matrix[0][3] * matrix[2][1] -
174 matrix[2][0] * matrix[0][1] * matrix[1][3] +
175 matrix[2][0] * matrix[0][3] * matrix[1][1];
176
177 inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] -
178 matrix[0][0] * matrix[1][2] * matrix[2][1] -
179 matrix[1][0] * matrix[0][1] * matrix[2][2] +
180 matrix[1][0] * matrix[0][2] * matrix[2][1] +
181 matrix[2][0] * matrix[0][1] * matrix[1][2] -
182 matrix[2][0] * matrix[0][2] * matrix[1][1];
183
184 K det = matrix[0][0] * inverse[0][0] + matrix[0][1] * inverse[1][0] +
185 matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0];
186
187 // return identity for singular or nearly singular matrices.
188 if (std::abs(det) < 1e-40) {
189 inverse = std::numeric_limits<K>::quiet_NaN();
190 throw NumericalProblem("Singular matrix");
191 } else
192 inverse *= 1.0 / det;
193
194 return det;
195}
196
197template<class K> using FMat4 = Dune::FieldMatrix<K,4,4>;
198
199template <typename K>
200static inline void invertMatrix(Dune::FieldMatrix<K,4,4>& matrix)
201{
202 FMat4<K> tmp(matrix);
203 invertMatrix4<FMat4>(tmp, matrix);
204}
205
206template <typename K>
207static inline void invertMatrix(Dune::DynamicMatrix<K>& matrix)
208{
209 // this function is only for 4 X 4 matrix
210 // for 4 X 4 matrix, using the invertMatrix() function above
211 // it is for temporary usage, mainly to reduce the huge burden of testing
212 // what algorithm should be used to invert 4 X 4 matrix will be handled
213 // as a seperate issue
214 if (matrix.rows() == 4) {
215 Dune::DynamicMatrix<K> A = matrix;
216 invertMatrix4(A, matrix);
217 return;
218 }
219
220 matrix.invert();
221}
222
223} // namespace detail
224
225template <class Scalar, int n, int m>
226class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
227{
228public:
229 using BaseType = Dune::FieldMatrix<Scalar, n, m> ;
230
231 using BaseType::operator= ;
232 using BaseType::rows;
233 using BaseType::cols;
234
235 MatrixBlock()
236 : BaseType(Scalar(0.0))
237 {}
238
239 explicit MatrixBlock(const Scalar value)
240 : BaseType(value)
241 {}
242
243 void invert()
244 { detail::invertMatrix(asBase()); }
245
246 const BaseType& asBase() const
247 { return static_cast<const BaseType&>(*this); }
248
249 BaseType& asBase()
250 { return static_cast<BaseType&>(*this); }
251};
252
253} // namespace Opm
254
255namespace Dune {
256
257template<class K, int n, int m>
258void print_row(std::ostream& s, const Opm::MatrixBlock<K, n, m>& A,
259 typename FieldMatrix<K, n, m>::size_type I,
260 typename FieldMatrix<K, n, m>::size_type J,
261 typename FieldMatrix<K, n, m>::size_type therow,
262 int width,
263 int precision)
264{ print_row(s, A.asBase(), I, J, therow, width, precision); }
265
266template <typename Scalar, int n, int m>
267struct MatrixDimension<Opm::MatrixBlock<Scalar, n, m> >
268 : public MatrixDimension<typename Opm::MatrixBlock<Scalar, n, m>::BaseType>
269{ };
270
271
272#if HAVE_UMFPACK
276template <typename T, typename A, int n, int m>
277class UMFPack<BCRSMatrix<Opm::MatrixBlock<T, n, m>, A> >
278 : public UMFPack<BCRSMatrix<FieldMatrix<T, n, m>, A> >
279{
280 using Base = UMFPack<BCRSMatrix<FieldMatrix<T, n, m>, A> >;
281 using Matrix = BCRSMatrix<FieldMatrix<T, n, m>, A>;
282
283public:
284 using RealMatrix = BCRSMatrix<Opm::MatrixBlock<T, n, m>, A>;
285
286 UMFPack(const RealMatrix& matrix, int verbose, bool)
287 : Base(reinterpret_cast<const Matrix&>(matrix), verbose)
288 {}
289};
290#endif
291
292#if HAVE_SUPERLU
296template <typename T, typename A, int n, int m>
297class SuperLU<BCRSMatrix<Opm::MatrixBlock<T, n, m>, A> >
298 : public SuperLU<BCRSMatrix<FieldMatrix<T, n, m>, A> >
299{
300 using Base = SuperLU<BCRSMatrix<FieldMatrix<T, n, m>, A> >;
301 using Matrix = BCRSMatrix<FieldMatrix<T, n, m>, A>;
302
303public:
304 using RealMatrix = BCRSMatrix<Opm::MatrixBlock<T, n, m>, A>;
305
306 SuperLU(const RealMatrix& matrix, int verb, bool reuse=true)
307 : Base(reinterpret_cast<const Matrix&>(matrix), verb, reuse)
308 {}
309};
310#endif
311
312template<typename T, int n, int m>
313struct IsNumber<Opm::MatrixBlock<T, n, m>>
314 : public IsNumber<Dune::FieldMatrix<T,n,m>>
315{};
316
317} // end namespace Dune
318
319
320#endif
Opm::MatrixBlock
Definition: matrixblock.hh:227
Opm::MatrixBlock::invert
void invert()
Definition: matrixblock.hh:243
Opm::MatrixBlock::asBase
BaseType & asBase()
Definition: matrixblock.hh:249
Opm::MatrixBlock::asBase
const BaseType & asBase() const
Definition: matrixblock.hh:246
Opm::MatrixBlock::MatrixBlock
MatrixBlock()
Definition: matrixblock.hh:235
Opm::MatrixBlock::BaseType
Dune::FieldMatrix< Scalar, n, m > BaseType
Definition: matrixblock.hh:229
Opm::MatrixBlock::MatrixBlock
MatrixBlock(const Scalar value)
Definition: matrixblock.hh:239
Dune
Definition: fvbaseprimaryvariables.hh:141
Dune::print_row
void print_row(std::ostream &s, const Opm::MatrixBlock< K, n, m > &A, typename FieldMatrix< K, n, m >::size_type I, typename FieldMatrix< K, n, m >::size_type J, typename FieldMatrix< K, n, m >::size_type therow, int width, int precision)
Definition: matrixblock.hh:258
Opm::detail::invertMatrix
static void invertMatrix(Dune::FieldMatrix< K, m, n > &matrix)
Definition: matrixblock.hh:42
Opm::detail::invertMatrix4
static K invertMatrix4(const Matrix< K > &matrix, Matrix< K > &inverse)
invert 4x4 Matrix without changing the original matrix
Definition: matrixblock.hh:70
Opm::detail::FMat4
Dune::FieldMatrix< K, 4, 4 > FMat4
Definition: matrixblock.hh:197
Opm::detail::invertMatrix
static void invertMatrix(Dune::DynamicMatrix< K > &matrix)
Definition: matrixblock.hh:207
Opm
Definition: blackoilboundaryratevector.hh:37