ENH: add 2D det() / inv() methods for Tensor/SymmTensor (#2724)

- for cases where a 3D tensor is being used to represent 2D content,
  the determinant is zero. Can use inv2D(excludeDirection) to compensate
  and invert as if it were only 2D.

ENH: consistent definitions for magSqr of symmTensors, diagSqr() norm

COMP: return scalar not component type for magSqr

- had inconsistent definitions with SymmTensor returning the component
  type and Tensor returning scalar. Only evident with complex.

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));
 }
 
 

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 = ",

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<complex> st1(SymmTensor<complex>::uniform({3, 4}));
+
+        Info<< "symmTensor: " << st1 << nl
+            << "  tr: " << tr(st1) << nl
+            << "  diagSqr: " << st1.diagSqr() << nl
+            << "  magSqr: " << magSqr(st1) << nl
+            << "  mag: " << mag(st1) << nl;
+
+        SymmTensor<scalar> st2(SymmTensor<scalar>::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<complex> dt1(SphericalTensor<complex>({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<complex>({3, 4});
+
+        Info<< "symmTensor: " << st1 << nl
+            << "  tr: " << tr(st1) << nl
+            << "  diagSqr: " << st1.diagSqr() << nl
+            << "  magSqr: " << magSqr(st1) << nl
+            << "  mag: " << mag(st1) << nl;
+    }
+
+    {
+        SymmTensor2D<complex> st1(SymmTensor2D<complex>::uniform({3, 4}));
+
+        Info<< "symmTensor: " << st1 << nl
+            << "  tr: " << tr(st1) << nl
+            << "  diagSqr: " << st1.diagSqr() << nl
+            << "  magSqr: " << magSqr(st1) << nl
+            << "  mag: " << mag(st1) << nl;
+    }
+
+    {
+        Tensor<complex> st1(Tensor<complex>::uniform({3, 4}));
+
+        Info<< "tensor: " << st1 << nl
+            << "  tr: " << tr(st1) << nl
+            << "  diagSqr: " << st1.diagSqr() << nl
+            << "  magSqr: " << magSqr(st1) << nl
+            << "  mag: " << mag(st1) << endl;
+
+        Tensor<scalar> st2(Tensor<scalar>::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);
 

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<class Cmpt2>
         inline DiagTensor(const VectorSpace<DiagTensor<Cmpt2>, Cmpt2, 3>&);
 
+        //- Construct given SphericalTensor
+        inline DiagTensor(const SphericalTensor<Cmpt>&);
+
         //- 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;
 };
 
 

src/OpenFOAM/primitives/DiagTensor/DiagTensorI.H

@@ -49,6 +49,15 @@ inline Foam::DiagTensor<Cmpt>::DiagTensor
 {}
 
 
+template<class Cmpt>
+inline Foam::DiagTensor<Cmpt>::DiagTensor(const SphericalTensor<Cmpt>& st)
+{
+    this->v_[XX] = st.ii();
+    this->v_[YY] = st.ii();
+    this->v_[ZZ] = st.ii();
+}
+
+
 template<class Cmpt>
 inline Foam::DiagTensor<Cmpt>::DiagTensor
 (
@@ -70,6 +79,20 @@ inline Foam::DiagTensor<Cmpt>::DiagTensor(Istream& is)
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Cmpt>
+inline Foam::scalar Foam::DiagTensor<Cmpt>::diagSqr() const
+{
+    return
+    (
+        Foam::magSqr(this->xx())
+      + Foam::magSqr(this->yy())
+      + Foam::magSqr(this->zz())
+    );
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -138,6 +161,13 @@ inline DiagTensor<Cmpt> diag(const SymmTensor<Cmpt>& st)
 }
 
 
+template<class Cmpt>
+inline Foam::scalar magSqr(const DiagTensor<Cmpt>& t)
+{
+    return t.diagSqr();
+}
+
+
 //- Linear interpolation of diagonal tensors a and b by factor t
 template<class Cmpt>
 inline DiagTensor<Cmpt> lerp

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<Cmpt>& T() const;
+        const SphericalTensor<Cmpt>& T() const noexcept { return *this; }
 };
 
 

src/OpenFOAM/primitives/SphericalTensor/SphericalTensorI.H

@@ -65,10 +65,9 @@ inline Foam::SphericalTensor<Cmpt>::SphericalTensor(Istream& is)
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Cmpt>
-inline const Foam::SphericalTensor<Cmpt>&
+inline Foam::scalar Foam::SphericalTensor<Cmpt>::diagSqr() const
-Foam::SphericalTensor<Cmpt>::T() const
 {
-    return *this;
+    return 3*Foam::magSqr(this->ii());
 }
 
 
@@ -111,8 +110,8 @@ inline SphericalTensor<Cmpt> inv(const SphericalTensor<Cmpt>& st)
     if (mag(st.ii()) < VSMALL)
     {
         FatalErrorInFunction
-            << "SphericalTensor is not invertible due to the zero determinant:"
+            << "SphericalTensor not invertible, determinant:"
-            << "det(SphericalTensor) = " << det(st)
+            << det(st) << " tensor:" << st << nl
             << abort(FatalError);
     }
     #endif
@@ -121,11 +120,11 @@ inline SphericalTensor<Cmpt> inv(const SphericalTensor<Cmpt>& st)
 }
 
 
-//- Return the square of Frobenius norm of a SphericalTensor as a Cmpt
+//- Return the square of Frobenius norm of a SphericalTensor
 template<class Cmpt>
-inline Cmpt magSqr(const SphericalTensor<Cmpt>& st)
+inline Foam::scalar magSqr(const SphericalTensor<Cmpt>& st)
 {
-    return Cmpt(3*mag(st.ii()*st.ii()));
+    return st.diagSqr();
 }
 
 

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<Cmpt>& T() const noexcept { return *this; }
 };
 
 

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<Cmpt>::SphericalTensor2D(Istream& is)
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Cmpt>
+inline Foam::scalar Foam::SphericalTensor2D<Cmpt>::diagSqr() const
+{
+    return 2*Foam::magSqr(this->ii());
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -100,8 +109,8 @@ inline SphericalTensor2D<Cmpt> inv(const SphericalTensor2D<Cmpt>& st)
     if (mag(st.ii()) < VSMALL)
     {
         FatalErrorInFunction
-            << "SphericalTensor2D is not invertible due to zero determinant:"
+            << "SphericalTensor2D not invertible, determinant:"
-            << "det(SphericalTensor2D) = " << det(st)
+            << det(st) << " tensor:" << st << nl
             << abort(FatalError);
     }
     #endif

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<Cmpt>& v);
 
+        //- Add to the diagonal
+        inline void addDiag(const Vector<Cmpt>& v);
+
+        //- Subtract from the diagonal
+        inline void subtractDiag(const Vector<Cmpt>& v);
+
+        //- The L2-norm squared of the diagonal
+        inline scalar diagSqr() const;
+
 
     // Tensor Operations
 
         //- Return non-Hermitian transpose
-        inline const SymmTensor<Cmpt>& T() const;
+        const SymmTensor<Cmpt>& 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<Cmpt> cof() const;
+
+        //- Return 2D cofactor matrix by excluding given direction
+        inline SymmTensor<Cmpt> cof2D(const direction excludeCmpt) const;
+
+        //- Return inverse
+        inline SymmTensor<Cmpt> inv() const;
+
+        //- Return inverse of 2D tensor (by excluding given direction)
+        inline SymmTensor<Cmpt> inv2D(const direction excludeCmpt) const;
 
 
     // Member Operators

src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H

@@ -213,9 +213,129 @@ inline void Foam::SymmTensor<Cmpt>::diag(const Vector<Cmpt>& v)
 
 
 template<class Cmpt>
-inline const Foam::SymmTensor<Cmpt>& Foam::SymmTensor<Cmpt>::T() const
+inline void Foam::SymmTensor<Cmpt>::addDiag(const Vector<Cmpt>& v)
 {
-    return *this;
+    this->v_[XX] += v.x(); this->v_[YY] += v.y(); this->v_[ZZ] += v.z();
+}
+
+
+template<class Cmpt>
+inline void Foam::SymmTensor<Cmpt>::subtractDiag(const Vector<Cmpt>& v)
+{
+    this->v_[XX] -= v.x(); this->v_[YY] -= v.y(); this->v_[ZZ] -= v.z();
+}
+
+
+template<class Cmpt>
+inline Foam::scalar Foam::SymmTensor<Cmpt>::diagSqr() const
+{
+    return
+    (
+        Foam::magSqr(this->xx())
+      + Foam::magSqr(this->yy())
+      + Foam::magSqr(this->zz())
+    );
+}
+
+
+template<class Cmpt>
+inline Cmpt Foam::SymmTensor<Cmpt>::det() const
+{
+    return
+    (
+        xx()*yy()*zz() + xy()*yz()*xz()
+      + xz()*xy()*yz() - xx()*yz()*yz()
+      - xy()*xy()*zz() - xz()*yy()*xz()
+    );
+}
+
+
+template<class Cmpt>
+inline Cmpt Foam::SymmTensor<Cmpt>::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<class Cmpt>
+inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::cof() const
+{
+    return SymmTensor<Cmpt>
+    (
+        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<class Cmpt>
+inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::cof2D
+(
+    const direction excludeCmpt
+) const
+{
+    switch (excludeCmpt)
+    {
+        case 0:  // Eliminate x
+        {
+            return SymmTensor<Cmpt>
+            (
+                Zero, Zero,  Zero,
+                      zz(), -yz(),
+                             yy()
+            );
+        }
+
+        case 1:  // Eliminate y
+        {
+            return SymmTensor<Cmpt>
+            (
+                zz(), Zero, -xz(),
+                      Zero,  Zero,
+                             xx()
+            );
+        }
+    }
+
+    // Fall-through: Eliminate z
+    return SymmTensor<Cmpt>
+    (
+        yy(), -xy(), Zero,
+               xx(), Zero,
+                     Zero
+    );
+}
+
+
+template<class Cmpt>
+inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::inv2D
+(
+    const direction excludeCmpt
+) const
+{
+    const Cmpt detval = this->det2D(excludeCmpt);
+
+    return this->cof2D(excludeCmpt).T()/detval;
 }
 
 
@@ -292,12 +412,7 @@ inline SymmTensor<Cmpt> dev2(const SymmTensor<Cmpt>& st)
 template<class Cmpt>
 inline Cmpt det(const SymmTensor<Cmpt>& st)
 {
-    return
+    return st.det();
-    (
-        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()
-    );
 }
 
 
@@ -305,35 +420,25 @@ inline Cmpt det(const SymmTensor<Cmpt>& st)
 template<class Cmpt>
 inline SymmTensor<Cmpt> cof(const SymmTensor<Cmpt>& st)
 {
-    return SymmTensor<Cmpt>
+    return st.cof();
-    (
-        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 the inverse of a SymmTensor by using the given determinant
+//- Return the inverse of a SymmTensor, using given determinant value
 template<class Cmpt>
-inline SymmTensor<Cmpt> inv(const SymmTensor<Cmpt>& st, const Cmpt detst)
+inline SymmTensor<Cmpt> inv(const SymmTensor<Cmpt>& st, const Cmpt detval)
 {
     #ifdef FULLDEBUG
-    if (mag(detst) < VSMALL)
+    if (mag(detval) < VSMALL)
     {
         FatalErrorInFunction
-            << "SymmTensor is not invertible due to the zero determinant:"
+            << "SymmTensor not properly invertible, determinant:"
-            << "det(symmTensor) = " << mag(detst)
+            << detval << " tensor:" << st << nl
             << abort(FatalError);
     }
     #endif
 
-    return cof(st).T()/detst;
+    return st.cof().T()/detval;
 }
 
 
@@ -341,7 +446,14 @@ inline SymmTensor<Cmpt> inv(const SymmTensor<Cmpt>& st, const Cmpt detst)
 template<class Cmpt>
 inline SymmTensor<Cmpt> inv(const SymmTensor<Cmpt>& st)
 {
-    return inv(st, det(st));
+    return inv(st, st.det());
+}
+
+
+template<class Cmpt>
+inline SymmTensor<Cmpt> SymmTensor<Cmpt>::inv() const
+{
+    return Foam::inv(*this, this->det());
 }
 
 
@@ -392,15 +504,15 @@ innerSqr(const SymmTensor<Cmpt>& st)
 }
 
 
-//- Return the square of Frobenius norm of a SymmTensor as a Cmpt
+//- Return the square of Frobenius norm of a SymmTensor
 template<class Cmpt>
-inline Cmpt magSqr(const SymmTensor<Cmpt>& st)
+inline Foam::scalar magSqr(const SymmTensor<Cmpt>& st)
 {
-    return Cmpt
+    return
     (
-        mag(st.xx()*st.xx()) + 2*mag(st.xy()*st.xy()) + 2*mag(st.xz()*st.xz())
+        magSqr(st.xx()) + 2*magSqr(st.xy()) + 2*magSqr(st.xz())
-                             +   mag(st.yy()*st.yy()) + 2*mag(st.yz()*st.yz())
+                        +   magSqr(st.yy()) + 2*magSqr(st.yz())
-                                                      +   mag(st.zz()*st.zz())
+                                            +   magSqr(st.zz())
     );
 }
 

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<Cmpt>& v);
 
+        //- The L2-norm squared of the diagonal
+        inline scalar diagSqr() const;
+
 
     // Tensor Operations
 
         //- Return non-Hermitian transpose
-        inline const SymmTensor2D<Cmpt>& T() const;
+        const SymmTensor2D<Cmpt>& T() const noexcept { return *this; }
+
+        //- The determinate
+        inline Cmpt det() const;
+
+        //- Return cofactor matrix
+        inline SymmTensor2D<Cmpt> cof() const;
+
+        //- Return inverse
+        inline SymmTensor2D<Cmpt> inv() const;
 
 
     // Member Operators

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<Cmpt>::diag(const Vector2D<Cmpt>& v)
 
 
 template<class Cmpt>
-inline const Foam::SymmTensor2D<Cmpt>& Foam::SymmTensor2D<Cmpt>::T() const
+inline Foam::scalar Foam::SymmTensor2D<Cmpt>::diagSqr() const
 {
-    return *this;
+    return
+    (
+        Foam::magSqr(this->xx())
+      + Foam::magSqr(this->yy())
+    );
+}
+
+
+template<class Cmpt>
+inline Cmpt Foam::SymmTensor2D<Cmpt>::det() const
+{
+    return (xx()*yy() - xy()*xy());
+}
+
+
+template<class Cmpt>
+inline Foam::SymmTensor2D<Cmpt> Foam::SymmTensor2D<Cmpt>::cof() const
+{
+    return SymmTensor2D<Cmpt>
+    (
+        yy(), -yx(),
+               xx()
+    );
 }
 
 
@@ -172,7 +194,7 @@ inline SymmTensor2D<Cmpt> dev2(const SymmTensor2D<Cmpt>& st)
 template<class Cmpt>
 inline Cmpt det(const SymmTensor2D<Cmpt>& st)
 {
-    return (st.xx()*st.yy() - st.xy()*st.xy());
+    return st.det();
 }
 
 
@@ -180,29 +202,25 @@ inline Cmpt det(const SymmTensor2D<Cmpt>& st)
 template<class Cmpt>
 inline SymmTensor2D<Cmpt> cof(const SymmTensor2D<Cmpt>& st)
 {
-    return SymmTensor2D<Cmpt>
+    return st.cof();
-    (
-        st.yy(), -st.yx(),
-                  st.xx()
-    );
 }
 
 
-//- Return the inverse of a SymmTensor2D using a given determinant
+//- Return the inverse of a SymmTensor2D, using given determinant value
 template<class Cmpt>
-inline SymmTensor2D<Cmpt> inv(const SymmTensor2D<Cmpt>& st, const Cmpt detst)
+inline SymmTensor2D<Cmpt> inv(const SymmTensor2D<Cmpt>& st, const Cmpt detval)
 {
     #ifdef FULLDEBUG
-    if (mag(detst) < SMALL)
+    if (mag(detval) < SMALL)
     {
         FatalErrorInFunction
-            << "SymmTensor2D is not invertible due to the zero determinant:"
+            << "SymmTensor2D not properly invertible, determinant:"
-            << "det(SymmTensor2D) = " << mag(detst)
+            << 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<Cmpt> inv(const SymmTensor2D<Cmpt>& st, const Cmpt detst)
 template<class Cmpt>
 inline SymmTensor2D<Cmpt> inv(const SymmTensor2D<Cmpt>& st)
 {
-    return inv(st, det(st));
+    return inv(st, st.det());
+}
+
+
+// The inverse of this SymmTensor2D
+template<class Cmpt>
+inline SymmTensor2D<Cmpt> SymmTensor2D<Cmpt>::inv() const
+{
+    return Foam::inv(*this, this->det());
 }
 
 
@@ -244,14 +270,14 @@ innerSqr(const SymmTensor2D<Cmpt>& st)
 }
 
 
-//- Return the square of Frobenius norm of a SymmTensor2D as a Cmpt
+//- Return the square of Frobenius norm of a SymmTensor2D
 template<class Cmpt>
-inline Cmpt magSqr(const SymmTensor2D<Cmpt>& st)
+inline Foam::scalar magSqr(const SymmTensor2D<Cmpt>& st)
 {
-    return Cmpt
+    return
     (
-        mag(st.xx()*st.xx()) + 2*mag(st.xy()*st.xy())
+        magSqr(st.xx()) + 2*magSqr(st.xy())
-                             +   mag(st.yy()*st.yy())
+                        +   magSqr(st.yy())
     );
 }
 

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<Cmpt>& v);
 
+        //- Add to the diagonal
+        inline void addDiag(const Vector<Cmpt>& v);
+
+        //- Subtract from the diagonal
+        inline void subtractDiag(const Vector<Cmpt>& v);
+
+        //- The L2-norm squared of the diagonal
+        inline scalar diagSqr() const;
+
 
     // Characteristics
 
@@ -272,9 +281,24 @@ public:
         //- Return non-Hermitian transpose
         inline Tensor<Cmpt> 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<Cmpt> cof() const;
+
+        //- Return 2D cofactor matrix by excluding given direction
+        inline Tensor<Cmpt> cof2D(const direction excludeCmpt) const;
+
         //- Return inverse
         inline Tensor<Cmpt> inv() const;
 
+        //- Return inverse of 2D tensor (by excluding given direction)
+        inline Tensor<Cmpt> inv2D(const direction excludeCmpt) const;
+
         //- Inner-product of this with another Tensor.
         inline Tensor<Cmpt> inner(const Tensor<Cmpt>& t2) const;
 

src/OpenFOAM/primitives/Tensor/TensorI.H

@@ -379,6 +379,32 @@ inline void Foam::Tensor<Cmpt>::diag(const Vector<Cmpt>& v)
 }
 
 
+template<class Cmpt>
+inline void Foam::Tensor<Cmpt>::addDiag(const Vector<Cmpt>& v)
+{
+    this->v_[XX] += v.x(); this->v_[YY] += v.y(); this->v_[ZZ] += v.z();
+}
+
+
+template<class Cmpt>
+inline void Foam::Tensor<Cmpt>::subtractDiag(const Vector<Cmpt>& v)
+{
+    this->v_[XX] -= v.x(); this->v_[YY] -= v.y(); this->v_[ZZ] -= v.z();
+}
+
+
+template<class Cmpt>
+inline Foam::scalar Foam::Tensor<Cmpt>::diagSqr() const
+{
+    return
+    (
+        Foam::magSqr(this->xx())
+      + Foam::magSqr(this->yy())
+      + Foam::magSqr(this->zz())
+    );
+}
+
+
 template<class Cmpt>
 inline bool Foam::Tensor<Cmpt>::is_identity(const scalar tol) const
 {
@@ -408,6 +434,110 @@ inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::T() const
 }
 
 
+template<class Cmpt>
+inline Cmpt Foam::Tensor<Cmpt>::det() const
+{
+    return
+    (
+        xx()*yy()*zz() + xy()*yz()*zx()
+      + xz()*yx()*zy() - xx()*yz()*zy()
+      - xy()*yx()*zz() - xz()*yy()*zx()
+    );
+}
+
+
+template<class Cmpt>
+inline Cmpt Foam::Tensor<Cmpt>::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<class Cmpt>
+inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::cof() const
+{
+    return Tensor<Cmpt>
+    (
+        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<class Cmpt>
+inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::cof2D
+(
+    const direction excludeCmpt
+) const
+{
+    switch (excludeCmpt)
+    {
+        case 0:  // Eliminate x
+        {
+            return Tensor<Cmpt>
+            (
+                Zero,  Zero,  Zero,
+                Zero,  zz(), -zy(),
+                Zero, -yz(),  yy()
+            );
+        }
+
+        case 1:  // Eliminate y
+        {
+            return Tensor<Cmpt>
+            (
+                zz(), Zero, -zx(),
+                Zero, Zero,  Zero,
+               -xz(), Zero,  xx()
+            );
+        }
+    }
+
+    // Fall-through: Eliminate z
+    return Tensor<Cmpt>
+    (
+        yy(), -yx(), Zero,
+       -xy(),  xx(), Zero,
+        Zero,  Zero, Zero
+    );
+}
+
+
+template<class Cmpt>
+inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::inv2D
+(
+    const direction excludeCmpt
+) const
+{
+    const Cmpt detval = this->det2D(excludeCmpt);
+
+    return this->cof2D(excludeCmpt).T()/detval;
+}
+
+
 template<class Cmpt>
 inline Foam::Tensor<Cmpt>
 Foam::Tensor<Cmpt>::inner(const Tensor<Cmpt>& t2) const
@@ -594,12 +724,7 @@ inline Tensor<Cmpt> dev2(const Tensor<Cmpt>& t)
 template<class Cmpt>
 inline Cmpt det(const Tensor<Cmpt>& t)
 {
-    return
+    return t.det();
-    (
-        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()
-    );
 }
 
 
@@ -607,39 +732,24 @@ inline Cmpt det(const Tensor<Cmpt>& t)
 template<class Cmpt>
 inline Tensor<Cmpt> cof(const Tensor<Cmpt>& t)
 {
-    return Tensor<Cmpt>
+    return t.cof();
-    (
-        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 the inverse of a Tensor by using the given determinant
+//- Return the inverse of a Tensor, using the given determinant value
 template<class Cmpt>
-inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t, const Cmpt dett)
+inline Tensor<Cmpt> inv(const Tensor<Cmpt>& 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 not properly invertible, determinant:"
-            << "    Tensor = " << t << nl
+            << detval << " tensor:" << t << nl
-            << "    det(Tensor) = " << dett
-            << abort(FatalError);
     }
     #endif
 
-    return cof(t).T()/dett;
+    return t.cof().T()/detval;
 }
 
 
@@ -647,15 +757,14 @@ inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t, const Cmpt dett)
 template<class Cmpt>
 inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t)
 {
-    return inv(t, det(t));
+    return inv(t, t.det());
 }
 
 
-//- Return the inverse of this Tensor
 template<class Cmpt>
 inline Tensor<Cmpt> Tensor<Cmpt>::inv() const
 {
-    return Foam::inv(*this);
+    return Foam::inv(*this, this->det());
 }
 
 

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<Cmpt>& v);
 
+        //- The L2-norm squared of the diagonal
+        inline scalar diagSqr() const;
+
 
     // Tensor Operations
 
         //- Return non-Hermitian transpose
         inline Tensor2D<Cmpt> T() const;
 
+        //- The determinate
+        inline Cmpt det() const;
+
+        //- Return cofactor matrix
+        inline Tensor2D<Cmpt> cof() const;
+
+        //- Return inverse
+        inline Tensor2D<Cmpt> inv() const;
+
         //- Inner-product of this with another Tensor2D.
         inline Tensor2D<Cmpt> inner(const Tensor2D<Cmpt>& t2) const;
 

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<Cmpt>::diag(const Vector2D<Cmpt>& v)
 }
 
 
+template<class Cmpt>
+inline Foam::scalar Foam::Tensor2D<Cmpt>::diagSqr() const
+{
+    return
+    (
+        Foam::magSqr(this->xx())
+      + Foam::magSqr(this->yy())
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Operations * * * * * * * * * * * * * //
 
 template<class Cmpt>
@@ -299,6 +310,24 @@ inline Foam::Tensor2D<Cmpt> Foam::Tensor2D<Cmpt>::T() const
 }
 
 
+template<class Cmpt>
+inline Cmpt Foam::Tensor2D<Cmpt>::det() const
+{
+    return (xx()*yy() - xy()*yx());
+}
+
+
+template<class Cmpt>
+inline Foam::Tensor2D<Cmpt> Foam::Tensor2D<Cmpt>::cof() const
+{
+    return Tensor2D<Cmpt>
+    (
+        yy(), -yx(),
+       -xy(),  xx()
+    );
+}
+
+
 template<class Cmpt>
 inline Foam::Tensor2D<Cmpt>
 Foam::Tensor2D<Cmpt>::inner(const Tensor2D<Cmpt>& t2) const
@@ -430,7 +459,7 @@ inline Tensor2D<Cmpt> dev2(const Tensor2D<Cmpt>& t)
 template<class Cmpt>
 inline Cmpt det(const Tensor2D<Cmpt>& t)
 {
-    return t.xx()*t.yy() - t.xy()*t.yx();
+    return t.det();
 }
 
 
@@ -438,29 +467,25 @@ inline Cmpt det(const Tensor2D<Cmpt>& t)
 template<class Cmpt>
 inline Tensor2D<Cmpt> cof(const Tensor2D<Cmpt>& t)
 {
-    return Tensor2D<Cmpt>
+    return t.cof();
-    (
-        t.yy(), -t.yx(),
-       -t.xy(),  t.xx()
-    );
 }
 
 
-//- Return the inverse of a Tensor2D using a given determinant
+//- Return the inverse of a Tensor2D, using given determinant value
 template<class Cmpt>
-inline Tensor2D<Cmpt> inv(const Tensor2D<Cmpt>& t, const Cmpt dett)
+inline Tensor2D<Cmpt> inv(const Tensor2D<Cmpt>& t, const Cmpt detval)
 {
     #ifdef FULLDEBUG
-    if (mag(dett) < SMALL)
+    if (mag(detval) < SMALL)
     {
         FatalErrorInFunction
-            << "Tensor2D is not invertible due to the zero determinant:"
+            << "Tensor2D not properly invertible, determinant:"
-            << "det(Tensor2D) = " << mag(dett)
+            << detval << " tensor:" << t << nl
             << abort(FatalError);
     }
     #endif
 
-    return cof(t).T()/dett;
+    return t.cof().T()/detval;
 }
 
 
@@ -468,7 +493,15 @@ inline Tensor2D<Cmpt> inv(const Tensor2D<Cmpt>& t, const Cmpt dett)
 template<class Cmpt>
 inline Tensor2D<Cmpt> inv(const Tensor2D<Cmpt>& t)
 {
-    return inv(t, det(t));
+    return inv(t, t.det());
+}
+
+
+// Return the inverse of this Tensor2D
+template<class Cmpt>
+inline Tensor2D<Cmpt> Tensor2D<Cmpt>::inv() const
+{
+    return Foam::inv(*this, this->det());
 }