OpenWalnut  1.4.0
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 #ifndef Q_MOC_RUN
36 #include <boost/array.hpp>
37 #endif
38 
39 #include <Eigen/Dense>
40 
41 #include "../WAssert.h"
42 #include "../WLimits.h"
43 #include "WCompileTimeFunctions.h"
44 #include "WTensor.h"
45 #include "WTensorSym.h"
46 #include "WMatrix.h"
47 
48 /**
49  * An eigensystem has all eigenvalues as well as its corresponding eigenvectors. A RealEigenSystem is an EigenSystem where all
50  * eigenvalues are real and not complex.
51  */
52 typedef boost::array< std::pair< double, WVector3d >, 3 > RealEigenSystem;
53 
54 /**
55  * An eigensystem has all eigenvalues as well its corresponding eigenvectors.
56  */
57 typedef boost::array< std::pair< std::complex< double >, WVector3d >, 3 > EigenSystem;
58 
59 std::ostream& operator<<( std::ostream& os, const RealEigenSystem& sys );
60 
61 /**
62  * Helper function to sort eigen values on their size. As each eigenvalues has an eigenvector they must be twiddled the same.
63  *
64  * \param es Eigensystem consisting of eigenvalues and eigenvectors.
65  */
66 void sortRealEigenSystem( RealEigenSystem* es )
67 {
68  if( ( *es )[0].first > ( *es )[2].first )
69  {
70  std::swap( ( *es )[0], ( *es )[2] );
71  }
72  if( ( *es )[0].first > ( *es )[1].first )
73  {
74  std::swap( ( *es )[0], ( *es )[1] );
75  }
76  if( ( *es )[1].first > ( *es )[2].first )
77  {
78  std::swap( ( *es )[1], ( *es )[2] );
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  * Compute all eigenvalues as well as the corresponding eigenvectors of a
196  * symmetric real Matrix.
197  *
198  * \note Data_T must be castable to double.
199  *
200  * \param[in] mat A real symmetric matrix.
201  * \param[out] RealEigenSystem A pointer to an RealEigenSystem.
202  */
203 template< typename Data_T >
204 void eigenEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat, RealEigenSystem* es )
205 {
206  RealEigenSystem& result = *es; // alias for the result
207  typedef Eigen::Matrix< Data_T, 3, 3 > MatrixType;
208 
209  MatrixType m;
210  m << mat( 0, 0 ), mat( 0, 1 ), mat( 0, 2 ),
211  mat( 1, 0 ), mat( 1, 1 ), mat( 1, 2 ),
212  mat( 2, 0 ), mat( 2, 1 ), mat( 2, 2 );
213 
214  Eigen::SelfAdjointEigenSolver< MatrixType > solv( m );
215 
216  for( std::size_t k = 0; k < 3; ++k )
217  {
218  result[ k ].first = solv.eigenvalues()[ k ];
219  for( std::size_t j = 0; j < 3; ++j )
220  {
221  result[ k ].second[ j ] = solv.eigenvectors()( j, k );
222  }
223  }
224 }
225 
226 /**
227  * Calculate eigenvectors via the characteristic polynomial. This is essentially the same
228  * function as in the GPU glyph shaders. This is for 3 dimensions only.
229  *
230  * \param m The symmetric matrix to calculate the eigenvalues from.
231  * \return A std::vector of 3 eigenvalues in descending order (of their magnitude).
232  */
233 std::vector< double > getEigenvaluesCardano( WTensorSym< 2, 3 > const& m );
234 
235 /**
236  * Multiply tensors of order 2. This is essentially a matrix-matrix multiplication.
237  *
238  * Tensors must have the same data type and dimension, however both symmetric and asymmetric
239  * tensors (in any combination) are allowed as operands. The returned tensor is always an asymmetric tensor.
240  *
241  * \param one The first operand, a WTensor or WTensorSym of order 2.
242  * \param other The second operand, a WTensor or WTensorSym of order 2.
243  *
244  * \return A WTensor that is the product of the operands.
245  */
246 template< template< std::size_t, std::size_t, typename > class TensorType1, // NOLINT brainlint can't find TensorType1
247  template< std::size_t, std::size_t, typename > class TensorType2, // NOLINT
248  std::size_t dim,
249  typename Data_T >
250 WTensor< 2, dim, Data_T > operator * ( TensorType1< 2, dim, Data_T > const& one,
251  TensorType2< 2, dim, Data_T > const& other )
252 {
254 
255  // calc res_ij = one_ik * other_kj
256  for( std::size_t i = 0; i < dim; ++i )
257  {
258  for( std::size_t j = 0; j < dim; ++j )
259  {
260  // components are initialized to zero, so there is no need to zero them here
261  for( std::size_t k = 0; k < dim; ++k )
262  {
263  res( i, j ) += one( i, k ) * other( k, j );
264  }
265  }
266  }
267 
268  return res;
269 }
270 
271 /**
272  * Evaluate a spherical function represented by a symmetric 4th-order tensor for a given gradient.
273  *
274  * \tparam Data_T The integral type used to store the tensor elements.
275  *
276  * \param tens The tensor representing the spherical function.
277  * \param gradient The normalized vector that represents the gradient direction.
278  *
279  * \note If the gradient is not normalized, the result is undefined.
280  */
281 template< typename Data_T >
282 double evaluateSphericalFunction( WTensorSym< 4, 3, Data_T > const& tens, WVector3d const& gradient )
283 {
284  // use symmetry to reduce computation overhead
285  // temporaries for some of the gradient element multiplications could further reduce
286  // computation time
287  return gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0, 0, 0 )
288  + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1, 1, 1 )
289  + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2, 2, 2 )
290  + static_cast< Data_T >( 4 ) *
291  ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * tens( 0, 0, 0, 1 )
292  + gradient[ 0 ] * gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * tens( 0, 0, 0, 2 )
293  + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 0 ] * tens( 1, 1, 1, 0 )
294  + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 0 ] * tens( 2, 2, 2, 0 )
295  + gradient[ 1 ] * gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * tens( 1, 1, 1, 2 )
296  + gradient[ 2 ] * gradient[ 2 ] * gradient[ 2 ] * gradient[ 1 ] * tens( 2, 2, 2, 1 ) )
297  + static_cast< Data_T >( 12 ) *
298  ( gradient[ 2 ] * gradient[ 1 ] * gradient[ 0 ] * gradient[ 0 ] * tens( 2, 1, 0, 0 )
299  + gradient[ 0 ] * gradient[ 2 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 2, 1, 1 )
300  + gradient[ 0 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 1, 2, 2 ) )
301  + static_cast< Data_T >( 6 ) *
302  ( gradient[ 0 ] * gradient[ 0 ] * gradient[ 1 ] * gradient[ 1 ] * tens( 0, 0, 1, 1 )
303  + gradient[ 0 ] * gradient[ 0 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 0, 0, 2, 2 )
304  + gradient[ 1 ] * gradient[ 1 ] * gradient[ 2 ] * gradient[ 2 ] * tens( 1, 1, 2, 2 ) );
305 }
306 
307 /**
308  * Evaluate a spherical function represented by a symmetric 2nd-order tensor for a given gradient.
309  *
310  * \tparam Data_T The integral type used to store the tensor elements.
311  *
312  * \param tens The tensor representing the spherical function.
313  * \param gradient The normalized vector that represents the gradient direction.
314  *
315  * \note If the gradient is not normalized, the result is undefined.
316  */
317 template< typename Data_T >
318 double evaluateSphericalFunction( WTensorSym< 2, 3, Data_T > const& tens, WVector3d const& gradient )
319 {
320  return gradient[ 0 ] * gradient[ 0 ] * tens( 0, 0 )
321  + gradient[ 1 ] * gradient[ 1 ] * tens( 1, 1 )
322  + gradient[ 2 ] * gradient[ 2 ] * tens( 2, 2 )
323  + static_cast< Data_T >( 2 ) *
324  ( gradient[ 0 ] * gradient[ 1 ] * tens( 0, 1 )
325  + gradient[ 0 ] * gradient[ 2 ] * tens( 0, 2 )
326  + gradient[ 1 ] * gradient[ 2 ] * tens( 1, 2 ) );
327 }
328 
329 /**
330  * Calculate the logarithm of the given real symmetric tensor.
331  *
332  * \param tens The tensor.
333  * \return The logarithm of the tensor.
334  */
335 template< typename T >
336 WTensorSym< 2, 3, T > tensorLog( WTensorSym< 2, 3, T > const& tens )
337 {
339 
340  // calculate eigenvalues and eigenvectors
341  RealEigenSystem sys;
342  eigenEigenvector3D( tens, &sys );
343 
344  // first, we check for negative or zero eigenvalues
345  if( sys[ 0 ].first <= 0.0 || sys[ 1 ].first <= 0.0 || sys[ 2 ].first <= 0.0 )
346  {
347  res( 0, 0 ) = res( 1, 1 ) = res( 2, 2 ) = 1.0;
348  return res;
349  }
350 
351  // this implements the matrix product U * log( E ) * U.transposed()
352  // note that u( i, j ) = jth value of the ith eigenvector
353  res( 0, 0 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 0 ] * log( sys[ 0 ].first )
354  + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 0 ] * log( sys[ 1 ].first )
355  + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 0 ] * log( sys[ 2 ].first );
356  res( 0, 1 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 1 ] * log( sys[ 0 ].first )
357  + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 1 ] * log( sys[ 1 ].first )
358  + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 1 ] * log( sys[ 2 ].first );
359  res( 0, 2 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
360  + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
361  + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
362  res( 1, 1 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 1 ] * log( sys[ 0 ].first )
363  + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 1 ] * log( sys[ 1 ].first )
364  + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 1 ] * log( sys[ 2 ].first );
365  res( 1, 2 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
366  + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
367  + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
368  res( 2, 2 ) = sys[ 0 ].second[ 2 ] * sys[ 0 ].second[ 2 ] * log( sys[ 0 ].first )
369  + sys[ 1 ].second[ 2 ] * sys[ 1 ].second[ 2 ] * log( sys[ 1 ].first )
370  + sys[ 2 ].second[ 2 ] * sys[ 2 ].second[ 2 ] * log( sys[ 2 ].first );
371  return res;
372 }
373 
374 /**
375  * Calculate the exponential of the given real symmetric tensor.
376  *
377  * \param tens The tensor.
378  * \return The exponential of the tensor.
379  */
380 template< typename T >
381 WTensorSym< 2, 3, T > tensorExp( WTensorSym< 2, 3, T > const& tens )
382 {
384 
385  // calculate eigenvalues and eigenvectors
386  RealEigenSystem sys;
387  eigenEigenvector3D( tens, &sys );
388 
389  // this implements the matrix product U * exp( E ) * U.transposed()
390  // note that u( i, j ) = jth value of the ith eigenvector
391  res( 0, 0 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 0 ] * exp( sys[ 0 ].first )
392  + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 0 ] * exp( sys[ 1 ].first )
393  + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 0 ] * exp( sys[ 2 ].first );
394  res( 0, 1 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 1 ] * exp( sys[ 0 ].first )
395  + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 1 ] * exp( sys[ 1 ].first )
396  + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 1 ] * exp( sys[ 2 ].first );
397  res( 0, 2 ) = sys[ 0 ].second[ 0 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
398  + sys[ 1 ].second[ 0 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
399  + sys[ 2 ].second[ 0 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
400  res( 1, 1 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 1 ] * exp( sys[ 0 ].first )
401  + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 1 ] * exp( sys[ 1 ].first )
402  + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 1 ] * exp( sys[ 2 ].first );
403  res( 1, 2 ) = sys[ 0 ].second[ 1 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
404  + sys[ 1 ].second[ 1 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
405  + sys[ 2 ].second[ 1 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
406  res( 2, 2 ) = sys[ 0 ].second[ 2 ] * sys[ 0 ].second[ 2 ] * exp( sys[ 0 ].first )
407  + sys[ 1 ].second[ 2 ] * sys[ 1 ].second[ 2 ] * exp( sys[ 1 ].first )
408  + sys[ 2 ].second[ 2 ] * sys[ 2 ].second[ 2 ] * exp( sys[ 2 ].first );
409  return res;
410 }
411 
412 /**
413  * Find the maxima of a spherical function represented by a symmetric tensor.
414  *
415  *
416  *
417  */
418 template< std::size_t order, typename Data_T >
419 void findSymmetricSphericalFunctionMaxima( WTensorSym< order, 3, Data_T > const& tensor,
420  double threshold, double minAngularSeparationCosine, double stepSize,
421  std::vector< WVector3d > const& startingPoints,
422  std::vector< WVector3d >& maxima ) // NOLINT
423 {
424  std::vector< double > values;
425 
426  std::vector< WVector3d >::const_iterator it;
427  for( it = startingPoints.begin(); it != startingPoints.end(); ++it )
428  {
429  WVector3d position = *it;
430  WVector3d gradient = approximateGradient( tensor, position );
431  int iter;
432 
433  for( iter = 0; iter < 100; ++iter )
434  {
435  WVector3d newPosition = position + stepSize * gradient;
436  normalize( newPosition );
437  WVector3d newGradient = approximateGradient( tensor, newPosition );
438 
439  double cos = dot( newGradient, gradient );
440  if( cos > 0.985 ) // less then 10 degree
441  {
442  stepSize *= 2.0;
443  }
444  else if( cos < 0.866 ) // more than 30 degree
445  {
446  stepSize *= 0.5;
447  if( stepSize < 0.000001 )
448  {
449  break;
450  }
451  }
452  if( cos > 0.886 ) // less than 30 degree
453  {
454  gradient = newGradient;
455  position = newPosition;
456  }
457  }
458  if( iter != 100 )
459  {
460  double newFuncValue = tensor.evaluateSphericalFunction( position );
461 
462  bool add = true;
463 
464  std::vector< int > m;
465  m.reserve( maxima.size() );
466 
467  for( std::size_t k = 0; k != maxima.size(); ++k )
468  {
469  // find all maxima that are to close to this one
470  if( dot( position, maxima[ k ] ) > minAngularSeparationCosine
471  || dot( maxima[ k ], -1.0 * position ) > minAngularSeparationCosine )
472  {
473  m.push_back( k );
474  if( values[ k ] >= newFuncValue )
475  {
476  add = false;
477  }
478  }
479  }
480  if( add )
481  {
482  maxima.push_back( position );
483  values.push_back( newFuncValue );
484 
485  for( int k = static_cast< int >( m.size() - 1 ); k > 0; --k )
486  {
487  maxima.erase( maxima.begin() + m[ k ] );
488  values.erase( values.begin() + m[ k ] );
489  }
490  }
491  }
492  }
493 
494  // remove maxima that are too small
495  double max = 0;
496  for( std::size_t k = 0; k != maxima.size(); ++k )
497  {
498  if( values[ k ] > max )
499  {
500  max = values[ k ];
501  }
502  }
503 
504  std::size_t k = 0;
505  while( k < maxima.size() )
506  {
507  if( values[ k ] < threshold * max )
508  {
509  maxima.erase( maxima.begin() + k );
510  values.erase( values.begin() + k );
511  }
512  else
513  {
514  ++k;
515  }
516  }
517 }
518 
519 template< std::size_t order, typename Data_T >
520 WVector3d approximateGradient( WTensorSym< order, 3, Data_T > const& tensor, WVector3d const& pos )
521 {
522  WVector3d eXY;
523  WVector3d eZ;
524 
525  if( fabs( fabs( pos[ 2 ] ) - 1.0 ) < 0.001 )
526  {
527  eXY = WVector3d( 1.0, 0.0, 0.0 );
528  eZ = WVector3d( 0.0, 1.0, 0.0 );
529  }
530  else
531  {
532  // this is vectorProduct( z, pos )
533  eXY = WVector3d( -pos[ 1 ], pos[ 0 ], 0.0 );
534  normalize( eXY );
535  // this is vectorProduct( eXY, pos )
536  eZ = WVector3d( eXY[ 1 ] * pos[ 2 ], -eXY[ 0 ] * pos[ 2 ], eXY[ 0 ] * pos[ 1 ] - eXY[ 1 ] * pos[ 0 ] );
537  normalize( eZ );
538  }
539 
540  double dXY = ( tensor.evaluateSphericalFunction( normalize( pos + eXY * 0.0001 ) )
541  - tensor.evaluateSphericalFunction( normalize( pos - eXY * 0.0001 ) ) )
542  / 0.0002;
543  double dZ = ( tensor.evaluateSphericalFunction( normalize( pos + eZ * 0.0001 ) )
544  - tensor.evaluateSphericalFunction( normalize( pos - eZ * 0.0001 ) ) )
545  / 0.0002;
546 
547  // std::sqrt( 1.0 - z² ) = sin( acos( z ) ) = sin( theta ) in spherical coordinates
548  double d = 1.0 - pos[ 2 ] * pos[ 2 ];
549  if( d < 0.0 ) // avoid possible numerical problems
550  {
551  d = 0.0;
552  }
553  WVector3d res = eZ * dZ + eXY * ( dXY / std::sqrt( d ) );
554  normalize( res );
555  return res;
556 }
557 
558 #endif // WTENSORFUNCTIONS_H
Implements a symmetric tensor that has the same number of components in every direction.
Definition: WTensorSym.h:74
A fixed size matrix class.
Definition: WMatrixFixed.h:153
const double DBL_EPS
Smallest double such: 1.0 + DBL_EPS == 1.0 is still true.
Definition: WLimits.cpp:36
Implements a tensor that has the same number of components in every direction.
Definition: WTensor.h:78
Data_T evaluateSphericalFunction(WValue< Data_T > const &gradient) const
Evaluate - for a given gradient - the spherical function represented by this symmetric tensor...
Definition: WTensorSym.h:154
const double MAX_DOUBLE
Maximum double value.
Definition: WLimits.cpp:31