1 /*! 2 * \file include/TFEL/Math/Stensor/Internals/StensorComputeEigenVectors.hxx 3 * \brief 4 * \author Thomas Helfer 5 * \date 13 Jul 2006 6 * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights 7 * reserved. 8 * This project is publicly released under either the GNU GPL Licence 9 * or the CECILL-A licence. A copy of thoses licences are delivered 10 * with the sources of TFEL. CEA or EDF may also distribute this 11 * project under specific licensing conditions. 12 */ 13 14 #ifndef LIB_TFEL_STENSORCOMPUTEEIGENVECTORS_HXX 15 #define LIB_TFEL_STENSORCOMPUTEEIGENVECTORS_HXX 16 17 #include <cmath> 18 #include <cassert> 19 #include <algorithm> 20 #include <type_traits> 21 22 #include "TFEL/Config/TFELConfig.hxx" 23 #include "TFEL/TypeTraits/IsReal.hxx" 24 #include "TFEL/TypeTraits/IsFundamentalNumericType.hxx" 25 #include "TFEL/Math/General/CubicRoots.hxx" 26 #include "TFEL/Math/tvector.hxx" 27 #include "TFEL/Math/tmatrix.hxx" 28 #include "TFEL/Math/stensor.hxx" 29 30 namespace tfel { 31 32 namespace math { 33 34 namespace internals { 35 36 template <unsigned short N> 37 struct StensorComputeEigenVectors; 38 39 template <> 40 struct StensorComputeEigenVectors<1u> { 41 template <typename T> testtfel::math::internals::StensorComputeEigenVectors42 static bool test(const T* const, 43 const tfel::math::tvector<3u, T>&, 44 const tfel::math::rotation_matrix<T>&) { 45 return true; 46 } 47 48 template <typename VectorType, typename T> computeEigenVectortfel::math::internals::StensorComputeEigenVectors49 static bool computeEigenVector(VectorType& v, 50 const T* const s, 51 const T vp) { 52 TFEL_STATIC_ASSERT(tfel::typetraits::IsFundamentalNumericType< 53 VectorNumType<VectorType>>::cond); 54 TFEL_STATIC_ASSERT( 55 (std::is_same<VectorNumType<VectorType>, T>::value)); 56 TFEL_STATIC_ASSERT( 57 tfel::typetraits::IsFundamentalNumericType<T>::cond); 58 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 59 constexpr const auto zero = T{0}; 60 constexpr const auto one = T{1}; 61 TFEL_CONSTEXPR const auto e = 10 * std::numeric_limits<T>::min(); 62 if (std::abs(s[0] - vp) < e) { 63 v = {one, zero, zero}; 64 return true; 65 } 66 if (std::abs(s[1] - vp) < e) { 67 v = {zero, one, zero}; 68 return true; 69 } 70 if (std::abs(s[2] - vp) < e) { 71 v = {zero, zero, one}; 72 return true; 73 } 74 return false; 75 } 76 77 template <typename T> exetfel::math::internals::StensorComputeEigenVectors78 static bool exe(const T* const s, 79 tfel::math::tvector<3u, T>& vp, 80 tfel::math::rotation_matrix<T>& m, 81 const bool) { 82 TFEL_STATIC_ASSERT( 83 tfel::typetraits::IsFundamentalNumericType<T>::cond); 84 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 85 vp = {s[0], s[1], s[2]}; 86 m = tfel::math::rotation_matrix<T>::Id(); 87 return true; 88 } 89 }; 90 91 template <> 92 struct StensorComputeEigenVectors<2u> { 93 template <typename T> testtfel::math::internals::StensorComputeEigenVectors94 static bool test(const T* const s, 95 const tfel::math::tvector<3u, T>& vp, 96 const tfel::math::rotation_matrix<T>& m) { 97 TFEL_STATIC_ASSERT( 98 tfel::typetraits::IsFundamentalNumericType<T>::cond); 99 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 100 using tfel::math::constexpr_fct::sqrt; 101 constexpr const auto icste = Cste<T>::isqrt2; 102 // first eigenvalue 103 auto y0 = s[0] * m(0, 0) + s[3] * icste * m(1, 0) - vp(0) * m(0, 0); 104 auto y1 = s[1] * m(1, 0) + s[3] * icste * m(0, 0) - vp(0) * m(1, 0); 105 auto tmp = std::sqrt(y0 * y0 + y1 * y1); 106 if (tmp > 100 * std::numeric_limits<T>::epsilon()) { 107 return false; 108 } 109 // second eigenvalue 110 y0 = s[0] * m(0, 1) + s[3] * icste * m(1, 1) - vp(1) * m(0, 1); 111 y1 = s[1] * m(1, 1) + s[3] * icste * m(0, 1) - vp(1) * m(1, 1); 112 tmp = std::sqrt(y0 * y0 + y1 * y1); 113 if (tmp > 100 * std::numeric_limits<T>::epsilon()) { 114 return false; 115 } 116 return true; 117 } 118 119 template <typename VectorType, typename T> computeEigenVectortfel::math::internals::StensorComputeEigenVectors120 static bool computeEigenVector(VectorType& v, 121 const T* const s, 122 const T vp) { 123 using namespace tfel::typetraits; 124 TFEL_STATIC_ASSERT( 125 IsFundamentalNumericType<VectorNumType<VectorType>>::cond); 126 TFEL_STATIC_ASSERT( 127 (std::is_same<VectorNumType<VectorType>, T>::value)); 128 TFEL_STATIC_ASSERT(IsFundamentalNumericType<T>::cond); 129 TFEL_STATIC_ASSERT(IsReal<T>::cond); 130 constexpr const auto zero = T{0}; 131 constexpr const auto one = T{1}; 132 if (std::abs(s[2] - vp) < 10 * std::numeric_limits<T>::min()) { 133 v(0) = zero; 134 v(1) = zero; 135 v(2) = one; 136 return true; 137 } 138 v(2) = zero; 139 return computeEigenVector(s[0], s[1], s[3], vp, v(0), v(1)); 140 } 141 142 template <typename T> exetfel::math::internals::StensorComputeEigenVectors143 static bool exe(const T* const s, 144 tfel::math::tvector<3u, T>& vp, 145 tfel::math::rotation_matrix<T>& vec, 146 const bool) { 147 TFEL_STATIC_ASSERT( 148 tfel::typetraits::IsFundamentalNumericType<T>::cond); 149 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 150 vp(2) = s[2]; 151 vec = tmatrix<3u, 3u, T>::Id(); 152 StensorComputeEigenVectors<2u>::computeEigenValues(s[0], s[1], s[3], 153 vp(0), vp(1)); 154 const auto prec = 155 10 * std::max(vp(0), vp(1)) * std::numeric_limits<T>::epsilon(); 156 if (std::abs(vp(0) - vp(1)) <= prec) { 157 return true; 158 } 159 if (StensorComputeEigenVectors<2u>::computeEigenVector( 160 s[0], s[1], s[3], vp(0), vec(0, 0), vec(1, 0)) == false) { 161 return false; 162 } 163 vec(0, 1) = -vec(1, 0); 164 vec(1, 1) = vec(0, 0); 165 #ifdef TFEL_PARANOIC_CHECK 166 return tfel::math::internals::StensorComputeEigenVectors<2u>::test( 167 s, vp, vec); 168 #else 169 return true; 170 #endif /* LIB_TFEL_STENSORCOMPUTEEIGENVECTORS_HXX */ 171 } 172 173 private: 174 template <typename T> computeEigenValuestfel::math::internals::StensorComputeEigenVectors175 static void computeEigenValues( 176 const T s0, const T s1, const T s3, T& vp1, T& vp2) { 177 TFEL_STATIC_ASSERT( 178 tfel::typetraits::IsFundamentalNumericType<T>::cond); 179 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 180 TFEL_CONSTEXPR const auto one_half = T{1} / T{2}; 181 const T tr = one_half * (s0 + s1); 182 const T tmp = s0 - s1; 183 const T tmp2 = std::sqrt(one_half * (tmp * tmp * one_half + s3 * s3)); 184 vp1 = tr + tmp2; 185 vp2 = tr - tmp2; 186 } 187 188 template <typename T> computeEigenVectortfel::math::internals::StensorComputeEigenVectors189 static bool computeEigenVector( 190 const T s_0, const T s_1, const T s3, const T& vp, T& x, T& y) { 191 TFEL_STATIC_ASSERT( 192 tfel::typetraits::IsFundamentalNumericType<T>::cond); 193 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 194 constexpr const auto zero = T{0}; 195 constexpr const auto one = T{1}; 196 constexpr const auto cste = Cste<T>::sqrt2; 197 auto s0 = s_0 - vp; 198 auto s1 = s_1 - vp; 199 if (std::abs(s3) < 200 std::max(std::min(s_0, s_1) * std::numeric_limits<T>::epsilon(), 201 10 * std::numeric_limits<T>::min())) { 202 if (std::abs(s0) > std::abs(s1)) { 203 x = zero; 204 y = one; 205 return true; 206 } else { 207 x = one; 208 y = zero; 209 return true; 210 } 211 } 212 if (std::abs(s0) > std::abs(s1)) { 213 if (std::abs(s0) < 100 * std::numeric_limits<T>::min()) { 214 return false; 215 } 216 y = one; 217 x = -s3 / (cste * s0); 218 } else { 219 if (std::abs(s1) < 100 * std::numeric_limits<T>::min()) { 220 return false; 221 } 222 x = one; 223 y = -s3 / (cste * s1); 224 } 225 s0 = std::sqrt(x * x + y * y); 226 if (s0 < 100 * std::numeric_limits<T>::min()) { 227 return false; 228 } 229 x /= s0; 230 y /= s0; 231 return true; 232 } 233 }; 234 235 template <> 236 struct StensorComputeEigenVectors<3u> { 237 template <typename T> testtfel::math::internals::StensorComputeEigenVectors238 static bool test(const T* const s, 239 const tfel::math::tvector<3u, T>& vp, 240 const tfel::math::rotation_matrix<T>& m) { 241 TFEL_STATIC_ASSERT( 242 tfel::typetraits::IsFundamentalNumericType<T>::cond); 243 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 244 constexpr const auto icste = Cste<T>::isqrt2; 245 // first eigenvalue 246 auto y0 = s[0] * m(0, 0) + s[3] * icste * m(1, 0) + 247 s[4] * icste * m(2, 0) - vp(0) * m(0, 0); 248 auto y1 = s[1] * m(1, 0) + s[3] * icste * m(0, 0) + 249 s[5] * icste * m(2, 0) - vp(0) * m(1, 0); 250 auto y2 = s[2] * m(2, 0) + s[4] * icste * m(0, 0) + 251 s[5] * icste * m(1, 0) - vp(0) * m(2, 0); 252 auto tmp = norm(y0, y1, y2); 253 if (tmp > 100 * std::numeric_limits<T>::epsilon()) { 254 return false; 255 } 256 // second eigenvalue 257 y0 = s[0] * m(0, 1) + s[3] * icste * m(1, 1) + 258 s[4] * icste * m(2, 1) - vp(1) * m(0, 1); 259 y1 = s[1] * m(1, 1) + s[3] * icste * m(0, 1) + 260 s[5] * icste * m(2, 1) - vp(1) * m(1, 1); 261 y2 = s[2] * m(2, 1) + s[4] * icste * m(0, 1) + 262 s[5] * icste * m(1, 1) - vp(1) * m(2, 1); 263 tmp = norm(y0, y1, y2); 264 if (tmp > 100 * std::numeric_limits<T>::epsilon()) { 265 return false; 266 } 267 // third eigenvalue 268 y0 = s[0] * m(0, 2) + s[3] * icste * m(1, 2) + 269 s[4] * icste * m(2, 2) - vp(2) * m(0, 2); 270 y1 = s[1] * m(1, 2) + s[3] * icste * m(0, 2) + 271 s[5] * icste * m(2, 2) - vp(2) * m(1, 2); 272 y2 = s[2] * m(2, 2) + s[4] * icste * m(0, 2) + 273 s[5] * icste * m(1, 2) - vp(2) * m(2, 2); 274 tmp = norm(y0, y1, y2); 275 if (tmp > 1000 * std::numeric_limits<T>::epsilon()) { 276 return false; 277 } 278 return true; 279 } 280 281 template <typename VectorType, typename T> computeEigenVectortfel::math::internals::StensorComputeEigenVectors282 static bool computeEigenVector(VectorType& v, 283 const T* const src, 284 const T vp) { 285 using namespace tfel::typetraits; 286 TFEL_STATIC_ASSERT( 287 IsFundamentalNumericType<VectorNumType<VectorType>>::cond); 288 TFEL_STATIC_ASSERT( 289 (std::is_same<VectorNumType<VectorType>, T>::value)); 290 TFEL_STATIC_ASSERT(IsFundamentalNumericType<T>::cond); 291 TFEL_STATIC_ASSERT(IsReal<T>::cond); 292 return computeEigenVector(src, vp, v(0), v(1), v(2)); 293 } 294 295 template <typename T> exetfel::math::internals::StensorComputeEigenVectors296 static bool exe(const T* const s, 297 tfel::math::tvector<3u, T>& vp, 298 tfel::math::rotation_matrix<T>& vec, 299 const bool b) { 300 using namespace std; 301 using tfel::math::tvector; 302 TFEL_STATIC_ASSERT( 303 tfel::typetraits::IsFundamentalNumericType<T>::cond); 304 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 305 TFEL_CONSTEXPR const auto rel_prec = 306 100 * std::numeric_limits<T>::epsilon(); 307 StensorComputeEigenValues<3u>::exe(s, vp(0), vp(1), vp(2), b); 308 const auto tr = (s[0] + s[1] + s[2]) / 3; 309 auto ms = T(0); 310 for (unsigned short i = 0; i != 6; ++i) { 311 ms = std::max(ms, std::abs(s[i])); 312 } 313 const bool bsmall = ((ms < 100 * std::numeric_limits<T>::min()) || 314 (ms * std::numeric_limits<T>::epsilon() < 315 100 * std::numeric_limits<T>::min())); 316 if (bsmall) { 317 // all eigenvalues are equal 318 vec = tmatrix<3u, 3u, T>::Id(); 319 #ifdef TFEL_PARANOIC_CHECK 320 return tfel::math::internals::StensorComputeEigenVectors<3>::test( 321 s, vp, vec); 322 #else 323 return true; 324 #endif /* LIB_TFEL_STENSORCOMPUTEEIGENVECTORS_HXX */ 325 } 326 const auto ims = T(1) / ms; 327 tvector<6u, T> s2(s); 328 tvector<3u, T> vp2(vp); 329 vp2(0) = (vp(0) - tr) * ims; 330 vp2(1) = (vp(1) - tr) * ims; 331 vp2(2) = (vp(2) - tr) * ims; 332 s2 *= ims; 333 s2(0) -= tr * ims; 334 s2(1) -= tr * ims; 335 s2(2) -= tr * ims; 336 const auto prec = 337 10 * std::max(rel_prec, 100 * std::numeric_limits<T>::min()); 338 if (((std::abs(vp2(0) - vp2(1)) <= prec) && 339 (std::abs(vp2(0) - vp2(2)) <= prec)) || 340 ((std::abs(vp2(1) - vp2(0)) <= prec) && 341 (std::abs(vp2(1) - vp2(2)) <= prec)) || 342 ((std::abs(vp2(2) - vp2(0)) <= prec) && 343 (std::abs(vp2(2) - vp2(1)) <= prec))) { 344 // all eigenvalues are equal 345 vec = tmatrix<3u, 3u, T>::Id(); 346 #ifdef TFEL_PARANOIC_CHECK 347 return tfel::math::internals::StensorComputeEigenVectors<3>::test( 348 s, vp, vec); 349 #else 350 return true; 351 #endif /* LIB_TFEL_STENSORCOMPUTEEIGENVECTORS_HXX */ 352 } 353 if ((std::abs(vp2(0) - vp2(1)) > prec) && 354 (std::abs(vp2(0) - vp2(2)) > prec)) { 355 // vp0 is single 356 if (computeEigenVector(s2.begin(), vp2(0), vec(0, 0), vec(1, 0), 357 vec(2, 0)) == false) { 358 return false; 359 } 360 if (std::abs(vp2(1) - vp2(2)) > prec) { 361 // vp1 is single 362 if (computeEigenVector(s2.begin(), vp2(1), vec(0, 1), vec(1, 1), 363 vec(2, 1)) == false) { 364 return false; 365 } 366 cross_product(vec(0, 2), vec(1, 2), vec(2, 2), vec(0, 0), 367 vec(1, 0), vec(2, 0), vec(0, 1), vec(1, 1), 368 vec(2, 1)); 369 } else { 370 // vp1 and vp2 are equal 371 find_perpendicular_vector(vec(0, 1), vec(1, 1), vec(2, 1), 372 vec(0, 0), vec(1, 0), vec(2, 0)); 373 cross_product(vec(0, 2), vec(1, 2), vec(2, 2), vec(0, 0), 374 vec(1, 0), vec(2, 0), vec(0, 1), vec(1, 1), 375 vec(2, 1)); 376 } 377 #ifdef TFEL_PARANOIC_CHECK 378 return tfel::math::internals::StensorComputeEigenVectors<3>::test( 379 s, vp, vec); 380 #else 381 return true; 382 #endif /* LIB_TFEL_STENSORCOMPUTEEIGENVECTORS_HXX */ 383 } else if ((std::abs(vp2(1) - vp2(0)) > prec) && 384 (std::abs(vp2(1) - vp2(2)) > prec)) { 385 // vp1 is single, vp0 and vp2 are equal 386 assert(std::abs(vp2(0) - vp2(2)) < prec); 387 if (computeEigenVector(s2.begin(), vp2(1), vec(0, 1), vec(1, 1), 388 vec(2, 1)) == false) { 389 return false; 390 } 391 find_perpendicular_vector(vec(0, 0), vec(1, 0), vec(2, 0), 392 vec(0, 1), vec(1, 1), vec(2, 1)); 393 cross_product(vec(0, 2), vec(1, 2), vec(2, 2), vec(0, 0), vec(1, 0), 394 vec(2, 0), vec(0, 1), vec(1, 1), vec(2, 1)); 395 #ifdef TFEL_PARANOIC_CHECK 396 return tfel::math::internals::StensorComputeEigenVectors<3>::test( 397 s, vp, vec); 398 #else 399 return true; 400 #endif /* LIB_TFEL_STENSORCOMPUTEEIGENVECTORS_HXX */ 401 } 402 assert((std::abs(vp2(2) - vp2(0)) > prec) && 403 (std::abs(vp2(2) - vp2(1)) > prec)); 404 assert(std::abs(vp2(0) - vp2(1)) < prec); 405 if (computeEigenVector(s2.begin(), vp2(2), vec(0, 2), vec(1, 2), 406 vec(2, 2)) == false) { 407 return false; 408 } 409 find_perpendicular_vector(vec(0, 0), vec(1, 0), vec(2, 0), vec(0, 2), 410 vec(1, 2), vec(2, 2)); 411 cross_product(vec(0, 1), vec(1, 1), vec(2, 1), vec(0, 2), vec(1, 2), 412 vec(2, 2), vec(0, 0), vec(1, 0), vec(2, 0)); 413 #ifdef TFEL_PARANOIC_CHECK 414 return tfel::math::internals::StensorComputeEigenVectors<3>::test( 415 s, vp, vec); 416 #else 417 return true; 418 #endif /* LIB_TFEL_STENSORCOMPUTEEIGENVECTORS_HXX */ 419 } 420 421 private: 422 template <typename T> conditionning_numbertfel::math::internals::StensorComputeEigenVectors423 static T conditionning_number(const T a, const T b, const T c) { 424 TFEL_STATIC_ASSERT( 425 tfel::typetraits::IsFundamentalNumericType<T>::cond); 426 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 427 auto tmp = a + b; 428 auto tmp2 = a - b; 429 const auto X = tmp / 2; 430 const auto Y = std::sqrt(tmp2 * tmp2 + 2 * c * c) / 2; 431 tmp = std::abs(X + Y); 432 tmp2 = std::abs(X - Y); 433 if (tmp > tmp2) { 434 if (tmp2 < 100 * std::numeric_limits<T>::min()) { 435 return std::numeric_limits<T>::max(); 436 } 437 return tmp / tmp2; 438 } else { 439 if (tmp < 100 * std::numeric_limits<T>::min()) { 440 return std::numeric_limits<T>::max(); 441 } 442 return tmp2 / tmp; 443 } 444 } 445 446 template <typename T> normtfel::math::internals::StensorComputeEigenVectors447 static T norm(const T& x0, const T& x1, const T& x2) { 448 TFEL_STATIC_ASSERT( 449 tfel::typetraits::IsFundamentalNumericType<T>::cond); 450 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 451 return std::sqrt(x0 * x0 + x1 * x1 + x2 * x2); 452 } 453 454 template <typename T> norm2tfel::math::internals::StensorComputeEigenVectors455 static T norm2(const T x0, const T x1, const T x2) { 456 TFEL_STATIC_ASSERT( 457 tfel::typetraits::IsFundamentalNumericType<T>::cond); 458 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 459 return x0 * x0 + x1 * x1 + x2 * x2; 460 } 461 462 template <typename T> cross_producttfel::math::internals::StensorComputeEigenVectors463 static void cross_product(T& z0, 464 T& z1, 465 T& z2, 466 const T x0, 467 const T x1, 468 const T x2, 469 const T y0, 470 const T y1, 471 const T y2) { 472 TFEL_STATIC_ASSERT( 473 tfel::typetraits::IsFundamentalNumericType<T>::cond); 474 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 475 z0 = x1 * y2 - x2 * y1; 476 z1 = x2 * y0 - x0 * y2; 477 z2 = x0 * y1 - x1 * y0; 478 } 479 480 template <typename T> find_perpendicular_vectortfel::math::internals::StensorComputeEigenVectors481 static void find_perpendicular_vector( 482 T& y0, T& y1, T& y2, const T x0, const T x1, const T x2) { 483 TFEL_STATIC_ASSERT( 484 tfel::typetraits::IsFundamentalNumericType<T>::cond); 485 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 486 constexpr const auto zero = T{0}; 487 constexpr const auto one = T{1}; 488 const auto norm2_x = norm2(x0, x1, x2); 489 if (norm2_x < 100 * std::numeric_limits<T>::min()) { 490 // x is null 491 y0 = one; 492 y1 = zero; 493 y2 = zero; 494 return; 495 } 496 if (std::abs(x0) < std::abs(x1)) { 497 if (std::abs(x0) < std::abs(x2)) { 498 //|x0| is min, (1 0 0) is a good choice 499 y0 = one - x0 * x0 / norm2_x; 500 y1 = -x0 * x1 / norm2_x; 501 y2 = -x0 * x2 / norm2_x; 502 } else { 503 //|x2| is min, (0 0 1) is a good choice 504 y0 = -x2 * x0 / norm2_x; 505 y1 = -x2 * x1 / norm2_x; 506 y2 = one - x2 * x2 / norm2_x; 507 } 508 } else if (std::abs(x1) < std::abs(x2)) { 509 // |x1| is min, (0 0 1) is a good choice 510 y0 = -x1 * x0 / norm2_x; 511 y1 = one - x1 * x1 / norm2_x; 512 y2 = -x1 * x2 / norm2_x; 513 } else { 514 // |x2| is min, (0 0 1) is a good choice 515 y0 = -x2 * x0 / norm2_x; 516 y1 = -x2 * x1 / norm2_x; 517 y2 = one - x2 * x2 / norm2_x; 518 } 519 const auto tmp = norm(y0, y1, y2); 520 y0 /= tmp; 521 y1 /= tmp; 522 y2 /= tmp; 523 } 524 525 template <typename T> computeEigenVectortfel::math::internals::StensorComputeEigenVectors526 static bool computeEigenVector( 527 const T* const src, const T vp, T& v0, T& v1, T& v2) { 528 using namespace tfel::math; 529 TFEL_STATIC_ASSERT( 530 tfel::typetraits::IsFundamentalNumericType<T>::cond); 531 TFEL_STATIC_ASSERT(tfel::typetraits::IsReal<T>::cond); 532 // Assume that vp is a single eigenvalue 533 constexpr const auto icste = Cste<T>::isqrt2; 534 const auto a = src[0] - vp; 535 const auto b = src[3] * icste; 536 const auto c = src[4] * icste; 537 const auto d = src[1] - vp; 538 const auto e = src[5] * icste; 539 const auto f = src[2] - vp; 540 const auto det3 = a * d - b * b; 541 const auto det2 = a * f - c * c; 542 const auto det1 = d * f - e * e; 543 const auto prec = 544 std::max(std::abs(vp) * std::numeric_limits<T>::epsilon(), 545 10 * std::numeric_limits<T>::min()); 546 if ((std::abs(det1) < prec) && (std::abs(det2) < prec) && 547 (std::abs(det3) < prec)) { 548 tvector<3u, T> v; 549 tmatrix<3u, 3u, T> m; 550 if (!StensorComputeEigenVectors<3u>::exe(src, v, m, true)) { 551 v0 = v1 = v2 = T(0); 552 return false; 553 } 554 const auto d0 = std::abs(v[0] - vp); 555 const auto d1 = std::abs(v[1] - vp); 556 const auto d2 = std::abs(v[2] - vp); 557 if ((d0 < d1) && (d0 < d2)) { 558 v0 = m(0, 0); 559 v1 = m(1, 0); 560 v2 = m(2, 0); 561 return true; 562 } else if ((d1 < d2) && (d1 < d0)) { 563 v0 = m(0, 1); 564 v1 = m(1, 1); 565 v2 = m(2, 1); 566 return true; 567 } 568 v0 = m(0, 2); 569 v1 = m(1, 2); 570 v2 = m(2, 2); 571 return true; 572 } 573 if ((std::abs(det3) >= std::abs(det1)) && 574 (std::abs(det3) >= std::abs(det2))) { 575 v0 = (b * e - c * d) / det3; 576 v1 = (b * c - a * e) / det3; 577 v2 = T(1); 578 } else if ((std::abs(det1) >= std::abs(det2)) && 579 (std::abs(det1) >= std::abs(det3))) { 580 v0 = T(1); 581 v1 = (c * e - b * f) / det1; 582 v2 = (b * e - c * d) / det1; 583 } else { 584 v0 = (c * e - b * f) / det2; 585 v1 = T(1); 586 v2 = (b * c - a * e) / det2; 587 } 588 const auto nr = norm(v0, v1, v2); 589 v0 /= nr; 590 v1 /= nr; 591 v2 /= nr; 592 return true; 593 } 594 }; 595 596 } // end of namespace internals 597 598 } // end of namespace math 599 600 } // end of namespace tfel 601 602 #endif /* LIB_TFEL_STENSORCOMPUTEEIGENVECTORS_HXX */ 603