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 // Note for maintainers: If updating this value, also update the value
189 // in gpuistl/detail/deviceBlockOperations.hpp
190 if (std::abs(det) < 1e-40) {
191 inverse = std::numeric_limits<K>::quiet_NaN();
192 throw NumericalProblem("Singular matrix");
193 } else
194 inverse *= 1.0 / det;
195
196 return det;
197}
198
199template<class K> using FMat4 = Dune::FieldMatrix<K,4,4>;
200
201template <typename K>
202static inline void invertMatrix(Dune::FieldMatrix<K,4,4>& matrix)
203{
204 FMat4<K> tmp(matrix);
205 invertMatrix4<FMat4>(tmp, matrix);
206}
207
208template <typename K>
209static inline void invertMatrix(Dune::DynamicMatrix<K>& matrix)
210{
211 // this function is only for 4 X 4 matrix
212 // for 4 X 4 matrix, using the invertMatrix() function above
213 // it is for temporary usage, mainly to reduce the huge burden of testing
214 // what algorithm should be used to invert 4 X 4 matrix will be handled
215 // as a seperate issue
216 if (matrix.rows() == 4) {
217 Dune::DynamicMatrix<K> A = matrix;
218 invertMatrix4(A, matrix);
219 return;
220 }
221
222 matrix.invert();
223}
224
225} // namespace detail
226
227template <class Scalar, int n, int m>
228class MatrixBlock : public Dune::FieldMatrix<Scalar, n, m>
229{
230public:
231 using BaseType = Dune::FieldMatrix<Scalar, n, m> ;
232
233 using BaseType::operator= ;
234 using BaseType::rows;
235 using BaseType::cols;
236
237 MatrixBlock()
238 : BaseType(Scalar(0.0))
239 {}
240
241 explicit MatrixBlock(const Scalar value)
242 : BaseType(value)
243 {}
244
245 void invert()
246 { detail::invertMatrix(asBase()); }
247
248 const BaseType& asBase() const
249 { return static_cast<const BaseType&>(*this); }
250
251 BaseType& asBase()
252 { return static_cast<BaseType&>(*this); }
253};
254
255} // namespace Opm
256
257namespace Dune {
258
259template<class K, int n, int m>
260void print_row(std::ostream& s, const Opm::MatrixBlock<K, n, m>& A,
261 typename FieldMatrix<K, n, m>::size_type I,
262 typename FieldMatrix<K, n, m>::size_type J,
263 typename FieldMatrix<K, n, m>::size_type therow,
264 int width,
265 int precision)
266{ print_row(s, A.asBase(), I, J, therow, width, precision); }
267
268template <typename Scalar, int n, int m>
269struct MatrixDimension<Opm::MatrixBlock<Scalar, n, m> >
270 : public MatrixDimension<typename Opm::MatrixBlock<Scalar, n, m>::BaseType>
271{ };
272
273
274#if HAVE_SUITESPARSE_UMFPACK
278template <typename T, typename A, int n, int m>
279class UMFPack<BCRSMatrix<Opm::MatrixBlock<T, n, m>, A> >
280 : public UMFPack<BCRSMatrix<FieldMatrix<T, n, m>, A> >
281{
282 using Base = UMFPack<BCRSMatrix<FieldMatrix<T, n, m>, A> >;
283 using Matrix = BCRSMatrix<FieldMatrix<T, n, m>, A>;
284
285public:
286 using RealMatrix = BCRSMatrix<Opm::MatrixBlock<T, n, m>, A>;
287
288 UMFPack(const RealMatrix& matrix, int verbose, bool)
289 : Base(reinterpret_cast<const Matrix&>(matrix), verbose)
290 {}
291};
292#endif
293
294#if HAVE_SUPERLU
298template <typename T, typename A, int n, int m>
299class SuperLU<BCRSMatrix<Opm::MatrixBlock<T, n, m>, A> >
300 : public SuperLU<BCRSMatrix<FieldMatrix<T, n, m>, A> >
301{
302 using Base = SuperLU<BCRSMatrix<FieldMatrix<T, n, m>, A> >;
303 using Matrix = BCRSMatrix<FieldMatrix<T, n, m>, A>;
304
305public:
306 using RealMatrix = BCRSMatrix<Opm::MatrixBlock<T, n, m>, A>;
307
308 SuperLU(const RealMatrix& matrix, int verb, bool reuse=true)
309 : Base(reinterpret_cast<const Matrix&>(matrix), verb, reuse)
310 {}
311};
312#endif
313
314template<typename T, int n, int m>
315struct IsNumber<Opm::MatrixBlock<T, n, m>>
316 : public IsNumber<Dune::FieldMatrix<T,n,m>>
317{};
318
319} // end namespace Dune
320
321
322#endif
Dune::UMFPack
Definition: MSWellHelpers.hpp:29
Opm::MatrixBlock
Definition: matrixblock.hh:229
Opm::MatrixBlock::invert
void invert()
Definition: matrixblock.hh:245
Opm::MatrixBlock::asBase
BaseType & asBase()
Definition: matrixblock.hh:251
Opm::MatrixBlock::asBase
const BaseType & asBase() const
Definition: matrixblock.hh:248
Opm::MatrixBlock::MatrixBlock
MatrixBlock()
Definition: matrixblock.hh:237
Opm::MatrixBlock::BaseType
Dune::FieldMatrix< Scalar, n, m > BaseType
Definition: matrixblock.hh:231
Opm::MatrixBlock::MatrixBlock
MatrixBlock(const Scalar value)
Definition: matrixblock.hh:241
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:260
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:199
Opm::detail::invertMatrix
static void invertMatrix(Dune::DynamicMatrix< K > &matrix)
Definition: matrixblock.hh:209
Opm
Definition: blackoilbioeffectsmodules.hh:43