00001
00002
00003
00004
00005
00006 #include<ctime>
00007
00008 #include <Intrepid_MiniTensor.h>
00009
00010 #include "Intrepid_FieldContainer.hpp"
00011 #include "Sacado.hpp"
00012 #include "Intrepid_MiniTensor.h"
00013 #include "Teuchos_UnitTestHarness.hpp"
00014
00015 using namespace std;
00016
00017 typedef double ScalarT;
00018
00019 namespace Intrepid
00020 {
00021
00022 TEUCHOS_UNIT_TEST(MiniTensor, Initialization)
00023 {
00024 FieldContainer<Real> FC(3, 3);
00025 FC(0, 0) = 1.0;
00026 FC(0, 1) = 2.0;
00027 FC(0, 2) = 3.0;
00028 FC(1, 0) = 4.0;
00029 FC(1, 1) = 5.0;
00030 FC(1, 2) = 6.0;
00031 FC(2, 0) = 7.0;
00032 FC(2, 1) = 8.0;
00033 FC(2, 2) = 9.0;
00034
00035 Real const * dataPtr0 = &FC(0, 0);
00036
00037 Index const N = 3;
00038 Vector<Real> u(N, dataPtr0);
00039
00040 TEST_COMPARE( u(0), ==, 1.0);
00041 TEST_COMPARE( u(1), ==, 2.0);
00042 TEST_COMPARE( u(2), ==, 3.0);
00043
00044 Real const * dataPtr1 = &FC(1, 0);
00045
00046 u = Vector<Real>(N, dataPtr1);
00047
00048 TEST_COMPARE( u(0), ==, 4.0);
00049 TEST_COMPARE( u(1), ==, 5.0);
00050 TEST_COMPARE( u(2), ==, 6.0);
00051
00052 Real const * dataPtr2 = &FC(2, 0);
00053
00054 u = Vector<Real>(N, dataPtr2);
00055
00056 TEST_COMPARE( u(0), ==, 7.0);
00057 TEST_COMPARE( u(1), ==, 8.0);
00058 TEST_COMPARE( u(2), ==, 9.0);
00059 }
00060
00061 TEUCHOS_UNIT_TEST(MiniTensor, VectorAddition)
00062 {
00063 Vector<Real> const u(1.0, 0.0, 0.0);
00064 Vector<Real> const v(0.0, 1.0, 0.0);
00065 Vector<Real> const w(1.0, 1.0, 0.0);
00066
00067 TEST_COMPARE( u + v == w, !=, 0);
00068 }
00069
00070 TEUCHOS_UNIT_TEST(MiniTensor, VectorSubtraction)
00071 {
00072 Vector<Real> u(3);
00073 Vector<Real> v(3);
00074 u(0) = 1.0;
00075 u(1) = 2.0;
00076 u(2) = 3.0;
00077
00078 v = u - u;
00079
00080 TEST_COMPARE(norm(v), <=, machine_epsilon<Real>());
00081 }
00082
00083 TEUCHOS_UNIT_TEST(MiniTensor, VectorScalarMultipliaction)
00084 {
00085 Vector<Real> u(3);
00086 Vector<Real> v(3);
00087 Vector<Real> w(3);
00088 u(0) = 1.0;
00089 u(1) = 2.0;
00090 u(2) = 3.0;
00091
00092 v(0) = -2.0;
00093 v(1) = -4.0;
00094 v(2) = -6.0;
00095
00096 w = 4.0 * u + 2.0 * v;
00097
00098 TEST_COMPARE( norm(w), <=, machine_epsilon<Real>());
00099 }
00100
00101 TEUCHOS_UNIT_TEST(MiniTensor, TensorInstantiation)
00102 {
00103 FieldContainer<Real> FC(2, 3, 3);
00104 FC(0, 0, 0) = 1.0;
00105 FC(0, 0, 1) = 2.0;
00106 FC(0, 0, 2) = 3.0;
00107 FC(0, 1, 0) = 4.0;
00108 FC(0, 1, 1) = 5.0;
00109 FC(0, 1, 2) = 6.0;
00110 FC(0, 2, 0) = 7.0;
00111 FC(0, 2, 1) = 8.0;
00112 FC(0, 2, 2) = 9.0;
00113 FC(1, 0, 0) = 10.0;
00114 FC(1, 0, 1) = 11.0;
00115 FC(1, 0, 2) = 12.0;
00116 FC(1, 1, 0) = 13.0;
00117 FC(1, 1, 1) = 14.0;
00118 FC(1, 1, 2) = 15.0;
00119 FC(1, 2, 0) = 16.0;
00120 FC(1, 2, 1) = 17.0;
00121 FC(1, 2, 2) = 18.0;
00122
00123 Real const * dataPtr0 = &FC(0, 0, 0);
00124
00125 Tensor<Real> const A(3, dataPtr0);
00126
00127 TEST_COMPARE( A(0,0), ==, 1.0);
00128 TEST_COMPARE( A(0,1), ==, 2.0);
00129 TEST_COMPARE( A(0,2), ==, 3.0);
00130 TEST_COMPARE( A(1,0), ==, 4.0);
00131 TEST_COMPARE( A(1,1), ==, 5.0);
00132 TEST_COMPARE( A(1,2), ==, 6.0);
00133 TEST_COMPARE( A(2,0), ==, 7.0);
00134 TEST_COMPARE( A(2,1), ==, 8.0);
00135 TEST_COMPARE( A(2,2), ==, 9.0);
00136
00137 Real const * dataPtr1 = &FC(1, 0, 0);
00138
00139 Tensor<Real> const B(3, dataPtr1);
00140
00141 TEST_COMPARE( B(0,0), ==, 10.0);
00142 TEST_COMPARE( B(0,1), ==, 11.0);
00143 TEST_COMPARE( B(0,2), ==, 12.0);
00144 TEST_COMPARE( B(1,0), ==, 13.0);
00145 TEST_COMPARE( B(1,1), ==, 14.0);
00146 TEST_COMPARE( B(1,2), ==, 15.0);
00147 TEST_COMPARE( B(2,0), ==, 16.0);
00148 TEST_COMPARE( B(2,1), ==, 17.0);
00149 TEST_COMPARE( B(2,2), ==, 18.0);
00150 }
00151
00152 TEUCHOS_UNIT_TEST(MiniTensor, TensorAddition)
00153 {
00154 Tensor<Real> const A(3, 1.0);
00155 Tensor<Real> const B(3, 2.0);
00156 Tensor<Real> const C(3, 3.0);
00157
00158 TEST_COMPARE( C == A + B, !=, 0);
00159 }
00160
00161 TEUCHOS_UNIT_TEST(MiniTensor, Inverse)
00162 {
00163 std::srand(std::time(NULL));
00164 Index const N = double(std::rand()) / double(RAND_MAX) * 7.0 + 3.0;
00165 Tensor<Real> A(N);
00166 Tensor<Real> B(N);
00167 Tensor<Real> C(N);
00168
00169 for (Index i = 0; i < N; ++i) {
00170 for (Index j = 0; j < N; ++j) {
00171 A(i, j) = double(std::rand()) / double(RAND_MAX) * 20.0 - 10.0;
00172 }
00173 }
00174
00175 B = inverse(A);
00176
00177 C = A * B;
00178
00179 Real const error = norm(C - eye<Real>(N)) / norm(A);
00180
00181 TEST_COMPARE(error, <=, 100.0 * machine_epsilon<Real>());
00182 }
00183
00184 TEUCHOS_UNIT_TEST(MiniTensor, TensorManipulation)
00185 {
00186 Tensor<Real> A = eye<Real>(3);
00187 Tensor<Real> B(3);
00188 Tensor<Real> C(3);
00189 Vector<Real> u(3);
00190
00191 A = 2.0 * A;
00192 A(1, 0) = A(0, 1) = 1.0;
00193 A(2, 1) = A(1, 2) = 1.0;
00194
00195 B = inverse(A);
00196
00197 C = A * B;
00198
00199 TEST_COMPARE(norm(C - eye<Real>(3)), <=, machine_epsilon<Real>());
00200
00201 Real I1_A = I1(A);
00202 Real I2_A = I2(A);
00203 Real I3_A = I3(A);
00204
00205 u(0) = I1_A - 6;
00206 u(1) = I2_A - 10;
00207 u(2) = I3_A - 4;
00208
00209 Real const error = norm(u);
00210
00211 TEST_COMPARE(error, <=, machine_epsilon<Real>());
00212 }
00213
00214 TEUCHOS_UNIT_TEST(MiniTensor, Exponential)
00215 {
00216 Tensor<Real> const A(1, 2, 3, 4, 5, 6, 7, 8, 9);
00217
00218 Tensor<Real> const B = exp_pade(A);
00219
00220 Tensor<Real> const C = exp_taylor(A);
00221
00222 Tensor<Real> const D = B - C;
00223
00224 Real const error = norm(D) / norm(B);
00225
00226 TEST_COMPARE( error, <=, 100.0 * machine_epsilon<Real>());
00227 }
00228
00229 TEUCHOS_UNIT_TEST(MiniTensor, SymmetricEigen)
00230 {
00231 Tensor<Real> A = eye<Real>(3);
00232 A(0, 1) = 0.1;
00233 A(1, 0) = 0.1;
00234
00235 Tensor<Real> V(3);
00236 Tensor<Real> D(3);
00237
00238 boost::tie(V, D) = eig_sym(A);
00239
00240 TEST_COMPARE(std::abs(D(0,0) - 1.1), <=, machine_epsilon<Real>());
00241 TEST_COMPARE(std::abs(D(1,1) - 1.0), <=, machine_epsilon<Real>());
00242 TEST_COMPARE(std::abs(D(2,2) - 0.9), <=, machine_epsilon<Real>());
00243 }
00244
00245 TEUCHOS_UNIT_TEST(MiniTensor, LeftPolarDecomposition)
00246 {
00247 Tensor<Real> V0(1.1, 0.2, 0.0, 0.2, 1.0, 0.0, 0.0, 0.0, 1.2);
00248
00249 Tensor<Real> R0(sqrt(2) / 2, -sqrt(2) / 2, 0.0, sqrt(2) / 2, sqrt(2) / 2,
00250 0.0, 0.0, 0.0, 1.0);
00251
00252 Tensor<Real> F = V0 * R0;
00253 Tensor<Real> V(3);
00254 Tensor<Real> R(3);
00255 boost::tie(V, R) = polar_left(F);
00256
00257 TEST_COMPARE(norm(V-V0), <=, 10.0*machine_epsilon<Real>());
00258 TEST_COMPARE(norm(R-R0), <=, machine_epsilon<Real>());
00259 }
00260
00261 TEUCHOS_UNIT_TEST(MiniTensor, LogRotation)
00262 {
00263 Tensor<Real> R = identity<Real>(3);
00264 Tensor<Real> R0(sqrt(2) / 2, -sqrt(2) / 2, 0.0, sqrt(2) / 2, sqrt(2) / 2,
00265 0.0, 0.0, 0.0, 1.0);
00266
00267 Tensor<Real> r = log_rotation(R);
00268 Tensor<Real> r0 = log_rotation(R0);
00269
00270 TEST_COMPARE(norm(r), <=, machine_epsilon<Real>());
00271
00272 TEST_COMPARE( std::abs(r0(0,1) + 0.785398163397448), <=,
00273 10.0*machine_epsilon<Real>());
00274
00275 TEST_COMPARE( std::abs(r0(0,1) + r0(1,0)), <=, machine_epsilon<Real>());
00276
00277 Real theta = std::acos(-1.0) + 10 * machine_epsilon<Real>();
00278
00279 R(0, 0) = cos(theta);
00280 R(1, 1) = cos(theta);
00281 R(0, 1) = sin(theta);
00282 R(1, 0) = -sin(theta);
00283 R(2, 2) = 1.0;
00284
00285 Tensor<Real> logR = log_rotation(R);
00286
00287 Tensor<Real> Rref(3, 0.0);
00288 Rref(0, 1) = -theta;
00289 Rref(1, 0) = theta;
00290
00291 TEST_COMPARE(norm(logR - Rref), <=, 100*machine_epsilon<Real>());
00292 }
00293
00294 TEUCHOS_UNIT_TEST(MiniTensor, BakerCampbellHausdorff)
00295 {
00296 Tensor<Real> F = 3.0 * identity<Real>(3);
00297 Tensor<Real> V(3), R(3), logV(3), logR(3);
00298
00299 boost::tie(V, R, logV) = polar_left_logV(F);
00300 logR = log_rotation(R);
00301
00302 Tensor<Real> f = bch(logV, logR);
00303
00304 TEST_COMPARE( std::abs(f(0,0) - std::log(3.0)), <=,
00305 machine_epsilon<Real>());
00306
00307 Vector<Real> u(3);
00308 u(0) = std::acos(-1.0) / std::sqrt(2.0);
00309 u(1) = u(0);
00310 u(2) = 0.0;
00311
00312 Tensor<Real> R1(3, 0.0);
00313 Tensor<Real> logR2(3, 0.0);
00314 logR2(0, 2) = u(1);
00315 logR2(1, 2) = -u(0);
00316 logR2(2, 0) = -u(1);
00317 logR2(2, 1) = u(0);
00318 logR2(0, 1) = -u(2);
00319 logR2(1, 0) = u(2);
00320
00321 R1 = exp_skew_symmetric(logR2);
00322 Tensor<Real> Rref = zero<Real>(3);
00323 Rref(0, 1) = 1.0;
00324 Rref(1, 0) = 1.0;
00325 Rref(2, 2) = -1.0;
00326
00327 TEST_COMPARE( norm(Rref-R1), <=, 100.0*machine_epsilon<Real>());
00328 TEST_COMPARE( norm(exp_skew_symmetric(logR) - R), <=,
00329 100.0*machine_epsilon<Real>());
00330 }
00331
00332 TEUCHOS_UNIT_TEST(MiniTensor, PolarLeftLog)
00333 {
00334 Tensor<Real> const F(3.60070151614402, 0.00545892068653966,
00335 0.144580850331452, -5.73345529510674, 0.176660910549112,
00336 1.39627497290058, 2.51510445213514, 0.453212159218359,
00337 -1.44616077859513);
00338
00339 Tensor<Real> const L(0.265620603957487, -1.066921781600734,
00340 -0.089540974250415, -1.066921781600734, 0.927394431410918,
00341 -0.942214085118614, -0.089540974250415, -0.942214085118613,
00342 0.105672693695746);
00343
00344 Tensor<Real> V(3), R(3), v(3), r(3);
00345
00346 boost::tie(V, R, v) = polar_left_logV(F);
00347
00348 Real const error = norm(v - L) / norm(L);
00349
00350 TEST_COMPARE( error, <=, 100*machine_epsilon<Real>());
00351 }
00352
00353 TEUCHOS_UNIT_TEST(MiniTensor, VolumetricDeviatoric)
00354 {
00355 Tensor<Real> A = 3.0 * eye<Real>(3);
00356
00357 TEST_COMPARE( norm(A - vol(A)), <=, 100.0*machine_epsilon<Real>());
00358
00359 Tensor<Real> B = dev(A);
00360
00361 A(0, 0) = 0.0;
00362 A(1, 1) = 0.0;
00363 A(2, 2) = 0.0;
00364
00365 TEST_COMPARE( norm(A - B), <=, 100.0*machine_epsilon<Real>());
00366 }
00367
00368 TEUCHOS_UNIT_TEST(MiniTensor, SVD2x2)
00369 {
00370 Real const phi = 1.0;
00371
00372 Real const psi = 2.0;
00373
00374 Real const s0 = sqrt(3.0);
00375
00376 Real const s1 = sqrt(2.0);
00377
00378 Real const cl = cos(phi);
00379
00380 Real const sl = sin(phi);
00381
00382 Real const cr = cos(psi);
00383
00384 Real const sr = sin(psi);
00385
00386 Tensor<Real> const X(cl, -sl, sl, cl);
00387
00388 Tensor<Real> const Y(cr, -sr, sr, cr);
00389
00390 Tensor<Real> const D(s0, 0.0, 0.0, s1);
00391
00392 Tensor<Real> const A = X * D * transpose(Y);
00393
00394 Tensor<Real> U(2), S(2), V(2);
00395
00396 boost::tie(U, S, V) = svd(A);
00397
00398 Tensor<Real> B = U * S * transpose(V);
00399
00400 Real const error = norm(A - B) / norm(A);
00401
00402 TEST_COMPARE(error, <=, 100.0*machine_epsilon<Real>());
00403 }
00404
00405 TEUCHOS_UNIT_TEST(MiniTensor, SVD3x3)
00406 {
00407 Tensor<Real> const A(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0);
00408
00409 Tensor<Real> U(3), S(3), V(3);
00410
00411 boost::tie(U, S, V) = svd(A);
00412
00413 Tensor<Real> const B = U * S * transpose(V);
00414
00415 Real const error = norm(A - B) / norm(A);
00416
00417 TEST_COMPARE(error, <=, 100.0*machine_epsilon<Real>());
00418 }
00419
00420 TEUCHOS_UNIT_TEST(MiniTensor, SymmetricEigen2x2)
00421 {
00422 Tensor<Real> const A(2.0, 1.0, 1.0, 2.0);
00423
00424 Tensor<Real> V(2), D(2);
00425
00426 boost::tie(V, D) = eig_sym(A);
00427
00428 Tensor<Real> const B = V * D * transpose(V);
00429
00430 Real const error = norm(A - B) / norm(A);
00431
00432 TEST_COMPARE(error, <=, 100.0*machine_epsilon<Real>());
00433 }
00434
00435 TEUCHOS_UNIT_TEST(MiniTensor, SymmetricEigen3x3)
00436 {
00437 Tensor<Real> const A(2.0, 1.0, 0.0, 1.0, 2.0, 1.0, 0.0, 1.0, 2.0);
00438
00439 Tensor<Real> V(3), D(3);
00440
00441 boost::tie(V, D) = eig_sym(A);
00442
00443 Tensor<Real> const B = V * D * transpose(V);
00444
00445 Real const error = norm(A - B) / norm(A);
00446
00447 TEST_COMPARE(error, <=, 100.0*machine_epsilon<Real>());
00448 }
00449
00450 TEUCHOS_UNIT_TEST(MiniTensor, Inverse4x4)
00451 {
00452 Tensor<Real> A = 2.0 * identity<Real>(4);
00453
00454 A(0, 1) = 1.0;
00455 A(1, 0) = 1.0;
00456
00457 A(1, 2) = 1.0;
00458 A(2, 1) = 1.0;
00459
00460 A(2, 3) = 1.0;
00461 A(3, 2) = 1.0;
00462
00463 Tensor<Real> const B = inverse(A);
00464
00465 Tensor<Real> const C = A * B;
00466
00467 Tensor<Real> const I = eye<Real>(4);
00468
00469 Real const error = norm(C - I) / norm(A);
00470
00471 TEST_COMPARE(error, <=, 100.0*machine_epsilon<Real>());
00472 }
00473
00474 TEUCHOS_UNIT_TEST(MiniTensor, Polar3x3)
00475 {
00476 Tensor<Real> A(2.0, 1.0, 0.0, 0.0, 2.0, 1.0, 0.0, 0.0, 2.0);
00477
00478 Tensor<Real> R(3), U(3);
00479
00480 boost::tie(R, U) = polar_right(A);
00481
00482 Tensor<Real> X(3), D(3), Y(3);
00483
00484 boost::tie(X, D, Y) = svd(A);
00485
00486 Tensor<Real> B = R - X * transpose(Y) + U - Y * D * transpose(Y);
00487
00488 Real const error = norm(B) / norm(A);
00489
00490 TEST_COMPARE(error, <=, 100.0*machine_epsilon<Real>());
00491 }
00492
00493 TEUCHOS_UNIT_TEST(MiniTensor, SVD3x3Fad)
00494 {
00495 Tensor<Sacado::Fad::DFad<double> > A(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0,
00496 9.0);
00497
00498 Tensor<Sacado::Fad::DFad<double> > U(3), S(3), V(3);
00499
00500 boost::tie(U, S, V) = svd(A);
00501
00502 Tensor<Sacado::Fad::DFad<double> > B = U * S * transpose(V);
00503
00504 Sacado::Fad::DFad<double> const error = norm(B - A) / norm(A);
00505
00506 TEST_COMPARE(error, <=, 100.0*machine_epsilon<Real>());
00507 }
00508
00509 TEUCHOS_UNIT_TEST(MiniTensor, Cholesky)
00510 {
00511 Tensor<Real> A(1.0, 1.0, 1.0, 1.0, 5.0, 3.0, 1.0, 3.0, 3.0);
00512
00513 Tensor<Real> G(3);
00514
00515 bool is_spd;
00516
00517 boost::tie(G, is_spd) = cholesky(A);
00518
00519 Tensor<Real> B(1.0, 0.0, 0.0, 1.0, 2.0, 0.0, 1.0, 1.0, 1.0);
00520
00521 Real const error = norm(G - B) / norm(A);
00522
00523 TEST_COMPARE(error, <=, 100.0*machine_epsilon<Real>());
00524 }
00525
00526 TEUCHOS_UNIT_TEST(MiniTensor, MechanicsTransforms)
00527 {
00528 Tensor<Real> F(0.0, -6.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 1.0 / 3.0);
00529
00530 Tensor<Real> sigma(0.0, 0.0, 0.0, 0.0, 50.0, 0.0, 0.0, 0.0, 0.0);
00531
00532 Tensor<Real> P = piola(F, sigma);
00533
00534 Real error = abs(P(1, 0) - 100.0) / 100.0;
00535
00536 TEST_COMPARE(error, <=, machine_epsilon<Real>());
00537
00538 sigma = piola_inverse(F, P);
00539
00540 error = abs(sigma(1, 1) - 50.0) / 50.0;
00541
00542 TEST_COMPARE(error, <=, machine_epsilon<Real>());
00543
00544 Tensor<Real> E = 0.5 * (t_dot(F, F) - eye<Real>(3));
00545
00546 Tensor<Real> e = 0.5 * (eye<Real>(3) - inverse(dot_t(F, F)));
00547
00548 Tensor<Real> g = push_forward_covariant(F, E);
00549
00550 error = norm(g - e) / norm(e);
00551
00552 TEST_COMPARE(error, <=, machine_epsilon<Real>());
00553
00554 Tensor<Real> G = pull_back_covariant(F, e);
00555
00556 error = norm(G - E) / norm(E);
00557
00558 TEST_COMPARE(error, <=, machine_epsilon<Real>());
00559 }
00560
00561 }
