TensorI.H 27.2 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
OpenFOAM bot's avatar
OpenFOAM bot committed
6 7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2016 OpenFOAM Foundation
9
    Copyright (C) 2016-2020 OpenCFD Ltd.
10 11 12 13
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

14 15 16 17
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
18 19 20 21 22 23 24

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
25
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
26 27

\*---------------------------------------------------------------------------*/
28
#include <type_traits>
29

30 31
#include "SymmTensor.H"

32 33
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

34
template<class Cmpt>
35
inline Foam::Tensor<Cmpt>::Tensor(const Foam::zero)
36
:
37
    Tensor::msType(Zero)
38 39 40 41 42 43 44 45 46 47 48 49 50 51
{}


template<class Cmpt>
template<class Cmpt2>
inline Foam::Tensor<Cmpt>::Tensor
(
    const MatrixSpace<Tensor<Cmpt2>, Cmpt2, 3, 3>& vs
)
:
    Tensor::msType(vs)
{}


Henry's avatar
Henry committed
52 53
template<class Cmpt>
template<class Cmpt2>
54 55 56 57
inline Foam::Tensor<Cmpt>::Tensor
(
    const VectorSpace<Tensor<Cmpt2>, Cmpt2, 9>& vs
)
58
:
59
    Tensor::msType(vs)
60 61 62
{}


Henry's avatar
Henry committed
63
template<class Cmpt>
64
inline Foam::Tensor<Cmpt>::Tensor(const SphericalTensor<Cmpt>& st)
65
{
66 67 68
    this->v_[XX] = st.ii();    this->v_[XY] = Zero;      this->v_[XZ] = Zero;
    this->v_[YX] = Zero;       this->v_[YY] = st.ii();   this->v_[YZ] = Zero;
    this->v_[ZX] = Zero;       this->v_[ZY] = Zero;      this->v_[ZZ] = st.ii();
69 70 71
}


Henry's avatar
Henry committed
72
template<class Cmpt>
73
inline Foam::Tensor<Cmpt>::Tensor(const SymmTensor<Cmpt>& st)
74 75 76 77 78 79 80
{
    this->v_[XX] = st.xx(); this->v_[XY] = st.xy(); this->v_[XZ] = st.xz();
    this->v_[YX] = st.xy(); this->v_[YY] = st.yy(); this->v_[YZ] = st.yz();
    this->v_[ZX] = st.xz(); this->v_[ZY] = st.yz(); this->v_[ZZ] = st.zz();
}


Henry's avatar
Henry committed
81
template<class Cmpt>
82 83 84 85 86 87 88 89
inline Foam::Tensor<Cmpt>::Tensor
(
    const Vector<Vector<Cmpt>>& vecs,
    const bool transposed
)
:
    Tensor<Cmpt>(vecs.x(), vecs.y(), vecs.z(), transposed)
{}
90 91


Henry's avatar
Henry committed
92
template<class Cmpt>
93
inline Foam::Tensor<Cmpt>::Tensor
94 95 96
(
    const Vector<Cmpt>& x,
    const Vector<Cmpt>& y,
97 98
    const Vector<Cmpt>& z,
    const bool transposed
99 100
)
{
101 102 103 104 105 106 107 108
    if (transposed)
    {
        this->cols(x, y, z);
    }
    else
    {
        this->rows(x, y, z);
    }
109 110 111
}


Henry's avatar
Henry committed
112
template<class Cmpt>
113
inline Foam::Tensor<Cmpt>::Tensor
114 115 116 117 118 119 120 121 122 123 124 125
(
    const Cmpt txx, const Cmpt txy, const Cmpt txz,
    const Cmpt tyx, const Cmpt tyy, const Cmpt tyz,
    const Cmpt tzx, const Cmpt tzy, const Cmpt tzz
)
{
    this->v_[XX] = txx; this->v_[XY] = txy; this->v_[XZ] = txz;
    this->v_[YX] = tyx; this->v_[YY] = tyy; this->v_[YZ] = tyz;
    this->v_[ZX] = tzx; this->v_[ZY] = tzy; this->v_[ZZ] = tzz;
}


Henry's avatar
Henry committed
126
template<class Cmpt>
127 128 129 130 131 132 133 134 135 136
template
<
    template<class, Foam::direction, Foam::direction> class Block2,
    Foam::direction BRowStart,
    Foam::direction BColStart
>
inline Foam::Tensor<Cmpt>::Tensor
(
    const Block2<Tensor<Cmpt>, BRowStart, BColStart>& block
)
137
:
138
    Tensor::msType(block)
139 140 141
{}


Henry's avatar
Henry committed
142
template<class Cmpt>
143 144 145 146
inline Foam::Tensor<Cmpt>::Tensor(Istream& is)
:
    Tensor::msType(is)
{}
147 148


149
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
150

Henry's avatar
Henry committed
151
template<class Cmpt>
152
inline const Cmpt& Foam::Tensor<Cmpt>::xx() const
153 154 155 156
{
    return this->v_[XX];
}

andy's avatar
andy committed
157

Henry's avatar
Henry committed
158
template<class Cmpt>
159
inline const Cmpt& Foam::Tensor<Cmpt>::xy() const
160 161 162 163
{
    return this->v_[XY];
}

andy's avatar
andy committed
164

Henry's avatar
Henry committed
165
template<class Cmpt>
166
inline const Cmpt& Foam::Tensor<Cmpt>::xz() const
167 168 169 170 171
{
    return this->v_[XZ];
}


Henry's avatar
Henry committed
172
template<class Cmpt>
173
inline const Cmpt& Foam::Tensor<Cmpt>::yx() const
174 175 176 177
{
    return this->v_[YX];
}

andy's avatar
andy committed
178

Henry's avatar
Henry committed
179
template<class Cmpt>
180
inline const Cmpt& Foam::Tensor<Cmpt>::yy() const
181 182 183 184
{
    return this->v_[YY];
}

andy's avatar
andy committed
185

Henry's avatar
Henry committed
186
template<class Cmpt>
187
inline const Cmpt& Foam::Tensor<Cmpt>::yz() const
188 189 190 191 192
{
    return this->v_[YZ];
}


Henry's avatar
Henry committed
193
template<class Cmpt>
194
inline const Cmpt& Foam::Tensor<Cmpt>::zx() const
195 196 197 198
{
    return this->v_[ZX];
}

andy's avatar
andy committed
199

Henry's avatar
Henry committed
200
template<class Cmpt>
201
inline const Cmpt& Foam::Tensor<Cmpt>::zy() const
202 203 204 205
{
    return this->v_[ZY];
}

andy's avatar
andy committed
206

Henry's avatar
Henry committed
207
template<class Cmpt>
208
inline const Cmpt& Foam::Tensor<Cmpt>::zz() const
209 210 211 212 213
{
    return this->v_[ZZ];
}


Henry's avatar
Henry committed
214
template<class Cmpt>
215
inline Cmpt& Foam::Tensor<Cmpt>::xx()
216 217 218 219
{
    return this->v_[XX];
}

andy's avatar
andy committed
220

Henry's avatar
Henry committed
221
template<class Cmpt>
222
inline Cmpt& Foam::Tensor<Cmpt>::xy()
223 224 225 226
{
    return this->v_[XY];
}

andy's avatar
andy committed
227

Henry's avatar
Henry committed
228
template<class Cmpt>
229
inline Cmpt& Foam::Tensor<Cmpt>::xz()
230 231 232 233 234
{
    return this->v_[XZ];
}


Henry's avatar
Henry committed
235
template<class Cmpt>
236
inline Cmpt& Foam::Tensor<Cmpt>::yx()
237 238 239 240
{
    return this->v_[YX];
}

andy's avatar
andy committed
241

Henry's avatar
Henry committed
242
template<class Cmpt>
243
inline Cmpt& Foam::Tensor<Cmpt>::yy()
244 245 246 247
{
    return this->v_[YY];
}

andy's avatar
andy committed
248

Henry's avatar
Henry committed
249
template<class Cmpt>
250
inline Cmpt& Foam::Tensor<Cmpt>::yz()
251 252 253 254 255
{
    return this->v_[YZ];
}


Henry's avatar
Henry committed
256
template<class Cmpt>
257
inline Cmpt& Foam::Tensor<Cmpt>::zx()
258 259 260 261
{
    return this->v_[ZX];
}

andy's avatar
andy committed
262

Henry's avatar
Henry committed
263
template<class Cmpt>
264
inline Cmpt& Foam::Tensor<Cmpt>::zy()
265 266 267 268
{
    return this->v_[ZY];
}

andy's avatar
andy committed
269

Henry's avatar
Henry committed
270
template<class Cmpt>
271
inline Cmpt& Foam::Tensor<Cmpt>::zz()
272 273 274 275 276
{
    return this->v_[ZZ];
}


277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298
template<class Cmpt>
inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::x() const
{
    return Vector<Cmpt>(this->v_[XX], this->v_[XY], this->v_[XZ]);
}


template<class Cmpt>
inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::y() const
{
    return Vector<Cmpt>(this->v_[YX], this->v_[YY], this->v_[YZ]);
}


template<class Cmpt>
inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::z() const
{
    return Vector<Cmpt>(this->v_[ZX], this->v_[ZY], this->v_[ZZ]);
}


template<class Cmpt>
299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321
inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::cx() const
{
    return Vector<Cmpt>(this->v_[XX], this->v_[YX], this->v_[ZX]);
}


template<class Cmpt>
inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::cy() const
{
    return Vector<Cmpt>(this->v_[XY], this->v_[YY], this->v_[ZY]);
}


template<class Cmpt>
inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::cz() const
{
    return Vector<Cmpt>(this->v_[XZ], this->v_[YZ], this->v_[ZZ]);
}


template<class Cmpt>
template<Foam::direction Col>
inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::col() const
322
{
323 324 325 326 327 328 329 330 331 332 333 334 335
    if (Col == 0) return cx();
    else if (Col == 1) return cy();
    else if (Col == 2) return cz();

    static_assert(Col < 3, "Invalid column access");
    return Zero;
}


template<class Cmpt>
inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::col(const direction c) const
{
    switch (c)
336
    {
337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426
        case 0: return cx(); break;
        case 1: return cy(); break;
        case 2: return cz(); break;
    }

    FatalErrorInFunction
        << "Invalid column access " << c << abort(FatalError);

    return Zero;
}


template<class Cmpt>
template<Foam::direction Row>
inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::row() const
{
    if (Row == 0) return x();
    else if (Row == 1) return y();
    else if (Row == 2) return z();

    static_assert(Row < 3, "Invalid row access");
    return Zero;
}


template<class Cmpt>
inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::row(const direction r) const
{
    switch (r)
    {
        case 0: return x(); break;
        case 1: return y(); break;
        case 2: return z(); break;
    }

    FatalErrorInFunction
        << "Invalid row access " << r << abort(FatalError);

    return Zero;
}


template<class Cmpt>
template<Foam::direction Col>
inline void Foam::Tensor<Cmpt>::col(const Vector<Cmpt>& v)
{
    if (Col == 0)
    {
        this->v_[XX] = v.x();
        this->v_[YX] = v.y();
        this->v_[ZX] = v.z();
    }
    else if (Col == 1)
    {
        this->v_[XY] = v.x();
        this->v_[YY] = v.y();
        this->v_[ZY] = v.z();
    }
    else if (Col == 2)
    {
        this->v_[XZ] = v.x();
        this->v_[YZ] = v.y();
        this->v_[ZZ] = v.z();
    }

    static_assert(Col < 3, "Invalid column access");
}


template<class Cmpt>
template<Foam::direction Row>
inline void Foam::Tensor<Cmpt>::row(const Vector<Cmpt>& v)
{
    if (Row == 0)
    {
        this->v_[XX] = v.x(); this->v_[XY] = v.y(); this->v_[XZ] = v.z();
    }
    else if (Row == 1)
    {
        this->v_[YX] = v.x(); this->v_[YY] = v.y(); this->v_[YZ] = v.z();
    }
    else if (Row == 2)
    {
        this->v_[ZX] = v.x(); this->v_[ZY] = v.y(); this->v_[ZZ] = v.z();
    }

    static_assert(Row < 3, "Invalid row access");
}


427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454
template<class Cmpt>
inline void Foam::Tensor<Cmpt>::cols
(
    const Vector<Cmpt>& x,
    const Vector<Cmpt>& y,
    const Vector<Cmpt>& z
)
{
    this->v_[XX] = x.x(); this->v_[XY] = y.x(); this->v_[XZ] = z.x();
    this->v_[YX] = x.y(); this->v_[YY] = y.y(); this->v_[YZ] = z.y();
    this->v_[ZX] = x.z(); this->v_[ZY] = y.z(); this->v_[ZZ] = z.z();
}


template<class Cmpt>
inline void Foam::Tensor<Cmpt>::rows
(
    const Vector<Cmpt>& x,
    const Vector<Cmpt>& y,
    const Vector<Cmpt>& z
)
{
    this->v_[XX] = x.x(); this->v_[XY] = x.y(); this->v_[XZ] = x.z();
    this->v_[YX] = y.x(); this->v_[YY] = y.y(); this->v_[YZ] = y.z();
    this->v_[ZX] = z.x(); this->v_[ZY] = z.y(); this->v_[ZZ] = z.z();
}


455 456 457 458 459 460 461 462 463 464 465
template<class Cmpt>
inline void Foam::Tensor<Cmpt>::col(const direction c, const Vector<Cmpt>& v)
{
    switch (c)
    {
        case 0: col<0>(v); break;
        case 1: col<1>(v); break;
        case 2: col<2>(v); break;
        default:
            FatalErrorInFunction
                << "Invalid column access " << c << abort(FatalError);
466
            break;
467 468 469 470 471 472 473 474 475 476 477 478 479 480 481
    }
}


template<class Cmpt>
inline void Foam::Tensor<Cmpt>::row(const direction r, const Vector<Cmpt>& v)
{
    switch (r)
    {
        case 0: row<0>(v); break;
        case 1: row<1>(v); break;
        case 2: row<2>(v); break;
        default:
            FatalErrorInFunction
                << "Invalid row access " << r << abort(FatalError);
482 483
            break;
    }
484
}
485 486


487
template<class Cmpt>
488 489 490 491 492 493 494 495
inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::diag() const
{
    return Vector<Cmpt>(this->v_[XX], this->v_[YY], this->v_[ZZ]);
}


template<class Cmpt>
inline void Foam::Tensor<Cmpt>::diag(const Vector<Cmpt>& v)
496
{
497
    this->v_[XX] = v.x(); this->v_[YY] = v.y(); this->v_[ZZ] = v.z();
498 499 500
}


501 502
// * * * * * * * * * * * * * * * Member Operations * * * * * * * * * * * * * //

Henry's avatar
Henry committed
503
template<class Cmpt>
504
inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::T() const
505 506 507 508 509 510 511 512 513 514
{
    return Tensor<Cmpt>
    (
        xx(), yx(), zx(),
        xy(), yy(), zy(),
        xz(), yz(), zz()
    );
}


515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536
template<class Cmpt>
inline Foam::Tensor<Cmpt>
Foam::Tensor<Cmpt>::inner(const Tensor<Cmpt>& t2) const
{
    const Tensor<Cmpt>& t1 = *this;

    return Tensor<Cmpt>
    (
        t1.xx()*t2.xx() + t1.xy()*t2.yx() + t1.xz()*t2.zx(),
        t1.xx()*t2.xy() + t1.xy()*t2.yy() + t1.xz()*t2.zy(),
        t1.xx()*t2.xz() + t1.xy()*t2.yz() + t1.xz()*t2.zz(),

        t1.yx()*t2.xx() + t1.yy()*t2.yx() + t1.yz()*t2.zx(),
        t1.yx()*t2.xy() + t1.yy()*t2.yy() + t1.yz()*t2.zy(),
        t1.yx()*t2.xz() + t1.yy()*t2.yz() + t1.yz()*t2.zz(),

        t1.zx()*t2.xx() + t1.zy()*t2.yx() + t1.zz()*t2.zx(),
        t1.zx()*t2.xy() + t1.zy()*t2.yy() + t1.zz()*t2.zy(),
        t1.zx()*t2.xz() + t1.zy()*t2.yz() + t1.zz()*t2.zz()
    );
}

537

538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561
template<class Cmpt>
inline Foam::Tensor<Cmpt>
Foam::Tensor<Cmpt>::schur(const Tensor<Cmpt>& t2) const
{
    const Tensor<Cmpt>& t1 = *this;

    return Tensor<Cmpt>
    (
        t1.xx()*t2.xx(), t1.xy()*t2.xy(), t1.xz()*t2.xz(),
        t1.yx()*t2.yx(), t1.yy()*t2.yy(), t1.yz()*t2.yz(),
        t1.zx()*t2.zx(), t1.zy()*t2.zy(), t1.zz()*t2.zz()
    );
}


// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //

template<class Cmpt>
inline void Foam::Tensor<Cmpt>::operator&=(const Tensor<Cmpt>& t)
{
    *this = this->inner(t);
}


562 563 564 565 566 567 568 569 570 571 572
template<class Cmpt>
template<class Cmpt2>
inline void Foam::Tensor<Cmpt>::operator=
(
    const VectorSpace<Tensor<Cmpt2>, Cmpt2, 9>& vs
)
{
    VectorSpace<Tensor<Cmpt>, Cmpt, 9>::operator=(vs);
}


Henry's avatar
Henry committed
573
template<class Cmpt>
574
inline void Foam::Tensor<Cmpt>::operator=(const SphericalTensor<Cmpt>& st)
575
{
576 577 578
    this->v_[XX] = st.ii();  this->v_[XY] = Zero;     this->v_[XZ] = Zero;
    this->v_[YX] = Zero;     this->v_[YY] = st.ii();  this->v_[YZ] = Zero;
    this->v_[ZX] = Zero;     this->v_[ZY] = Zero;     this->v_[ZZ] = st.ii();
579 580 581
}


Henry's avatar
Henry committed
582
template<class Cmpt>
583
inline void Foam::Tensor<Cmpt>::operator=(const SymmTensor<Cmpt>& st)
584 585 586 587 588 589 590
{
    this->v_[XX] = st.xx(); this->v_[XY] = st.xy(); this->v_[XZ] = st.xz();
    this->v_[YX] = st.xy(); this->v_[YY] = st.yy(); this->v_[YZ] = st.yz();
    this->v_[ZX] = st.xz(); this->v_[ZY] = st.yz(); this->v_[ZZ] = st.zz();
}


Henry's avatar
Henry committed
591
template<class Cmpt>
592
inline void Foam::Tensor<Cmpt>::operator=(const Vector<Vector<Cmpt>>& tr)
593 594 595 596 597 598 599 600 601 602 603 604 605 606 607
{
    this->v_[XX] = tr.x().x();
    this->v_[XY] = tr.x().y();
    this->v_[XZ] = tr.x().z();

    this->v_[YX] = tr.y().x();
    this->v_[YY] = tr.y().y();
    this->v_[YZ] = tr.y().z();

    this->v_[ZX] = tr.z().x();
    this->v_[ZY] = tr.z().y();
    this->v_[ZZ] = tr.z().z();
}


608 609 610 611 612
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

613 614
// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //

615
//- Return the trace of a Tensor
Henry's avatar
Henry committed
616
template<class Cmpt>
617 618 619 620 621 622
inline Cmpt tr(const Tensor<Cmpt>& t)
{
    return t.xx() + t.yy() + t.zz();
}


623
//- Return the spherical part of a Tensor
Henry's avatar
Henry committed
624
template<class Cmpt>
625 626
inline SphericalTensor<Cmpt> sph(const Tensor<Cmpt>& t)
{
627 628 629 630
    return SphericalTensor<Cmpt>
    (
        (1.0/3.0)*tr(t)
    );
631 632 633
}


634
//- Return the symmetric part of a Tensor
Henry's avatar
Henry committed
635
template<class Cmpt>
636 637 638 639 640 641 642 643 644 645 646
inline SymmTensor<Cmpt> symm(const Tensor<Cmpt>& t)
{
    return SymmTensor<Cmpt>
    (
        t.xx(), 0.5*(t.xy() + t.yx()), 0.5*(t.xz() + t.zx()),
                t.yy(),                0.5*(t.yz() + t.zy()),
                                       t.zz()
    );
}


647
//- Return twice the symmetric part of a Tensor
Henry's avatar
Henry committed
648
template<class Cmpt>
649 650 651 652 653 654 655 656 657 658 659
inline SymmTensor<Cmpt> twoSymm(const Tensor<Cmpt>& t)
{
    return SymmTensor<Cmpt>
    (
        2*t.xx(), (t.xy() + t.yx()), (t.xz() + t.zx()),
                  2*t.yy(),          (t.yz() + t.zy()),
                                     2*t.zz()
    );
}


660
//- Return the skew-symmetric part of a Tensor
Henry's avatar
Henry committed
661
template<class Cmpt>
662 663 664 665
inline Tensor<Cmpt> skew(const Tensor<Cmpt>& t)
{
    return Tensor<Cmpt>
    (
666 667 668
        Zero,                  0.5*(t.xy() - t.yx()), 0.5*(t.xz() - t.zx()),
        0.5*(t.yx() - t.xy()), Zero,                  0.5*(t.yz() - t.zy()),
        0.5*(t.zx() - t.xz()), 0.5*(t.zy() - t.yz()), Zero
669 670 671 672
    );
}


673
//- Return the skew-symmetric part of a SymmTensor as a Tensor
Henry's avatar
Henry committed
674
template<class Cmpt>
675 676 677 678 679 680
inline const Tensor<Cmpt>& skew(const SymmTensor<Cmpt>& st)
{
    return Tensor<Cmpt>::zero;
}


681
//- Return the deviatoric part of a Tensor
Henry's avatar
Henry committed
682
template<class Cmpt>
683 684
inline Tensor<Cmpt> dev(const Tensor<Cmpt>& t)
{
685
    return t - sph(t);
686 687 688
}


689
//- Return the two-third deviatoric part of a Tensor
Henry's avatar
Henry committed
690
template<class Cmpt>
691 692
inline Tensor<Cmpt> dev2(const Tensor<Cmpt>& t)
{
693
    return t - 2*sph(t);
694 695 696
}


697
//- Return the determinant of a Tensor
Henry's avatar
Henry committed
698
template<class Cmpt>
699 700 701 702 703 704 705 706 707 708 709
inline Cmpt det(const Tensor<Cmpt>& t)
{
    return
    (
        t.xx()*t.yy()*t.zz() + t.xy()*t.yz()*t.zx()
      + t.xz()*t.yx()*t.zy() - t.xx()*t.yz()*t.zy()
      - t.xy()*t.yx()*t.zz() - t.xz()*t.yy()*t.zx()
    );
}


710
//- Return the cofactor Tensor of a Tensor
Henry's avatar
Henry committed
711
template<class Cmpt>
712
inline Tensor<Cmpt> cof(const Tensor<Cmpt>& t)
713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730
{
    return Tensor<Cmpt>
    (
        t.yy()*t.zz() - t.zy()*t.yz(),
        t.zx()*t.yz() - t.yx()*t.zz(),
        t.yx()*t.zy() - t.yy()*t.zx(),

        t.xz()*t.zy() - t.xy()*t.zz(),
        t.xx()*t.zz() - t.xz()*t.zx(),
        t.xy()*t.zx() - t.xx()*t.zy(),

        t.xy()*t.yz() - t.xz()*t.yy(),
        t.yx()*t.xz() - t.xx()*t.yz(),
        t.xx()*t.yy() - t.yx()*t.xy()
    );
}


731
//- Return the inverse of a Tensor by using the given determinant
Henry's avatar
Henry committed
732
template<class Cmpt>
733 734
inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t, const Cmpt dett)
{
735 736 737 738 739 740 741 742 743
    #ifdef FULLDEBUG
    if (mag(dett) < SMALL)
    {
        FatalErrorInFunction
            << "Tensor is not invertible due to the zero determinant:"
            << "det(Tensor) = " << mag(dett)
            << abort(FatalError);
    }
    #endif
744

745
    return cof(t).T()/dett;
746 747 748
}


749
//- Return the inverse of a Tensor
Henry's avatar
Henry committed
750
template<class Cmpt>
751 752 753 754 755 756
inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t)
{
    return inv(t, det(t));
}


757
//- Return the inverse of this Tensor
758 759 760 761 762 763 764
template<class Cmpt>
inline Tensor<Cmpt> Tensor<Cmpt>::inv() const
{
    return Foam::inv(*this);
}


765
//- Return the 1st invariant of a Tensor
Henry's avatar
Henry committed
766
template<class Cmpt>
767 768 769 770 771 772
inline Cmpt invariantI(const Tensor<Cmpt>& t)
{
    return tr(t);
}


773
//- Return the 2nd invariant of a Tensor
Henry's avatar
Henry committed
774
template<class Cmpt>
775 776 777 778
inline Cmpt invariantII(const Tensor<Cmpt>& t)
{
    return
    (
779 780
        t.xx()*t.yy() + t.yy()*t.zz() + t.xx()*t.zz()
      - t.xy()*t.yx() - t.yz()*t.zy() - t.xz()*t.zx()
781 782 783 784
    );
}


785
//- Return the 3rd invariant of a Tensor
Henry's avatar
Henry committed
786
template<class Cmpt>
787 788 789 790 791 792
inline Cmpt invariantIII(const Tensor<Cmpt>& t)
{
    return det(t);
}


793
// * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
794

795
//- Sum of a SphericalTensor and a Tensor
Henry's avatar
Henry committed
796
template<class Cmpt>
797 798 799 800 801 802 803 804 805 806 807 808
inline Tensor<Cmpt>
operator+(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
{
    return Tensor<Cmpt>
    (
        st1.ii() + t2.xx(), t2.xy(),            t2.xz(),
        t2.yx(),            st1.ii() + t2.yy(), t2.yz(),
        t2.zx(),            t2.zy(),            st1.ii() + t2.zz()
    );
}


809
//- Sum of a Tensor and a SphericalTensor
Henry's avatar
Henry committed
810
template<class Cmpt>
811 812 813 814 815 816 817 818 819 820 821 822
inline Tensor<Cmpt>
operator+(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
{
    return Tensor<Cmpt>
    (
        t1.xx() + st2.ii(), t1.xy(),            t1.xz(),
        t1.yx(),            t1.yy() + st2.ii(), t1.yz(),
        t1.zx(),            t1.zy(),            t1.zz() + st2.ii()
    );
}


823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850
//- Sum of a SymmTensor and a Tensor
template<class Cmpt>
inline Tensor<Cmpt>
operator+(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
{
    return Tensor<Cmpt>
    (
        st1.xx() + t2.xx(), st1.xy() + t2.xy(), st1.xz() + t2.xz(),
        st1.xy() + t2.yx(), st1.yy() + t2.yy(), st1.yz() + t2.yz(),
        st1.xz() + t2.zx(), st1.yz() + t2.zy(), st1.zz() + t2.zz()
    );
}


//- Sum of a Tensor and a SymmTensor
template<class Cmpt>
inline Tensor<Cmpt>
operator+(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
{
    return Tensor<Cmpt>
    (
        t1.xx() + st2.xx(), t1.xy() + st2.xy(), t1.xz() + st2.xz(),
        t1.yx() + st2.xy(), t1.yy() + st2.yy(), t1.yz() + st2.yz(),
        t1.zx() + st2.xz(), t1.zy() + st2.yz(), t1.zz() + st2.zz()
    );
}


851
//- Subtract a Tensor from a SphericalTensor
Henry's avatar
Henry committed
852
template<class Cmpt>
853 854 855 856 857 858 859 860 861 862 863 864
inline Tensor<Cmpt>
operator-(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
{
    return Tensor<Cmpt>
    (
        st1.ii() - t2.xx(), -t2.xy(),            -t2.xz(),
       -t2.yx(),             st1.ii() - t2.yy(), -t2.yz(),
       -t2.zx(),            -t2.zy(),             st1.ii() - t2.zz()
    );
}


865
//- Subtract a SphericalTensor from a Tensor
Henry's avatar
Henry committed
866
template<class Cmpt>
867 868 869 870 871 872 873 874 875 876 877 878
inline Tensor<Cmpt>
operator-(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
{
    return Tensor<Cmpt>
    (
        t1.xx() - st2.ii(), t1.xy(),            t1.xz(),
        t1.yx(),            t1.yy() - st2.ii(), t1.yz(),
        t1.zx(),            t1.zy(),            t1.zz() - st2.ii()
    );
}


879
//- Subtract a Tensor from a SymmTensor
Henry's avatar
Henry committed
880
template<class Cmpt>
881
inline Tensor<Cmpt>
882
operator-(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
883 884 885
{
    return Tensor<Cmpt>
    (
886 887 888
        st1.xx() - t2.xx(), st1.xy() - t2.xy(), st1.xz() - t2.xz(),
        st1.xy() - t2.yx(), st1.yy() - t2.yy(), st1.yz() - t2.yz(),
        st1.xz() - t2.zx(), st1.yz() - t2.zy(), st1.zz() - t2.zz()
889 890 891 892
    );
}


893
//- Subtract a SymmTensor from a Tensor
Henry's avatar
Henry committed
894
template<class Cmpt>
895
inline Tensor<Cmpt>
896
operator-(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
897 898 899
{
    return Tensor<Cmpt>
    (
900 901 902
        t1.xx() - st2.xx(), t1.xy() - st2.xy(), t1.xz() - st2.xz(),
        t1.yx() - st2.xy(), t1.yy() - st2.yy(), t1.yz() - st2.yz(),
        t1.zx() - st2.xz(), t1.zy() - st2.yz(), t1.zz() - st2.zz()
903 904 905 906
    );
}


907
//- Return the Hodge dual of a Tensor as a Vector
Henry's avatar
Henry committed
908
template<class Cmpt>
909
inline Vector<Cmpt> operator*(const Tensor<Cmpt>& t)
910
{
911
    return Vector<Cmpt>(t.yz(), -t.xz(), t.xy());
912 913 914
}


915
//- Return the Hodge dual of a Vector as a Tensor
Henry's avatar
Henry committed
916
template<class Cmpt>
917
inline Tensor<Cmpt> operator*(const Vector<Cmpt>& v)
918
{
919 920 921 922 923 924
    return Tensor<Cmpt>
    (
         Zero,  -v.z(),  v.y(),
         v.z(),  Zero,  -v.x(),
        -v.y(),  v.x(),  Zero
    );
925 926
}

andy's avatar
andy committed
927

928
//- Division of a Vector by a Tensor
929
template<class Cmpt>
930 931
inline typename innerProduct<Vector<Cmpt>, Tensor<Cmpt>>::type
operator/(const Vector<Cmpt>& v, const Tensor<Cmpt>& t)
932
{
933 934
    return inv(t) & v;
}
935 936


937
//- Division of a Tensor by a Cmpt
Henry's avatar
Henry committed
938
template<class Cmpt>
939
inline Tensor<Cmpt>
940
operator/(const Tensor<Cmpt>& t, const Cmpt s)
941
{
942 943 944 945 946 947 948 949 950 951 952
    #ifdef FULLDEBUG
    if (mag(s) < VSMALL)
    {
        FatalErrorInFunction
            << "Tensor = " << t
            << " is not divisible due to a zero value in Cmpt:"
            << "Cmpt = " << s
            << abort(FatalError);
    }
    #endif

953 954
    return Tensor<Cmpt>
    (
955 956 957
        t.xx()/s, t.xy()/s, t.xz()/s,
        t.yx()/s, t.yy()/s, t.yz()/s,
        t.zx()/s, t.zy()/s, t.zz()/s
958 959 960 961
    );
}


962
//- Inner-product of a Tensor and a Tensor
Henry's avatar
Henry committed
963
template<class Cmpt>
964 965
inline typename innerProduct<Tensor<Cmpt>, Tensor<Cmpt>>::type
operator&(const Tensor<Cmpt>& t1, const Tensor<Cmpt>& t2)
966
{
967
    return t1.inner(t2);
968 969 970
}


971
//- Inner-product of a SphericalTensor and a Tensor
Henry's avatar
Henry committed
972
template<class Cmpt>
973
inline Tensor<Cmpt>
974
operator&(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
975 976 977
{
    return Tensor<Cmpt>
    (
978 979 980
        st1.ii()*t2.xx(), st1.ii()*t2.xy(), st1.ii()*t2.xz(),
        st1.ii()*t2.yx(), st1.ii()*t2.yy(), st1.ii()*t2.yz(),
        st1.ii()*t2.zx(), st1.ii()*t2.zy(), st1.ii()*t2.zz()
981 982 983 984
    );
}


985
//- Inner-product of a Tensor and a SphericalTensor
Henry's avatar
Henry committed
986
template<class Cmpt>
987
inline Tensor<Cmpt>
988
operator&(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
989 990 991
{
    return Tensor<Cmpt>
    (
992 993 994
        t1.xx()*st2.ii(), t1.xy()*st2.ii(), t1.xz()*st2.ii(),
        t1.yx()*st2.ii(), t1.yy()*st2.ii(), t1.yz()*st2.ii(),
        t1.zx()*st2.ii(), t1.zy()*st2.ii(), t1.zz()*st2.ii()
995 996 997 998
    );
}


999
//- Inner-product of a SymmTensor and a Tensor
Henry's avatar
Henry committed
1000
template<class Cmpt>
1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020
inline Tensor<Cmpt>
operator&(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
{
    return Tensor<Cmpt>
    (
        st1.xx()*t2.xx() + st1.xy()*t2.yx() + st1.xz()*t2.zx(),
        st1.xx()*t2.xy() + st1.xy()*t2.yy() + st1.xz()*t2.zy(),
        st1.xx()*t2.xz() + st1.xy()*t2.yz() + st1.xz()*t2.zz(),

        st1.xy()*t2.xx() + st1.yy()*t2.yx() + st1.yz()*t2.zx(),
        st1.xy()*t2.xy() + st1.yy()*t2.yy() + st1.yz()*t2.zy(),
        st1.xy()*t2.xz() + st1.yy()*t2.yz() + st1.yz()*t2.zz(),

        st1.xz()*t2.xx() + st1.yz()*t2.yx() + st1.zz()*t2.zx(),
        st1.xz()*t2.xy() + st1.yz()*t2.yy() + st1.zz()*t2.zy(),
        st1.xz()*t2.xz() + st1.yz()*t2.yz() + st1.zz()*t2.zz()
    );
}


1021
//- Inner-product of a Tensor and a SymmTensor
Henry's avatar
Henry committed
1022
template<class Cmpt>
1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042
inline Tensor<Cmpt>
operator&(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
{
    return Tensor<Cmpt>
    (
        t1.xx()*st2.xx() + t1.xy()*st2.xy() + t1.xz()*st2.xz(),
        t1.xx()*st2.xy() + t1.xy()*st2.yy() + t1.xz()*st2.yz(),
        t1.xx()*st2.xz() + t1.xy()*st2.yz() + t1.xz()*st2.zz(),

        t1.yx()*st2.xx() + t1.yy()*st2.xy() + t1.yz()*st2.xz(),
        t1.yx()*st2.xy() + t1.yy()*st2.yy() + t1.yz()*st2.yz(),
        t1.yx()*st2.xz() + t1.yy()*st2.yz() + t1.yz()*st2.zz(),

        t1.zx()*st2.xx() + t1.zy()*st2.xy() + t1.zz()*st2.xz(),
        t1.zx()*st2.xy() + t1.zy()*st2.yy() + t1.zz()*st2.yz(),
        t1.zx()*st2.xz() + t1.zy()*st2.yz() + t1.zz()*st2.zz()
    );
}


1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088
//- Inner-product of a Tensor and a Vector
template<class Cmpt>
inline typename innerProduct<Tensor<Cmpt>, Vector<Cmpt>>::type
operator&(const Tensor<Cmpt>& t, const Vector<Cmpt>& v)
{
    return Vector<Cmpt>
    (
        t.xx()*v.x() + t.xy()*v.y() + t.xz()*v.z(),
        t.yx()*v.x() + t.yy()*v.y() + t.yz()*v.z(),
        t.zx()*v.x() + t.zy()*v.y() + t.zz()*v.z()
    );
}


//- Inner-product of a Vector and a Tensor
template<class Cmpt>
inline typename innerProduct<Vector<Cmpt>, Tensor<Cmpt>>::type
operator&(const Vector<Cmpt>& v, const Tensor<Cmpt>& t)
{
    return Vector<Cmpt>
    (
        v.x()*t.xx() + v.y()*t.yx() + v.z()*t.zx(),
        v.x()*t.xy() + v.y()*t.yy() + v.z()*t.zy(),
        v.x()*t.xz() + v.y()*t.yz() + v.z()*t.zz()
    );
}


//- Double-inner-product of a SphericalTensor and a Tensor
template<class Cmpt>
inline Cmpt
operator&&(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
{
    return (st1.ii()*t2.xx() + st1.ii()*t2.yy() + st1.ii()*t2.zz());
}


//- Double-inner-product of a Tensor and a SphericalTensor
template<class Cmpt>
inline Cmpt
operator&&(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
{
    return (t1.xx()*st2.ii() + t1.yy()*st2.ii() + t1.zz()*st2.ii());
}


1089
//- Double-inner-product of a SymmTensor and a Tensor
Henry's avatar
Henry committed
1090
template<class Cmpt>
1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102
inline Cmpt
operator&&(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
{
    return
    (
        st1.xx()*t2.xx() + st1.xy()*t2.xy() + st1.xz()*t2.xz() +
        st1.xy()*t2.yx() + st1.yy()*t2.yy() + st1.yz()*t2.yz() +
        st1.xz()*t2.zx() + st1.yz()*t2.zy() + st1.zz()*t2.zz()
    );
}


1103
//- Double-inner-produ