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_TEST_H
00026 #define WTENSORFUNCTIONS_TEST_H
00027
00028 #include <string>
00029 #include <vector>
00030 #include <algorithm>
00031
00032 #include <cxxtest/TestSuite.h>
00033
00034 #include "../../WAssert.h"
00035 #include "../../WLimits.h"
00036 #include "../../WStringUtils.h"
00037 #include "../WTensorFunctions.h"
00038
00039
00040
00041
00042 class WTensorFunctionsTest : public CxxTest::TestSuite
00043 {
00044 public:
00045
00046
00047
00048
00049
00050 void testSpecialSymMatrixEigenvalueTestCaseNumericalStability( void )
00051 {
00052 WTensorSym< 2, 3 > t;
00053 t( 0, 0 ) = 0.000179516;
00054 t( 0, 1 ) = 2.09569e-05;
00055 t( 0, 2 ) = 2.76557e-06;
00056 t( 1, 1 ) = 0.000170189;
00057 t( 1, 2 ) = -5.52619e-07;
00058 t( 2, 2 ) = 0.00015239;
00059 RealEigenSystem sys;
00060 jacobiEigenvector3D( t, &sys );
00061
00062 TS_ASSERT_DELTA( sys[0].first, 1.5062467240725114e-04, 1e-9 );
00063 TS_ASSERT_DELTA( sys[1].first, 1.5507354000104679e-04, 1e-9 );
00064 TS_ASSERT_DELTA( sys[2].first, 1.9639678759170208e-04, 1e-9 );
00065 }
00066
00067
00068
00069
00070 void testJacobiEigenvectors()
00071 {
00072
00073 WTensorSym< 2, 3 > t;
00074
00075
00076
00077 t( 0, 0 ) = 1.0;
00078 t( 1, 1 ) = 2.0;
00079 t( 2, 2 ) = 3.0;
00080
00081 RealEigenSystem sys;
00082 jacobiEigenvector3D( t, &sys );
00083
00084 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
00085 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
00086 TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 );
00087
00088 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00089 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00090 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00091 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00092 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00093 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00094
00095
00096 t( 0, 0 ) = 1.0;
00097 t( 1, 1 ) = 2.0;
00098 t( 2, 2 ) = -3.0;
00099
00100 jacobiEigenvector3D( t, &sys );
00101
00102 TS_ASSERT_DELTA( sys[0].first, -3.0, 1e-6 );
00103 TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 );
00104 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 );
00105 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00106 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00107 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00108 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00109 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00110 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00111
00112
00113 t( 0, 0 ) = 1.0;
00114 t( 1, 1 ) = 2.0;
00115 t( 2, 2 ) = 2.0;
00116
00117 jacobiEigenvector3D( t, &sys );
00118
00119 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
00120 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
00121 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 );
00122 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00123 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00124 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00125 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00126 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00127 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00128
00129
00130 t( 0, 0 ) = -1.0;
00131 t( 1, 1 ) = -1.0;
00132 t( 2, 2 ) = -1.0;
00133
00134 jacobiEigenvector3D( t, &sys );
00135
00136 TS_ASSERT_DELTA( sys[0].first, -1.0, 1e-6 );
00137 TS_ASSERT_DELTA( sys[1].first, -1.0, 1e-6 );
00138 TS_ASSERT_DELTA( sys[2].first, -1.0, 1e-6 );
00139 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00140 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00141 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00142 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00143 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00144 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00145
00146
00147 t( 0, 0 ) = 1.0;
00148 t( 1, 1 ) = 0.0;
00149 t( 2, 2 ) = 1.0;
00150
00151 jacobiEigenvector3D( t, &sys );
00152
00153 TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 );
00154 TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 );
00155 TS_ASSERT_DELTA( sys[2].first, 1.0, 1e-6 );
00156 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00157 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00158 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00159 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00160 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00161 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00162
00163
00164 t( 0, 0 ) = 0.0;
00165 t( 1, 1 ) = 0.0;
00166 t( 2, 2 ) = 0.0;
00167
00168 jacobiEigenvector3D( t, &sys );
00169
00170 TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 );
00171 TS_ASSERT_DELTA( sys[1].first, 0.0, 1e-6 );
00172 TS_ASSERT_DELTA( sys[2].first, 0.0, 1e-6 );
00173 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00174 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00175 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00176 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00177 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00178 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00179
00180
00181
00182 t( 0, 0 ) = 2.000001;
00183 t( 1, 1 ) = 0.0;
00184 t( 2, 2 ) = 1.999998;
00185
00186 jacobiEigenvector3D( t, &sys );
00187
00188 TS_ASSERT_DELTA( sys[0].first, 0.0, 1e-6 );
00189 TS_ASSERT_DELTA( sys[1].first, 1.999998, 1e-6 );
00190 TS_ASSERT_DELTA( sys[2].first, 2.000001, 1e-6 );
00191 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00192 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00193 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00194 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00195 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00196 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00197
00198
00199
00200 t( 0, 0 ) = 3.824572321236e30;
00201 t( 1, 1 ) = 1.0;
00202 t( 2, 2 ) = 2.0;
00203
00204 jacobiEigenvector3D( t, &sys );
00205
00206 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
00207 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
00208 TS_ASSERT_DELTA( sys[2].first, 3.824572321236e30, 1e-6 );
00209 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00210 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00211 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00212 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00213 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00214 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00215
00216
00217
00218 t( 0, 0 ) = 3.824572321236e-30;
00219 t( 1, 1 ) = 1.0;
00220 t( 2, 2 ) = 2.0;
00221
00222 jacobiEigenvector3D( t, &sys );
00223
00224 TS_ASSERT_DELTA( sys[0].first, 3.824572321236e-30, 1e-6 );
00225 TS_ASSERT_DELTA( sys[1].first, 1.0, 1e-6 );
00226 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 );
00227 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00228 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00229 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00230 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00231 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00232 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00233
00234
00235
00236 t( 0, 0 ) = 1;
00237 t( 1, 1 ) = 2;
00238 t( 2, 2 ) = 3;
00239
00240
00241
00242 t = similarity_rotate_givens( t, 0, 2, 2.78 );
00243
00244 jacobiEigenvector3D( t, &sys );
00245
00246 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
00247 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
00248 TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 );
00249 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00250 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00251 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00252 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00253 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00254 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00255 compare_results( t, sys );
00256
00257 t = WTensorSym< 2, 3 >();
00258 t( 0, 0 ) = 2;
00259 t( 1, 1 ) = 1;
00260 t( 2, 2 ) = 3;
00261
00262 t = similarity_rotate_givens( t, 0, 2, 2.79 );
00263 t = similarity_rotate_givens( t, 1, 2, -3.44 );
00264 t = similarity_rotate_givens( t, 1, 0, -0.46 );
00265 t = similarity_rotate_givens( t, 2, 1, 5.98 );
00266
00267 jacobiEigenvector3D( t, &sys );
00268
00269 TS_ASSERT_DELTA( sys[0].first, 1.0, 1e-6 );
00270 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
00271 TS_ASSERT_DELTA( sys[2].first, 3.0, 1e-6 );
00272 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00273 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00274 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00275 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00276 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00277 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00278 compare_results( t, sys );
00279
00280 t = WTensorSym< 2, 3 >();
00281 t( 0, 0 ) = 2;
00282 t( 1, 1 ) = 2;
00283 t( 2, 2 ) = 2;
00284
00285 t = similarity_rotate_givens( t, 1, 2, -3.44 );
00286 t = similarity_rotate_givens( t, 2, 1, 5.98 );
00287 t = similarity_rotate_givens( t, 0, 2, 2.79 );
00288 t = similarity_rotate_givens( t, 1, 0, -0.46 );
00289
00290 jacobiEigenvector3D( t, &sys );
00291
00292 TS_ASSERT_DELTA( sys[0].first, 2.0, 1e-6 );
00293 TS_ASSERT_DELTA( sys[1].first, 2.0, 1e-6 );
00294 TS_ASSERT_DELTA( sys[2].first, 2.0, 1e-6 );
00295 TS_ASSERT_DELTA( dot( sys[0].second, sys[1].second ), 0.0, 1e-9 );
00296 TS_ASSERT_DELTA( dot( sys[1].second, sys[2].second ), 0.0, 1e-9 );
00297 TS_ASSERT_DELTA( dot( sys[2].second, sys[0].second ), 0.0, 1e-9 );
00298 TS_ASSERT_DELTA( length( sys[0].second ), 1.0, 1e-9 );
00299 TS_ASSERT_DELTA( length( sys[1].second ), 1.0, 1e-9 );
00300 TS_ASSERT_DELTA( length( sys[2].second ), 1.0, 1e-9 );
00301 compare_results( t, sys );
00302 }
00303
00304
00305
00306
00307 void testCardanoEigenvalues()
00308 {
00309
00310 WTensorSym< 2, 3 > t;
00311
00312 std::vector< double > d( 3 );
00313
00314
00315
00316
00317 t( 0, 0 ) = 1.0;
00318 t( 1, 1 ) = 2.0;
00319 t( 2, 2 ) = 3.0;
00320
00321 d = getEigenvaluesCardano( t );
00322
00323 TS_ASSERT_DELTA( d[ 0 ], 3.0, 1e-6 );
00324 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 );
00325 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 );
00326
00327
00328 t( 0, 0 ) = 1.0;
00329 t( 1, 1 ) = 2.0;
00330 t( 2, 2 ) = -3.0;
00331
00332 d = getEigenvaluesCardano( t );
00333
00334 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 );
00335 TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 );
00336 TS_ASSERT_DELTA( d[ 2 ], -3.0, 1e-6 );
00337
00338
00339 t( 0, 0 ) = 1.0;
00340 t( 1, 1 ) = 2.0;
00341 t( 2, 2 ) = 2.0;
00342
00343 d = getEigenvaluesCardano( t );
00344
00345 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 );
00346 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 );
00347 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 );
00348
00349
00350 t( 0, 0 ) = -1.0;
00351 t( 1, 1 ) = -1.0;
00352 t( 2, 2 ) = -1.0;
00353
00354 d = getEigenvaluesCardano( t );
00355
00356 TS_ASSERT_DELTA( d[ 0 ], -1.0, 1e-6 );
00357 TS_ASSERT_DELTA( d[ 1 ], -1.0, 1e-6 );
00358 TS_ASSERT_DELTA( d[ 2 ], -1.0, 1e-6 );
00359
00360
00361 t( 0, 0 ) = 1.0;
00362 t( 1, 1 ) = 0.0;
00363 t( 2, 2 ) = 1.0;
00364
00365 d = getEigenvaluesCardano( t );
00366
00367 TS_ASSERT_DELTA( d[ 0 ], 1.0, 1e-6 );
00368 TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 );
00369 TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 );
00370
00371
00372 t( 0, 0 ) = 0.0;
00373 t( 1, 1 ) = 0.0;
00374 t( 2, 2 ) = 0.0;
00375
00376 d = getEigenvaluesCardano( t );
00377
00378 TS_ASSERT_DELTA( d[ 0 ], 0.0, 1e-6 );
00379 TS_ASSERT_DELTA( d[ 1 ], 0.0, 1e-6 );
00380 TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 );
00381
00382
00383
00384 t( 0, 0 ) = 2.000001;
00385 t( 1, 1 ) = 0.0;
00386 t( 2, 2 ) = 1.999998;
00387
00388 d = getEigenvaluesCardano( t );
00389
00390 TS_ASSERT_DELTA( d[ 0 ], 2.000001, 1e-6 );
00391 TS_ASSERT_DELTA( d[ 1 ], 1.999998, 1e-6 );
00392 TS_ASSERT_DELTA( d[ 2 ], 0.0, 1e-6 );
00393
00394
00395
00396 t( 0, 0 ) = 3.824572321236e10;
00397 t( 1, 1 ) = 1.0;
00398 t( 2, 2 ) = 2.0;
00399
00400 d = getEigenvaluesCardano( t );
00401
00402 TS_ASSERT_DELTA( d[ 0 ], 3.824572321236e10, 1e-6 );
00403 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 );
00404 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 );
00405
00406
00407
00408 t( 0, 0 ) = 3.824572321236e-30;
00409 t( 1, 1 ) = 1.0;
00410 t( 2, 2 ) = 2.0;
00411
00412 d = getEigenvaluesCardano( t );
00413
00414 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 );
00415 TS_ASSERT_DELTA( d[ 1 ], 1.0, 1e-6 );
00416 TS_ASSERT_DELTA( d[ 2 ], 3.824572321236e-30, 1e-6 );
00417
00418
00419
00420 t( 0, 0 ) = 1;
00421 t( 1, 1 ) = 2;
00422 t( 2, 2 ) = 3;
00423
00424
00425
00426 t = similarity_rotate_givens( t, 0, 2, 2.78 );
00427
00428 d = getEigenvaluesCardano( t );
00429
00430 TS_ASSERT_DELTA( d[ 0 ], 3.0, 1e-6 );
00431 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 );
00432 TS_ASSERT_DELTA( d[ 2 ], 1.0, 1e-6 );
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451 t = WTensorSym< 2, 3 >();
00452 t( 0, 0 ) = 2;
00453 t( 1, 1 ) = 2;
00454 t( 2, 2 ) = 2;
00455
00456 t = similarity_rotate_givens( t, 1, 2, -3.44 );
00457 t = similarity_rotate_givens( t, 2, 1, 5.98 );
00458 t = similarity_rotate_givens( t, 0, 2, 2.79 );
00459 t = similarity_rotate_givens( t, 1, 0, -0.46 );
00460
00461 d = getEigenvaluesCardano( t );
00462
00463 TS_ASSERT_DELTA( d[ 0 ], 2.0, 1e-6 );
00464 TS_ASSERT_DELTA( d[ 1 ], 2.0, 1e-6 );
00465 TS_ASSERT_DELTA( d[ 2 ], 2.0, 1e-6 );
00466 }
00467
00468 private:
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481 template< std::size_t dim, typename Data_T >
00482 WTensorSym< 2, dim, Data_T > similarity_rotate_givens( WTensorSym< 2, dim, Data_T > const& m,
00483 std::size_t i,
00484 std::size_t j,
00485 double angle )
00486 {
00487 WAssert( i != j, "" );
00488
00489 double s = sin( angle );
00490 double c = cos( angle );
00491 WTensor< 2, dim, Data_T > r;
00492 for( std::size_t k = 0; k < dim; ++k )
00493 {
00494 if( k == i || k == j )
00495 {
00496 r( k, k ) = c;
00497 }
00498 else
00499 {
00500 r( k, k ) = 1.0;
00501 }
00502 }
00503 r( i, j ) = s;
00504 r( j, i ) = -s;
00505 WTensor< 2, dim, Data_T > t( r );
00506 std::swap( t( i, j ), t( j, i ) );
00507
00508 r = t * m * r;
00509 WTensorSym< 2, dim, Data_T > res;
00510 res( 0, 0 ) = r( 0, 0 );
00511 res( 1, 0 ) = r( 1, 0 );
00512 res( 2, 0 ) = r( 2, 0 );
00513 res( 1, 1 ) = r( 1, 1 );
00514 res( 2, 1 ) = r( 2, 1 );
00515 res( 2, 2 ) = r( 2, 2 );
00516
00517 return res;
00518 }
00519
00520
00521
00522
00523
00524
00525
00526
00527 template< std::size_t dim, typename Data_T >
00528 void compare_results( WTensorSym< 2, dim, Data_T > const& m, RealEigenSystem const& sys )
00529 {
00530 WTensor< 2, dim, Data_T > t, r;
00531 for( std::size_t i = 0; i < dim; ++i )
00532 {
00533 for( std::size_t j = 0; j < dim; ++j )
00534 {
00535 t( j, i ) = sys[i].second[j];
00536 r( j, i ) = sys[i].first * sys[i].second[j];
00537 }
00538 }
00539
00540 t = m * t;
00541 for( std::size_t i = 0; i < dim; ++i )
00542 {
00543 for( std::size_t j = 0; j < dim; ++j )
00544 {
00545 TS_ASSERT_DELTA( t( i, j ), r( i, j ), 1e-15 );
00546 }
00547 }
00548 }
00549 };
00550
00551
00552
00553
00554 class WTensorOperatorsTest : public CxxTest::TestSuite
00555 {
00556 public:
00557
00558
00559
00560 void testMultiplyTensorsOperator()
00561 {
00562
00563 TS_ASSERT_EQUALS( zero * zero, zero );
00564 TS_ASSERT_EQUALS( zero * one, zero );
00565 TS_ASSERT_EQUALS( one * zero, zero );
00566 TS_ASSERT_EQUALS( one * one, one );
00567 TS_ASSERT_EQUALS( one * rdm1, rdm1 );
00568 TS_ASSERT_EQUALS( one * rdm2, rdm2 );
00569 TS_ASSERT_EQUALS( rdm1 * one, rdm1 );
00570 TS_ASSERT_EQUALS( rdm2 * one, rdm2 );
00571
00572
00573 TS_ASSERT_EQUALS( rdm1 * rdm2, res1 );
00574
00575
00576 TS_ASSERT_EQUALS( szero * szero, zero );
00577 TS_ASSERT_EQUALS( szero * sone, zero );
00578 TS_ASSERT_EQUALS( sone * szero, zero );
00579 TS_ASSERT_EQUALS( sone * sone, one );
00580 TS_ASSERT_EQUALS( sone * srdm1, res2 );
00581 TS_ASSERT_EQUALS( srdm1 * sone, res2 );
00582 TS_ASSERT_EQUALS( srdm1 * srdm2, res3 );
00583
00584
00585 TS_ASSERT_EQUALS( zero * szero, zero );
00586 TS_ASSERT_EQUALS( one * sone, one );
00587 TS_ASSERT_EQUALS( zero * sone, zero );
00588 TS_ASSERT_EQUALS( one * szero, zero );
00589
00590
00591 TS_ASSERT_EQUALS( szero * zero, zero );
00592 TS_ASSERT_EQUALS( sone * one, one );
00593 TS_ASSERT_EQUALS( szero * one, zero );
00594 TS_ASSERT_EQUALS( sone * zero, zero );
00595 TS_ASSERT_EQUALS( srdm1 * rdm1, res4 );
00596 }
00597
00598
00599
00600
00601 void testEvaluateSphericalFunction()
00602 {
00603 WTensorSym< 4, 3, double > t;
00604
00605 t( 0, 0, 0, 0 ) = 2.5476;
00606 t( 1, 1, 1, 1 ) = 3.5476;
00607 t( 2, 2, 2, 2 ) = 4.5476;
00608 t( 0, 0, 0, 1 ) = 5.5476;
00609 t( 0, 0, 0, 2 ) = 6.5476;
00610 t( 1, 1, 1, 0 ) = 7.5476;
00611 t( 1, 1, 1, 2 ) = 8.5476;
00612 t( 2, 2, 2, 0 ) = 9.5476;
00613 t( 2, 2, 2, 1 ) = 10.5476;
00614 t( 0, 0, 1, 2 ) = 11.5476;
00615 t( 1, 1, 0, 2 ) = 12.5476;
00616 t( 2, 2, 0, 1 ) = 13.5476;
00617 t( 0, 0, 1, 1 ) = 14.5476;
00618 t( 0, 0, 2, 2 ) = 15.5476;
00619 t( 1, 1, 2, 2 ) = 16.5476;
00620
00621
00622 std::vector< WVector3d > gradients;
00623 gradients.push_back( WVector3d( 1.0, 0.0, 0.0 ) );
00624 gradients.push_back( WVector3d( 0.0, 1.0, 0.0 ) );
00625 gradients.push_back( normalize( WVector3d( 1.0, 1.0, 0.0 ) ) );
00626 gradients.push_back( normalize( WVector3d( 0.3, 0.4, 0.5 ) ) );
00627 gradients.push_back( normalize( WVector3d( -7.0, 3.0, -1.0 ) ) );
00628
00629 for( int k = 0; k < 5; ++k )
00630 {
00631 double res = calcTens( t, gradients[ k ] );
00632 TS_ASSERT_DELTA( res, evaluateSphericalFunction( t, gradients[ k ] ), 0.001 );
00633 }
00634 }
00635
00636 private:
00637
00638
00639
00640 void setUp()
00641 {
00642 one( 0, 0 ) = one( 1, 1 ) = one( 2, 2 ) = 1.0;
00643 sone( 0, 0 ) = sone( 1, 1 ) = sone( 2, 2 ) = 1.0;
00644
00645 rdm1( 0, 0 ) = 2;
00646 rdm1( 0, 1 ) = 3;
00647 rdm1( 0, 2 ) = 1;
00648 rdm1( 1, 0 ) = 0;
00649 rdm1( 1, 1 ) = 4;
00650 rdm1( 1, 2 ) = 0;
00651 rdm1( 2, 0 ) = 5;
00652 rdm1( 2, 1 ) = -3;
00653 rdm1( 2, 2 ) = -7;
00654
00655 rdm2( 0, 0 ) = -1;
00656 rdm2( 0, 1 ) = -1;
00657 rdm2( 0, 2 ) = -1;
00658 rdm2( 1, 0 ) = 3;
00659 rdm2( 1, 1 ) = 0;
00660 rdm2( 1, 2 ) = -2;
00661 rdm2( 2, 0 ) = -1;
00662 rdm2( 2, 1 ) = -2;
00663 rdm2( 2, 2 ) = -3;
00664
00665 res1( 0, 0 ) = 6;
00666 res1( 0, 1 ) = -4;
00667 res1( 0, 2 ) = -11;
00668 res1( 1, 0 ) = 12;
00669 res1( 1, 1 ) = 0;
00670 res1( 1, 2 ) = -8;
00671 res1( 2, 0 ) = -7;
00672 res1( 2, 1 ) = 9;
00673 res1( 2, 2 ) = 22;
00674
00675 srdm1( 0, 0 ) = 2;
00676 srdm1( 0, 1 ) = 3;
00677 srdm1( 0, 2 ) = 1;
00678 srdm1( 1, 1 ) = -1;
00679 srdm1( 1, 2 ) = 4;
00680 srdm1( 2, 2 ) = -3;
00681
00682 srdm2( 0, 0 ) = -2;
00683 srdm2( 0, 1 ) = 5;
00684 srdm2( 0, 2 ) = 2;
00685 srdm2( 1, 1 ) = 0;
00686 srdm2( 1, 2 ) = -3;
00687 srdm2( 2, 2 ) = 1;
00688
00689
00690 res2 = srdm1;
00691
00692 res3( 0, 0 ) = 13;
00693 res3( 0, 1 ) = 7;
00694 res3( 0, 2 ) = -4;
00695 res3( 1, 0 ) = -3;
00696 res3( 1, 1 ) = 3;
00697 res3( 1, 2 ) = 13;
00698 res3( 2, 0 ) = 12;
00699 res3( 2, 1 ) = 14;
00700 res3( 2, 2 ) = -13;
00701
00702 res4( 0, 0 ) = 9;
00703 res4( 0, 1 ) = 15;
00704 res4( 0, 2 ) = -5;
00705 res4( 1, 0 ) = 26;
00706 res4( 1, 1 ) = -7;
00707 res4( 1, 2 ) = -25;
00708 res4( 2, 0 ) = -13;
00709 res4( 2, 1 ) = 28;
00710 res4( 2, 2 ) = 22;
00711 }
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721 double calcTens( WTensorSym< 4, 3, double > const& t, WVector3d const& v )
00722 {
00723 double res = 0.0;
00724 for( int a = 0; a < 3; ++a )
00725 {
00726 for( int b = 0; b < 3; ++b )
00727 {
00728 for( int c = 0; c < 3; ++c )
00729 {
00730 for( int d = 0; d < 3; ++d )
00731 {
00732 res += v[ a ] * v[ b ] * v[ c ] * v[ d ] * t( a, b, c, d );
00733 }
00734 }
00735 }
00736 }
00737 return res;
00738 }
00739
00740
00741 WTensor< 2, 3, int > one;
00742
00743 WTensor< 2, 3, int > zero;
00744
00745 WTensor< 2, 3, int > rdm1;
00746
00747 WTensor< 2, 3, int > rdm2;
00748
00749 WTensor< 2, 3, int > res1;
00750
00751 WTensor< 2, 3, int > res2;
00752
00753 WTensor< 2, 3, int > res3;
00754
00755 WTensor< 2, 3, int > res4;
00756
00757 WTensorSym< 2, 3, int > sone;
00758
00759 WTensorSym< 2, 3, int > szero;
00760
00761 WTensorSym< 2, 3, int > srdm1;
00762
00763 WTensorSym< 2, 3, int > srdm2;
00764 };
00765
00766 #endif // WTENSORFUNCTIONS_TEST_H