OpenWalnut  1.2.5
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
WTensorFunctions.cpp
1 //---------------------------------------------------------------------------
2 //
3 // Project: OpenWalnut ( http://www.openwalnut.org )
4 //
5 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
6 // For more information see http://www.openwalnut.org/copying
7 //
8 // This file is part of OpenWalnut.
9 //
10 // OpenWalnut is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // OpenWalnut is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public License
21 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
22 //
23 //---------------------------------------------------------------------------
24 
25 #include <vector>
26 
27 #include "WTensorFunctions.h"
28 
29 std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m )
30 {
31  // this is copied from the gpu glyph shader
32  // src/graphicsEngine/shaders/tensorTools.fs
33  // originally implemented by Mario Hlawitschka
34  const double M_SQRT3 = 1.73205080756887729352744634151;
35  double de = m( 1, 2 ) * m( 1, 0 );
36  double dd = m( 1, 2 ) * m( 1, 2 );
37  double ee = m( 1, 0 ) * m( 1, 0 );
38  double ff = m( 2, 0 ) * m( 2, 0 );
39  double m0 = m( 0, 0 ) + m( 1, 1 ) + m( 2, 2 );
40  double c1 = m( 0, 0 ) * m( 1, 1 ) + m( 0, 0 ) * m( 2, 2 ) + m( 1, 1 ) * m( 2, 2 )
41  - ( dd + ee + ff );
42  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;
43 
44  double p, sqrt_p, q, c, s, phi;
45  p = m0 * m0 - 3. * c1;
46  q = m0 * ( p - ( 3. / 2. ) * c1 ) - ( 27. / 2. ) * c0;
47  sqrt_p = sqrt( fabs( p ) );
48 
49  phi = 27. * ( 0.25 * c1 * c1 * ( p - c1 ) + c0 * ( q + 27. / 4. * c0 ) );
50  phi = ( 1. / 3. ) * atan2( sqrt( fabs( phi ) ), q );
51 
52  c = sqrt_p * cos( phi );
53  s = ( 1. / M_SQRT3 ) * sqrt_p * sin( phi );
54 
55  std::vector< double > w( 3 );
56  // fix: swapped w[ 2 ] and w[ 1 ]
57  w[ 2 ] = ( 1. / 3. ) * ( m0 - c );
58  w[ 1 ] = w[ 2 ] + s;
59  w[ 0 ] = w[ 2 ] + c;
60  w[ 2 ] -= s;
61  return w;
62 }