OpenWalnut 1.2.5

WTensorFunctions.h

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 #ifndef WTENSORFUNCTIONS_H
00026 #define WTENSORFUNCTIONS_H
00027 
00028 #include <algorithm>
00029 #include <cmath>
00030 #include <complex>
00031 #include <iostream>
00032 #include <utility>
00033 #include <vector>
00034 
00035 #include <boost/array.hpp>
00036 
00037 #include "../WAssert.h"
00038 #include "../WLimits.h"
00039 #include "WCompileTimeFunctions.h"
00040 #include "WTensor.h"
00041 #include "WTensorSym.h"
00042 #include "linearAlgebra/WLinearAlgebra.h"
00043 
00044 /**
00045  * An eigensystem has all eigenvalues as well as its corresponding eigenvectors. A RealEigenSystem is an EigenSystem where all
00046  * eigenvalues are real and not complex.
00047  */
00048 typedef boost::array< std::pair< double, WVector3d >, 3 > RealEigenSystem;
00049 
00050 /**
00051  * An eigensystem has all eigenvalues as well its corresponding eigenvectors.
00052  */
00053 typedef boost::array< std::pair< std::complex< double >, WVector3d >, 3 > EigenSystem;
00054 
00055 std::ostream& operator<<( std::ostream& os, const RealEigenSystem& sys )
00056 {
00057     os << sys[0].first << ", " << sys[0].second << std::endl;
00058     os << sys[1].first << ", " << sys[1].second << std::endl;
00059     os << sys[2].first << ", " << sys[2].second << std::endl;
00060     return os;
00061 }
00062 
00063 namespace
00064 {
00065     void sortRealEigenSystem( RealEigenSystem* es )
00066     {
00067         if( ( *es )[0].first > ( *es )[2].first )
00068         {
00069             std::swap( ( *es )[0], ( *es )[2] );
00070         }
00071         if( ( *es )[0].first > ( *es )[1].first )
00072         {
00073             std::swap( ( *es )[0], ( *es )[1] );
00074         }
00075         if( ( *es )[1].first > ( *es )[2].first )
00076         {
00077             std::swap( ( *es )[1], ( *es )[2] );
00078         }
00079     }
00080 }
00081 
00082 /**
00083  * Compute all eigenvalues as well as the corresponding eigenvectors of a
00084  * symmetric real Matrix.
00085  *
00086  * \note Data_T must be castable to double.
00087  *
00088  * \param[in] mat A real symmetric matrix.
00089  * \param[out] RealEigenSystem A pointer to an RealEigenSystem.
00090  */
00091 template< typename Data_T >
00092 void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat, RealEigenSystem* es )
00093 {
00094     RealEigenSystem& result = *es; // alias for the result
00095     WTensorSym< 2, 3, Data_T > in( mat );
00096     WTensor< 2, 3, Data_T > ev;
00097     ev( 0, 0 ) = ev( 1, 1 ) = ev( 2, 2 ) = 1.0;
00098 
00099     int iter = 50;
00100     Data_T evp[ 3 ];
00101     Data_T evq[ 3 ];
00102 
00103     while( iter >= 0 )
00104     {
00105         int p = 1;
00106         int q = 0;
00107 
00108         for( int i = 0; i < 2; ++i )
00109         {
00110             if( fabs( in( 2, i ) ) > fabs( in( p, q ) ) )
00111             {
00112                 p = 2;
00113                 q = i;
00114             }
00115         }
00116 
00117         // Note: If all non diagonal elements sum up to nearly zero, we may quit already!
00118         // Thereby the chosen threshold 1.0e-50 was taken arbitrarily and is just a guess.
00119         if( std::abs( in( 0, 1 ) ) + std::abs( in( 0, 2 ) ) + std::abs( in( 1, 2 ) ) < 1.0e-50 )
00120         {
00121             for( int i = 0; i < 3; ++i )
00122             {
00123                 result[i].first = in( i, i );
00124                 for( int j = 0; j < 3; ++j )
00125                 {
00126                     result[i].second[j] = static_cast< double >( ev( j, i ) );
00127                 }
00128             }
00129             sortRealEigenSystem( es );
00130             return;
00131         }
00132 
00133         Data_T r = in( q, q ) - in( p, p );
00134         Data_T o = r / ( 2.0 * in( p, q ) );
00135 
00136         Data_T t;
00137         Data_T signofo = ( o < 0.0 ? -1.0 : 1.0 );
00138         if( o * o > wlimits::MAX_DOUBLE )
00139         {
00140             t = signofo / ( 2.0 * fabs( o ) );
00141         }
00142         else
00143         {
00144             t = signofo / ( fabs( o ) + sqrt( o * o + 1.0 ) );
00145         }
00146 
00147         Data_T c;
00148 
00149         if( t * t < wlimits::DBL_EPS )
00150         {
00151             c = 1.0;
00152         }
00153         else
00154         {
00155             c = 1.0 / sqrt( t * t + 1.0 );
00156         }
00157 
00158         Data_T s = c * t;
00159 
00160         int k = 0;
00161         while( k == q || k == p )
00162         {
00163             ++k;
00164         }
00165 
00166         Data_T u = ( 1.0 - c ) / s;
00167 
00168         Data_T x = in( k, p );
00169         Data_T y = in( k, q );
00170         in( p, k ) = in( k, p ) = x - s * ( y + u * x );
00171         in( q, k ) = in( k, q ) = y + s * ( x - u * y );
00172         x = in( p, p );
00173         y = in( q, q );
00174         in( p, p ) = x - t * in( p, q );
00175         in( q, q ) = y + t * in( p, q );
00176         in( q, p ) = in( p, q ) = 0.0;
00177 
00178         for( int l = 0; l < 3; ++l )
00179         {
00180             evp[ l ] = ev( l, p );
00181             evq[ l ] = ev( l, q );
00182         }
00183         for( int l = 0; l < 3; ++l )
00184         {
00185             ev( l, p ) = c * evp[ l ] - s * evq[ l ];
00186             ev( l, q ) = s * evp[ l ] + c * evq[ l ];
00187         }
00188 
00189         --iter;
00190     }
00191     WAssert( iter >= 0, "Jacobi eigenvector iteration did not converge." );
00192 }
00193 
00194 /**
00195  * Calculate eigenvectors via the characteristic polynomial. This is essentially the same
00196  * function as in the GPU glyph shaders. This is for 3 dimensions only.
00197  *
00198  * \param m The symmetric matrix to calculate the eigenvalues from.
00199  * \return A std::vector of 3 eigenvalues in descending order (of their magnitude).
00200  */
00201 std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m );
00202 
00203 /**
00204  * Multiply tensors of order 2. This is essentially a matrix-matrix multiplication.
00205  *
00206  * Tensors must have the same data type and dimension, however both symmetric and asymmetric
00207  * tensors (in any combination) are allowed as operands. The returned tensor is always an asymmetric tensor.
00208  *
00209  * \param one The first operand, a WTensor or WTensorSym of order 2.
00210  * \param other The second operand, a WTensor or WTensorSym of order 2.
00211  *
00212  * \return A WTensor that is the product of the operands.
00213  */
00214 template< template< std::size_t, std::size_t, typename > class TensorType1, // NOLINT brainlint can't find TensorType1
00215           template< std::size_t, std::size_t, typename > class TensorType2, // NOLINT
00216           std::size_t dim,
00217           typename Data_T >
00218 WTensor< 2, dim, Data_T > operator * ( TensorType1< 2, dim, Data_T > const& one,
00219                                        TensorType2< 2, dim, Data_T > const& other )
00220 {
00221     WTensor< 2, dim, Data_T > res;
00222 
00223     // calc res_ij = one_ik * other_kj
00224     for( std::size_t i = 0; i < dim; ++i )
00225     {
00226         for( std::size_t j = 0; j < dim; ++j )
00227         {
00228             // components are initialized to zero, so there is no need to zero them here
00229             for( std::size_t k = 0; k < dim; ++k )
00230             {
00231                 res( i, j ) += one( i, k ) * other( k, j );
00232             }
00233         }
00234     }
00235 
00236     return res;
00237 }
00238 
00239 /**
00240  * Evaluate a spherical function represented by a symmetric 4th-order tensor for a given gradient.
00241  *
00242  * \tparam Data_T The integral type used to store the tensor elements.
00243  *
00244  * \param tens The tensor representing the spherical function.
00245  * \param gradient The normalized vector that represents the gradient direction.
00246  *
00247  * \note If the gradient is not normalized, the result is undefined.
00248  */
00249 template< typename Data_T >
00250 double evaluateSphericalFunction( WTensorSym< 4, 3, Data_T > const& tens, WVector3d const& gradient )
00251 {
00252     // use symmetry to reduce computation overhead
00253     // temporaries for some of the gradient element multiplications could further reduce
00254     // computation time
00255     return gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0, 0, 0 )
00256          + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1, 1, 1 )
00257          + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2, 2, 2 )
00258          + static_cast< Data_T >( 4 ) *
00259          ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * tens( 0, 0, 0, 1 )
00260          + gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * tens( 0, 0, 0, 2 )
00261          + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 0 ] * tens( 1, 1, 1, 0 )
00262          + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 0 ] * tens( 2, 2, 2, 0 )
00263          + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * tens( 1, 1, 1, 2 )
00264          + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 1 ] * tens( 2, 2, 2, 1 ) )
00265          + static_cast< Data_T >( 12 ) *
00266          ( gradient[ 2 ] * gradient[ 1 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 2, 1, 0, 0 )
00267          + gradient[ 0 ] * gradient[ 2 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 2, 1, 1 )
00268          + gradient[ 0 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 1, 2, 2 ) )
00269          + static_cast< Data_T >( 6 ) *
00270          ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 0, 1, 1 )
00271          + gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 0, 2, 2 )
00272          + gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 1, 1, 2, 2 ) );
00273 }
00274 
00275 /**
00276  * Evaluate a spherical function represented by a symmetric 2nd-order tensor for a given gradient.
00277  *
00278  * \tparam Data_T The integral type used to store the tensor elements.
00279  *
00280  * \param tens The tensor representing the spherical function.
00281  * \param gradient The normalized vector that represents the gradient direction.
00282  *
00283  * \note If the gradient is not normalized, the result is undefined.
00284  */
00285 template< typename Data_T >
00286 double evaluateSphericalFunction( WTensorSym< 2, 3, Data_T > const& tens, WVector3d const& gradient )
00287 {
00288     return gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0 )
00289          + gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1 )
00290          + gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2 )
00291          + static_cast< Data_T >( 2 ) *
00292          ( gradient[ 0 ] * gradient[ 1 ] * tens( 0, 1 )
00293          + gradient[ 0 ] * gradient[ 2 ] * tens( 0, 2 )
00294          + gradient[ 1 ] * gradient[ 2 ] * tens( 1, 2 ) );
00295 }
00296 
00297 #endif  // WTENSORFUNCTIONS_H
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends