diff --git a/applications/test/SymmTensor/Test-SymmTensor.C b/applications/test/SymmTensor/Test-SymmTensor.C index d1ea6a0b1d..0d73f4d0d7 100644 --- a/applications/test/SymmTensor/Test-SymmTensor.C +++ b/applications/test/SymmTensor/Test-SymmTensor.C @@ -370,7 +370,7 @@ void test_global_funcs(Type) Type(126) ) ); - cmp(" Square of Frobenius norm = ", magSqr(sT), Type(205)); + cmp(" Square of Frobenius norm = ", magSqr(sT), scalar(205)); } diff --git a/applications/test/SymmTensor2D/Test-SymmTensor2D.C b/applications/test/SymmTensor2D/Test-SymmTensor2D.C index da9b102b5f..315e9e808e 100644 --- a/applications/test/SymmTensor2D/Test-SymmTensor2D.C +++ b/applications/test/SymmTensor2D/Test-SymmTensor2D.C @@ -325,7 +325,7 @@ void test_global_funcs(Type) Type(13) ) ); - cmp(" Square of Frobenius norm = ", magSqr(sT), Type(17.999999999999996)); + cmp(" Square of Frobenius norm = ", magSqr(sT), scalar(17.999999999999996)); cmp ( " Outer-product of a Vector2D with itself = ", diff --git a/applications/test/complex/Test-complex.C b/applications/test/complex/Test-complex.C index e8940e07a7..f4875f42ec 100644 --- a/applications/test/complex/Test-complex.C +++ b/applications/test/complex/Test-complex.C @@ -34,6 +34,9 @@ Description #include "complex.H" #include "complexFields.H" #include "scalarField.H" +#include "diagTensor.H" +#include "symmTensor.H" +#include "symmTensor2D.H" #include "ListOps.H" #include "ops.H" @@ -164,6 +167,72 @@ int main(int argc, char *argv[]) << " => " << imags << nl; } + { + SymmTensor st1(SymmTensor::uniform({3, 4})); + + Info<< "symmTensor: " << st1 << nl + << " tr: " << tr(st1) << nl + << " diagSqr: " << st1.diagSqr() << nl + << " magSqr: " << magSqr(st1) << nl + << " mag: " << mag(st1) << nl; + + SymmTensor st2(SymmTensor::uniform(5)); + + Info<< "symmTensor: " << st2 << nl + << " tr: " << tr(st2) << nl + << " diagSqr: " << st2.diagSqr() << nl + << " magSqr: " << magSqr(st2) << nl + << " mag: " << mag(st2) << nl; + + st2 = Zero; + + DiagTensor dt1(SphericalTensor({3, 4})); + + Info<< "diagTensor: " << dt1 << nl + << " tr: " << tr(dt1) << nl + << " diagSqr: " << dt1.diagSqr() << nl + << " magSqr: " << magSqr(dt1) << nl + << " mag: " << mag(dt1) << nl; + + + // A bit ugly... + st1 = SphericalTensor({3, 4}); + + Info<< "symmTensor: " << st1 << nl + << " tr: " << tr(st1) << nl + << " diagSqr: " << st1.diagSqr() << nl + << " magSqr: " << magSqr(st1) << nl + << " mag: " << mag(st1) << nl; + } + + { + SymmTensor2D st1(SymmTensor2D::uniform({3, 4})); + + Info<< "symmTensor: " << st1 << nl + << " tr: " << tr(st1) << nl + << " diagSqr: " << st1.diagSqr() << nl + << " magSqr: " << magSqr(st1) << nl + << " mag: " << mag(st1) << nl; + } + + { + Tensor st1(Tensor::uniform({3, 4})); + + Info<< "tensor: " << st1 << nl + << " tr: " << tr(st1) << nl + << " diagSqr: " << st1.diagSqr() << nl + << " magSqr: " << magSqr(st1) << nl + << " mag: " << mag(st1) << endl; + + Tensor st2(Tensor::uniform(5)); + + Info<< "Tensor: " << st2 << nl + << " tr: " << tr(st2) << nl + << " diagSqr: " << st2.diagSqr() << nl + << " magSqr: " << magSqr(st2) << nl + << " mag: " << mag(st2) << endl; + } + complexField fld1(3, complex(2.0, 1.0)); complexField fld2(fld1); diff --git a/src/OpenFOAM/primitives/DiagTensor/DiagTensor.H b/src/OpenFOAM/primitives/DiagTensor/DiagTensor.H index 2f5817037a..c4d141d667 100644 --- a/src/OpenFOAM/primitives/DiagTensor/DiagTensor.H +++ b/src/OpenFOAM/primitives/DiagTensor/DiagTensor.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019-2022 OpenCFD Ltd. + Copyright (C) 2019-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -43,6 +43,7 @@ SourceFiles #define Foam_DiagTensor_H #include "Tensor.H" +#include "SphericalTensor.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -97,6 +98,9 @@ public: template inline DiagTensor(const VectorSpace, Cmpt2, 3>&); + //- Construct given SphericalTensor + inline DiagTensor(const SphericalTensor&); + //- Construct given three components inline DiagTensor(const Cmpt& txx, const Cmpt& tyy, const Cmpt& tzz); @@ -115,6 +119,12 @@ public: Cmpt& xx() noexcept { return this->v_[XX]; } Cmpt& yy() noexcept { return this->v_[YY]; } Cmpt& zz() noexcept { return this->v_[ZZ]; } + + + // Diagonal access and manipulation + + //- The L2-norm squared of the diagonal + inline scalar diagSqr() const; }; diff --git a/src/OpenFOAM/primitives/DiagTensor/DiagTensorI.H b/src/OpenFOAM/primitives/DiagTensor/DiagTensorI.H index 2432ac4d06..e898e8181d 100644 --- a/src/OpenFOAM/primitives/DiagTensor/DiagTensorI.H +++ b/src/OpenFOAM/primitives/DiagTensor/DiagTensorI.H @@ -49,6 +49,15 @@ inline Foam::DiagTensor::DiagTensor {} +template +inline Foam::DiagTensor::DiagTensor(const SphericalTensor& st) +{ + this->v_[XX] = st.ii(); + this->v_[YY] = st.ii(); + this->v_[ZZ] = st.ii(); +} + + template inline Foam::DiagTensor::DiagTensor ( @@ -70,6 +79,20 @@ inline Foam::DiagTensor::DiagTensor(Istream& is) {} +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +inline Foam::scalar Foam::DiagTensor::diagSqr() const +{ + return + ( + Foam::magSqr(this->xx()) + + Foam::magSqr(this->yy()) + + Foam::magSqr(this->zz()) + ); +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -138,6 +161,13 @@ inline DiagTensor diag(const SymmTensor& st) } +template +inline Foam::scalar magSqr(const DiagTensor& t) +{ + return t.diagSqr(); +} + + //- Linear interpolation of diagonal tensors a and b by factor t template inline DiagTensor lerp diff --git a/src/OpenFOAM/primitives/SphericalTensor/SphericalTensor.H b/src/OpenFOAM/primitives/SphericalTensor/SphericalTensor.H index 768f239e8d..0dba5685bd 100644 --- a/src/OpenFOAM/primitives/SphericalTensor/SphericalTensor.H +++ b/src/OpenFOAM/primitives/SphericalTensor/SphericalTensor.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019-2022 OpenCFD Ltd. + Copyright (C) 2019-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -124,10 +124,16 @@ public: Cmpt& ii() noexcept { return this->v_[II]; } + // Diagonal access and manipulation + + //- The L2-norm squared of the diagonal + inline scalar diagSqr() const; + + // Tensor Operations //- Return non-Hermitian transpose (no-op) - inline const SphericalTensor& T() const; + const SphericalTensor& T() const noexcept { return *this; } }; diff --git a/src/OpenFOAM/primitives/SphericalTensor/SphericalTensorI.H b/src/OpenFOAM/primitives/SphericalTensor/SphericalTensorI.H index f7d2f5b387..3f0ed117d2 100644 --- a/src/OpenFOAM/primitives/SphericalTensor/SphericalTensorI.H +++ b/src/OpenFOAM/primitives/SphericalTensor/SphericalTensorI.H @@ -65,10 +65,9 @@ inline Foam::SphericalTensor::SphericalTensor(Istream& is) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -inline const Foam::SphericalTensor& -Foam::SphericalTensor::T() const +inline Foam::scalar Foam::SphericalTensor::diagSqr() const { - return *this; + return 3*Foam::magSqr(this->ii()); } @@ -111,8 +110,8 @@ inline SphericalTensor inv(const SphericalTensor& st) if (mag(st.ii()) < VSMALL) { FatalErrorInFunction - << "SphericalTensor is not invertible due to the zero determinant:" - << "det(SphericalTensor) = " << det(st) + << "SphericalTensor not invertible, determinant:" + << det(st) << " tensor:" << st << nl << abort(FatalError); } #endif @@ -121,11 +120,11 @@ inline SphericalTensor inv(const SphericalTensor& st) } -//- Return the square of Frobenius norm of a SphericalTensor as a Cmpt +//- Return the square of Frobenius norm of a SphericalTensor template -inline Cmpt magSqr(const SphericalTensor& st) +inline Foam::scalar magSqr(const SphericalTensor& st) { - return Cmpt(3*mag(st.ii()*st.ii())); + return st.diagSqr(); } diff --git a/src/OpenFOAM/primitives/SphericalTensor2D/SphericalTensor2D.H b/src/OpenFOAM/primitives/SphericalTensor2D/SphericalTensor2D.H index e40a4542c2..72936c4a9c 100644 --- a/src/OpenFOAM/primitives/SphericalTensor2D/SphericalTensor2D.H +++ b/src/OpenFOAM/primitives/SphericalTensor2D/SphericalTensor2D.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019-2022 OpenCFD Ltd. + Copyright (C) 2019-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -121,6 +121,18 @@ public: const Cmpt& ii() const noexcept { return this->v_[II]; } Cmpt& ii() noexcept { return this->v_[II]; } + + + // Diagonal access and manipulation + + //- The L2-norm squared of the diagonal + inline scalar diagSqr() const; + + + // Tensor Operations + + //- Return non-Hermitian transpose (no-op) + const SphericalTensor2D& T() const noexcept { return *this; } }; diff --git a/src/OpenFOAM/primitives/SphericalTensor2D/SphericalTensor2DI.H b/src/OpenFOAM/primitives/SphericalTensor2D/SphericalTensor2DI.H index 3f3297054a..389637527d 100644 --- a/src/OpenFOAM/primitives/SphericalTensor2D/SphericalTensor2DI.H +++ b/src/OpenFOAM/primitives/SphericalTensor2D/SphericalTensor2DI.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2020 OpenCFD Ltd. + Copyright (C) 2020-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -61,6 +61,15 @@ inline Foam::SphericalTensor2D::SphericalTensor2D(Istream& is) {} +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +inline Foam::scalar Foam::SphericalTensor2D::diagSqr() const +{ + return 2*Foam::magSqr(this->ii()); +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -100,8 +109,8 @@ inline SphericalTensor2D inv(const SphericalTensor2D& st) if (mag(st.ii()) < VSMALL) { FatalErrorInFunction - << "SphericalTensor2D is not invertible due to zero determinant:" - << "det(SphericalTensor2D) = " << det(st) + << "SphericalTensor2D not invertible, determinant:" + << det(st) << " tensor:" << st << nl << abort(FatalError); } #endif diff --git a/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H b/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H index 45ff7aabac..d3f3f5f88f 100644 --- a/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H +++ b/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019-2022 OpenCFD Ltd. + Copyright (C) 2019-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -224,11 +224,38 @@ public: //- Set values of the diagonal inline void diag(const Vector& v); + //- Add to the diagonal + inline void addDiag(const Vector& v); + + //- Subtract from the diagonal + inline void subtractDiag(const Vector& v); + + //- The L2-norm squared of the diagonal + inline scalar diagSqr() const; + // Tensor Operations //- Return non-Hermitian transpose - inline const SymmTensor& T() const; + const SymmTensor& T() const noexcept { return *this; } + + //- The determinate + inline Cmpt det() const; + + //- The 2D determinant by excluding given direction + inline Cmpt det2D(const direction excludeCmpt) const; + + //- Return cofactor matrix + inline SymmTensor cof() const; + + //- Return 2D cofactor matrix by excluding given direction + inline SymmTensor cof2D(const direction excludeCmpt) const; + + //- Return inverse + inline SymmTensor inv() const; + + //- Return inverse of 2D tensor (by excluding given direction) + inline SymmTensor inv2D(const direction excludeCmpt) const; // Member Operators diff --git a/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H b/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H index ed665450e7..56ad9eb4c6 100644 --- a/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H +++ b/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H @@ -213,9 +213,129 @@ inline void Foam::SymmTensor::diag(const Vector& v) template -inline const Foam::SymmTensor& Foam::SymmTensor::T() const +inline void Foam::SymmTensor::addDiag(const Vector& v) { - return *this; + this->v_[XX] += v.x(); this->v_[YY] += v.y(); this->v_[ZZ] += v.z(); +} + + +template +inline void Foam::SymmTensor::subtractDiag(const Vector& v) +{ + this->v_[XX] -= v.x(); this->v_[YY] -= v.y(); this->v_[ZZ] -= v.z(); +} + + +template +inline Foam::scalar Foam::SymmTensor::diagSqr() const +{ + return + ( + Foam::magSqr(this->xx()) + + Foam::magSqr(this->yy()) + + Foam::magSqr(this->zz()) + ); +} + + +template +inline Cmpt Foam::SymmTensor::det() const +{ + return + ( + xx()*yy()*zz() + xy()*yz()*xz() + + xz()*xy()*yz() - xx()*yz()*yz() + - xy()*xy()*zz() - xz()*yy()*xz() + ); +} + + +template +inline Cmpt Foam::SymmTensor::det2D(const direction excludeCmpt) const +{ + switch (excludeCmpt) + { + case 0: // Eliminate x + { + return (yy()*zz() - yz()*yz()); + } + + case 1: // Eliminate y + { + return (xx()*zz() - xz()*xz()); + } + } + + // Fall-through: Eliminate z + return (xx()*yy() - xy()*xy()); +} + + +template +inline Foam::SymmTensor Foam::SymmTensor::cof() const +{ + return SymmTensor + ( + yy()*zz() - yz()*yz(), + xz()*yz() - xy()*zz(), + xy()*yz() - xz()*yy(), + + xx()*zz() - xz()*xz(), + xy()*xz() - xx()*yz(), + + xx()*yy() - xy()*xy() + ); +} + + +template +inline Foam::SymmTensor Foam::SymmTensor::cof2D +( + const direction excludeCmpt +) const +{ + switch (excludeCmpt) + { + case 0: // Eliminate x + { + return SymmTensor + ( + Zero, Zero, Zero, + zz(), -yz(), + yy() + ); + } + + case 1: // Eliminate y + { + return SymmTensor + ( + zz(), Zero, -xz(), + Zero, Zero, + xx() + ); + } + } + + // Fall-through: Eliminate z + return SymmTensor + ( + yy(), -xy(), Zero, + xx(), Zero, + Zero + ); +} + + +template +inline Foam::SymmTensor Foam::SymmTensor::inv2D +( + const direction excludeCmpt +) const +{ + const Cmpt detval = this->det2D(excludeCmpt); + + return this->cof2D(excludeCmpt).T()/detval; } @@ -292,12 +412,7 @@ inline SymmTensor dev2(const SymmTensor& st) template inline Cmpt det(const SymmTensor& st) { - return - ( - st.xx()*st.yy()*st.zz() + st.xy()*st.yz()*st.xz() - + st.xz()*st.xy()*st.yz() - st.xx()*st.yz()*st.yz() - - st.xy()*st.xy()*st.zz() - st.xz()*st.yy()*st.xz() - ); + return st.det(); } @@ -305,35 +420,25 @@ inline Cmpt det(const SymmTensor& st) template inline SymmTensor cof(const SymmTensor& st) { - return SymmTensor - ( - st.yy()*st.zz() - st.yz()*st.yz(), - st.xz()*st.yz() - st.xy()*st.zz(), - st.xy()*st.yz() - st.xz()*st.yy(), - - st.xx()*st.zz() - st.xz()*st.xz(), - st.xy()*st.xz() - st.xx()*st.yz(), - - st.xx()*st.yy() - st.xy()*st.xy() - ); + return st.cof(); } -//- Return the inverse of a SymmTensor by using the given determinant +//- Return the inverse of a SymmTensor, using given determinant value template -inline SymmTensor inv(const SymmTensor& st, const Cmpt detst) +inline SymmTensor inv(const SymmTensor& st, const Cmpt detval) { #ifdef FULLDEBUG - if (mag(detst) < VSMALL) + if (mag(detval) < VSMALL) { FatalErrorInFunction - << "SymmTensor is not invertible due to the zero determinant:" - << "det(symmTensor) = " << mag(detst) + << "SymmTensor not properly invertible, determinant:" + << detval << " tensor:" << st << nl << abort(FatalError); } #endif - return cof(st).T()/detst; + return st.cof().T()/detval; } @@ -341,7 +446,14 @@ inline SymmTensor inv(const SymmTensor& st, const Cmpt detst) template inline SymmTensor inv(const SymmTensor& st) { - return inv(st, det(st)); + return inv(st, st.det()); +} + + +template +inline SymmTensor SymmTensor::inv() const +{ + return Foam::inv(*this, this->det()); } @@ -392,15 +504,15 @@ innerSqr(const SymmTensor& st) } -//- Return the square of Frobenius norm of a SymmTensor as a Cmpt +//- Return the square of Frobenius norm of a SymmTensor template -inline Cmpt magSqr(const SymmTensor& st) +inline Foam::scalar magSqr(const SymmTensor& st) { - return Cmpt + return ( - mag(st.xx()*st.xx()) + 2*mag(st.xy()*st.xy()) + 2*mag(st.xz()*st.xz()) - + mag(st.yy()*st.yy()) + 2*mag(st.yz()*st.yz()) - + mag(st.zz()*st.zz()) + magSqr(st.xx()) + 2*magSqr(st.xy()) + 2*magSqr(st.xz()) + + magSqr(st.yy()) + 2*magSqr(st.yz()) + + magSqr(st.zz()) ); } diff --git a/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2D.H b/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2D.H index a10e93ad4b..1308b015e9 100644 --- a/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2D.H +++ b/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2D.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019-2022 OpenCFD Ltd. + Copyright (C) 2019-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -140,11 +140,23 @@ public: //- Set values of the diagonal inline void diag(const Vector2D& v); + //- The L2-norm squared of the diagonal + inline scalar diagSqr() const; + // Tensor Operations //- Return non-Hermitian transpose - inline const SymmTensor2D& T() const; + const SymmTensor2D& T() const noexcept { return *this; } + + //- The determinate + inline Cmpt det() const; + + //- Return cofactor matrix + inline SymmTensor2D cof() const; + + //- Return inverse + inline SymmTensor2D inv() const; // Member Operators diff --git a/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2DI.H b/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2DI.H index 5f0ab5fdf0..b583b669c5 100644 --- a/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2DI.H +++ b/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2DI.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2019-2020 OpenCFD Ltd. + Copyright (C) 2019-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -91,9 +91,31 @@ inline void Foam::SymmTensor2D::diag(const Vector2D& v) template -inline const Foam::SymmTensor2D& Foam::SymmTensor2D::T() const +inline Foam::scalar Foam::SymmTensor2D::diagSqr() const { - return *this; + return + ( + Foam::magSqr(this->xx()) + + Foam::magSqr(this->yy()) + ); +} + + +template +inline Cmpt Foam::SymmTensor2D::det() const +{ + return (xx()*yy() - xy()*xy()); +} + + +template +inline Foam::SymmTensor2D Foam::SymmTensor2D::cof() const +{ + return SymmTensor2D + ( + yy(), -yx(), + xx() + ); } @@ -172,7 +194,7 @@ inline SymmTensor2D dev2(const SymmTensor2D& st) template inline Cmpt det(const SymmTensor2D& st) { - return (st.xx()*st.yy() - st.xy()*st.xy()); + return st.det(); } @@ -180,29 +202,25 @@ inline Cmpt det(const SymmTensor2D& st) template inline SymmTensor2D cof(const SymmTensor2D& st) { - return SymmTensor2D - ( - st.yy(), -st.yx(), - st.xx() - ); + return st.cof(); } -//- Return the inverse of a SymmTensor2D using a given determinant +//- Return the inverse of a SymmTensor2D, using given determinant value template -inline SymmTensor2D inv(const SymmTensor2D& st, const Cmpt detst) +inline SymmTensor2D inv(const SymmTensor2D& st, const Cmpt detval) { #ifdef FULLDEBUG - if (mag(detst) < SMALL) + if (mag(detval) < SMALL) { FatalErrorInFunction - << "SymmTensor2D is not invertible due to the zero determinant:" - << "det(SymmTensor2D) = " << mag(detst) + << "SymmTensor2D not properly invertible, determinant:" + << detval << " tensor:" << st << nl << abort(FatalError); } #endif - // cofactor and adjugate matrices of SymmTensor are the same, hence no .T() - return cof(st)/detst; + + return st.cof().T()/detval; } @@ -210,7 +228,15 @@ inline SymmTensor2D inv(const SymmTensor2D& st, const Cmpt detst) template inline SymmTensor2D inv(const SymmTensor2D& st) { - return inv(st, det(st)); + return inv(st, st.det()); +} + + +// The inverse of this SymmTensor2D +template +inline SymmTensor2D SymmTensor2D::inv() const +{ + return Foam::inv(*this, this->det()); } @@ -244,14 +270,14 @@ innerSqr(const SymmTensor2D& st) } -//- Return the square of Frobenius norm of a SymmTensor2D as a Cmpt +//- Return the square of Frobenius norm of a SymmTensor2D template -inline Cmpt magSqr(const SymmTensor2D& st) +inline Foam::scalar magSqr(const SymmTensor2D& st) { - return Cmpt + return ( - mag(st.xx()*st.xx()) + 2*mag(st.xy()*st.xy()) - + mag(st.yy()*st.yy()) + magSqr(st.xx()) + 2*magSqr(st.xy()) + + magSqr(st.yy()) ); } diff --git a/src/OpenFOAM/primitives/Tensor/Tensor.H b/src/OpenFOAM/primitives/Tensor/Tensor.H index dfcb5b1c10..eac4435f94 100644 --- a/src/OpenFOAM/primitives/Tensor/Tensor.H +++ b/src/OpenFOAM/primitives/Tensor/Tensor.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2018-2022 OpenCFD Ltd. + Copyright (C) 2018-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -260,6 +260,15 @@ public: //- Set values of the diagonal inline void diag(const Vector& v); + //- Add to the diagonal + inline void addDiag(const Vector& v); + + //- Subtract from the diagonal + inline void subtractDiag(const Vector& v); + + //- The L2-norm squared of the diagonal + inline scalar diagSqr() const; + // Characteristics @@ -272,9 +281,24 @@ public: //- Return non-Hermitian transpose inline Tensor T() const; + //- The determinate + inline Cmpt det() const; + + //- The 2D determinant by excluding given direction + inline Cmpt det2D(const direction excludeCmpt) const; + + //- Return cofactor matrix + inline Tensor cof() const; + + //- Return 2D cofactor matrix by excluding given direction + inline Tensor cof2D(const direction excludeCmpt) const; + //- Return inverse inline Tensor inv() const; + //- Return inverse of 2D tensor (by excluding given direction) + inline Tensor inv2D(const direction excludeCmpt) const; + //- Inner-product of this with another Tensor. inline Tensor inner(const Tensor& t2) const; diff --git a/src/OpenFOAM/primitives/Tensor/TensorI.H b/src/OpenFOAM/primitives/Tensor/TensorI.H index 19173ee725..a0254112bb 100644 --- a/src/OpenFOAM/primitives/Tensor/TensorI.H +++ b/src/OpenFOAM/primitives/Tensor/TensorI.H @@ -379,6 +379,32 @@ inline void Foam::Tensor::diag(const Vector& v) } +template +inline void Foam::Tensor::addDiag(const Vector& v) +{ + this->v_[XX] += v.x(); this->v_[YY] += v.y(); this->v_[ZZ] += v.z(); +} + + +template +inline void Foam::Tensor::subtractDiag(const Vector& v) +{ + this->v_[XX] -= v.x(); this->v_[YY] -= v.y(); this->v_[ZZ] -= v.z(); +} + + +template +inline Foam::scalar Foam::Tensor::diagSqr() const +{ + return + ( + Foam::magSqr(this->xx()) + + Foam::magSqr(this->yy()) + + Foam::magSqr(this->zz()) + ); +} + + template inline bool Foam::Tensor::is_identity(const scalar tol) const { @@ -408,6 +434,110 @@ inline Foam::Tensor Foam::Tensor::T() const } +template +inline Cmpt Foam::Tensor::det() const +{ + return + ( + xx()*yy()*zz() + xy()*yz()*zx() + + xz()*yx()*zy() - xx()*yz()*zy() + - xy()*yx()*zz() - xz()*yy()*zx() + ); +} + + +template +inline Cmpt Foam::Tensor::det2D(const direction excludeCmpt) const +{ + switch (excludeCmpt) + { + case 0: // Eliminate x + { + return (yy()*zz() - yz()*zy()); + } + + case 1: // Eliminate y + { + return (xx()*zz() - xz()*zx()); + } + } + + // Fall-through: Eliminate z + return (xx()*yy() - xy()*yx()); +} + + +template +inline Foam::Tensor Foam::Tensor::cof() const +{ + return Tensor + ( + yy()*zz() - zy()*yz(), + zx()*yz() - yx()*zz(), + yx()*zy() - yy()*zx(), + + xz()*zy() - xy()*zz(), + xx()*zz() - xz()*zx(), + xy()*zx() - xx()*zy(), + + xy()*yz() - xz()*yy(), + yx()*xz() - xx()*yz(), + xx()*yy() - yx()*xy() + ); +} + + +template +inline Foam::Tensor Foam::Tensor::cof2D +( + const direction excludeCmpt +) const +{ + switch (excludeCmpt) + { + case 0: // Eliminate x + { + return Tensor + ( + Zero, Zero, Zero, + Zero, zz(), -zy(), + Zero, -yz(), yy() + ); + } + + case 1: // Eliminate y + { + return Tensor + ( + zz(), Zero, -zx(), + Zero, Zero, Zero, + -xz(), Zero, xx() + ); + } + } + + // Fall-through: Eliminate z + return Tensor + ( + yy(), -yx(), Zero, + -xy(), xx(), Zero, + Zero, Zero, Zero + ); +} + + +template +inline Foam::Tensor Foam::Tensor::inv2D +( + const direction excludeCmpt +) const +{ + const Cmpt detval = this->det2D(excludeCmpt); + + return this->cof2D(excludeCmpt).T()/detval; +} + + template inline Foam::Tensor Foam::Tensor::inner(const Tensor& t2) const @@ -594,12 +724,7 @@ inline Tensor dev2(const Tensor& t) template inline Cmpt det(const Tensor& 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() - ); + return t.det(); } @@ -607,39 +732,24 @@ inline Cmpt det(const Tensor& t) template inline Tensor cof(const Tensor& t) { - return Tensor - ( - 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() - ); + return t.cof(); } -//- Return the inverse of a Tensor by using the given determinant +//- Return the inverse of a Tensor, using the given determinant value template -inline Tensor inv(const Tensor& t, const Cmpt dett) +inline Tensor inv(const Tensor& t, const Cmpt detval) { #ifdef FULLDEBUG - if (mag(dett) < VSMALL) + if (mag(detval) < VSMALL) { FatalErrorInFunction - << "Tensor is not invertible due to the (almost) zero determinant:" - << " Tensor = " << t << nl - << " det(Tensor) = " << dett - << abort(FatalError); + << "Tensor not properly invertible, determinant:" + << detval << " tensor:" << t << nl } #endif - return cof(t).T()/dett; + return t.cof().T()/detval; } @@ -647,15 +757,14 @@ inline Tensor inv(const Tensor& t, const Cmpt dett) template inline Tensor inv(const Tensor& t) { - return inv(t, det(t)); + return inv(t, t.det()); } -//- Return the inverse of this Tensor template inline Tensor Tensor::inv() const { - return Foam::inv(*this); + return Foam::inv(*this, this->det()); } diff --git a/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H b/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H index 079aa9501a..f12b68104f 100644 --- a/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H +++ b/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2018-2022 OpenCFD Ltd. + Copyright (C) 2018-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -200,12 +200,24 @@ public: //- Set values of the diagonal inline void diag(const Vector2D& v); + //- The L2-norm squared of the diagonal + inline scalar diagSqr() const; + // Tensor Operations //- Return non-Hermitian transpose inline Tensor2D T() const; + //- The determinate + inline Cmpt det() const; + + //- Return cofactor matrix + inline Tensor2D cof() const; + + //- Return inverse + inline Tensor2D inv() const; + //- Inner-product of this with another Tensor2D. inline Tensor2D inner(const Tensor2D& t2) const; diff --git a/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H b/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H index 178ad48496..9367f4c82b 100644 --- a/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H +++ b/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2018-2022 OpenCFD Ltd. + Copyright (C) 2018-2023 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -286,6 +286,17 @@ inline void Foam::Tensor2D::diag(const Vector2D& v) } +template +inline Foam::scalar Foam::Tensor2D::diagSqr() const +{ + return + ( + Foam::magSqr(this->xx()) + + Foam::magSqr(this->yy()) + ); +} + + // * * * * * * * * * * * * * * * Member Operations * * * * * * * * * * * * * // template @@ -299,6 +310,24 @@ inline Foam::Tensor2D Foam::Tensor2D::T() const } +template +inline Cmpt Foam::Tensor2D::det() const +{ + return (xx()*yy() - xy()*yx()); +} + + +template +inline Foam::Tensor2D Foam::Tensor2D::cof() const +{ + return Tensor2D + ( + yy(), -yx(), + -xy(), xx() + ); +} + + template inline Foam::Tensor2D Foam::Tensor2D::inner(const Tensor2D& t2) const @@ -430,7 +459,7 @@ inline Tensor2D dev2(const Tensor2D& t) template inline Cmpt det(const Tensor2D& t) { - return t.xx()*t.yy() - t.xy()*t.yx(); + return t.det(); } @@ -438,29 +467,25 @@ inline Cmpt det(const Tensor2D& t) template inline Tensor2D cof(const Tensor2D& t) { - return Tensor2D - ( - t.yy(), -t.yx(), - -t.xy(), t.xx() - ); + return t.cof(); } -//- Return the inverse of a Tensor2D using a given determinant +//- Return the inverse of a Tensor2D, using given determinant value template -inline Tensor2D inv(const Tensor2D& t, const Cmpt dett) +inline Tensor2D inv(const Tensor2D& t, const Cmpt detval) { #ifdef FULLDEBUG - if (mag(dett) < SMALL) + if (mag(detval) < SMALL) { FatalErrorInFunction - << "Tensor2D is not invertible due to the zero determinant:" - << "det(Tensor2D) = " << mag(dett) + << "Tensor2D not properly invertible, determinant:" + << detval << " tensor:" << t << nl << abort(FatalError); } #endif - return cof(t).T()/dett; + return t.cof().T()/detval; } @@ -468,7 +493,15 @@ inline Tensor2D inv(const Tensor2D& t, const Cmpt dett) template inline Tensor2D inv(const Tensor2D& t) { - return inv(t, det(t)); + return inv(t, t.det()); +} + + +// Return the inverse of this Tensor2D +template +inline Tensor2D Tensor2D::inv() const +{ + return Foam::inv(*this, this->det()); }