Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
#define DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
/** \file
* \brief max eigenvalue computations for the FieldMatrix class
*/
#include <iostream>
#include <cmath>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fmatrixev.hh>
template <typename K>
static K maxEigenvalue(const Dune::FieldMatrix<K, 1, 1>& matrix) {
return matrix[0][0];
}
/** \brief calculates the eigenvalues of a symmetric field matrix
\param[in] matrix matrix eigenvalues are calculated for
\param[out] max eigenvalue
*/
template <typename K>
static K maxEigenvalue(const Dune::FieldMatrix<K, 2, 2>& matrix) {
const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
const K p2 = p - matrix[1][1];
K q = p2 * p2 + matrix[1][0] * matrix[0][1];
if( q < 0 && q > -1e-14 ) q = 0;
if (q < 0) {
std::cout << matrix << std::endl;
// Complex eigenvalues are either caused by non-symmetric matrices or by round-off errors
DUNE_THROW(Dune::MathError, "Complex eigenvalue detected (which this implementation cannot handle).");
}
// get square root
q = std::sqrt(q);
return p + q;
}
/** \brief Calculates the eigenvalues of a symmetric 3x3 field matrix
\param[in] matrix eigenvalues are calculated for
\param[out] max eigenvalue
\note If the input matrix is not symmetric the behavior of this method is undefined.
*/
template <typename K>
static K maxEigenvalue(const Dune::FieldMatrix<K, 3, 3>& matrix) {
Dune::FieldVector<K, 3> eigenvalues;
Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
return eigenvalues[0];
}
/** @} end documentation */
#endif