OpenWalnut  1.5.0dev
WTensorFunctions_test.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_TEST_H
26 #define WTENSORFUNCTIONS_TEST_H
27 
28 #include <string>
29 #include <vector>
30 #include <algorithm>
31 
32 #include <cxxtest/TestSuite.h>
33 
34 #include "../../WAssert.h"
35 #include "../../WLimits.h"
36 #include "../../WStringUtils.h"
37 #include "../WTensorFunctions.h"
38 
39 /**
40  * Test class for some tensor functions.
41  */
42 class WTensorFunctionsTest : public CxxTest::TestSuite
43 {
44 public:
45  /**
46  * The eigenvalue of the symmetrical matrix:
47  * 0.000179516, 2.09569e-05, 2.76557e-06, 0.000170189, -5.52619e-07, 0.00015239
48  * (0.000196397;0.000155074;0.000150625)
49  */
50  void testSpecialSymMatrixEigenvalueTestCaseNumericalStability( void )
51  {
52  WTensorSym< 2, 3 > t;
53  t( 0, 0 ) = 0.000179516;
54  t( 0, 1 ) = 2.09569e-05;
55  t( 0, 2 ) = 2.76557e-06;
56  t( 1, 1 ) = 0.000170189;
57  t( 1, 2 ) = -5.52619e-07;
58  t( 2, 2 ) = 0.00015239;
59  RealEigenSystem sys;
60  jacobiEigenvector3D( t, &sys );
61 
62  TS_ASSERT_DELTA( sys[0].first, 1.5062467240725114e-04, 1e-9 );
63  TS_ASSERT_DELTA( sys[1].first, 1.5507354000104679e-04, 1e-9 );
64  TS_ASSERT_DELTA( sys[2].first, 1.9639678759170208e-04, 1e-9 );
65  }
66 
67  /**
68  * Test the jacobi eigenvector calculation.
69  */
70  void testJacobiEigenvectors()
71  {
72  // the test matrix
73  WTensorSym< 2, 3 > t;
74 
75  // simple diagonal matrices should be decomposed correctly
76  // 1 2 3
77  t( 0, 0 ) = 1.0;
78  t( 1, 1 ) = 2.0;
79  t( 2, 2 ) = 3.0;
80 
81  RealEigenSystem sys;
82  jacobiEigenvector3D( t, &sys );
83 
84  TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
85  TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
86  TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 );
87  // eigenvectors should be perpendicular
88  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
89  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
90  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
91  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
92  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
93  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
94 
95  // 1 2 -3
96  t( 0, 0 ) = 1.0;
97  t( 1, 1 ) = 2.0;
98  t( 2, 2 ) = -3.0;
99 
100  jacobiEigenvector3D( t, &sys );
101 
102  TS_ASSERT_DELTA( sys[0].first, -3.0, 1e-6 );
103  TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 );
104  TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 );
105  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
106  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
107  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
108  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
109  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
110  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
111 
112  // 1 2 2
113  t( 0, 0 ) = 1.0;
114  t( 1, 1 ) = 2.0;
115  t( 2, 2 ) = 2.0;
116 
117  jacobiEigenvector3D( t, &sys );
118 
119  TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
120  TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
121  TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 );
122  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
123  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
124  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
125  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
126  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
127  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
128 
129  // -1 -1 -1
130  t( 0, 0 ) = -1.0;
131  t( 1, 1 ) = -1.0;
132  t( 2, 2 ) = -1.0;
133 
134  jacobiEigenvector3D( t, &sys );
135 
136  TS_ASSERT_DELTA( sys[0].first, -1.0, 1e-6 );
137  TS_ASSERT_DELTA( sys[1].first, -1.0, 1e-6 );
138  TS_ASSERT_DELTA( sys[2].first, -1.0, 1e-6 );
139  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
140  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
141  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
142  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
143  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
144  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
145 
146  // 1 0 1
147  t( 0, 0 ) = 1.0;
148  t( 1, 1 ) = 0.0;
149  t( 2, 2 ) = 1.0;
150 
151  jacobiEigenvector3D( t, &sys );
152 
153  TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 );
154  TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 );
155  TS_ASSERT_DELTA( sys[2].first, 1.0, 1e-6 );
156  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
157  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
158  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
159  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
160  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
161  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
162 
163  // 0 0 0
164  t( 0, 0 ) = 0.0;
165  t( 1, 1 ) = 0.0;
166  t( 2, 2 ) = 0.0;
167 
168  jacobiEigenvector3D( t, &sys );
169 
170  TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 );
171  TS_ASSERT_DELTA( sys[1].first, 0.0, 1e-6 );
172  TS_ASSERT_DELTA( sys[2].first, 0.0, 1e-6 );
173  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
174  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
175  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
176  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
177  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
178  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
179 
180  // similar eigenvalues
181  // 2.000001 0.0 1.999998
182  t( 0, 0 ) = 2.000001;
183  t( 1, 1 ) = 0.0;
184  t( 2, 2 ) = 1.999998;
185 
186  jacobiEigenvector3D( t, &sys );
187 
188  TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 );
189  TS_ASSERT_DELTA( sys[1].first, 1.999998, 1e-6 );
190  TS_ASSERT_DELTA( sys[2].first, 2.000001, 1e-6 );
191  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
192  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
193  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
194  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
195  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
196  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
197 
198  // very large eigenvalues
199  // 3.824572321236e1000 1 2
200  t( 0, 0 ) = 3.824572321236e30;
201  t( 1, 1 ) = 1.0;
202  t( 2, 2 ) = 2.0;
203 
204  jacobiEigenvector3D( t, &sys );
205 
206  TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
207  TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
208  TS_ASSERT_DELTA( sys[2].first, 3.824572321236e30, 1e-6 );
209  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
210  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
211  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
212  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
213  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
214  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
215 
216  // very small eigenvalues
217  // 3.824572321236e-1000 1 2
218  t( 0, 0 ) = 3.824572321236e-30;
219  t( 1, 1 ) = 1.0;
220  t( 2, 2 ) = 2.0;
221 
222  jacobiEigenvector3D( t, &sys );
223 
224  TS_ASSERT_DELTA( sys[0].first, 3.824572321236e-30, 1e-6 );
225  TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 );
226  TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 );
227  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
228  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
229  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
230  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
231  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
232  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
233 
234  // some more sophisticated tests
235  // (using similarity transformations on diagonal matrices to create test cases)
236  t( 0, 0 ) = 1;
237  t( 1, 1 ) = 2;
238  t( 2, 2 ) = 3;
239 
240  // note that the jacobi eigenvector functions does not sort the eigenvalues and
241  // eigenvectors that were found
242  t = similarity_rotate_givens( t, 0, 2, 2.78 );
243 
244  jacobiEigenvector3D( t, &sys );
245 
246  TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
247  TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
248  TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 );
249  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
250  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
251  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
252  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
253  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
254  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
255  compare_results( t, sys );
256 
257  t = WTensorSym< 2, 3 >();
258  t( 0, 0 ) = 2;
259  t( 1, 1 ) = 1;
260  t( 2, 2 ) = 3;
261 
262  t = similarity_rotate_givens( t, 0, 2, 2.79 );
263  t = similarity_rotate_givens( t, 1, 2, -3.44 );
264  t = similarity_rotate_givens( t, 1, 0, -0.46 );
265  t = similarity_rotate_givens( t, 2, 1, 5.98 );
266 
267  jacobiEigenvector3D( t, &sys );
268 
269  TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
270  TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
271  TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 );
272  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
273  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
274  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
275  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
276  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
277  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
278  compare_results( t, sys );
279 
280  t = WTensorSym< 2, 3 >();
281  t( 0, 0 ) = 2;
282  t( 1, 1 ) = 2;
283  t( 2, 2 ) = 2;
284 
285  t = similarity_rotate_givens( t, 1, 2, -3.44 );
286  t = similarity_rotate_givens( t, 2, 1, 5.98 );
287  t = similarity_rotate_givens( t, 0, 2, 2.79 );
288  t = similarity_rotate_givens( t, 1, 0, -0.46 );
289 
290  jacobiEigenvector3D( t, &sys );
291 
292  TS_ASSERT_DELTA( sys[0].first, 2.0, 1e-6 );
293  TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
294  TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 );
295  TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
296  TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
297  TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
298  TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
299  TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
300  TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
301  compare_results( t, sys );
302  }
303 
304  /**
305  * Test the cardano eigenvalue calculation.
306  */
307  void testCardanoEigenvalues()
308  {
309  // the test matrix
310  WTensorSym< 2, 3 > t;
311  // variables for the output values
312  std::vector< double > d( 3 );
313 
314  // simple diagonal matrices should be decomposed correctly
315 
316  // 1 2 3
317  t( 0, 0 ) = 1.0;
318  t( 1, 1 ) = 2.0;
319  t( 2, 2 ) = 3.0;
320 
321  d = getEigenvaluesCardano( t );
322 
323  TS_ASSERT_DELTA( d[ 0 ], 3.0, 1e-6 );
324  TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 );
325  TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 );
326 
327  // 1 2 -3
328  t( 0, 0 ) = 1.0;
329  t( 1, 1 ) = 2.0;
330  t( 2, 2 ) = -3.0;
331 
332  d = getEigenvaluesCardano( t );
333 
334  TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 );
335  TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 );
336  TS_ASSERT_DELTA( d[ 2 ], -3.0, 1e-6 );
337 
338  // 1 2 2
339  t( 0, 0 ) = 1.0;
340  t( 1, 1 ) = 2.0;
341  t( 2, 2 ) = 2.0;
342 
343  d = getEigenvaluesCardano( t );
344 
345  TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 );
346  TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 );
347  TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 );
348 
349  // -1 -1 -1
350  t( 0, 0 ) = -1.0;
351  t( 1, 1 ) = -1.0;
352  t( 2, 2 ) = -1.0;
353 
354  d = getEigenvaluesCardano( t );
355 
356  TS_ASSERT_DELTA( d[ 0 ], -1.0, 1e-6 );
357  TS_ASSERT_DELTA( d[ 1 ], -1.0, 1e-6 );
358  TS_ASSERT_DELTA( d[ 2 ], -1.0, 1e-6 );
359 
360  // 1 0 1
361  t( 0, 0 ) = 1.0;
362  t( 1, 1 ) = 0.0;
363  t( 2, 2 ) = 1.0;
364 
365  d = getEigenvaluesCardano( t );
366 
367  TS_ASSERT_DELTA( d[ 0 ], 1.0, 1e-6 );
368  TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 );
369  TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 );
370 
371  // 0 0 0
372  t( 0, 0 ) = 0.0;
373  t( 1, 1 ) = 0.0;
374  t( 2, 2 ) = 0.0;
375 
376  d = getEigenvaluesCardano( t );
377 
378  TS_ASSERT_DELTA( d[ 0 ], 0.0, 1e-6 );
379  TS_ASSERT_DELTA( d[ 1 ], 0.0, 1e-6 );
380  TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 );
381 
382  // similar eigenvalues
383  // 2.000001 0.0 1.999998
384  t( 0, 0 ) = 2.000001;
385  t( 1, 1 ) = 0.0;
386  t( 2, 2 ) = 1.999998;
387 
388  d = getEigenvaluesCardano( t );
389 
390  TS_ASSERT_DELTA( d[ 0 ], 2.000001, 1e-6 );
391  TS_ASSERT_DELTA( d[ 1 ], 1.999998, 1e-6 );
392  TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 );
393 
394  // very large eigenvalues
395  // 3.824572321236e10 1 2
396  t( 0, 0 ) = 3.824572321236e10;
397  t( 1, 1 ) = 1.0;
398  t( 2, 2 ) = 2.0;
399 
400  d = getEigenvaluesCardano( t );
401 
402  TS_ASSERT_DELTA( d[ 0 ], 3.824572321236e10, 1e-6 );
403  TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 );
404  TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 );
405 
406  // very small eigenvalues
407  // 3.824572321236e-10 1 2
408  t( 0, 0 ) = 3.824572321236e-30;
409  t( 1, 1 ) = 1.0;
410  t( 2, 2 ) = 2.0;
411 
412  d = getEigenvaluesCardano( t );
413 
414  TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 );
415  TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 );
416  TS_ASSERT_DELTA( d[ 2 ], 3.824572321236e-30, 1e-6 );
417 
418  // some more sophisticated tests
419  // (using similarity transformations on diagonal matrices to create test cases)
420  t( 0, 0 ) = 1;
421  t( 1, 1 ) = 2;
422  t( 2, 2 ) = 3;
423 
424  // note that the jacobi eigenvector functions does not sort the eigenvalues and
425  // eigenvectors that were found
426  t = similarity_rotate_givens( t, 0, 2, 2.78 );
427 
428  d = getEigenvaluesCardano( t );
429 
430  TS_ASSERT_DELTA( d[ 0 ], 3.0, 1e-6 );
431  TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 );
432  TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 );
433 
434  // the eigenvalues calculated for this test are a bit off
435  //t = WTensorSym< 2, 3 >();
436  //t( 0, 0 ) = 2;
437  //t( 1, 1 ) = 1;
438  //t( 2, 2 ) = 3;
439 
440  //t = similarity_rotate_givens( t, 0, 2, 2.79 );
441  //t = similarity_rotate_givens( t, 1, 2, -3.44 );
442  //t = similarity_rotate_givens( t, 1, 0, -0.46 );
443  //t = similarity_rotate_givens( t, 2, 1, 5.98 );
444 
445  //d = getEigenvaluesCardano( t );
446 
447  //TS_ASSERT_DELTA( d[ 0 ], 3.0, 1e-6 );
448  //TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 );
449  //TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 );
450 
451  t = WTensorSym< 2, 3 >();
452  t( 0, 0 ) = 2;
453  t( 1, 1 ) = 2;
454  t( 2, 2 ) = 2;
455 
456  t = similarity_rotate_givens( t, 1, 2, -3.44 );
457  t = similarity_rotate_givens( t, 2, 1, 5.98 );
458  t = similarity_rotate_givens( t, 0, 2, 2.79 );
459  t = similarity_rotate_givens( t, 1, 0, -0.46 );
460 
461  d = getEigenvaluesCardano( t );
462 
463  TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 );
464  TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 );
465  TS_ASSERT_DELTA( d[ 2 ], 2.0, 1e-6 );
466  }
467 
468  /**
469  * Test if tensor log and exp functions behave correctly.
470  */
471  void testLogAndExp()
472  {
473  // init some tensor
474  WTensorSym< 2, 3, double > tens;
475  tens( 0, 0 ) = 8.0;
476  tens( 0, 1 ) = 10.0;
477  tens( 0, 2 ) = 4.0;
478  tens( 1, 1 ) = 17.0;
479  tens( 1, 2 ) = 14.0;
480  tens( 2, 2 ) = 20.0;
481 
482  WTensorSym< 2, 3, double > s = tensorExp( tensorLog( tens ) );
483 
484  TS_ASSERT_DELTA( s( 0, 0 ), tens( 0, 0 ), 1e-4 );
485  TS_ASSERT_DELTA( s( 0, 1 ), tens( 0, 1 ), 1e-4 );
486  TS_ASSERT_DELTA( s( 0, 2 ), tens( 0, 2 ), 1e-4 );
487  TS_ASSERT_DELTA( s( 1, 1 ), tens( 1, 1 ), 1e-4 );
488  TS_ASSERT_DELTA( s( 1, 2 ), tens( 1, 2 ), 1e-4 );
489  TS_ASSERT_DELTA( s( 2, 2 ), tens( 2, 2 ), 1e-4 );
490  }
491 
492 private:
493  /**
494  * A helper function performing a similarity transform using a givens rotation.
495  *
496  * \param m The symmetric tensor to transform.
497  * \param i A row index.
498  * \param j A column index.
499  * \param angle The rotation angle (in radians).
500  *
501  * \note i must not have the same values as j
502  *
503  * \return The new tensor
504  */
505  template< std::size_t dim, typename Data_T >
506  WTensorSym< 2, dim, Data_T > similarity_rotate_givens( WTensorSym< 2, dim, Data_T > const& m,
507  std::size_t i,
508  std::size_t j,
509  double angle )
510  {
511  WAssert( i != j, "" );
512 
513  double s = sin( angle );
514  double c = cos( angle );
515  WTensor< 2, dim, Data_T > r;
516  for( std::size_t k = 0; k < dim; ++k )
517  {
518  if( k == i || k == j )
519  {
520  r( k, k ) = c;
521  }
522  else
523  {
524  r( k, k ) = 1.0;
525  }
526  }
527  r( i, j ) = s;
528  r( j, i ) = -s;
529  WTensor< 2, dim, Data_T > t( r );
530  std::swap( t( i, j ), t( j, i ) );
531 
532  r = t * m * r;
533  WTensorSym< 2, dim, Data_T > res;
534  res( 0, 0 ) = r( 0, 0 );
535  res( 1, 0 ) = r( 1, 0 );
536  res( 2, 0 ) = r( 2, 0 );
537  res( 1, 1 ) = r( 1, 1 );
538  res( 2, 1 ) = r( 2, 1 );
539  res( 2, 2 ) = r( 2, 2 );
540 
541  return res;
542  }
543 
544  /**
545  * Test if the given vectors are eigenvectors to the given eigenvalues of a
546  * symmetric matrix.
547  *
548  * \param m A symmetric matrix.
549  * \param sys The eigen system ( eigenvalues and eigenvectors )
550  */
551  template< std::size_t dim, typename Data_T >
552  void compare_results( WTensorSym< 2, dim, Data_T > const& m, RealEigenSystem const& sys )
553  {
554  WTensor< 2, dim, Data_T > t, r;
555  for( std::size_t i = 0; i < dim; ++i )
556  {
557  for( std::size_t j = 0; j < dim; ++j )
558  {
559  t( j, i ) = sys[i].second[j];
560  r( j, i ) = sys[i].first * sys[i].second[j];
561  }
562  }
563 
564  t = m * t;
565  for( std::size_t i = 0; i < dim; ++i )
566  {
567  for( std::size_t j = 0; j < dim; ++j )
568  {
569  TS_ASSERT_DELTA( t( i, j ), r( i, j ), 1e-15 );
570  }
571  }
572  }
573 };
574 
575 /**
576  * Test class for all tensor operators.
577  */
578 class WTensorOperatorsTest : public CxxTest::TestSuite
579 {
580 public:
581  /**
582  * Test order 2 tensor multiplication.
583  */
584  void testMultiplyTensorsOperator()
585  {
586  // some special cases, WTensor * WTensor
587  TS_ASSERT_EQUALS( zero * zero, zero );
588  TS_ASSERT_EQUALS( zero * one, zero );
589  TS_ASSERT_EQUALS( one * zero, zero );
590  TS_ASSERT_EQUALS( one * one, one );
591  TS_ASSERT_EQUALS( one * rdm1, rdm1 );
592  TS_ASSERT_EQUALS( one * rdm2, rdm2 );
593  TS_ASSERT_EQUALS( rdm1 * one, rdm1 );
594  TS_ASSERT_EQUALS( rdm2 * one, rdm2 );
595 
596  // a normal case
597  TS_ASSERT_EQUALS( rdm1 * rdm2, res1 );
598 
599  // some special cases, WTensorSym * WTensorSym
600  TS_ASSERT_EQUALS( szero * szero, zero );
601  TS_ASSERT_EQUALS( szero * sone, zero );
602  TS_ASSERT_EQUALS( sone * szero, zero );
603  TS_ASSERT_EQUALS( sone * sone, one );
604  TS_ASSERT_EQUALS( sone * srdm1, res2 );
605  TS_ASSERT_EQUALS( srdm1 * sone, res2 );
606  TS_ASSERT_EQUALS( srdm1 * srdm2, res3 );
607 
608  // WTensor * WTensorSym
609  TS_ASSERT_EQUALS( zero * szero, zero );
610  TS_ASSERT_EQUALS( one * sone, one );
611  TS_ASSERT_EQUALS( zero * sone, zero );
612  TS_ASSERT_EQUALS( one * szero, zero );
613 
614  // WTensorSym * WTensor
615  TS_ASSERT_EQUALS( szero * zero, zero );
616  TS_ASSERT_EQUALS( sone * one, one );
617  TS_ASSERT_EQUALS( szero * one, zero );
618  TS_ASSERT_EQUALS( sone * zero, zero );
619  TS_ASSERT_EQUALS( srdm1 * rdm1, res4 );
620  }
621 
622  /**
623  * The optimizations for symmetric tensors should not corrupt the result.
624  */
625  void testEvaluateSphericalFunction()
626  {
627  WTensorSym< 4, 3, double > t;
628  // the tensor
629  t( 0, 0, 0, 0 ) = 2.5476;
630  t( 1, 1, 1, 1 ) = 3.5476;
631  t( 2, 2, 2, 2 ) = 4.5476;
632  t( 0, 0, 0, 1 ) = 5.5476;
633  t( 0, 0, 0, 2 ) = 6.5476;
634  t( 1, 1, 1, 0 ) = 7.5476;
635  t( 1, 1, 1, 2 ) = 8.5476;
636  t( 2, 2, 2, 0 ) = 9.5476;
637  t( 2, 2, 2, 1 ) = 10.5476;
638  t( 0, 0, 1, 2 ) = 11.5476;
639  t( 1, 1, 0, 2 ) = 12.5476;
640  t( 2, 2, 0, 1 ) = 13.5476;
641  t( 0, 0, 1, 1 ) = 14.5476;
642  t( 0, 0, 2, 2 ) = 15.5476;
643  t( 1, 1, 2, 2 ) = 16.5476;
644 
645  // the gradients
646  std::vector< WVector3d > gradients;
647  gradients.push_back( WVector3d( 1.0, 0.0, 0.0 ) );
648  gradients.push_back( WVector3d( 0.0, 1.0, 0.0 ) );
649  gradients.push_back( normalize( WVector3d( 1.0, 1.0, 0.0 ) ) );
650  gradients.push_back( normalize( WVector3d( 0.3, 0.4, 0.5 ) ) );
651  gradients.push_back( normalize( WVector3d( -7.0, 3.0, -1.0 ) ) );
652 
653  for( int k = 0; k < 5; ++k )
654  {
655  double res = calcTens( t, gradients[ k ] );
656  TS_ASSERT_DELTA( res, evaluateSphericalFunction( t, gradients[ k ] ), 0.001 );
657  }
658  }
659 
660 private:
661  /**
662  * Initialize a lot of tensors.
663  */
664  void setUp()
665  {
666  one( 0, 0 ) = one( 1, 1 ) = one( 2, 2 ) = 1.0;
667  sone( 0, 0 ) = sone( 1, 1 ) = sone( 2, 2 ) = 1.0;
668 
669  rdm1( 0, 0 ) = 2;
670  rdm1( 0, 1 ) = 3;
671  rdm1( 0, 2 ) = 1;
672  rdm1( 1, 0 ) = 0;
673  rdm1( 1, 1 ) = 4;
674  rdm1( 1, 2 ) = 0;
675  rdm1( 2, 0 ) = 5;
676  rdm1( 2, 1 ) = -3;
677  rdm1( 2, 2 ) = -7;
678 
679  rdm2( 0, 0 ) = -1;
680  rdm2( 0, 1 ) = -1;
681  rdm2( 0, 2 ) = -1;
682  rdm2( 1, 0 ) = 3;
683  rdm2( 1, 1 ) = 0;
684  rdm2( 1, 2 ) = -2;
685  rdm2( 2, 0 ) = -1;
686  rdm2( 2, 1 ) = -2;
687  rdm2( 2, 2 ) = -3;
688 
689  res1( 0, 0 ) = 6;
690  res1( 0, 1 ) = -4;
691  res1( 0, 2 ) = -11;
692  res1( 1, 0 ) = 12;
693  res1( 1, 1 ) = 0;
694  res1( 1, 2 ) = -8;
695  res1( 2, 0 ) = -7;
696  res1( 2, 1 ) = 9;
697  res1( 2, 2 ) = 22;
698 
699  srdm1( 0, 0 ) = 2;
700  srdm1( 0, 1 ) = 3;
701  srdm1( 0, 2 ) = 1;
702  srdm1( 1, 1 ) = -1;
703  srdm1( 1, 2 ) = 4;
704  srdm1( 2, 2 ) = -3;
705 
706  srdm2( 0, 0 ) = -2;
707  srdm2( 0, 1 ) = 5;
708  srdm2( 0, 2 ) = 2;
709  srdm2( 1, 1 ) = 0;
710  srdm2( 1, 2 ) = -3;
711  srdm2( 2, 2 ) = 1;
712 
713  // copy symmetric tensor srdm1 to an asymmetric tensor
714  res2 = srdm1;
715 
716  res3( 0, 0 ) = 13;
717  res3( 0, 1 ) = 7;
718  res3( 0, 2 ) = -4;
719  res3( 1, 0 ) = -3;
720  res3( 1, 1 ) = 3;
721  res3( 1, 2 ) = 13;
722  res3( 2, 0 ) = 12;
723  res3( 2, 1 ) = 14;
724  res3( 2, 2 ) = -13;
725 
726  res4( 0, 0 ) = 9;
727  res4( 0, 1 ) = 15;
728  res4( 0, 2 ) = -5;
729  res4( 1, 0 ) = 26;
730  res4( 1, 1 ) = -7;
731  res4( 1, 2 ) = -25;
732  res4( 2, 0 ) = -13;
733  res4( 2, 1 ) = 28;
734  res4( 2, 2 ) = 22;
735  }
736 
737  /**
738  * A helper function that implements the simple approach to tensor evaluation.
739  *
740  * \param t The tensor.
741  * \param v The gradient.
742  *
743  * \return value
744  */
745  double calcTens( WTensorSym< 4, 3, double > const& t, WVector3d const& v )
746  {
747  double res = 0.0;
748  for( int a = 0; a < 3; ++a )
749  {
750  for( int b = 0; b < 3; ++b )
751  {
752  for( int c = 0; c < 3; ++c )
753  {
754  for( int d = 0; d < 3; ++d )
755  {
756  res += v[ a ] * v[ b ] * v[ c ] * v[ d ] * t( a, b, c, d );
757  }
758  }
759  }
760  }
761  return res;
762  }
763 
764  //! a test tensor
765  WTensor< 2, 3, int > one;
766  //! a test tensor
767  WTensor< 2, 3, int > zero;
768  //! a test tensor
769  WTensor< 2, 3, int > rdm1;
770  //! a test tensor
771  WTensor< 2, 3, int > rdm2;
772  //! a test tensor
773  WTensor< 2, 3, int > res1;
774  //! a test tensor
775  WTensor< 2, 3, int > res2;
776  //! a test tensor
777  WTensor< 2, 3, int > res3;
778  //! a test tensor
779  WTensor< 2, 3, int > res4;
780  //! a test tensor
781  WTensorSym< 2, 3, int > sone;
782  //! a test tensor
783  WTensorSym< 2, 3, int > szero;
784  //! a test tensor
785  WTensorSym< 2, 3, int > srdm1;
786  //! a test tensor
787  WTensorSym< 2, 3, int > srdm2;
788 };
789 
790 #endif // WTENSORFUNCTIONS_TEST_H
WTensorFunctionsTest
Test class for some tensor functions.
Definition: WTensorFunctions_test.h:43
WTensorFunctionsTest::similarity_rotate_givens
WTensorSym< 2, dim, Data_T > similarity_rotate_givens(WTensorSym< 2, dim, Data_T > const &m, std::size_t i, std::size_t j, double angle)
A helper function performing a similarity transform using a givens rotation.
Definition: WTensorFunctions_test.h:506
WTensorFunctionsTest::testJacobiEigenvectors
void testJacobiEigenvectors()
Test the jacobi eigenvector calculation.
Definition: WTensorFunctions_test.h:70
WTensorFunctionsTest::testCardanoEigenvalues
void testCardanoEigenvalues()
Test the cardano eigenvalue calculation.
Definition: WTensorFunctions_test.h:307
WTensorFunctionsTest::testSpecialSymMatrixEigenvalueTestCaseNumericalStability
void testSpecialSymMatrixEigenvalueTestCaseNumericalStability(void)
The eigenvalue of the symmetrical matrix: 0.000179516, 2.09569e-05, 2.76557e-06, 0....
Definition: WTensorFunctions_test.h:50
WTensorFunctionsTest::compare_results
void compare_results(WTensorSym< 2, dim, Data_T > const &m, RealEigenSystem const &sys)
Test if the given vectors are eigenvectors to the given eigenvalues of a symmetric matrix.
Definition: WTensorFunctions_test.h:552
WTensorFunctionsTest::testLogAndExp
void testLogAndExp()
Test if tensor log and exp functions behave correctly.
Definition: WTensorFunctions_test.h:471
WTensorOperatorsTest
Test class for all tensor operators.
Definition: WTensorFunctions_test.h:579
WTensorOperatorsTest::zero
WTensor< 2, 3, int > zero
a test tensor
Definition: WTensorFunctions_test.h:767
WTensorOperatorsTest::res1
WTensor< 2, 3, int > res1
a test tensor
Definition: WTensorFunctions_test.h:773
WTensorOperatorsTest::rdm1
WTensor< 2, 3, int > rdm1
a test tensor
Definition: WTensorFunctions_test.h:769
WTensorOperatorsTest::res2
WTensor< 2, 3, int > res2
a test tensor
Definition: WTensorFunctions_test.h:775
WTensorOperatorsTest::szero
WTensorSym< 2, 3, int > szero
a test tensor
Definition: WTensorFunctions_test.h:783
WTensorOperatorsTest::calcTens
double calcTens(WTensorSym< 4, 3, double > const &t, WVector3d const &v)
A helper function that implements the simple approach to tensor evaluation.
Definition: WTensorFunctions_test.h:745
WTensorOperatorsTest::srdm1
WTensorSym< 2, 3, int > srdm1
a test tensor
Definition: WTensorFunctions_test.h:785
WTensorOperatorsTest::res4
WTensor< 2, 3, int > res4
a test tensor
Definition: WTensorFunctions_test.h:779
WTensorOperatorsTest::rdm2
WTensor< 2, 3, int > rdm2
a test tensor
Definition: WTensorFunctions_test.h:771
WTensorOperatorsTest::testMultiplyTensorsOperator
void testMultiplyTensorsOperator()
Test order 2 tensor multiplication.
Definition: WTensorFunctions_test.h:584
WTensorOperatorsTest::one
WTensor< 2, 3, int > one
a test tensor
Definition: WTensorFunctions_test.h:765
WTensorOperatorsTest::testEvaluateSphericalFunction
void testEvaluateSphericalFunction()
The optimizations for symmetric tensors should not corrupt the result.
Definition: WTensorFunctions_test.h:625
WTensorOperatorsTest::setUp
void setUp()
Initialize a lot of tensors.
Definition: WTensorFunctions_test.h:664
WTensorOperatorsTest::res3
WTensor< 2, 3, int > res3
a test tensor
Definition: WTensorFunctions_test.h:777
WTensorOperatorsTest::sone
WTensorSym< 2, 3, int > sone
a test tensor
Definition: WTensorFunctions_test.h:781
WTensorOperatorsTest::srdm2
WTensorSym< 2, 3, int > srdm2
a test tensor
Definition: WTensorFunctions_test.h:787
WTensorSym
Implements a symmetric tensor that has the same number of components in every direction.
Definition: WTensorSym.h:73
WTensor
Implements a tensor that has the same number of components in every direction.
Definition: WTensor.h:79
Generated by doxygen 1.9.1