25 #ifndef WTENSORFUNCTIONS_H
26 #define WTENSORFUNCTIONS_H
35 #include <boost/array.hpp>
37 #include "../WAssert.h"
38 #include "../WLimits.h"
39 #include "WCompileTimeFunctions.h"
41 #include "WTensorSym.h"
42 #include "linearAlgebra/WLinearAlgebra.h"
48 typedef boost::array< std::pair< double, WVector3d >, 3 > RealEigenSystem;
53 typedef boost::array< std::pair< std::complex< double >,
WVector3d >, 3 > EigenSystem;
55 std::ostream& operator<<( std::ostream& os,
const RealEigenSystem& sys )
57 os << sys[0].first <<
", " << sys[0].second << std::endl;
58 os << sys[1].first <<
", " << sys[1].second << std::endl;
59 os << sys[2].first <<
", " << sys[2].second << std::endl;
65 void sortRealEigenSystem( RealEigenSystem* es )
67 if( ( *es )[0].first > ( *es )[2].first )
69 std::swap( ( *es )[0], ( *es )[2] );
71 if( ( *es )[0].first > ( *es )[1].first )
73 std::swap( ( *es )[0], ( *es )[1] );
75 if( ( *es )[1].first > ( *es )[2].first )
77 std::swap( ( *es )[1], ( *es )[2] );
91 template<
typename Data_T >
94 RealEigenSystem& result = *es;
97 ev( 0, 0 ) = ev( 1, 1 ) = ev( 2, 2 ) = 1.0;
108 for(
int i = 0; i < 2; ++i )
110 if( fabs( in( 2, i ) ) > fabs( in( p, q ) ) )
119 if( std::abs( in( 0, 1 ) ) + std::abs( in( 0, 2 ) ) + std::abs( in( 1, 2 ) ) < 1.0e-50 )
121 for(
int i = 0; i < 3; ++i )
123 result[i].first = in( i, i );
124 for(
int j = 0; j < 3; ++j )
126 result[i].second[j] =
static_cast< double >( ev( j, i ) );
129 sortRealEigenSystem( es );
133 Data_T r = in( q, q ) - in( p, p );
134 Data_T o = r / ( 2.0 * in( p, q ) );
137 Data_T signofo = ( o < 0.0 ? -1.0 : 1.0 );
140 t = signofo / ( 2.0 * fabs( o ) );
144 t = signofo / ( fabs( o ) + sqrt( o * o + 1.0 ) );
155 c = 1.0 / sqrt( t * t + 1.0 );
161 while( k == q || k == p )
166 Data_T u = ( 1.0 - c ) / s;
168 Data_T x = in( k, p );
169 Data_T y = in( k, q );
170 in( p, k ) = in( k, p ) = x - s * ( y + u * x );
171 in( q, k ) = in( k, q ) = y + s * ( x - u * y );
174 in( p, p ) = x - t * in( p, q );
175 in( q, q ) = y + t * in( p, q );
176 in( q, p ) = in( p, q ) = 0.0;
178 for(
int l = 0; l < 3; ++l )
180 evp[ l ] = ev( l, p );
181 evq[ l ] = ev( l, q );
183 for(
int l = 0; l < 3; ++l )
185 ev( l, p ) = c * evp[ l ] - s * evq[ l ];
186 ev( l, q ) = s * evp[ l ] + c * evq[ l ];
191 WAssert( iter >= 0,
"Jacobi eigenvector iteration did not converge." );
214 template<
template< std::
size_t, std::
size_t,
typename >
class TensorType1,
215 template< std::
size_t, std::
size_t,
typename >
class TensorType2,
219 TensorType2< 2, dim, Data_T >
const& other )
224 for( std::size_t i = 0; i < dim; ++i )
226 for( std::size_t j = 0; j < dim; ++j )
229 for( std::size_t k = 0; k < dim; ++k )
231 res( i, j ) += one( i, k ) * other( k, j );
249 template<
typename Data_T >
255 return gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0, 0, 0 )
256 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1, 1, 1 )
257 + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2, 2, 2 )
258 +
static_cast< Data_T
>( 4 ) *
259 ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * tens( 0, 0, 0, 1 )
260 + gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * tens( 0, 0, 0, 2 )
261 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 0 ] * tens( 1, 1, 1, 0 )
262 + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 0 ] * tens( 2, 2, 2, 0 )
263 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * tens( 1, 1, 1, 2 )
264 + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 1 ] * tens( 2, 2, 2, 1 ) )
265 + static_cast< Data_T >( 12 ) *
266 ( gradient[ 2 ] * gradient[ 1 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 2, 1, 0, 0 )
267 + gradient[ 0 ] * gradient[ 2 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 2, 1, 1 )
268 + gradient[ 0 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 1, 2, 2 ) )
269 + static_cast< Data_T >( 6 ) *
270 ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 0, 1, 1 )
271 + gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 0, 2, 2 )
272 + gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 1, 1, 2, 2 ) );
285 template<
typename Data_T >
288 return gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0 )
289 + gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1 )
290 + gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2 )
291 +
static_cast< Data_T
>( 2 ) *
292 ( gradient[ 0 ] * gradient[ 1 ] * tens( 0, 1 )
293 + gradient[ 0 ] * gradient[ 2 ] * tens( 0, 2 )
294 + gradient[ 1 ] * gradient[ 2 ] * tens( 1, 2 ) );
297 #endif // WTENSORFUNCTIONS_H