OpenWalnut  1.2.5
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
WTensorFunctions.h
1 //---------------------------------------------------------------------------
2 //
3 // Project: OpenWalnut ( http://www.openwalnut.org )
4 //
5 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
6 // For more information see http://www.openwalnut.org/copying
7 //
8 // This file is part of OpenWalnut.
9 //
10 // OpenWalnut is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // OpenWalnut is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public License
21 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
22 //
23 //---------------------------------------------------------------------------
24 
25 #ifndef WTENSORFUNCTIONS_H
26 #define WTENSORFUNCTIONS_H
27 
28 #include <algorithm>
29 #include <cmath>
30 #include <complex>
31 #include <iostream>
32 #include <utility>
33 #include <vector>
34 
35 #include <boost/array.hpp>
36 
37 #include "../WAssert.h"
38 #include "../WLimits.h"
39 #include "WCompileTimeFunctions.h"
40 #include "WTensor.h"
41 #include "WTensorSym.h"
42 #include "linearAlgebra/WLinearAlgebra.h"
43 
44 /**
45  * An eigensystem has all eigenvalues as well as its corresponding eigenvectors. A RealEigenSystem is an EigenSystem where all
46  * eigenvalues are real and not complex.
47  */
48 typedef boost::array< std::pair< double, WVector3d >, 3 > RealEigenSystem;
49 
50 /**
51  * An eigensystem has all eigenvalues as well its corresponding eigenvectors.
52  */
53 typedef boost::array< std::pair< std::complex< double >, WVector3d >, 3 > EigenSystem;
54 
55 std::ostream& operator<<( std::ostream& os, const RealEigenSystem& sys )
56 {
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;
60  return os;
61 }
62 
63 namespace
64 {
65  void sortRealEigenSystem( RealEigenSystem* es )
66  {
67  if( ( *es )[0].first > ( *es )[2].first )
68  {
69  std::swap( ( *es )[0], ( *es )[2] );
70  }
71  if( ( *es )[0].first > ( *es )[1].first )
72  {
73  std::swap( ( *es )[0], ( *es )[1] );
74  }
75  if( ( *es )[1].first > ( *es )[2].first )
76  {
77  std::swap( ( *es )[1], ( *es )[2] );
78  }
79  }
80 }
81 
82 /**
83  * Compute all eigenvalues as well as the corresponding eigenvectors of a
84  * symmetric real Matrix.
85  *
86  * \note Data_T must be castable to double.
87  *
88  * \param[in] mat A real symmetric matrix.
89  * \param[out] RealEigenSystem A pointer to an RealEigenSystem.
90  */
91 template< typename Data_T >
92 void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat, RealEigenSystem* es )
93 {
94  RealEigenSystem& result = *es; // alias for the result
97  ev( 0, 0 ) = ev( 1, 1 ) = ev( 2, 2 ) = 1.0;
98 
99  int iter = 50;
100  Data_T evp[ 3 ];
101  Data_T evq[ 3 ];
102 
103  while( iter >= 0 )
104  {
105  int p = 1;
106  int q = 0;
107 
108  for( int i = 0; i < 2; ++i )
109  {
110  if( fabs( in( 2, i ) ) > fabs( in( p, q ) ) )
111  {
112  p = 2;
113  q = i;
114  }
115  }
116 
117  // Note: If all non diagonal elements sum up to nearly zero, we may quit already!
118  // Thereby the chosen threshold 1.0e-50 was taken arbitrarily and is just a guess.
119  if( std::abs( in( 0, 1 ) ) + std::abs( in( 0, 2 ) ) + std::abs( in( 1, 2 ) ) < 1.0e-50 )
120  {
121  for( int i = 0; i < 3; ++i )
122  {
123  result[i].first = in( i, i );
124  for( int j = 0; j < 3; ++j )
125  {
126  result[i].second[j] = static_cast< double >( ev( j, i ) );
127  }
128  }
129  sortRealEigenSystem( es );
130  return;
131  }
132 
133  Data_T r = in( q, q ) - in( p, p );
134  Data_T o = r / ( 2.0 * in( p, q ) );
135 
136  Data_T t;
137  Data_T signofo = ( o < 0.0 ? -1.0 : 1.0 );
138  if( o * o > wlimits::MAX_DOUBLE )
139  {
140  t = signofo / ( 2.0 * fabs( o ) );
141  }
142  else
143  {
144  t = signofo / ( fabs( o ) + sqrt( o * o + 1.0 ) );
145  }
146 
147  Data_T c;
148 
149  if( t * t < wlimits::DBL_EPS )
150  {
151  c = 1.0;
152  }
153  else
154  {
155  c = 1.0 / sqrt( t * t + 1.0 );
156  }
157 
158  Data_T s = c * t;
159 
160  int k = 0;
161  while( k == q || k == p )
162  {
163  ++k;
164  }
165 
166  Data_T u = ( 1.0 - c ) / s;
167 
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 );
172  x = in( p, p );
173  y = in( q, q );
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;
177 
178  for( int l = 0; l < 3; ++l )
179  {
180  evp[ l ] = ev( l, p );
181  evq[ l ] = ev( l, q );
182  }
183  for( int l = 0; l < 3; ++l )
184  {
185  ev( l, p ) = c * evp[ l ] - s * evq[ l ];
186  ev( l, q ) = s * evp[ l ] + c * evq[ l ];
187  }
188 
189  --iter;
190  }
191  WAssert( iter >= 0, "Jacobi eigenvector iteration did not converge." );
192 }
193 
194 /**
195  * Calculate eigenvectors via the characteristic polynomial. This is essentially the same
196  * function as in the GPU glyph shaders. This is for 3 dimensions only.
197  *
198  * \param m The symmetric matrix to calculate the eigenvalues from.
199  * \return A std::vector of 3 eigenvalues in descending order (of their magnitude).
200  */
201 std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m );
202 
203 /**
204  * Multiply tensors of order 2. This is essentially a matrix-matrix multiplication.
205  *
206  * Tensors must have the same data type and dimension, however both symmetric and asymmetric
207  * tensors (in any combination) are allowed as operands. The returned tensor is always an asymmetric tensor.
208  *
209  * \param one The first operand, a WTensor or WTensorSym of order 2.
210  * \param other The second operand, a WTensor or WTensorSym of order 2.
211  *
212  * \return A WTensor that is the product of the operands.
213  */
214 template< template< std::size_t, std::size_t, typename > class TensorType1, // NOLINT brainlint can't find TensorType1
215  template< std::size_t, std::size_t, typename > class TensorType2, // NOLINT
216  std::size_t dim,
217  typename Data_T >
218 WTensor< 2, dim, Data_T > operator * ( TensorType1< 2, dim, Data_T > const& one,
219  TensorType2< 2, dim, Data_T > const& other )
220 {
222 
223  // calc res_ij = one_ik * other_kj
224  for( std::size_t i = 0; i < dim; ++i )
225  {
226  for( std::size_t j = 0; j < dim; ++j )
227  {
228  // components are initialized to zero, so there is no need to zero them here
229  for( std::size_t k = 0; k < dim; ++k )
230  {
231  res( i, j ) += one( i, k ) * other( k, j );
232  }
233  }
234  }
235 
236  return res;
237 }
238 
239 /**
240  * Evaluate a spherical function represented by a symmetric 4th-order tensor for a given gradient.
241  *
242  * \tparam Data_T The integral type used to store the tensor elements.
243  *
244  * \param tens The tensor representing the spherical function.
245  * \param gradient The normalized vector that represents the gradient direction.
246  *
247  * \note If the gradient is not normalized, the result is undefined.
248  */
249 template< typename Data_T >
250 double evaluateSphericalFunction( WTensorSym< 4, 3, Data_T > const& tens, WVector3d const& gradient )
251 {
252  // use symmetry to reduce computation overhead
253  // temporaries for some of the gradient element multiplications could further reduce
254  // computation time
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 ) );
273 }
274 
275 /**
276  * Evaluate a spherical function represented by a symmetric 2nd-order tensor for a given gradient.
277  *
278  * \tparam Data_T The integral type used to store the tensor elements.
279  *
280  * \param tens The tensor representing the spherical function.
281  * \param gradient The normalized vector that represents the gradient direction.
282  *
283  * \note If the gradient is not normalized, the result is undefined.
284  */
285 template< typename Data_T >
286 double evaluateSphericalFunction( WTensorSym< 2, 3, Data_T > const& tens, WVector3d const& gradient )
287 {
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 ) );
295 }
296 
297 #endif // WTENSORFUNCTIONS_H