23#ifndef EWOMS_MATRIX_BLOCK_HH
24#define EWOMS_MATRIX_BLOCK_HH
26#include <dune/common/dynmatrix.hh>
27#include <dune/common/fmatrix.hh>
28#include <dune/common/typetraits.hh>
30#include <dune/istl/superlu.hh>
31#include <dune/istl/umfpack.hh>
32#include <dune/istl/matrixutils.hh>
34#include <opm/common/Exceptions.hpp>
41template <
typename K,
int m,
int n>
42static inline void invertMatrix(Dune::FieldMatrix<K,m,n>& matrix)
48static inline void invertMatrix(Dune::FieldMatrix<K,1,1>& matrix)
50 Dune::FieldMatrix<K,1,1> tmp(matrix);
51 Dune::FMatrixHelp::invertMatrix(tmp,matrix);
55static inline void invertMatrix(Dune::FieldMatrix<K,2,2>& matrix)
57 Dune::FieldMatrix<K,2,2> tmp(matrix);
58 Dune::FMatrixHelp::invertMatrix(tmp,matrix);
62static inline void invertMatrix(Dune::FieldMatrix<K,3,3>& matrix)
64 Dune::FieldMatrix<K,3,3> tmp(matrix);
65 Dune::FMatrixHelp::invertMatrix(tmp,matrix);
69template <
template<
class K>
class Matrix,
typename K>
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
188 if (std::abs(
det) < 1
e-40) {
189 inverse = std::numeric_limits<K>::quiet_NaN();
197template<
class K>
using FMat4 = Dune::FieldMatrix<K,4,4>;
200static inline void invertMatrix(Dune::FieldMatrix<K,4,4>& matrix)
207static inline void invertMatrix(Dune::DynamicMatrix<K>& matrix)
214 if (matrix.rows() == 4) {
215 Dune::DynamicMatrix<K>
A = matrix;
216 invertMatrix4(
A, matrix);
225template <
class Scalar,
int n,
int m>
229 using BaseType = Dune::FieldMatrix<Scalar, n, m> ;
231 using BaseType::operator= ;
232 using BaseType::rows;
233 using BaseType::cols;
236 : BaseType(Scalar(0.0))
244 { detail::invertMatrix(asBase()); }
246 const BaseType& asBase()
const
247 {
return static_cast<const BaseType&
>(*this); }
250 {
return static_cast<BaseType&
>(*this); }
257template<
class K,
int n,
int m>
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,
264{ print_row(s, A.asBase(), I, J, therow, width, precision); }
266template <
typename Scalar,
int n,
int m>
267struct MatrixDimension<
Opm::MatrixBlock<Scalar, n, m> >
268 :
public MatrixDimension<typename Opm::MatrixBlock<Scalar, n, m>::BaseType>
272#if HAVE_SUITESPARSE_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> >
281 using Matrix = BCRSMatrix<FieldMatrix<T, n, m>, A>;
284 using RealMatrix = BCRSMatrix<Opm::MatrixBlock<T, n, m>, A>;
286 UMFPack(
const RealMatrix& matrix,
int verbose,
bool)
287 : Base(reinterpret_cast<const Matrix&>(matrix), verbose)
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> >
300 using Base = SuperLU<BCRSMatrix<FieldMatrix<T, n, m>, A> >;
301 using Matrix = BCRSMatrix<FieldMatrix<T, n, m>, A>;
304 using RealMatrix = BCRSMatrix<Opm::MatrixBlock<T, n, m>, A>;
306 SuperLU(
const RealMatrix& matrix,
int verb,
bool reuse=
true)
307 : Base(reinterpret_cast<const Matrix&>(matrix), verb, reuse)
312template<
typename T,
int n,
int m>
313struct IsNumber<
Opm::MatrixBlock<T, n, m>>
314 :
public IsNumber<Dune::FieldMatrix<T,n,m>>
Definition MSWellHelpers.hpp:29
Definition matrixblock.hh:227
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242