OpenWalnut  1.4.0
WTensorFunctions.cpp
00001 //---------------------------------------------------------------------------
00002 //
00003 // Project: OpenWalnut ( http://www.openwalnut.org )
00004 //
00005 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
00006 // For more information see http://www.openwalnut.org/copying
00007 //
00008 // This file is part of OpenWalnut.
00009 //
00010 // OpenWalnut is free software: you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as published by
00012 // the Free Software Foundation, either version 3 of the License, or
00013 // (at your option) any later version.
00014 //
00015 // OpenWalnut is distributed in the hope that it will be useful,
00016 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018 // GNU Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public License
00021 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
00022 //
00023 //---------------------------------------------------------------------------
00024 
00025 #include <vector>
00026 
00027 #include "WTensorFunctions.h"
00028 
00029 std::ostream& operator<<( std::ostream& os, const RealEigenSystem& sys )
00030 {
00031     os << sys[0].first << ", " << sys[0].second << std::endl;
00032     os << sys[1].first << ", " << sys[1].second << std::endl;
00033     os << sys[2].first << ", " << sys[2].second << std::endl;
00034     return os;
00035 }
00036 
00037 std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m )
00038 {
00039     // this is copied from the gpu glyph shader
00040     // src/graphicsEngine/shaders/tensorTools.fs
00041     // originally implemented by Mario Hlawitschka
00042     const double M_SQRT3 = 1.73205080756887729352744634151;
00043     double de = m( 1, 2 ) * m( 1, 0 );
00044     double dd = m( 1, 2 ) * m( 1, 2 );
00045     double ee = m( 1, 0 ) * m( 1, 0 );
00046     double ff = m( 2, 0 ) * m( 2, 0 );
00047     double m0 = m( 0, 0 ) + m( 1, 1 ) + m( 2, 2 );
00048     double c1 = m( 0, 0 ) * m( 1, 1 ) + m( 0, 0 ) * m( 2, 2 ) + m( 1, 1 ) * m( 2, 2 )
00049              - ( dd + ee + ff );
00050     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;
00051 
00052     double p, sqrt_p, q, c, s, phi;
00053     p = m0 * m0 - 3. * c1;
00054     q = m0 * ( p - ( 3. / 2. ) * c1 ) - ( 27. / 2. ) * c0;
00055     sqrt_p = sqrt( fabs( p ) );
00056 
00057     phi = 27. * ( 0.25 * c1 * c1 * ( p - c1 ) + c0 * ( q + 27. / 4. * c0 ) );
00058     phi = ( 1. / 3. ) * atan2( sqrt( fabs( phi ) ), q );
00059 
00060     c = sqrt_p * cos( phi );
00061     s = ( 1. / M_SQRT3 ) * sqrt_p * sin( phi );
00062 
00063     std::vector< double > w( 3 );
00064     // fix: swapped w[ 2 ] and w[ 1 ]
00065     w[ 2 ] = ( 1. / 3. ) * ( m0 - c );
00066     w[ 1 ] = w[ 2 ] + s;
00067     w[ 0 ] = w[ 2 ] + c;
00068     w[ 2 ] -= s;
00069     return w;
00070 }