00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
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
00046
00047
00048 typedef boost::array< std::pair< double, WVector3d >, 3 > RealEigenSystem;
00049
00050
00051
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
00084
00085
00086
00087
00088
00089
00090
00091 template< typename Data_T >
00092 void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat, RealEigenSystem* es )
00093 {
00094 RealEigenSystem& result = *es;
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
00118
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
00196
00197
00198
00199
00200
00201 std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m );
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214 template< template< std::size_t, std::size_t, typename > class TensorType1,
00215 template< std::size_t, std::size_t, typename > class TensorType2,
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
00224 for( std::size_t i = 0; i < dim; ++i )
00225 {
00226 for( std::size_t j = 0; j < dim; ++j )
00227 {
00228
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
00241
00242
00243
00244
00245
00246
00247
00248
00249 template< typename Data_T >
00250 double evaluateSphericalFunction( WTensorSym< 4, 3, Data_T > const& tens, WVector3d const& gradient )
00251 {
00252
00253
00254
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
00277
00278
00279
00280
00281
00282
00283
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