1 2/*! 3 * \file include/TFEL/Math/Stensor/DecompositionInPositiveAndNegativeParts.ixx 4 * \brief 5 * \author Thomas Helfer 6 * \date 22 nov. 2014 7 * \copyright Copyright (C) 2006-2018 CEA/DEN, EDF R&D. All rights 8 * reserved. 9 * This project is publicly released under either the GNU GPL Licence 10 * or the CECILL-A licence. A copy of thoses licences are delivered 11 * with the sources of TFEL. CEA or EDF may also distribute this 12 * project under specific licensing conditions. 13 */ 14 15#ifndef LIB_TFEL_MATH_STENSOR_DECOMPOSITIONINPOSITIVEANDNEGATIVEPARTSIXX 16#define LIB_TFEL_MATH_STENSOR_DECOMPOSITIONINPOSITIVEANDNEGATIVEPARTSIXX 17 18#include"TFEL/Math/General/MathConstants.hxx" 19 20namespace tfel{ 21 22 namespace math{ 23 24 namespace internals{ 25 26 template<typename NumType> 27 TFEL_MATH_INLINE NumType stensor_ppos(const NumType x) 28 { 29 return x >= NumType(0) ? x : NumType(0); 30 } 31 32 template<typename NumType> 33 TFEL_MATH_INLINE NumType stensor_pneg(const NumType x) 34 { 35 return x <= NumType(0) ? x : NumType(0); 36 } 37 38 } 39 40 template<typename DPPType,typename PPType,typename StensorType> 41 typename std::enable_if< 42 tfel::meta::Implements<DPPType,ST2toST2Concept>::cond&& 43 tfel::meta::Implements<PPType,StensorConcept>::cond&& 44 tfel::meta::Implements<StensorType,StensorConcept>::cond&& 45 ST2toST2Traits<DPPType>::dime==StensorTraits<StensorType>::dime && 46 StensorTraits<PPType>::dime==StensorTraits<StensorType>::dime && 47 StensorTraits<StensorType>::dime==1u && 48 tfel::typetraits::IsAssignableTo<StensorNumType<StensorType>, 49 StensorNumType<PPType>>::cond && 50 tfel::typetraits::IsAssignableTo<tfel::typetraits::base_type<StensorNumType<StensorType>>, 51 ST2toST2NumType<DPPType>>::cond, 52 void>::type 53 computeStensorPositivePartAndDerivative(DPPType& dpp, 54 PPType& pp, 55 const StensorType& s, 56 const StensorNumType<StensorType> eps) 57 { 58 using std::abs; 59 typedef StensorNumType<StensorType> NumType; 60 typedef tfel::typetraits::base_type<NumType> real; 61 TFEL_CONSTEXPR const auto one_half = real(1)/(real(2)); 62 dpp(0,1) = dpp(0,2) = dpp(1,0) = dpp(1,2) = dpp(2,0) = dpp(2,1) = real(0); 63 if(abs(s(0))<eps){ 64 dpp(0,0) = one_half; 65 pp(0)=real(0); 66 } else if (s(0)>NumType(0)){ 67 dpp(0,0) = real(1); 68 pp(0)=s(0); 69 } else { 70 dpp(0,0) = real(0); 71 pp(0)=NumType(0); 72 } 73 if(abs(s(1))<eps){ 74 dpp(1,1) = one_half; 75 pp(1)=real(0); 76 } else if (s(1)>NumType(0)){ 77 dpp(1,1) = real(1); 78 pp(1)=s(1); 79 } else { 80 dpp(1,1) = real(0); 81 pp(1)=NumType(0); 82 } 83 if(abs(s(2))<eps){ 84 dpp(2,2) = one_half; 85 pp(2)=real(0); 86 } else if (s(2)>NumType(0)){ 87 dpp(2,2) = real(1); 88 pp(2)=s(2); 89 } else { 90 dpp(2,2) = real(0); 91 pp(2)=NumType(0); 92 } 93 } 94 95 template<typename DPPType, 96 typename PPType, 97 typename StensorType> 98 typename std::enable_if< 99 tfel::meta::Implements<DPPType,ST2toST2Concept>::cond&& 100 tfel::meta::Implements<PPType,StensorConcept>::cond&& 101 tfel::meta::Implements<StensorType,StensorConcept>::cond&& 102 ST2toST2Traits<DPPType>::dime==StensorTraits<StensorType>::dime && 103 StensorTraits<PPType>::dime==StensorTraits<StensorType>::dime && 104 StensorTraits<StensorType>::dime==2u && 105 tfel::typetraits::IsAssignableTo<StensorNumType<StensorType>, 106 StensorNumType<PPType>>::cond && 107 tfel::typetraits::IsAssignableTo<tfel::typetraits::base_type<StensorNumType<StensorType>>, 108 ST2toST2NumType<DPPType>>::cond, 109 void>::type 110 computeStensorPositivePartAndDerivative(DPPType& dpp, 111 PPType& pp, 112 const StensorType& s, 113 const StensorNumType<StensorType> eps) 114 { 115 using std::abs; 116 using tfel::math::internals::stensor_ppos; 117 using tfel::math::internals::stensor_pneg; 118 typedef StensorNumType<StensorType> NumType; 119 typedef tfel::typetraits::base_type<NumType> real; 120 constexpr const auto cste = Cste<real>::sqrt2; 121 TFEL_CONSTEXPR const auto one_half = real(1)/(real(2)); 122 stensor<2u,NumType> ls(s); // local copy of s 123 stensor<2u,real> n0; 124 stensor<2u,real> n1; 125 stensor<2u,real> n2; 126 tvector<3u,NumType> vp; 127 tmatrix<3u,3u,real> m; 128 ls.computeEigenVectors(vp,m); 129 stensor<2u,real>::computeEigenTensors(n0,n1,n2,m); 130 const tvector<3u,real> v0 = m.template column_view<0u>(); 131 const tvector<3u,real> v1 = m.template column_view<1u>(); 132 const stensor<2u,real> n01 = stensor<2u,real>::buildFromVectorsSymmetricDiadicProduct(v0,v1)/cste; 133 if(abs(vp(0)-vp(1))<eps){ 134 const NumType vpm = (vp(0)+vp(1))*one_half; 135 if(abs(vpm)<eps){ 136 dpp = ((n0^n0)+(n1^n1)+(n01^n01))*one_half; 137 pp = stensor<2u,NumType>(real(0)); 138 } else if(vpm>NumType(0)){ 139 dpp = (n0^n0)+(n1^n1)+(n01^n01); 140 pp = vpm*(n0+n1); 141 } else { 142 dpp = st2tost2<2u,real>(real(0)); 143 pp = stensor<2u,NumType>(real(0)); 144 } 145 } else { 146 if(abs(vp(0))<eps){ 147 dpp = (n0^n0)*one_half; 148 pp = stensor<2u,NumType>(real(0)); 149 } else if(vp(0)>NumType(0)){ 150 dpp = (n0^n0); 151 pp = vp(0)*n0; 152 } else { 153 dpp = st2tost2<2u,real>(real(0)); 154 pp = stensor<2u,NumType>(real(0)); 155 } 156 if(abs(vp(1))<eps){ 157 dpp += (n1^n1)*one_half; 158 } else if(vp(1)>NumType(0)){ 159 dpp += (n1^n1); 160 pp += vp(1)*n1; 161 } 162 dpp += (stensor_ppos(vp(1))-stensor_ppos(vp(0)))/(vp(1)-vp(0))*(n01^n01); 163 } 164 if(abs(vp(2))<eps){ 165 dpp += (n2^n2)*one_half; 166 } else if (vp(2)>0){ 167 dpp += (n2^n2); 168 pp(2) += vp(2); 169 } 170 } 171 172 template<typename DPPType, 173 typename PPType, 174 typename StensorType> 175 typename std::enable_if< 176 tfel::meta::Implements<DPPType,ST2toST2Concept>::cond&& 177 tfel::meta::Implements<PPType,StensorConcept>::cond&& 178 tfel::meta::Implements<StensorType,StensorConcept>::cond&& 179 ST2toST2Traits<DPPType>::dime==StensorTraits<StensorType>::dime && 180 StensorTraits<PPType>::dime==StensorTraits<StensorType>::dime && 181 StensorTraits<StensorType>::dime==3u && 182 tfel::typetraits::IsAssignableTo<StensorNumType<StensorType>, 183 StensorNumType<PPType>>::cond && 184 tfel::typetraits::IsAssignableTo<tfel::typetraits::base_type<StensorNumType<StensorType>>, 185 ST2toST2NumType<DPPType>>::cond, 186 void>::type 187 computeStensorPositivePartAndDerivative(DPPType& dpp, 188 PPType& pp, 189 const StensorType& s, 190 const StensorNumType<StensorType> eps) 191 { 192 using tfel::math::internals::stensor_ppos; 193 using tfel::math::internals::stensor_pneg; 194 typedef StensorNumType<StensorType> NumType; 195 typedef tfel::typetraits::base_type<NumType> real; 196 constexpr const auto cste = Cste<real>::sqrt2; 197 TFEL_CONSTEXPR const auto one_half = real(1)/(real(2)); 198 stensor<3u,NumType> ls(s); // local copy of s 199 stensor<3u,real> n0; 200 stensor<3u,real> n1; 201 stensor<3u,real> n2; 202 tvector<3u,NumType> vp; 203 tmatrix<3u,3u,real> m; 204 ls.computeEigenVectors(vp,m); 205 stensor<3u,real>::computeEigenTensors(n0,n1,n2,m); 206 if((abs(vp(0)-vp(1))<eps)&&(abs(vp(0)-vp(2))<eps)){ 207 NumType vpm = (vp(0)+vp(1)+vp(2))/3; 208 if(abs(vpm)<eps){ 209 dpp = st2tost2<3u,real>::Id()*one_half; 210 pp = stensor<3u,NumType>(NumType(0)); 211 } else if(vpm>NumType(0)){ 212 dpp = st2tost2<3u,real>::Id(); 213 pp = s; 214 } else { 215 dpp = st2tost2<3u,real>(real(0)); 216 pp = stensor<3u,real>(real(0)); 217 } 218 } else { 219 const tvector<3u,real> v0 = m.template column_view<0u>(); 220 const tvector<3u,real> v1 = m.template column_view<1u>(); 221 const tvector<3u,real> v2 = m.template column_view<2u>(); 222 const stensor<3u,real> n01 = stensor<3u,real>::buildFromVectorsSymmetricDiadicProduct(v0,v1)/cste; 223 const stensor<3u,real> n02 = stensor<3u,real>::buildFromVectorsSymmetricDiadicProduct(v0,v2)/cste; 224 const stensor<3u,real> n12 = stensor<3u,real>::buildFromVectorsSymmetricDiadicProduct(v1,v2)/cste; 225 if(abs(vp(0)-vp(1))<eps){ 226 NumType vpm = (vp(0)+vp(1))*one_half; 227 if(abs(vpm)<eps){ 228 dpp = ((n0^n0)+(n1^n1)+(n01^n01))*one_half; 229 pp = stensor<3u,NumType>(NumType(0)); 230 } else if(vpm>NumType(0)){ 231 dpp = (n0^n0)+(n1^n1)+(n01^n01); 232 pp = vpm*(n0+n1); 233 } else { 234 dpp = st2tost2<3u,real>(real(0)); 235 pp = stensor<3u,NumType>(NumType(0)); 236 } 237 if(abs(vp(2))<eps){ 238 dpp += (n2^n2)*one_half; 239 } else if(vp(2)>NumType(0)){ 240 dpp += n2^n2; 241 pp += vp(2)*n2; 242 } 243 dpp += (stensor_ppos(vpm)-stensor_ppos(vp(2)))/(vpm-vp(2))*((n02^n02)+(n12^n12)); 244 } else if(abs(vp(0)-vp(2))<eps){ 245 NumType vpm = (vp(0)+vp(2))*one_half; 246 if(abs(vpm)<eps){ 247 dpp = ((n0^n0)+(n2^n2)+(n02^n02))*one_half; 248 pp = stensor<3u,NumType>(NumType(0)); 249 } else if(vpm>NumType(0)){ 250 dpp = (n0^n0)+(n2^n2)+(n02^n02); 251 pp = vpm*(n0+n2); 252 } else { 253 dpp = st2tost2<3u,real>(real(0)); 254 pp = stensor<3u,NumType>(NumType(0)); 255 } 256 if(abs(vp(1))<eps){ 257 dpp += (n1^n1)*one_half; 258 } else if(vp(1)>NumType(0)){ 259 dpp += (n1^n1); 260 pp += vp(1)*n1; 261 } 262 dpp += (stensor_ppos(vpm)-stensor_ppos(vp(1)))/(vpm-vp(1))*((n01^n01)+(n12^n12)); 263 } else if(abs(vp(1)-vp(2))<eps){ 264 NumType vpm = (vp(1)+vp(2))*one_half; 265 if(abs(vpm)<eps){ 266 dpp = ((n1^n1)+(n2^n2)+(n12^n12))*one_half; 267 pp = stensor<3u,NumType>(NumType(0)); 268 } else if(vpm>NumType(0)){ 269 dpp = (n1^n1)+(n2^n2)+(n12^n12); 270 pp = vpm*(n1+n2); 271 } else { 272 dpp = st2tost2<3u,real>(real(0)); 273 pp = stensor<3u,NumType>(NumType(0)); 274 } 275 if(abs(vp(0))<eps){ 276 dpp += (n0^n0)*one_half; 277 } else if(vp(0)>NumType(0)){ 278 dpp += (n0^n0); 279 pp += vp(0)*n0; 280 } 281 dpp += (stensor_ppos(vpm)-stensor_ppos(vp(0)))/(vpm-vp(0))*((n01^n01)+(n02^n02)); 282 } else { 283 // all eigenvalues are distincts 284 if(abs(vp(0))<eps){ 285 dpp = (n0^n0)*one_half; 286 pp = stensor<3u,NumType>(NumType(0)); 287 } else if (vp(0)>NumType(0)){ 288 dpp = (n0^n0)+vp(0)/(vp(0)-vp(1))*(n01^n01)+vp(0)/(vp(0)-vp(2))*(n02^n02); 289 pp = vp(0)*n0; 290 } else { 291 dpp = st2tost2<3u,real>(real(0)); 292 pp = stensor<3u,NumType>(NumType(0)); 293 } 294 if(abs(vp(1))<eps){ 295 dpp += (n1^n1)*one_half; 296 } else if (vp(1)>NumType(0)){ 297 dpp += (n1^n1)+vp(1)/(vp(1)-vp(0))*(n01^n01)+vp(1)/(vp(1)-vp(2))*(n12^n12); 298 pp += vp(1)*n1; 299 } else { 300 dpp += st2tost2<3u,real>(real(0)); 301 pp += stensor<3u,NumType>(NumType(0)); 302 } 303 if(abs(vp(2))<eps){ 304 dpp += (n2^n2)*one_half; 305 } else if (vp(2)>NumType(0)){ 306 dpp += (n2^n2)+vp(2)/(vp(2)-vp(0))*(n02^n02)+vp(2)/(vp(2)-vp(1))*(n12^n12); 307 pp += vp(2)*n2; 308 } else { 309 dpp += st2tost2<3u,real>(real(0)); 310 pp += stensor<3u,NumType>(NumType(0)); 311 } 312 } 313 } 314 } 315 316 template<typename DPPType, 317 typename DNPType, 318 typename PPType, 319 typename NPType, 320 typename StensorType> 321 typename std::enable_if< 322 tfel::meta::Implements<DPPType,ST2toST2Concept>::cond&& 323 tfel::meta::Implements<DNPType,ST2toST2Concept>::cond&& 324 tfel::meta::Implements<PPType,StensorConcept>::cond&& 325 tfel::meta::Implements<NPType,StensorConcept>::cond&& 326 tfel::meta::Implements<StensorType,StensorConcept>::cond&& 327 ST2toST2Traits<DPPType>::dime==StensorTraits<StensorType>::dime && 328 ST2toST2Traits<DNPType>::dime==StensorTraits<StensorType>::dime && 329 StensorTraits<PPType>::dime==StensorTraits<StensorType>::dime && 330 StensorTraits<NPType>::dime==StensorTraits<StensorType>::dime && 331 StensorTraits<StensorType>::dime==1u && 332 tfel::typetraits::IsAssignableTo<StensorNumType<StensorType>, 333 StensorNumType<PPType>>::cond && 334 tfel::typetraits::IsAssignableTo<StensorNumType<StensorType>, 335 StensorNumType<NPType>>::cond && 336 tfel::typetraits::IsAssignableTo<tfel::typetraits::base_type<StensorNumType<StensorType>>, 337 ST2toST2NumType<DPPType>>::cond && 338 tfel::typetraits::IsAssignableTo<tfel::typetraits::base_type<StensorNumType<StensorType>>, 339 ST2toST2NumType<DNPType>>::cond, 340 void>::type 341 computeStensorDecompositionInPositiveAndNegativeParts(DPPType& dpp, 342 DNPType& dnp, 343 PPType& pp, 344 NPType& np, 345 const StensorType& s, 346 const StensorNumType<StensorType> eps) 347 { 348 using std::abs; 349 typedef StensorNumType<StensorType> NumType; 350 typedef tfel::typetraits::base_type<NumType> real; 351 const real one_half = real(1)/(real(2)); 352 dpp(0,1) = dpp(0,2) = dpp(1,0) = dpp(1,2) = dpp(2,0) = dpp(2,1) = real(0); 353 dnp(0,1) = dnp(0,2) = dnp(1,0) = dnp(1,2) = dnp(2,0) = dnp(2,1) = real(0); 354 if(abs(s(0))<eps){ 355 dpp(0,0) = one_half; 356 dnp(0,0) = one_half; 357 pp(0)=np(0)=real(0); 358 } else if (s(0)>NumType(0)){ 359 dpp(0,0) = real(1); 360 dnp(0,0) = real(0); 361 pp(0)=s(0); 362 np(0)=NumType(0); 363 } else { 364 dpp(0,0) = real(0); 365 dnp(0,0) = real(1); 366 pp(0)=NumType(0); 367 np(0)=s(0); 368 } 369 if(abs(s(1))<eps){ 370 dpp(1,1) = one_half; 371 dnp(1,1) = one_half; 372 pp(1)=np(1)=real(0); 373 } else if (s(1)>NumType(0)){ 374 dpp(1,1) = real(1); 375 dnp(1,1) = real(0); 376 pp(1)=s(1); 377 np(1)=NumType(0); 378 } else { 379 dpp(1,1) = real(0); 380 dnp(1,1) = real(1); 381 pp(1)=NumType(0); 382 np(1)=s(1); 383 } 384 if(abs(s(2))<eps){ 385 dpp(2,2) = one_half; 386 dnp(2,2) = one_half; 387 pp(2)=np(2)=real(0); 388 } else if (s(2)>NumType(0)){ 389 dpp(2,2) = real(1); 390 dnp(2,2) = real(0); 391 pp(2)=s(2); 392 np(2)=NumType(0); 393 } else { 394 dpp(2,2) = real(0); 395 dnp(2,2) = real(1); 396 pp(2)=NumType(0); 397 np(2)=s(2); 398 } 399 } // end of computeStensorDecompositionInPositiveAndNegativeParts 400 401 template<typename DPPType, 402 typename DNPType, 403 typename PPType, 404 typename NPType, 405 typename StensorType> 406 typename std::enable_if< 407 tfel::meta::Implements<DPPType,ST2toST2Concept>::cond&& 408 tfel::meta::Implements<DNPType,ST2toST2Concept>::cond&& 409 tfel::meta::Implements<PPType,StensorConcept>::cond&& 410 tfel::meta::Implements<NPType,StensorConcept>::cond&& 411 tfel::meta::Implements<StensorType,StensorConcept>::cond&& 412 ST2toST2Traits<DPPType>::dime==StensorTraits<StensorType>::dime && 413 ST2toST2Traits<DNPType>::dime==StensorTraits<StensorType>::dime && 414 StensorTraits<PPType>::dime==StensorTraits<StensorType>::dime && 415 StensorTraits<NPType>::dime==StensorTraits<StensorType>::dime && 416 StensorTraits<StensorType>::dime==2u && 417 tfel::typetraits::IsAssignableTo<StensorNumType<StensorType>, 418 StensorNumType<PPType>>::cond && 419 tfel::typetraits::IsAssignableTo<StensorNumType<StensorType>, 420 StensorNumType<NPType>>::cond && 421 tfel::typetraits::IsAssignableTo<tfel::typetraits::base_type<StensorNumType<StensorType>>, 422 ST2toST2NumType<DPPType>>::cond && 423 tfel::typetraits::IsAssignableTo<tfel::typetraits::base_type<StensorNumType<StensorType>>, 424 ST2toST2NumType<DNPType>>::cond, 425 void>::type 426 computeStensorDecompositionInPositiveAndNegativeParts(DPPType& dpp, 427 DNPType& dnp, 428 PPType& pp, 429 NPType& np, 430 const StensorType& s, 431 const StensorNumType<StensorType> eps) 432 { 433 using std::abs; 434 using tfel::math::internals::stensor_ppos; 435 using tfel::math::internals::stensor_pneg; 436 typedef StensorNumType<StensorType> NumType; 437 typedef tfel::typetraits::base_type<NumType> real; 438 constexpr const auto cste = Cste<real>::sqrt2; 439 TFEL_CONSTEXPR const auto one_half = real(1)/(real(2)); 440 stensor<2u,NumType> ls(s); // local copy of s 441 stensor<2u,real> n0; 442 stensor<2u,real> n1; 443 stensor<2u,real> n2; 444 tvector<3u,NumType> vp; 445 tmatrix<3u,3u,real> m; 446 ls.computeEigenVectors(vp,m); 447 stensor<2u,real>::computeEigenTensors(n0,n1,n2,m); 448 const tvector<3u,real> v0 = m.template column_view<0u>(); 449 const tvector<3u,real> v1 = m.template column_view<1u>(); 450 const stensor<2u,real> n01 = stensor<2u,real>::buildFromVectorsSymmetricDiadicProduct(v0,v1)/cste; 451 if(abs(vp(0)-vp(1))<eps){ 452 const NumType vpm = (vp(0)+vp(1))*one_half; 453 if(abs(vpm)<eps){ 454 dpp = dnp = ((n0^n0)+(n1^n1)+(n01^n01))*one_half; 455 pp = np = stensor<2u,NumType>(real(0)); 456 } else if(vpm>NumType(0)){ 457 dpp = (n0^n0)+(n1^n1)+(n01^n01); 458 dnp = st2tost2<2u,real>(real(0)); 459 pp = vpm*(n0+n1); 460 np = stensor<2u,NumType>(real(0)); 461 } else { 462 dpp = st2tost2<2u,real>(real(0)); 463 dnp = (n0^n0)+(n1^n1)+(n01^n01); 464 pp = stensor<2u,NumType>(real(0)); 465 np = vpm*(n0+n1); 466 } 467 } else { 468 if(abs(vp(0))<eps){ 469 dpp = dnp = (n0^n0)*one_half; 470 pp = np = stensor<2u,NumType>(real(0)); 471 } else if(vp(0)>NumType(0)){ 472 dpp = (n0^n0); 473 dnp = st2tost2<2u,real>(real(0)); 474 pp = vp(0)*n0; 475 np = stensor<2u,NumType>(real(0)); 476 } else { 477 dpp = st2tost2<2u,real>(real(0)); 478 dnp = (n0^n0); 479 pp = stensor<2u,NumType>(real(0)); 480 np = vp(0)*n0; 481 } 482 if(abs(vp(1))<eps){ 483 dpp += (n1^n1)*one_half; 484 dnp += (n1^n1)*one_half; 485 } else if(vp(1)>NumType(0)){ 486 dpp += (n1^n1); 487 pp += vp(1)*n1; 488 } else { 489 dnp += (n1^n1); 490 np += vp(1)*n1; 491 } 492 dpp += (stensor_ppos(vp(1))-stensor_ppos(vp(0)))/(vp(1)-vp(0))*(n01^n01); 493 dnp += (stensor_pneg(vp(1))-stensor_pneg(vp(0)))/(vp(1)-vp(0))*(n01^n01); 494 } 495 if(abs(vp(2))<eps){ 496 dpp += (n2^n2)*one_half; 497 dnp += (n2^n2)*one_half; 498 } else if (vp(2)>0){ 499 dpp += (n2^n2); 500 pp(2) += vp(2); 501 } else { 502 dnp += (n2^n2); 503 np(2) += vp(2); 504 } 505 } // end of computeStensorDecompositionInPositiveAndNegativeParts 506 507 template<typename DPPType, 508 typename DNPType, 509 typename PPType, 510 typename NPType, 511 typename StensorType> 512 typename std::enable_if< 513 tfel::meta::Implements<DPPType,ST2toST2Concept>::cond&& 514 tfel::meta::Implements<DNPType,ST2toST2Concept>::cond&& 515 tfel::meta::Implements<PPType,StensorConcept>::cond&& 516 tfel::meta::Implements<NPType,StensorConcept>::cond&& 517 tfel::meta::Implements<StensorType,StensorConcept>::cond&& 518 ST2toST2Traits<DPPType>::dime==StensorTraits<StensorType>::dime && 519 ST2toST2Traits<DNPType>::dime==StensorTraits<StensorType>::dime && 520 StensorTraits<PPType>::dime==StensorTraits<StensorType>::dime && 521 StensorTraits<NPType>::dime==StensorTraits<StensorType>::dime && 522 StensorTraits<StensorType>::dime==3u && 523 tfel::typetraits::IsAssignableTo<StensorNumType<StensorType>, 524 StensorNumType<PPType>>::cond && 525 tfel::typetraits::IsAssignableTo<StensorNumType<StensorType>, 526 StensorNumType<NPType>>::cond && 527 tfel::typetraits::IsAssignableTo<tfel::typetraits::base_type<StensorNumType<StensorType>>, 528 ST2toST2NumType<DPPType>>::cond && 529 tfel::typetraits::IsAssignableTo<tfel::typetraits::base_type<StensorNumType<StensorType>>, 530 ST2toST2NumType<DNPType>>::cond, 531 void>::type 532 computeStensorDecompositionInPositiveAndNegativeParts(DPPType& dpp, 533 DNPType& dnp, 534 PPType& pp, 535 NPType& np, 536 const StensorType& s, 537 const StensorNumType<StensorType> eps) 538 { 539 using tfel::math::internals::stensor_ppos; 540 using tfel::math::internals::stensor_pneg; 541 typedef StensorNumType<StensorType> NumType; 542 typedef tfel::typetraits::base_type<NumType> real; 543 constexpr const auto cste = Cste<real>::sqrt2; 544 TFEL_CONSTEXPR const auto one_half = real(1)/(real(2)); 545 stensor<3u,NumType> ls(s); // local copy of s 546 stensor<3u,real> n0; 547 stensor<3u,real> n1; 548 stensor<3u,real> n2; 549 tvector<3u,NumType> vp; 550 tmatrix<3u,3u,real> m; 551 ls.computeEigenVectors(vp,m); 552 stensor<3u,real>::computeEigenTensors(n0,n1,n2,m); 553 if((abs(vp(0)-vp(1))<eps)&&(abs(vp(0)-vp(2))<eps)){ 554 NumType vpm = (vp(0)+vp(1)+vp(2))/3; 555 if(abs(vpm)<eps){ 556 dpp = dnp = st2tost2<3u,real>::Id()*one_half; 557 pp = np = stensor<3u,NumType>(NumType(0)); 558 } else if(vpm>NumType(0)){ 559 dpp = st2tost2<3u,real>::Id(); 560 pp = s; 561 dnp = st2tost2<3u,real>(real(0)); 562 np = stensor<3u,real>(real(0)); 563 } else { 564 dpp = st2tost2<3u,real>(real(0)); 565 pp = stensor<3u,real>(real(0)); 566 dnp = st2tost2<3u,real>::Id(); 567 np = s; 568 } 569 } else { 570 const tvector<3u,real> v0 = m.template column_view<0u>(); 571 const tvector<3u,real> v1 = m.template column_view<1u>(); 572 const tvector<3u,real> v2 = m.template column_view<2u>(); 573 const stensor<3u,real> n01 = stensor<3u,real>::buildFromVectorsSymmetricDiadicProduct(v0,v1)/cste; 574 const stensor<3u,real> n02 = stensor<3u,real>::buildFromVectorsSymmetricDiadicProduct(v0,v2)/cste; 575 const stensor<3u,real> n12 = stensor<3u,real>::buildFromVectorsSymmetricDiadicProduct(v1,v2)/cste; 576 if(abs(vp(0)-vp(1))<eps){ 577 NumType vpm = (vp(0)+vp(1))*one_half; 578 if(abs(vpm)<eps){ 579 dpp = dnp = ((n0^n0)+(n1^n1)+(n01^n01))*one_half; 580 pp = np = stensor<3u,NumType>(NumType(0)); 581 } else if(vpm>NumType(0)){ 582 dpp = (n0^n0)+(n1^n1)+(n01^n01); 583 dnp = st2tost2<3u,real>(real(0)); 584 pp = vpm*(n0+n1); 585 np = stensor<3u,NumType>(NumType(0)); 586 } else { 587 dpp = st2tost2<3u,real>(real(0)); 588 dnp = (n0^n0)+(n1^n1)+(n01^n01); 589 pp = stensor<3u,NumType>(NumType(0)); 590 np = vpm*(n0+n1); 591 } 592 if(abs(vp(2))<eps){ 593 dpp += (n2^n2)*one_half; 594 dnp += (n2^n2)*one_half; 595 } else if(vp(2)>NumType(0)){ 596 dpp += n2^n2; 597 pp += vp(2)*n2; 598 } else { 599 dnp += n2^n2; 600 np += vp(2)*n2; 601 } 602 dpp += (stensor_ppos(vpm)-stensor_ppos(vp(2)))/(vpm-vp(2))*((n02^n02)+(n12^n12)); 603 dnp += (stensor_pneg(vpm)-stensor_pneg(vp(2)))/(vpm-vp(2))*((n02^n02)+(n12^n12)); 604 } else if(abs(vp(0)-vp(2))<eps){ 605 NumType vpm = (vp(0)+vp(2))*one_half; 606 if(abs(vpm)<eps){ 607 dpp = dnp = ((n0^n0)+(n2^n2)+(n02^n02))*one_half; 608 pp = np = stensor<3u,NumType>(NumType(0)); 609 } else if(vpm>NumType(0)){ 610 dpp = (n0^n0)+(n2^n2)+(n02^n02); 611 dnp = st2tost2<3u,real>(real(0)); 612 pp = vpm*(n0+n2); 613 np = stensor<3u,NumType>(NumType(0)); 614 } else { 615 dpp = st2tost2<3u,real>(real(0)); 616 dnp = (n0^n0)+(n2^n2)+(n02^n02); 617 pp = stensor<3u,NumType>(NumType(0)); 618 np = vpm*(n0+n2); 619 } 620 if(abs(vp(1))<eps){ 621 dpp += (n1^n1)*one_half; 622 dnp += (n1^n1)*one_half; 623 } else if(vp(1)>NumType(0)){ 624 dpp += (n1^n1); 625 pp += vp(1)*n1; 626 } else { 627 dnp += (n1^n1); 628 np += vp(1)*n1; 629 } 630 dpp += (stensor_ppos(vpm)-stensor_ppos(vp(1)))/(vpm-vp(1))*((n01^n01)+(n12^n12)); 631 dnp += (stensor_pneg(vpm)-stensor_pneg(vp(1)))/(vpm-vp(1))*((n01^n01)+(n12^n12)); 632 } else if(abs(vp(1)-vp(2))<eps){ 633 NumType vpm = (vp(1)+vp(2))*one_half; 634 if(abs(vpm)<eps){ 635 dpp = dnp = ((n1^n1)+(n2^n2)+(n12^n12))*one_half; 636 pp = np = stensor<3u,NumType>(NumType(0)); 637 } else if(vpm>NumType(0)){ 638 dpp = (n1^n1)+(n2^n2)+(n12^n12); 639 dnp = st2tost2<3u,real>(real(0)); 640 pp = vpm*(n1+n2); 641 np = stensor<3u,NumType>(NumType(0)); 642 } else { 643 dpp = st2tost2<3u,real>(real(0)); 644 dnp = (n1^n1)+(n2^n2)+(n12^n12); 645 pp = stensor<3u,NumType>(NumType(0)); 646 np = vpm*(n1+n2); 647 } 648 if(abs(vp(0))<eps){ 649 dpp += (n0^n0)*one_half; 650 dnp += (n0^n0)*one_half; 651 } else if(vp(0)>NumType(0)){ 652 dpp += (n0^n0); 653 pp += vp(0)*n0; 654 } else { 655 dnp += (n0^n0); 656 np += vp(0)*n0; 657 } 658 dpp += (stensor_ppos(vpm)-stensor_ppos(vp(0)))/(vpm-vp(0))*((n01^n01)+(n02^n02)); 659 dnp += (stensor_pneg(vpm)-stensor_pneg(vp(0)))/(vpm-vp(0))*((n01^n01)+(n02^n02)); 660 } else { 661 // all eigenvalues are distincts 662 if(abs(vp(0))<eps){ 663 dpp = (n0^n0)*one_half; 664 dnp = (n0^n0)*one_half; 665 pp = stensor<3u,NumType>(NumType(0)); 666 np = stensor<3u,NumType>(NumType(0)); 667 } else if (vp(0)>NumType(0)){ 668 dpp = (n0^n0)+vp(0)/(vp(0)-vp(1))*(n01^n01)+vp(0)/(vp(0)-vp(2))*(n02^n02); 669 dnp = st2tost2<3u,real>(real(0)); 670 pp = vp(0)*n0; 671 np = stensor<3u,NumType>(NumType(0)); 672 } else { 673 dpp = st2tost2<3u,real>(real(0)); 674 dnp = (n0^n0)+vp(0)/(vp(0)-vp(1))*(n01^n01)+vp(0)/(vp(0)-vp(2))*(n02^n02); 675 pp = stensor<3u,NumType>(NumType(0)); 676 np = vp(0)*n0; 677 } 678 if(abs(vp(1))<eps){ 679 dpp += (n1^n1)*one_half; 680 dnp += (n1^n1)*one_half; 681 } else if (vp(1)>NumType(0)){ 682 dpp += (n1^n1)+vp(1)/(vp(1)-vp(0))*(n01^n01)+vp(1)/(vp(1)-vp(2))*(n12^n12); 683 dnp += st2tost2<3u,real>(real(0)); 684 pp += vp(1)*n1; 685 np += stensor<3u,NumType>(NumType(0)); 686 } else { 687 dpp += st2tost2<3u,real>(real(0)); 688 dnp += (n1^n1)+vp(1)/(vp(1)-vp(0))*(n01^n01)+vp(1)/(vp(1)-vp(2))*(n12^n12); 689 pp += stensor<3u,NumType>(NumType(0)); 690 np += vp(1)*n1; 691 } 692 if(abs(vp(2))<eps){ 693 dpp += (n2^n2)*one_half; 694 dnp += (n2^n2)*one_half; 695 } else if (vp(2)>NumType(0)){ 696 dpp += (n2^n2)+vp(2)/(vp(2)-vp(0))*(n02^n02)+vp(2)/(vp(2)-vp(1))*(n12^n12); 697 dnp += st2tost2<3u,real>(real(0)); 698 pp += vp(2)*n2; 699 np += stensor<3u,NumType>(NumType(0)); 700 } else { 701 dpp += st2tost2<3u,real>(real(0)); 702 dnp += (n2^n2)+vp(2)/(vp(2)-vp(1))*(n02^n02)+vp(2)/(vp(2)-vp(1))*(n12^n12); 703 pp += stensor<3u,NumType>(NumType(0)); 704 np += vp(2)*n2; 705 } 706 } 707 } 708 } // end of computeStensorDecompositionInPositiveAndNegativeParts 709 710 } // end of namespace math 711} // end of namespace tfel 712 713#endif /* LIB_TFEL_MATH_STENSOR_DECOMPOSITIONINPOSITIVEANDNEGATIVEPARTSIXX */ 714