src/foam/primitives/VectorN/TensorNI.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  27
  28 namespace Foam
  29 {
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 template <class Cmpt, int length>
  34 const char* const TensorN<Cmpt, length>::typeName =
  35     ("tensor" + name(length)).c_str();
  36
  37 template <class Cmpt, int length>
  38 const TensorN<Cmpt, length> TensorN<Cmpt, length>::zero(0);
  39
  40 template <class Cmpt, int length>
  41 const TensorN<Cmpt, length> TensorN<Cmpt, length>::one(1);
  42
  43
  44 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  45
  46 //- Construct null
  47 template <class Cmpt, int length>
  48 inline TensorN<Cmpt, length>::TensorN()
  49 {}
  50
  51
  52 //- Construct given VectorSpace
  53 template <class Cmpt, int length>
  54 inline TensorN<Cmpt, length>::TensorN
  55 (
  56     const VectorSpace<TensorN<Cmpt, length>, Cmpt, length*length>& vs
  57 )
  58 :
  59     VectorSpace<TensorN<Cmpt, length>, Cmpt, length*length>(vs)
  60 {}
  61
  62
  63 //- Construct from component
  64 template <class Cmpt, int length>
  65 inline TensorN<Cmpt, length>::TensorN(const Cmpt& tx)
  66 {
  67     VectorSpaceOps<TensorN<Cmpt, length>::nComponents,0>::eqOpS
  68     (
  69         *this,
  70         tx,
  71         eqOp<Cmpt>()
  72     );
  73 }
  74
  75
  76 //- Construct from Istream
  77 template <class Cmpt, int length>
  78 inline TensorN<Cmpt, length>::TensorN(Istream& is)
  79 :
  80     VectorSpace<TensorN<Cmpt, length>, Cmpt, length*length>(is)
  81 {}
  82
  83
  84 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  85
  86 //- Return tensor transpose
  87 template <class Cmpt, int length>
  88 inline TensorN<Cmpt, length> TensorN<Cmpt, length>::T() const
  89 {
  90     TensorN<Cmpt, length> transpose;
  91
  92     label i = 0;
  93     for (label row = 0; row < TensorN<Cmpt, length>::rowLength; row++)
  94     {
  95         label j = row;
  96         for (label col = 0; col < TensorN<Cmpt, length>::rowLength; col++)
  97         {
  98             transpose.v_[i] = this->v_[j];
  99             i++;
 100             j += TensorN<Cmpt, length>::rowLength;
 101         }
 102     }
 103
 104     return transpose;
 105 }
 106
 107 //- Return tensor diagonal
 108 template <class Cmpt, int length>
 109 inline DiagTensorN<Cmpt, length> TensorN<Cmpt, length>::diag() const
 110 {
 111     DiagTensorN<Cmpt, length> dt;
 112
 113     label diagI=0;
 114     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 115     {
 116         dt[i] = this->v_[diagI];
 117         diagI += TensorN<Cmpt, length>::rowLength + 1;
 118     }
 119
 120     return dt;
 121 }
 122
 123 //- Negative sum the vertical off-diagonal components
 124 template <class Cmpt, int length>
 125 inline TensorN<Cmpt, length> TensorN<Cmpt, length>::negSumDiag() const
 126 {
 127     TensorN<Cmpt, length> negsumdiag;
 128
 129     // Zero main diagonal
 130     label diagI = 0;
 131     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 132     {
 133         negsumdiag.v_[diagI] = 0.0;
 134         diagI += TensorN<Cmpt, length>::rowLength + 1;
 135     }
 136
 137     label k = 0;
 138     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 139     {
 140         diagI = 0;
 141         for (label j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
 142         {
 143             if (k != diagI)
 144             {
 145                 negsumdiag.v_[k] = this->v_[k];
 146                 negsumdiag.v_[diagI] -= this->v_[k];
 147             }
 148             k++;
 149             diagI += TensorN<Cmpt, length>::rowLength + 1;
 150         }
 151     }
 152
 153     return negsumdiag;
 154 }
 155
 156 //- Assign to a SphericalTensorN
 157 template <class Cmpt, int length>
 158 inline void
 159 TensorN<Cmpt, length>::operator=(const SphericalTensorN<Cmpt, length>& st)
 160 {
 161     label diag = 0;
 162     for (label i = 0; i < TensorN<Cmpt, length>::nComponents; i++)
 163     {
 164         if (i == diag)
 165         {
 166             this->v_[i] = st[0];
 167             diag += TensorN<Cmpt, length>::rowLength + 1;
 168         }
 169         else
 170         {
 171             this->v_[i] = pTraits<Cmpt>::zero;
 172         }
 173     }
 174 }
 175
 176
 177 //- Assign to a DiagTensorN
 178 template <class Cmpt, int length>
 179 inline void
 180 TensorN<Cmpt, length>::operator=(const DiagTensorN<Cmpt, length>& dt)
 181 {
 182     label diag = 0;
 183     label k = 0;
 184     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 185     {
 186         for (label j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
 187         {
 188             if (j == diag)
 189             {
 190                 this->v_[k] = dt[i];
 191             }
 192             else
 193             {
 194                 this->v_[k] = pTraits<Cmpt>::zero;
 195             }
 196             k++;
 197         }
 198         diag++;
 199     }
 200 }
 201
 202
 203 //- Transform the tensor
 204 //- The components are assumed to be individual scalars
 205 //- i.e. transform has no effect
 206 template<class Cmpt, int length>
 207 inline TensorN<Cmpt, length> transform
 208 (
 209     const tensor& tt,
 210     const TensorN<Cmpt, length>& v
 211 )
 212 {
 213     return v;
 214 }
 215
 216
 217 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 218
 219 template <class Cmpt, int length>
 220 inline direction TensorN<Cmpt, length>::cmpt
 221 (
 222     const direction i,
 223     const direction j
 224 )
 225 {
 226     if
 227     (
 228         i > TensorN<Cmpt, length>::rowLength - 1
 229      || j > TensorN<Cmpt, length>::rowLength - 1
 230     )
 231     {
 232         FatalErrorIn
 233         (
 234             "direction TensorN<Cmpt, length>::cmpt()\n"
 235             "(\n"
 236             "    const direction i,\n"
 237             "    const direction j\n"
 238             ") const"
 239         )   << "Direction out of range (0 "
 240             << TensorN<Cmpt, length>::rowLength - 1 << ")"
 241             << " and (i j) = ("  << i << " " << j << ")"
 242             << abort(FatalError);
 243     }
 244
 245     return i*TensorN<Cmpt, length>::rowLength + j;
 246 }
 247
 248
 249 template <class Cmpt, int length>
 250 inline const Cmpt& TensorN<Cmpt, length>::operator()
 251 (
 252     const direction i,
 253     const direction j
 254 ) const
 255 {
 256     return this->operator[](i*TensorN<Cmpt, length>::rowLength + j);
 257 }
 258
 259
 260 template <class Cmpt, int length>
 261 inline Cmpt& TensorN<Cmpt, length>::operator()
 262 (
 263     const direction i,
 264     const direction j
 265 )
 266 {
 267     return this->operator[](i*TensorN<Cmpt, length>::rowLength + j);
 268 }
 269
 270
 271 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 272
 273 //- Inner-product between two tensors
 274 template <class Cmpt, int length>
 275 inline typename
 276 innerProduct<TensorN<Cmpt, length>, TensorN<Cmpt, length> >::type
 277 operator&(const TensorN<Cmpt, length>& t1, const TensorN<Cmpt, length>& t2)
 278 {
 279     TensorN<Cmpt, length> result(TensorN<Cmpt, length>::zero);
 280
 281     label i = 0;
 282     label j = 0;
 283
 284     label m, n;
 285
 286     for (label row = 0; row < TensorN<Cmpt, length>::rowLength; row++)
 287     {
 288         for (label col = 0; col < TensorN<Cmpt, length>::rowLength; col++)
 289         {
 290             Cmpt& r = result.v_[i];
 291             m = j;
 292             n = col;
 293
 294             for
 295             (
 296                 label row2 = 0;
 297                 row2 < TensorN<Cmpt, length>::rowLength;
 298                 row2++
 299             )
 300             {
 301                 r += t1.v_[m]*t2.v_[n];
 302                 m++;
 303                 n += TensorN<Cmpt, length>::rowLength;
 304             }
 305             i++;
 306         }
 307         j += TensorN<Cmpt, length>::rowLength;
 308     }
 309
 310     return result;
 311 }
 312
 313 //- Inner-product between diagonal tensor and tensor
 314 template <class Cmpt, int length>
 315 inline typename
 316 innerProduct<DiagTensorN<Cmpt, length>, TensorN<Cmpt, length> >::type
 317 operator&
 318 (
 319     const DiagTensorN<Cmpt, length>& dt1,
 320     const TensorN<Cmpt, length>& t2
 321 )
 322 {
 323     TensorN<Cmpt, length> result;
 324
 325     label k = 0;
 326     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 327     {
 328         const Cmpt& xx = dt1.v_[i];
 329
 330         for (label j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
 331         {
 332             result.v_[k] = xx*t2.v_[k];
 333             k++;
 334         }
 335     }
 336
 337     return result;
 338 }
 339
 340 //- Inner-product between tensor and diagonal tensor
 341 template <class Cmpt, int length>
 342 inline typename
 343 innerProduct<TensorN<Cmpt, length>, DiagTensorN<Cmpt, length> >::type
 344 operator&
 345 (
 346     const TensorN<Cmpt, length>& t1,
 347     const DiagTensorN<Cmpt, length>& dt2
 348 )
 349 {
 350     TensorN<Cmpt, length> result;
 351
 352     label k = 0;
 353     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 354     {
 355         for (label j = 0; j < TensorN<Cmpt, length>::rowLength; j++)
 356         {
 357             result.v_[k] = t1.v_[k]*dt2.v_[j];
 358             k++;
 359         }
 360     }
 361
 362     return result;
 363 }
 364
 365
 366 //- Inner-product between spherical tensor and tensor
 367 template <class Cmpt, int length>
 368 inline typename
 369 innerProduct<SphericalTensorN<Cmpt, length>, TensorN<Cmpt, length> >::type
 370 operator&
 371 (
 372     const SphericalTensorN<Cmpt, length>& st1,
 373     const TensorN<Cmpt, length>& t2
 374 )
 375 {
 376     const Cmpt& s = st1.v_[0];
 377     TensorN<Cmpt, length> res;
 378     VectorSpaceOps<TensorN<Cmpt, length>::nComponents,0>::opSV
 379     (
 380         res,
 381         s,
 382         t2,
 383         multiplyOp<Cmpt>()
 384     );
 385
 386     return res;
 387 }
 388
 389
 390 //- Inner-product between tensor and spherical tensor
 391 template <class Cmpt, int length>
 392 inline typename
 393 innerProduct<TensorN<Cmpt, length>, SphericalTensorN<Cmpt, length> >::type
 394 operator&
 395 (
 396     const TensorN<Cmpt, length>& t1,
 397     const SphericalTensorN<Cmpt, length>& st2
 398 )
 399 {
 400     const Cmpt& s = st2.v_[0];
 401     TensorN<Cmpt, length> res;
 402     VectorSpaceOps<TensorN<Cmpt, length>::nComponents,0>::opVS
 403     (
 404         res,
 405         t1,
 406         s,
 407         multiplyOp<Cmpt>()
 408     );
 409
 410     return res;
 411 }
 412
 413
 414 //- Inner-product between a tensor and a vector
 415 template <class Cmpt, int length>
 416 inline typename
 417 innerProduct<TensorN<Cmpt, length>, VectorN<Cmpt, length> >::type
 418 operator&(const TensorN<Cmpt, length>& t, const VectorN<Cmpt, length>& v)
 419 {
 420     VectorN<Cmpt, length> result(VectorN<Cmpt, length>::zero);
 421
 422     label i = 0;
 423     for (label row = 0; row < TensorN<Cmpt, length>::rowLength; row++)
 424     {
 425         Cmpt& r = result.v_[row];
 426
 427         for (label col = 0; col < TensorN<Cmpt, length>::rowLength; col++)
 428         {
 429             r += t.v_[i]*v.v_[col];
 430             i++;
 431         }
 432     }
 433
 434     return result;
 435 }
 436
 437
 438 //- Inner-product between a vector and a tensor
 439 template <class Cmpt, int length>
 440 inline typename
 441 innerProduct<VectorN<Cmpt, length>, TensorN<Cmpt, length> >::type
 442 operator&(const VectorN<Cmpt, length>& v, const TensorN<Cmpt, length>& t)
 443 {
 444     VectorN<Cmpt, length> result(VectorN<Cmpt, length>::zero);
 445
 446     for (label col = 0; col < TensorN<Cmpt, length>::rowLength; col++)
 447     {
 448         label j = col;
 449         Cmpt& r = result.v_[col];
 450
 451         for (label row = 0; row < TensorN<Cmpt, length>::rowLength; row++)
 452         {
 453             r += v.v_[row]*t.v_[j];
 454             j += TensorN<Cmpt, length>::rowLength;
 455         }
 456     }
 457
 458     return result;
 459 }
 460
 461
 462 //- Outer-product between two vectors
 463 template <class Cmpt, int length>
 464 inline typename
 465 outerProduct<VectorN<Cmpt, length>, VectorN<Cmpt, length> >::type
 466 operator*(const VectorN<Cmpt, length>& v1, const VectorN<Cmpt, length>& v2)
 467 {
 468     TensorN<Cmpt, length> result(TensorN<Cmpt, length>::zero);
 469
 470     label i = 0;
 471     for (label row = 0; row < TensorN<Cmpt, length>::rowLength; row++)
 472     {
 473         for (label col = 0; col < TensorN<Cmpt, length>::rowLength; col++)
 474         {
 475             result.v_[i] = v1.v_[row]*v2.v_[col];
 476             i++;
 477         }
 478     }
 479
 480     return result;
 481 }
 482
 483
 484 //- Addition of TensorN and TensorN
 485 template <class Cmpt, int length>
 486 inline TensorN<Cmpt, length>
 487 operator+(const TensorN<Cmpt, length>& t1, const TensorN<Cmpt, length>& t2)
 488 {
 489     TensorN<Cmpt, length> res;
 490     VectorSpaceOps<TensorN<Cmpt, length>::nComponents,0>::op
 491     (
 492         res,
 493         t1,
 494         t2,
 495         plusOp<Cmpt>()
 496     );
 497
 498     return res;
 499 }
 500
 501
 502 //- Addition of TensorN and DiagTensorN
 503 template <class Cmpt, int length>
 504 inline TensorN<Cmpt, length>
 505 operator+(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
 506 {
 507     TensorN<Cmpt, length> result(t1);
 508
 509     label diag = 0;
 510     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 511     {
 512         result.v_[diag] += dt2.v_[i];
 513         diag += TensorN<Cmpt, length>::rowLength + 1;
 514     }
 515
 516     return result;
 517 }
 518
 519
 520 //- Addition of DiagTensorN and TensorN
 521 template <class Cmpt, int length>
 522 inline TensorN<Cmpt, length>
 523 operator+(const DiagTensorN<Cmpt, length>& dt1, const TensorN<Cmpt, length>& t2)
 524 {
 525     TensorN<Cmpt, length> result(t2);
 526
 527     label diag = 0;
 528     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 529     {
 530         result.v_[diag] += dt1.v_[i];
 531         diag += TensorN<Cmpt, length>::rowLength + 1;
 532     }
 533
 534     return result;
 535 }
 536
 537
 538 //- Addition of TensorN and SphericalTensorN
 539 template <class Cmpt, int length>
 540 inline TensorN<Cmpt, length>
 541 operator+
 542 (
 543     const TensorN<Cmpt, length>& t1,
 544     const SphericalTensorN<Cmpt, length>& st2
 545 )
 546 {
 547     TensorN<Cmpt, length> result(t1);
 548
 549     const Cmpt& s = st2.v_[0];
 550     label diag = 0;
 551     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 552     {
 553         result.v_[diag] += s;
 554         diag += TensorN<Cmpt, length>::rowLength + 1;
 555     }
 556
 557     return result;
 558 }
 559
 560
 561 //- Addition of SphericalTensorN and TensorN
 562 template <class Cmpt, int length>
 563 inline TensorN<Cmpt, length>
 564 operator+
 565 (
 566     const SphericalTensorN<Cmpt, length>& st1,
 567     const TensorN<Cmpt, length>& t2
 568 )
 569 {
 570     TensorN<Cmpt, length> result(t2);
 571
 572     const Cmpt& s = st1.v_[0];
 573     label diag = 0;
 574     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 575     {
 576         result.v_[diag] += s;
 577         diag += TensorN<Cmpt, length>::rowLength + 1;
 578     }
 579
 580     return result;
 581 }
 582
 583
 584 //- Subtraction of TensorN and TensorN
 585 template <class Cmpt, int length>
 586 inline TensorN<Cmpt, length>
 587 operator-(const TensorN<Cmpt, length>& t1, const TensorN<Cmpt, length>& t2)
 588 {
 589     TensorN<Cmpt, length> res;
 590     VectorSpaceOps<TensorN<Cmpt, length>::nComponents,0>::op
 591     (
 592         res,
 593         t1,
 594         t2,
 595         minusOp<Cmpt>()
 596     );
 597
 598     return res;
 599 }
 600
 601
 602 //- Subtraction of TensorN and DiagTensorN
 603 template <class Cmpt, int length>
 604 inline TensorN<Cmpt, length>
 605 operator-(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
 606 {
 607     TensorN<Cmpt, length> result(t1);
 608
 609     label diag = 0;
 610     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 611     {
 612         result.v_[diag] -= dt2.v_[i];
 613         diag += TensorN<Cmpt, length>::rowLength + 1;
 614     }
 615
 616     return result;
 617 }
 618
 619
 620 //- Subtraction of DiagTensorN and TensorN
 621 template <class Cmpt, int length>
 622 inline TensorN<Cmpt, length>
 623 operator-(const DiagTensorN<Cmpt, length>& dt1, const TensorN<Cmpt, length>& t2)
 624 {
 625     TensorN<Cmpt, length> result(-t2);
 626
 627     label diag = 0;
 628     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 629     {
 630         result.v_[diag] += dt1.v_[i];
 631         diag += TensorN<Cmpt, length>::rowLength + 1;
 632     }
 633
 634     return result;
 635 }
 636
 637
 638 //- Subtraction of TensorN and SphericalTensorN
 639 template <class Cmpt, int length>
 640 inline TensorN<Cmpt, length>
 641 operator-
 642 (
 643     const TensorN<Cmpt, length>& t1,
 644     const SphericalTensorN<Cmpt, length>& st2
 645 )
 646 {
 647     TensorN<Cmpt, length> result(t1);
 648
 649     const Cmpt& s = st2.v_[0];
 650     label diag = 0;
 651     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 652     {
 653         result.v_[diag] -= s;
 654         diag += TensorN<Cmpt, length>::rowLength + 1;
 655     }
 656
 657     return result;
 658 }
 659
 660
 661 //- Subtraction of SphericalTensorN and TensorN
 662 template <class Cmpt, int length>
 663 inline TensorN<Cmpt, length>
 664 operator-
 665 (
 666     const SphericalTensorN<Cmpt, length>& st1,
 667     const TensorN<Cmpt, length>& t2
 668 )
 669 {
 670     TensorN<Cmpt, length> result(-t2);
 671
 672     const Cmpt& s = st1.v_[0];
 673     label diag = 0;
 674     for (label i = 0; i < TensorN<Cmpt, length>::rowLength; i++)
 675     {
 676         result.v_[diag] += s;
 677         diag += TensorN<Cmpt, length>::rowLength + 1;
 678     }
 679
 680     return result;
 681 }
 682
 683
 684 //- Division of a scalar by a TensorN
 685 template <class Cmpt, int length>
 686 inline TensorN<Cmpt, length>
 687 operator/(const scalar s, const TensorN<Cmpt, length>& t)
 688 {
 689     return s*inv(t);
 690 }
 691
 692 //- Inner Product of a VectorN by an inverse TensorN
 693 template <class Cmpt, int length>
 694 inline VectorN<Cmpt, length>
 695 operator/(const VectorN<Cmpt, length>& v, const TensorN<Cmpt, length>& t)
 696 {
 697     return v & inv(t);
 698 }
 699
 700 //- Inner Product of a TensorN by an inverse TensorN
 701 template <class Cmpt, int length>
 702 inline TensorN<Cmpt, length>
 703 operator/(const TensorN<Cmpt, length>& t1, const TensorN<Cmpt, length>& t2)
 704 {
 705     return t1 & inv(t2);
 706 }
 707
 708
 709 //- Inner Product of a DiagTensorN and an inverse TensorN
 710 template <class Cmpt, int length>
 711 inline TensorN<Cmpt, length>
 712 operator/(const DiagTensorN<Cmpt, length>& dt1, const TensorN<Cmpt, length>& t2)
 713 {
 714     return dt1 & inv(t2);
 715 }
 716
 717
 718 //- Inner Product of a TensorN and an inverse DiagTensorN
 719 template <class Cmpt, int length>
 720 inline TensorN<Cmpt, length>
 721 operator/(const TensorN<Cmpt, length>& t1, const DiagTensorN<Cmpt, length>& dt2)
 722 {
 723     return t1 & inv(dt2);
 724 }
 725
 726
 727 //- Inner Product of a SphericalTensorN and an inverse TensorN
 728 template <class Cmpt, int length>
 729 inline TensorN<Cmpt, length>
 730 operator/
 731 (
 732     const SphericalTensorN<Cmpt, length>& st1,
 733     const TensorN<Cmpt, length>& t2
 734 )
 735 {
 736     return st1.v_[0] * inv(t2);
 737 }
 738
 739
 740 //- Inner Product of a TensorN and an inverse SphericalTensorN
 741 template <class Cmpt, int length>
 742 inline TensorN<Cmpt, length>
 743 operator/
 744 (
 745     const TensorN<Cmpt, length>& t1,
 746     const SphericalTensorN<Cmpt, length>& st2
 747 )
 748 {
 749     TensorN<Cmpt, length> result;
 750
 751     const Cmpt& s = st2[0];
 752     for (label i = 0; i < TensorN<Cmpt, length>::nComponents; i++)
 753     {
 754         result.v_[i] = t1.v_[i]/s;
 755     }
 756
 757     return result;
 758 }
 759
 760
 761 template <class Cmpt, int length>
 762 inline Cmpt det(const TensorN<Cmpt, length>& t)
 763 {
 764     // Calculate determinant via sub-determinants
 765     Cmpt result = pTraits<Cmpt>::zero;
 766
 767     TensorN<Cmpt, length - 1> subMatrix;
 768
 769     for (label i = 0; i < length; i++)
 770     {
 771         label nj = 0;
 772
 773         // Build sub-matrix, skipping the
 774         for (label j = 0; j < length; j++)
 775         {
 776             // Skip i-th column
 777             if (j == i) continue;
 778
 779             for (label k = 1; k < length; k++)
 780             {
 781                 subMatrix(nj, k) = t(j, k);
 782             }
 783
 784             nj++;
 785         }
 786
 787         // Handle +/- sign switch
 788         result += pow(-1, i)*t(i, 0)*det(subMatrix);
 789     }
 790
 791     return result;
 792 }
 793
 794
 795 // Determinant: partial specialisation for rank 1
 796 template<class Cmpt>
 797 inline Cmpt det(const TensorN<Cmpt, 1>& t)
 798 {
 799     return t(0, 0);
 800 }
 801
 802
 803 // Determinant: partial specialisation for rank 2
 804 template<class Cmpt>
 805 inline Cmpt det(const TensorN<Cmpt, 2>& t)
 806 {
 807     return t(0, 0)*t(1, 1) - t(0, 1)*t(1, 0);
 808 }
 809
 810
 811 //- Return the inverse of a tensor give the determinant
 812 //  Uses Gauss-Jordan Elimination with full pivoting
 813 template <class Cmpt, int length>
 814 inline TensorN<Cmpt, length> inv(const TensorN<Cmpt, length>& t)
 815 {
 816     TensorN<Cmpt, length> result(t);
 817
 818     label iRow = 0, iCol = 0;
 819     Cmpt largestCoeff, temp;
 820     Cmpt* __restrict__ srcIter;
 821     Cmpt* __restrict__ destIter;
 822
 823     // Lists used for bookkeeping on the pivoting
 824     List<label> indexCol(length), indexRow(length), iPivot(length);
 825
 826     iPivot = 0;
 827
 828     label curRowOffset;
 829
 830     // Main loop over columns to be reduced
 831     for (label i = 0; i < length; i++)
 832     {
 833         largestCoeff = pTraits<Cmpt>::zero;
 834
 835         // Search for pivot element
 836         curRowOffset = 0;
 837
 838         for (label j = 0; j < length; j++)
 839         {
 840             if (iPivot[j] != 1)
 841             {
 842                 for (int k = 0; k < length; k++)
 843                 {
 844                     if (iPivot[k] == 0)
 845                     {
 846                         if
 847                         (
 848                             (temp = Foam::mag(result[curRowOffset+k]))
 849                          >= largestCoeff
 850                         )
 851                         {
 852                             largestCoeff = temp;
 853                             iRow = j;
 854                             iCol = k;
 855                         }
 856                     }
 857                 }
 858             }
 859             curRowOffset += length;
 860         }
 861         ++(iPivot[iCol]);
 862
 863         // We now have the pivot element
 864         // Interchange rows if needed
 865         if (iRow != iCol)
 866         {
 867             srcIter = &result(iRow, 0);
 868             destIter = &result(iCol, 0);
 869
 870             for (label j = 0; j < length; j++)
 871             {
 872                 Swap((*srcIter), (*destIter));
 873                 srcIter++;
 874                 destIter++;
 875             }
 876         }
 877         indexRow[i] = iRow;
 878         indexCol[i] = iCol;
 879
 880         //Check for singularity
 881         srcIter = &result(iCol, iCol);  // Dummy pointer to reduce indexing
 882         if ((*srcIter) == Cmpt(0))
 883         {
 884             FatalErrorIn
 885             (
 886                 "inline TensorN<Cmpt, length> inv("
 887                 "const TensorN<Cmpt, length>& t)"
 888             )   << "Singular tensor" << length << Foam::abort(FatalError);
 889         }
 890
 891         // Divide the pivot row by pivot element
 892         temp = pTraits<Cmpt>::one/(*srcIter);
 893         (*srcIter) = pTraits<Cmpt>::one;
 894
 895         srcIter = &result(iCol, 0);
 896         for (label j = 0; j < length; j++)
 897         {
 898             (*srcIter) *= temp;
 899             srcIter++;
 900         }
 901
 902         // Reduce the rows, excluding the pivot row
 903         for (label j = 0; j < length; j++)
 904         {
 905             if (j != iCol)
 906             {
 907                 destIter = &result(j, 0);
 908                 srcIter = &result(iCol, 0);
 909
 910                 temp=destIter[iCol];
 911                 destIter[iCol] = pTraits<Cmpt>::zero;
 912
 913                 for (label k = 0; k < length; k++)
 914                 {
 915                     (*destIter) -= (*srcIter)*temp;
 916                     srcIter++;
 917                     destIter++;
 918                 }
 919             }
 920         }
 921     }
 922
 923     // Unscamble the solution
 924     for (label i = length - 1; i >= 0; i--)
 925     {
 926         if (indexRow[i] != indexCol[i])
 927         {
 928             srcIter = &result[indexRow[i]];
 929             destIter = &result[indexCol[i]];
 930             for (label j = 0; j < length; j++)
 931             {
 932                 Swap((*srcIter), (*destIter));
 933                 srcIter += length;
 934                 destIter += length;
 935             }
 936         }
 937     }
 938
 939     return result;
 940 }
 941
 942
 943 //- Inverse: partial template specialisation for rank 1
 944 template <class Cmpt>
 945 inline TensorN<Cmpt, 1> inv(const TensorN<Cmpt, 1>& t)
 946 {
 947     return TensorN<Cmpt, 1>(1/t(0, 0));
 948 }
 949
 950
 951 //- Inverse: partial template specialisation for rank 2
 952 template <class Cmpt>
 953 inline TensorN<Cmpt, 2> inv(const TensorN<Cmpt, 2>& t)
 954 {
 955     TensorN<Cmpt, 2> result(t);
 956
 957     Cmpt oneOverDet = 1/(t(0, 0)*t(1, 1) - t(0, 1)*t(1, 0));
 958
 959     result(0, 0) = oneOverDet*t(1, 1);
 960     result(0, 1) = -oneOverDet*t(0, 1);
 961     result(1, 0) = -oneOverDet*t(1, 0);
 962     result(1, 1) = oneOverDet*t(0, 0);
 963
 964     return result;
 965 }
 966
 967
 968 //- Return tensor diagonal
 969 template <class Cmpt, int length>
 970 inline DiagTensorN<Cmpt, length> diag(const TensorN<Cmpt, length>& t)
 971 {
 972     return t.diag();
 973 }
 974
 975 //- Return tensor diagonal
 976 template <class Cmpt, int length>
 977 inline TensorN<Cmpt, length> negSumDiag(const TensorN<Cmpt, length>& t)
 978 {
 979     return t.negSumDiag();
 980 }
 981
 982 //- Return the component sum
 983 // template <class Cmpt, int length>
 984 // inline Cmpt sum(const TensorN<Cmpt, length>& t)
 985 // {
 986 //     Cmpt result=Cmpt::zero;
 987 //     for(int i=0; i<TensorN<Cmpt, length>::nComponents; i++)
 988 //     {
 989 //         result += t[i];
 990 //     }
 991 //     return result;
 992 // }
 993
 994 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 995
 996 template<class Cmpt, int length>
 997 class outerProduct<TensorN<Cmpt, length>, Cmpt>
 998 {
 999 public:
1000
1001     typedef TensorN<Cmpt, length> type;
1002 };
1003
1004 template<class Cmpt, int length>
1005 class outerProduct<Cmpt, TensorN<Cmpt, length> >
1006 {
1007 public:
1008
1009     typedef TensorN<Cmpt, length> type;
1010 };
1011
1012 template<class Cmpt, int length>
1013 class innerProduct<DiagTensorN<Cmpt, length>, TensorN<Cmpt, length> >
1014 {
1015 public:
1016
1017     typedef TensorN<Cmpt, length> type;
1018 };
1019
1020
1021 template<class Cmpt, int length>
1022 class innerProduct<TensorN<Cmpt, length>, DiagTensorN<Cmpt, length> >
1023 {
1024 public:
1025
1026     typedef TensorN<Cmpt, length> type;
1027 };
1028
1029
1030 template<class Cmpt, int length>
1031 class innerProduct<SphericalTensorN<Cmpt, length>, TensorN<Cmpt, length> >
1032 {
1033 public:
1034
1035     typedef TensorN<Cmpt, length> type;
1036 };
1037
1038
1039 template<class Cmpt, int length>
1040 class innerProduct<TensorN<Cmpt, length>, SphericalTensorN<Cmpt, length> >
1041 {
1042 public:
1043
1044     typedef TensorN<Cmpt, length> type;
1045 };
1046
1047
1048 template<class Cmpt, int length>
1049 class innerProduct<VectorN<Cmpt, length>, TensorN<Cmpt, length> >
1050 {
1051 public:
1052
1053     typedef VectorN<Cmpt, length> type;
1054 };
1055
1056
1057 template<class Cmpt, int length>
1058 class innerProduct<TensorN<Cmpt, length>, VectorN<Cmpt, length> >
1059 {
1060 public:
1061
1062     typedef VectorN<Cmpt, length> type;
1063 };
1064
1065
1066 template<class Cmpt, int length>
1067 class innerProduct<TensorN<Cmpt, length>, TensorN<Cmpt, length> >
1068 {
1069 public:
1070
1071     typedef TensorN<Cmpt, length> type;
1072 };
1073
1074
1075 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1076
1077 } // End namespace Foam
1078
1079 // ************************************************************************* //