diff --git a/opm/simulators/linalg/matrixblock.hh b/opm/simulators/linalg/matrixblock.hh index 1f84bcfe0..a2cb8df58 100644 --- a/opm/simulators/linalg/matrixblock.hh +++ b/opm/simulators/linalg/matrixblock.hh @@ -36,191 +36,6 @@ #include namespace Opm { -namespace detail { - -template -static inline void invertMatrix(Dune::FieldMatrix& matrix) -{ - matrix.invert(); -} - -template -static inline void invertMatrix(Dune::FieldMatrix& matrix) -{ - Dune::FieldMatrix tmp(matrix); - Dune::FMatrixHelp::invertMatrix(tmp,matrix); -} - -template -static inline void invertMatrix(Dune::FieldMatrix& matrix) -{ - Dune::FieldMatrix tmp(matrix); - Dune::FMatrixHelp::invertMatrix(tmp,matrix); -} - -template -static inline void invertMatrix(Dune::FieldMatrix& matrix) -{ - Dune::FieldMatrix tmp(matrix); - Dune::FMatrixHelp::invertMatrix(tmp,matrix); -} - -//! invert 4x4 Matrix without changing the original matrix -template class Matrix, typename K> -static inline K invertMatrix4(const Matrix& matrix, Matrix& inverse) -{ - inverse[0][0] = matrix[1][1] * matrix[2][2] * matrix[3][3] - - matrix[1][1] * matrix[2][3] * matrix[3][2] - - matrix[2][1] * matrix[1][2] * matrix[3][3] + - matrix[2][1] * matrix[1][3] * matrix[3][2] + - matrix[3][1] * matrix[1][2] * matrix[2][3] - - matrix[3][1] * matrix[1][3] * matrix[2][2]; - - inverse[1][0] = -matrix[1][0] * matrix[2][2] * matrix[3][3] + - matrix[1][0] * matrix[2][3] * matrix[3][2] + - matrix[2][0] * matrix[1][2] * matrix[3][3] - - matrix[2][0] * matrix[1][3] * matrix[3][2] - - matrix[3][0] * matrix[1][2] * matrix[2][3] + - matrix[3][0] * matrix[1][3] * matrix[2][2]; - - inverse[2][0] = matrix[1][0] * matrix[2][1] * matrix[3][3] - - matrix[1][0] * matrix[2][3] * matrix[3][1] - - matrix[2][0] * matrix[1][1] * matrix[3][3] + - matrix[2][0] * matrix[1][3] * matrix[3][1] + - matrix[3][0] * matrix[1][1] * matrix[2][3] - - matrix[3][0] * matrix[1][3] * matrix[2][1]; - - inverse[3][0] = -matrix[1][0] * matrix[2][1] * matrix[3][2] + - matrix[1][0] * matrix[2][2] * matrix[3][1] + - matrix[2][0] * matrix[1][1] * matrix[3][2] - - matrix[2][0] * matrix[1][2] * matrix[3][1] - - matrix[3][0] * matrix[1][1] * matrix[2][2] + - matrix[3][0] * matrix[1][2] * matrix[2][1]; - - inverse[0][1]= -matrix[0][1] * matrix[2][2] * matrix[3][3] + - matrix[0][1] * matrix[2][3] * matrix[3][2] + - matrix[2][1] * matrix[0][2] * matrix[3][3] - - matrix[2][1] * matrix[0][3] * matrix[3][2] - - matrix[3][1] * matrix[0][2] * matrix[2][3] + - matrix[3][1] * matrix[0][3] * matrix[2][2]; - - inverse[1][1] = matrix[0][0] * matrix[2][2] * matrix[3][3] - - matrix[0][0] * matrix[2][3] * matrix[3][2] - - matrix[2][0] * matrix[0][2] * matrix[3][3] + - matrix[2][0] * matrix[0][3] * matrix[3][2] + - matrix[3][0] * matrix[0][2] * matrix[2][3] - - matrix[3][0] * matrix[0][3] * matrix[2][2]; - - inverse[2][1] = -matrix[0][0] * matrix[2][1] * matrix[3][3] + - matrix[0][0] * matrix[2][3] * matrix[3][1] + - matrix[2][0] * matrix[0][1] * matrix[3][3] - - matrix[2][0] * matrix[0][3] * matrix[3][1] - - matrix[3][0] * matrix[0][1] * matrix[2][3] + - matrix[3][0] * matrix[0][3] * matrix[2][1]; - - inverse[3][1] = matrix[0][0] * matrix[2][1] * matrix[3][2] - - matrix[0][0] * matrix[2][2] * matrix[3][1] - - matrix[2][0] * matrix[0][1] * matrix[3][2] + - matrix[2][0] * matrix[0][2] * matrix[3][1] + - matrix[3][0] * matrix[0][1] * matrix[2][2] - - matrix[3][0] * matrix[0][2] * matrix[2][1]; - - inverse[0][2] = matrix[0][1] * matrix[1][2] * matrix[3][3] - - matrix[0][1] * matrix[1][3] * matrix[3][2] - - matrix[1][1] * matrix[0][2] * matrix[3][3] + - matrix[1][1] * matrix[0][3] * matrix[3][2] + - matrix[3][1] * matrix[0][2] * matrix[1][3] - - matrix[3][1] * matrix[0][3] * matrix[1][2]; - - inverse[1][2] = -matrix[0][0] * matrix[1][2] * matrix[3][3] + - matrix[0][0] * matrix[1][3] * matrix[3][2] + - matrix[1][0] * matrix[0][2] * matrix[3][3] - - matrix[1][0] * matrix[0][3] * matrix[3][2] - - matrix[3][0] * matrix[0][2] * matrix[1][3] + - matrix[3][0] * matrix[0][3] * matrix[1][2]; - - inverse[2][2] = matrix[0][0] * matrix[1][1] * matrix[3][3] - - matrix[0][0] * matrix[1][3] * matrix[3][1] - - matrix[1][0] * matrix[0][1] * matrix[3][3] + - matrix[1][0] * matrix[0][3] * matrix[3][1] + - matrix[3][0] * matrix[0][1] * matrix[1][3] - - matrix[3][0] * matrix[0][3] * matrix[1][1]; - - inverse[3][2] = -matrix[0][0] * matrix[1][1] * matrix[3][2] + - matrix[0][0] * matrix[1][2] * matrix[3][1] + - matrix[1][0] * matrix[0][1] * matrix[3][2] - - matrix[1][0] * matrix[0][2] * matrix[3][1] - - matrix[3][0] * matrix[0][1] * matrix[1][2] + - matrix[3][0] * matrix[0][2] * matrix[1][1]; - - inverse[0][3] = -matrix[0][1] * matrix[1][2] * matrix[2][3] + - matrix[0][1] * matrix[1][3] * matrix[2][2] + - matrix[1][1] * matrix[0][2] * matrix[2][3] - - matrix[1][1] * matrix[0][3] * matrix[2][2] - - matrix[2][1] * matrix[0][2] * matrix[1][3] + - matrix[2][1] * matrix[0][3] * matrix[1][2]; - - inverse[1][3] = matrix[0][0] * matrix[1][2] * matrix[2][3] - - matrix[0][0] * matrix[1][3] * matrix[2][2] - - matrix[1][0] * matrix[0][2] * matrix[2][3] + - matrix[1][0] * matrix[0][3] * matrix[2][2] + - matrix[2][0] * matrix[0][2] * matrix[1][3] - - matrix[2][0] * matrix[0][3] * matrix[1][2]; - - inverse[2][3] = -matrix[0][0] * matrix[1][1] * matrix[2][3] + - matrix[0][0] * matrix[1][3] * matrix[2][1] + - matrix[1][0] * matrix[0][1] * matrix[2][3] - - matrix[1][0] * matrix[0][3] * matrix[2][1] - - matrix[2][0] * matrix[0][1] * matrix[1][3] + - matrix[2][0] * matrix[0][3] * matrix[1][1]; - - inverse[3][3] = matrix[0][0] * matrix[1][1] * matrix[2][2] - - matrix[0][0] * matrix[1][2] * matrix[2][1] - - matrix[1][0] * matrix[0][1] * matrix[2][2] + - matrix[1][0] * matrix[0][2] * matrix[2][1] + - matrix[2][0] * matrix[0][1] * matrix[1][2] - - matrix[2][0] * matrix[0][2] * matrix[1][1]; - - K det = matrix[0][0] * inverse[0][0] + matrix[0][1] * inverse[1][0] + - matrix[0][2] * inverse[2][0] + matrix[0][3] * inverse[3][0]; - - // return identity for singular or nearly singular matrices. - if (std::abs(det) < 1e-40) { - inverse = std::numeric_limits::quiet_NaN(); - throw NumericalProblem("Singular matrix"); - } else - inverse *= 1.0 / det; - - return det; -} - -template using FMat4 = Dune::FieldMatrix; - -template -static inline void invertMatrix(Dune::FieldMatrix& matrix) -{ - FMat4 tmp(matrix); - invertMatrix4(tmp, matrix); -} - -template -static inline void invertMatrix(Dune::DynamicMatrix& matrix) -{ - // this function is only for 4 X 4 matrix - // for 4 X 4 matrix, using the invertMatrix() function above - // it is for temporary usage, mainly to reduce the huge burden of testing - // what algorithm should be used to invert 4 X 4 matrix will be handled - // as a seperate issue - if (matrix.rows() == 4) { - Dune::DynamicMatrix A = matrix; - invertMatrix4(A, matrix); - return; - } - - matrix.invert(); -} - -} // namespace detail template class MatrixBlock : public Dune::FieldMatrix @@ -240,9 +55,6 @@ public: : BaseType(value) {} - void invert() - { detail::invertMatrix(asBase()); } - const BaseType& asBase() const { return static_cast(*this); }