Commit 14b509ea authored by Sebastian Eichelbaum's avatar Sebastian Eichelbaum

[MERGE]

parents e8bdeed5 0ca83ad8
......@@ -22,8 +22,44 @@
//
//---------------------------------------------------------------------------
#include <vector>
#include "WTensorFunctions.h"
namespace wmath
{
std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m )
{
// this is copied from the gpu glyph shader
// src/graphicsEngine/shaders/tensorTools.fs
// originally implemented by Mario Hlawitschka
const double M_SQRT3 = 1.73205080756887729352744634151;
double de = m( 1, 2 ) * m( 1, 0 );
double dd = m( 1, 2 ) * m( 1, 2 );
double ee = m( 1, 0 ) * m( 1, 0 );
double ff = m( 2, 0 ) * m( 2, 0 );
double m0 = m( 0, 0 ) + m( 1, 1 ) + m( 2, 2 );
double c1 = m( 0, 0 ) * m( 1, 1 ) + m( 0, 0 ) * m( 2, 2 ) + m( 1, 1 ) * m( 2, 2 )
- ( dd + ee + ff );
double c0 = m( 2, 2 ) * dd + m( 0, 0 ) * ee + m( 1, 1 ) * ff - m( 0, 0 ) * m( 1, 1 ) * m( 2, 2 ) - 2. * m( 2, 0 ) * de;
double p, sqrt_p, q, c, s, phi;
p = m0 * m0 - 3. * c1;
q = m0 * ( p - ( 3. / 2. ) * c1 ) - ( 27. / 2. ) * c0;
sqrt_p = sqrt( fabs( p ) );
phi = 27. * ( 0.25 * c1 * c1 * ( p - c1 ) + c0 * ( q + 27. / 4. * c0 ) );
phi = ( 1. / 3. ) * atan2( sqrt( fabs( phi ) ), q );
c = sqrt_p * cos( phi );
s = ( 1. / M_SQRT3 ) * sqrt_p * sin( phi );
std::vector< double > w( 3 );
// fix: swapped w[ 2 ] and w[ 1 ]
w[ 2 ] = ( 1. / 3. ) * ( m0 - c );
w[ 1 ] = w[ 2 ] + s;
w[ 0 ] = w[ 2 ] + c;
w[ 2 ] -= s;
return w;
}
} // namespace wmath
......@@ -37,12 +37,11 @@
namespace wmath
{
// TODO(reichenbach): needs to be tested!
/**
* Compute all eigenvalues as well as the corresponding eigenvectors of a
* symmetric real Matrix.
*
* \pre Data_T must be castable to double.
* \note Data_T must be castable to double.
*
* \param[in] mat A real symmetric matrix.
* \param[out] eigenValues A pointer to a vector of eigenvalues.
......@@ -55,6 +54,7 @@ void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat,
{
WTensorSym< 2, 3, Data_T > in( mat );
WTensor< 2, 3, Data_T > ev;
ev( 0, 0 ) = ev( 1, 1 ) = ev( 2, 2 ) = 1.0;
int iter = 50;
Data_T evp[ 3 ];
......@@ -81,7 +81,7 @@ void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat,
eigenValues->at( i ) = in( i, i );
for( int j = 0; j < 3; ++j )
{
eigenVectors->at( i )[ j ] = static_cast< double >( ev( i, j ) );
eigenVectors->at( i )[ j ] = static_cast< double >( ev( j, i ) );
}
}
return;
......@@ -149,6 +149,15 @@ void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat,
WAssert( iter >= 0, "Jacobi eigenvector iteration did not converge." );
}
/**
* Calculate eigenvectors via the characteristic polynomial. This is essentially the same
* function as in the gpu glyph shaders. This is for 3 dimensions only.
*
* \param m The symmetric matrix to calculate the eigenvalues from.
* \return A std::vector of 3 eigenvalues in descending order.
*/
std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m );
/**
* Multiply tensors of order 2. This is essentially a matrix-matrix multiplication.
*
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment