OpenWalnut
1.2.5
|
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