src/VectorN/OpenFOAM/primitives/TensorNI.H

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