ENH: improve funcs and opers in Tensor types

- ensures each Tensor-container operates for the following base types:
    - floatScalar
    - doubleScalar
    - complex

  - adds/improves test applications for each container and base type:
    - constructors
    - member functions
    - global functions
    - global operators

  - misc:
    - silently removes `invariantIII()` for `tensor2D` and `symmTensor2D`
      since the 3rd invariant does not exist for 2x2 matrices
    - fixes `invariantII()` algorithm for `tensor2D` and `symmTensor2D`
    - adds `Cmpt` multiplication to `Vector2D` and `Vector`
    - adds missing access funcs for symmetric containers
    - improves func/header documentations

applications/test/DiagTensor/Make/files

@@ -0,0 +1,3 @@
+Test-DiagTensor.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-DiagTensor

applications/test/DiagTensor/Test-DiagTensor.C

@@ -0,0 +1,439 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    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.
+
+    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
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-DiagTensor
+
+Description
+    Tests for \c DiagTensor constructors, member functions and operators
+    using \c floatScalar, \c doubleScalar, and \c complex base types.
+
+    Cross-checks were obtained from 'NumPy 1.15.1' and 'SciPy 1.1.0' if no
+    theoretical cross-check exists (like eigendecomposition relations), and
+    were hard-coded for elementwise comparisons.
+
+    For \c complex base type, the cross-checks do only involve zero imag part.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Tensor.H"
+#include "SymmTensor.H"
+#include "SphericalTensor.H"
+#include "DiagTensor.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Total number of unit tests
+unsigned nTest_ = 0;
+
+
+// Total number of failed unit tests
+unsigned nFail_ = 0;
+
+
+// Compare two floating point types, and print output.
+// Do ++nFail_ if values of two objects are not equal within a given tolerance.
+// The function is converted from PEP-485.
+template<class Type>
+typename std::enable_if
+<
+    std::is_same<floatScalar, Type>::value ||
+    std::is_same<doubleScalar, Type>::value ||
+    std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,  //<! are values the same within 8 decimals
+    const scalar absTol = 0      //<! useful for cmps near zero
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
+    {
+        ++nFail;
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Compare two containers elementwise, and print output.
+// Do ++nFail_ if two components are not equal within a given tolerance.
+// The function is converted from PEP-485
+template<class Type>
+typename std::enable_if
+<
+    !std::is_same<floatScalar, Type>::value &&
+    !std::is_same<doubleScalar, Type>::value &&
+    !std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,
+    const scalar absTol = 0
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    for (label i = 0; i < pTraits<Type>::nComponents; ++i)
+    {
+        if (max(absTol, relTol*max(mag(x[i]), mag(y[i]))) < mag(x[i] - y[i]))
+        {
+            ++nFail;
+        }
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Create each constructor of DiagTensor<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+    {
+        Info<< "# Construct initialized to zero:" << nl;
+        const DiagTensor<Type> dT(Zero);
+        Info<< dT << endl;
+    }
+    {
+        Info<< "# Construct given VectorSpace of the same rank:" << nl;
+        const VectorSpace<DiagTensor<Type>, Type, 3> V(Zero);
+        const DiagTensor<Type> dT(V);
+        Info<< dT << endl;
+    }
+    {
+        Info<< "# Construct given the three components:" << nl;
+        const DiagTensor<Type> dT
+        (
+            Type(1),
+                    Type(5),
+                            Type(-9)
+        );
+        Info<< dT << endl;
+    }
+    {
+        Info<< "# Copy construct:" << nl;
+        const DiagTensor<Type> dT(Zero);
+        const DiagTensor<Type> copydT(dT);
+        Info<< dT << tab << copydT << endl;
+    }
+}
+
+
+// Execute each member function of DiagTensor<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+    DiagTensor<Type> dT(Type(1), Type(5), Type(-9));
+    const DiagTensor<Type> cdT(Type(-9), Type(5), Type(1));
+
+    Info<< "# Operand: " << nl
+        << "  DiagTensor = " << dT << endl;
+
+
+    {
+        Info<< "# Component access:" << nl;
+
+        DiagTensor<Type> cpdT(dT.xx(), dT.yy(), dT.zz());
+        cmp("  'DiagTensor' access:", dT, cpdT);
+
+        const DiagTensor<Type> cpcdT(cdT.xx(), cdT.yy(), cdT.zz());
+        cmp("  'const DiagTensor' access:", cdT, cpcdT);
+    }
+}
+
+
+// Execute each global function of DiagTensor<Type>, and print output
+template<class Type>
+void test_global_funcs(Type)
+{
+    const Tensor<Type> T
+    (
+        Type(-1), Type(2), Type(-3),
+        Type(4),  Type(5), Type(-6),
+        Type(7),  Type(8), Type(-9)
+    );
+    const SymmTensor<Type> sT
+    (
+        Type(-1), Type(2), Type(-3),
+                  Type(5), Type(-6),
+                           Type(-9)
+    );
+    const DiagTensor<Type> dT(Type(1), Type(5), Type(-9));
+
+    Info<< "# Operands: " << nl
+        << "  Tensor = " << T << nl
+        << "  SymmTensor = " << sT << nl
+        << "  DiagTensor = " << dT << endl;
+
+
+    cmp("  Trace = ", tr(dT), Type(-3));
+    cmp("  Spherical part = ", sph(dT), SphericalTensor<Type>(tr(dT)/Type(3)));
+    cmp("  Determinant = ", det(dT), Type(-44.99999999999999));
+    cmp
+    (
+        "  Inverse = ",
+        inv(dT),
+        DiagTensor<Type>(Type(1), Type(0.2), Type(-0.11111111))
+    );
+    cmp
+    (
+        "  Diagonal of Tensor = ",
+        diag(T),
+        DiagTensor<Type>(Type(-1), Type(5), Type(-9))
+    );
+    cmp
+    (
+        "  Diagonal of SymmTensor = ",
+        diag(sT),
+        DiagTensor<Type>(Type(-1), Type(5), Type(-9))
+    );
+}
+
+
+// Execute each global operator of DiagTensor<Type>, and print output
+template<class Type>
+void test_global_opers(Type)
+{
+    const Tensor<Type> T
+    (
+        Type(-1), Type(2), Type(-3),
+        Type(4),  Type(5), Type(-6),
+        Type(7),  Type(8), Type(-9)
+    );
+    const SymmTensor<Type> sT
+    (
+        Type(-1), Type(2), Type(-3),
+                  Type(5), Type(-6),
+                           Type(-9)
+    );
+    const DiagTensor<Type> dT(Type(1), Type(5), Type(-9));
+    const SphericalTensor<Type> spT(Type(1));
+    const Vector<Type> v(Type(3), Type(2), Type(1));
+    const Type x(4);
+
+    Info<< "# Operands:" << nl
+        << "  Tensor = " << T << nl
+        << "  SymmTensor = " << sT << nl
+        << "  DiagTensor = " << dT << nl
+        << "  SphericalTensor = " << spT << nl
+        << "  Vector = " << v << nl
+        << "  Type = " << x << endl;
+
+
+    cmp
+    (
+        "  Sum of DiagTensor-Tensor = ",
+        (dT + T),
+        Tensor<Type>
+        (
+            Type(0), Type(2),  Type(-3),
+            Type(4), Type(10), Type(-6),
+            Type(7), Type(8),  Type(-18)
+        )
+    );
+    cmp
+    (
+        "  Sum of Tensor-DiagTensor = ",
+        (T + dT),
+        Tensor<Type>
+        (
+            Type(0), Type(2),  Type(-3),
+            Type(4), Type(10), Type(-6),
+            Type(7), Type(8),  Type(-18)
+        )
+    );
+    cmp
+    (
+        "  Subtract Tensor from DiagTensor = ",
+        (dT - T),
+        Tensor<Type>
+        (
+            Type(2),  Type(-2), Type(3),
+            Type(-4), Type(0),  Type(6),
+            Type(-7), Type(-8), Type(0)
+        )
+    );
+    cmp
+    (
+        "  Subtract DiagTensor from Tensor = ",
+        (T - dT),
+        Tensor<Type>
+        (
+            Type(-2), Type(2), Type(-3),
+            Type(4),  Type(0), Type(-6),
+            Type(7),  Type(8), Type(0)
+        )
+    );
+    cmp
+    (
+        "  Division of Type by DiagTensor = ",
+        (x/dT),
+        DiagTensor<Type>(Type(4), Type(0.8), Type(-0.44444444))
+    );
+    cmp
+    (
+        "  Division of DiagTensor by Type = ",
+        (dT/x),
+        DiagTensor<Type>(Type(0.25), Type(1.25), Type(-2.25))
+    );
+    cmp
+    (
+        "  Division of Vector by DiagTensor = ",
+        (v/dT),
+        Vector<Type>(Type(3), Type(0.4), Type(-0.11111111))
+    );
+    cmp
+    (
+        "  Inner-product of DiagTensor-DiagTensor = ",
+        (dT & dT),
+        DiagTensor<Type>(Type(1), Type(25), Type(81))
+    );
+    cmp
+    (
+        "  Inner-product of DiagTensor-Tensor = ",
+        (dT & T),
+        Tensor<Type>
+        (
+            Type(-1),  Type(2),   Type(-3),
+            Type(20),  Type(25),  Type(-30),
+            Type(-63), Type(-72), Type(81)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of Tensor-DiagTensor = ",
+        (T & dT),
+        Tensor<Type>
+        (
+            Type(-1), Type(10), Type(27),
+            Type(4),  Type(25), Type(54),
+            Type(7),  Type(40), Type(81)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of DiagTensor-Vector = ",
+        (dT & v),
+        Vector<Type>(Type(3), Type(10), Type(-9))
+    );
+    cmp
+    (
+        "  Inner-product of Vector-DiagTensor = ",
+        (v & dT),
+        Vector<Type>(Type(3), Type(10), Type(-9))
+    );
+}
+
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
+    test_member_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global functions: "<< typeID[I] << " ##" << nl;
+    test_global_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global operators: "<< typeID[I] <<" ##" << nl;
+    test_global_opers(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main()
+{
+    const std::tuple<floatScalar, doubleScalar, complex> types
+    (
+        std::make_tuple(Zero, Zero, Zero)
+    );
+
+    const List<word> typeID
+    ({
+        "DiagTensor<floatScalar>",
+        "DiagTensor<doubleScalar>",
+        "DiagTensor<complex>"
+    });
+
+    run_tests(types, typeID);
+
+
+    if (nFail_)
+    {
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
+    }
+
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
+    return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/test/SphericalTensor/Make/files

@@ -0,0 +1,3 @@
+Test-SphericalTensor.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-SphericalTensor

applications/test/SphericalTensor/Make/options

@@ -0,0 +1,2 @@
+/* EXE_INC = */
+/* EXE_LIBS = */

applications/test/SphericalTensor/Test-SphericalTensor.C

@@ -0,0 +1,348 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    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.
+
+    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
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-SphericalTensor
+
+Description
+    Tests for \c SphericalTensor constructors, member functions and operators
+    using \c floatScalar, \c doubleScalar, and \c complex base types.
+
+    Cross-checks were obtained from 'NumPy 1.15.1' and 'SciPy 1.1.0' if no
+    theoretical cross-check exists (like eigendecomposition relations), and
+    were hard-coded for elementwise comparisons.
+
+    For \c complex base type, the cross-checks do only involve zero imag part.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Tensor.H"
+#include "SymmTensor.H"
+#include "SphericalTensor.H"
+#include "DiagTensor.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Total number of unit tests
+unsigned nTest_ = 0;
+
+
+// Total number of failed unit tests
+unsigned nFail_ = 0;
+
+
+// Compare two floating point types, and print output.
+// Do ++nFail_ if values of two objects are not equal within a given tolerance.
+// The function is converted from PEP-485.
+template<class Type>
+typename std::enable_if
+<
+    std::is_same<floatScalar, Type>::value ||
+    std::is_same<doubleScalar, Type>::value ||
+    std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,  //<! are values the same within 8 decimals
+    const scalar absTol = 0      //<! useful for cmps near zero
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
+    {
+        ++nFail;
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Compare two containers elementwise, and print output.
+// Do ++nFail_ if two components are not equal within a given tolerance.
+// The function is converted from PEP-485
+template<class Type>
+typename std::enable_if
+<
+    !std::is_same<floatScalar, Type>::value &&
+    !std::is_same<doubleScalar, Type>::value &&
+    !std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,
+    const scalar absTol = 0
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    for (label i = 0; i < pTraits<Type>::nComponents; ++i)
+    {
+        if (max(absTol, relTol*max(mag(x[i]), mag(y[i]))) < mag(x[i] - y[i]))
+        {
+            ++nFail;
+        }
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Create each constructor of SphericalTensor<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+    {
+        Info<< "# Construct initialized to zero:" << nl;
+        const SphericalTensor<Type> spT(Zero);
+        Info<< spT << endl;
+    }
+    {
+        Info<< "# Construct given VectorSpace of the same rank:" << nl;
+        const VectorSpace<SphericalTensor<Type>, Type, 1> V(Zero);
+        const SphericalTensor<Type> spT(V);
+        Info<< spT << endl;
+    }
+    {
+        Info<< "# Construct given the component:" << nl;
+        const SphericalTensor<Type> spT(Type(1));
+        Info<< spT << endl;
+    }
+    {
+        Info<< "# Copy construct:" << nl;
+        const SphericalTensor<Type> spT(Zero);
+        const SphericalTensor<Type> copyspT(spT);
+        Info<< spT << tab << copyspT << endl;
+    }
+}
+
+
+// Execute each member function of SphericalTensor<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+    SphericalTensor<Type> spT(Type(1));
+    const SphericalTensor<Type> cspT(Type(-9));
+
+    Info<< "# Operand: " << nl
+        << "  SphericalTensor = " << spT << endl;
+
+
+    {
+        Info<< "# Component access:" << nl;
+
+        SphericalTensor<Type> cpspT(spT.ii());
+        cmp("  'SphericalTensor' access:", spT, cpspT);
+
+        const SphericalTensor<Type> cpcspT(cspT.ii());
+        cmp("  'const SphericalTensor' access:", cspT, cpcspT);
+    }
+    {
+        Info<< "# SphericalTensor operations:" << nl;
+
+        Info<< "  Transpose:" << nl;
+        cmp("  'SphericalTensor'.T():", spT.T(), spT);
+    }
+}
+
+
+// Execute each global function of SphericalTensor<Type>, and print output
+template<class Type>
+void test_global_funcs(Type)
+{
+    const SphericalTensor<Type> spT(Type(5));
+
+    Info<< "# Operand: " << nl
+        << "  SphericalTensor = " << spT << endl;
+
+
+    cmp("  Trace = ", tr(spT), Type(15));
+    cmp("  Spherical part = ", sph(spT), spT);
+    cmp("  Determinant = ", det(spT), Type(124.99999999999994));
+    cmp
+    (
+        "  Inverse = ",
+        inv(spT),
+        SphericalTensor<Type>(Type(0.2))
+    );
+    cmp("  Square of Frobenius norm = ", magSqr(spT), Type(75));
+    cmp("  Max component = ", cmptMax(spT), Type(5));
+    cmp("  Min component = ", cmptMax(spT), Type(5));
+    cmp("  Sum of components = ", cmptSum(spT), Type(15));
+    cmp("  Arithmetic average of components = ", cmptAv(spT), Type(5));
+}
+
+
+// Execute each global operator of SphericalTensor<Type>, and print output
+template<class Type>
+void test_global_opers(Type)
+{
+    const Tensor<Type> T
+    (
+        Type(-1), Type(2), Type(-3),
+        Type(4),  Type(5), Type(-6),
+        Type(7),  Type(8), Type(-9)
+    );
+    const SymmTensor<Type> sT
+    (
+        Type(-1), Type(2), Type(-3),
+                  Type(5), Type(-6),
+                           Type(-9)
+    );
+    const DiagTensor<Type> dT(Type(1), Type(5), Type(-9));
+    const SphericalTensor<Type> spT(Type(-2));
+    const Vector<Type> v(Type(3), Type(2), Type(1));
+    const Type x(4);
+
+    Info<< "# Operands:" << nl
+        << "  Tensor = " << T << nl
+        << "  SymmTensor = " << sT << nl
+        << "  DiagTensor = " << dT << nl
+        << "  SphericalTensor = " << spT << nl
+        << "  Vector = " << v << nl
+        << "  Type = " << x << endl;
+
+
+    cmp
+    (
+        "  Division of Type by SpTensor = ",
+        (x/spT),
+        SphericalTensor<Type>(Type(-2))
+    );
+    cmp
+    (
+        "  Division of SpTensor by Type = ",
+        (spT/x),
+        SphericalTensor<Type>(Type(-0.5))
+    );
+    cmp
+    (
+        "  Inner-product of SpTensor-SpTensor = ",
+        (spT & spT),
+        SphericalTensor<Type>(Type(4))
+    );
+    cmp
+    (
+        "  Inner-product of SpTensor-Vector = ",
+        (spT & v),
+        Vector<Type>(Type(-6), Type(-4), Type(-2)) // Column-vector
+    );
+    cmp
+    (
+        "  Inner-product of Vector-SpTensor = ",
+        (v & spT),
+        Vector<Type>(Type(-6), Type(-4), Type(-2)) // Row-vector
+    );
+    cmp("  D-inner-product of SpTensor-SpTensor = ", (spT && spT), Type(12));
+}
+
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
+    test_member_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global functions: "<< typeID[I] << " ##" << nl;
+    test_global_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global operators: "<< typeID[I] <<" ##" << nl;
+    test_global_opers(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main()
+{
+    const std::tuple<floatScalar, doubleScalar, complex> types
+    (
+        std::make_tuple(Zero, Zero, Zero)
+    );
+
+    const List<word> typeID
+    ({
+        "SphericalTensor<floatScalar>",
+        "SphericalTensor<doubleScalar>",
+        "SphericalTensor<complex>"
+    });
+
+    run_tests(types, typeID);
+
+
+    if (nFail_)
+    {
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
+    }
+
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
+    return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/test/SphericalTensor2D/Make/files

@@ -0,0 +1,3 @@
+Test-SphericalTensor2D.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-SphericalTensor2D

applications/test/SphericalTensor2D/Make/options

@@ -0,0 +1,2 @@
+/* EXE_INC = */
+/* EXE_LIBS = */

applications/test/SphericalTensor2D/Test-SphericalTensor2D.C

@@ -0,0 +1,331 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    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.
+
+    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
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-SphericalTensor2D
+
+Description
+    Tests for \c SphericalTensor2D constructors, member functions and operators
+    using \c floatScalar, \c doubleScalar, and \c complex base types.
+
+    Cross-checks were obtained from 'NumPy 1.15.1' and 'SciPy 1.1.0' if no
+    theoretical cross-check exists (like eigendecomposition relations), and
+    were hard-coded for elementwise comparisons.
+
+    For \c complex base type, the cross-checks do only involve zero imag part.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Tensor2D.H"
+#include "SymmTensor2D.H"
+#include "SphericalTensor2D.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Total number of unit tests
+unsigned nTest_ = 0;
+
+
+// Total number of failed unit tests
+unsigned nFail_ = 0;
+
+
+// Compare two floating point types, and print output.
+// Do ++nFail_ if values of two objects are not equal within a given tolerance.
+// The function is converted from PEP-485.
+template<class Type>
+typename std::enable_if
+<
+    std::is_same<floatScalar, Type>::value ||
+    std::is_same<doubleScalar, Type>::value ||
+    std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,  //<! are values the same within 8 decimals
+    const scalar absTol = 0      //<! useful for cmps near zero
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
+    {
+        ++nFail;
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Compare two containers elementwise, and print output.
+// Do ++nFail_ if two components are not equal within a given tolerance.
+// The function is converted from PEP-485
+template<class Type>
+typename std::enable_if
+<
+    !std::is_same<floatScalar, Type>::value &&
+    !std::is_same<doubleScalar, Type>::value &&
+    !std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,
+    const scalar absTol = 0
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    for (label i = 0; i < pTraits<Type>::nComponents; ++i)
+    {
+        if (max(absTol, relTol*max(mag(x[i]), mag(y[i]))) < mag(x[i] - y[i]))
+        {
+            ++nFail;
+        }
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Create each constructor of SphericalTensor2D<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+    {
+        Info<< "# Construct initialized to zero:" << nl;
+        const SphericalTensor2D<Type> spT(Zero);
+        Info<< spT << endl;
+    }
+    {
+        Info<< "# Construct given VectorSpace of the same rank:" << nl;
+        const VectorSpace<SphericalTensor2D<Type>, Type, 1> V(Zero);
+        const SphericalTensor2D<Type> spT(V);
+        Info<< spT << endl;
+    }
+    {
+        Info<< "# Construct given the component:" << nl;
+        const SphericalTensor2D<Type> spT(Type(1));
+        Info<< spT << endl;
+    }
+    {
+        Info<< "# Copy construct:" << nl;
+        const SphericalTensor2D<Type> spT(Zero);
+        const SphericalTensor2D<Type> copyspT(spT);
+        Info<< spT << tab << copyspT << endl;
+    }
+}
+
+
+// Execute each member function of SphericalTensor2D<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+    SphericalTensor2D<Type> spT(Type(1));
+    const SphericalTensor2D<Type> cspT(Type(-9));
+
+    Info<< "# Operand: " << nl
+        << "  SphericalTensor2D = " << spT << endl;
+
+
+    {
+        Info<< "# Component access:" << nl;
+
+        SphericalTensor2D<Type> cpspT(spT.ii());
+        cmp("  'SphericalTensor2D' access:", spT, cpspT);
+
+        const SphericalTensor2D<Type> cpcspT(cspT.ii());
+        cmp("  'const SphericalTensor2D' access:", cspT, cpcspT);
+    }
+}
+
+
+// Execute each global function of SphericalTensor2D<Type>, and print output
+template<class Type>
+void test_global_funcs(Type)
+{
+    const SphericalTensor2D<Type> spT(Type(5));
+
+    Info<< "# Operand: " << nl
+        << "  SphericalTensor2D = " << spT << endl;
+
+
+    cmp("  Trace = ", tr(spT), Type(10));
+    cmp("  Spherical part = ", sph(spT), spT);
+    cmp("  Determinant = ", det(spT), Type(24.99999999999994));
+    cmp
+    (
+        "  Inverse = ",
+        inv(spT),
+        SphericalTensor2D<Type>(Type(0.2))
+    );
+}
+
+
+// Execute each global operator of SphericalTensor2D<Type>, and print output
+template<class Type>
+void test_global_opers(Type)
+{
+    const Tensor2D<Type> T
+    (
+        Type(-1), Type(2),
+        Type(4),  Type(5)
+    );
+    const SymmTensor2D<Type> sT
+    (
+        Type(-1), Type(2),
+                  Type(5)
+    );
+    const SphericalTensor2D<Type> spT(Type(-2));
+    const Vector2D<Type> v(Type(3), Type(2));
+    const Type x(4);
+
+    Info<< "# Operands:" << nl
+        << "  Tensor2D = " << T << nl
+        << "  SymmTensor2D = " << sT << nl
+        << "  SphericalTensor2D = " << spT << nl
+        << "  Vector2D = " << v << nl
+        << "  Type = " << x << endl;
+
+
+    cmp
+    (
+        "  Division of Type by SpTensor2D = ",
+        (x/spT),
+        SphericalTensor2D<Type>(Type(-2))
+    );
+    cmp
+    (
+        "  Division of SpTensor2D by Type = ",
+        (spT/x),
+        SphericalTensor2D<Type>(Type(-0.5))
+    );
+    cmp
+    (
+        "  Inner-product of SpTensor2D-SpTensor2D = ",
+        (spT & spT),
+        SphericalTensor2D<Type>(Type(4))
+    );
+    cmp
+    (
+        "  Inner-product of SpTensor2D-Vector2D = ",
+        (spT & v),
+        Vector2D<Type>(Type(-6), Type(-4)) // Column-vector
+    );
+    cmp
+    (
+        "  Inner-product of Vector2D-SpTensor2D = ",
+        (v & spT),
+        Vector2D<Type>(Type(-6), Type(-4)) // Row-vector
+    );
+}
+
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
+    test_member_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global functions: "<< typeID[I] << " ##" << nl;
+    test_global_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global operators: "<< typeID[I] <<" ##" << nl;
+    test_global_opers(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main()
+{
+    const std::tuple<floatScalar, doubleScalar, complex> types
+    (
+        std::make_tuple(Zero, Zero, Zero)
+    );
+
+    const List<word> typeID
+    ({
+        "SphericalTensor2D<floatScalar>",
+        "SphericalTensor2D<doubleScalar>",
+        "SphericalTensor2D<complex>"
+    });
+
+    run_tests(types, typeID);
+
+
+    if (nFail_)
+    {
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
+    }
+
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
+    return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/test/SymmTensor/Make/files

@@ -0,0 +1,3 @@
+Test-SymmTensor.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-SymmTensor

applications/test/SymmTensor/Make/options

@@ -0,0 +1,2 @@
+/* EXE_INC = */
+/* EXE_LIBS = */

applications/test/SymmTensor/Test-SymmTensor.C

@@ -0,0 +1,576 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    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.
+
+    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
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-SymmTensor
+
+Description
+    Tests for \c SymmTensor constructors, member functions and operators
+    using \c floatScalar, \c doubleScalar, and \c complex base types.
+
+    Eigen decomposition tests for \c symmTensor, i.e. SymmTensor<scalar>.
+
+    Cross-checks were obtained from 'NumPy 1.15.1' and 'SciPy 1.1.0' if no
+    theoretical cross-check exists (like eigendecomposition relations), and
+    were hard-coded for elementwise comparisons.
+
+    For \c complex base type, the cross-checks do only involve zero imag part.
+
+\*---------------------------------------------------------------------------*/
+
+#include "symmTensor.H"
+#include "transform.H"
+#include "Random.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Total number of unit tests
+unsigned nTest_ = 0;
+
+
+// Total number of failed unit tests
+unsigned nFail_ = 0;
+
+
+// Create a random symmTensor
+symmTensor makeRandomContainer(Random& rnd)
+{
+    symmTensor A(Zero);
+    std::generate(A.begin(), A.end(), [&]{ return rnd.GaussNormal<scalar>(); });
+    return A;
+}
+
+
+// Compare two floating point types, and print output.
+// Do ++nFail_ if values of two objects are not equal within a given tolerance.
+// The function is converted from PEP-485.
+template<class Type>
+typename std::enable_if
+<
+    std::is_same<floatScalar, Type>::value ||
+    std::is_same<doubleScalar, Type>::value ||
+    std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,  //<! are values the same within 8 decimals
+    const scalar absTol = 0      //<! useful for cmps near zero
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
+    {
+        ++nFail;
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Compare two containers elementwise, and print output.
+// Do ++nFail_ if two components are not equal within a given tolerance.
+// The function is converted from PEP-485
+template<class Type>
+typename std::enable_if
+<
+    !std::is_same<floatScalar, Type>::value &&
+    !std::is_same<doubleScalar, Type>::value &&
+    !std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,
+    const scalar absTol = 0
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    for (label i = 0; i < pTraits<Type>::nComponents; ++i)
+    {
+        if (max(absTol, relTol*max(mag(x[i]), mag(y[i]))) < mag(x[i] - y[i]))
+        {
+            ++nFail;
+        }
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Create each constructor of SymmTensor<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+    {
+        Info<< "# Construct initialized to zero:" << nl;
+        const SymmTensor<Type> sT(Zero);
+        Info<< sT << endl;
+    }
+    {
+        Info<< "# Construct given VectorSpace of the same rank:" << nl;
+        const VectorSpace<SymmTensor<Type>, Type, 6> M(Zero);
+        const SymmTensor<Type> sT(M);
+        Info<< sT << endl;
+    }
+    {
+        Info<< "# Construct given SphericalTensor:" << nl;
+        const SphericalTensor<Type> Sp(Type(5));
+        const SymmTensor<Type> sT(Sp);
+        Info<< sT << endl;
+    }
+    {
+        Info<< "# Construct given the six components:" << nl;
+        const SymmTensor<Type> sT
+        (
+            Type(1), Type(2), Type(-3),
+                     Type(5), Type(-6),
+                              Type(-9)
+        );
+        Info<< sT << endl;
+    }
+    {
+        Info<< "# Copy construct:" << nl;
+        const SymmTensor<Type> sT(Zero);
+        const SymmTensor<Type> copysT(sT);
+        Info<< sT << tab << copysT << endl;
+    }
+}
+
+
+// Execute each member function of SymmTensor<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+    SymmTensor<Type> sT
+    (
+        Type(1), Type(2), Type(-3),
+                 Type(5), Type(-6),
+                          Type(-9)
+    );
+    const SymmTensor<Type> csT
+    (
+        Type(1), Type(2), Type(-3),
+                 Type(5), Type(-6),
+                          Type(-9)
+    );
+
+    Info<< "# Operand: " << nl
+        << "  SymmTensor = " << sT << endl;
+
+
+    {
+        Info<< "# Component access:" << nl;
+
+        SymmTensor<Type> cpsT
+        (
+            sT.xx(), sT.xy(), sT.xz(),
+                     sT.yy(), sT.yz(),
+                              sT.zz()
+        );
+        cmp("  'SymmTensor' access:", sT, cpsT);
+        cmp("  xy()=yx():", sT.xy(), sT.yx());
+        cmp("  xz()=zx():", sT.xz(), sT.zx());
+        cmp("  yz()=zy():", sT.yz(), sT.zy());
+
+        const SymmTensor<Type> cpcsT
+        (
+            csT.xx(), csT.xy(), csT.xz(),
+                      csT.yy(), csT.yz(),
+                                csT.zz()
+        );
+        cmp("  'const SymmTensor' access:", csT, cpcsT);
+        cmp("  xy()=yx():", sT.xy(), sT.yx());
+        cmp("  xz()=zx():", sT.xz(), sT.zx());
+        cmp("  yz()=zy():", sT.yz(), sT.zy());
+    }
+    {
+        Info<< "# Diagonal access:" << nl;
+
+        cmp
+        (
+            "  'SymmTensor'.diag():",
+            sT.diag(),
+            Vector<Type>(Type(1), Type(5), Type(-9))
+        );
+        cmp
+        (
+            "  'const SymmTensor'.diag():",
+            csT.diag(),
+            Vector<Type>(Type(1), Type(5), Type(-9))
+        );
+
+
+        Info<< "# Diagonal manipulation:" << nl;
+
+        sT.diag(Vector<Type>(Type(-10), Type(-15), Type(-20)));
+        cmp
+        (
+            "  'SymmTensor'.diag('Vector'):",
+            sT.diag(),
+            Vector<Type>(Type(-10), Type(-15), Type(-20))
+        );
+    }
+    {
+        Info<< "# Tensor operations:" << nl;
+
+        Info<< "  Transpose:" << nl;
+        cmp("  'SymmTensor'.T():", sT.T(), sT);
+    }
+    {
+        Info<< "# Member operators:" << nl;
+
+        sT = SphericalTensor<Type>(Type(5));
+        cmp
+        (
+            "  Assign to a SphericalTensor:",
+            sT,
+            SymmTensor<Type>
+            (
+                Type(5),  Zero,    Zero,
+                          Type(5), Zero,
+                                   Type(5)
+            )
+        );
+    }
+}
+
+
+// Execute each global function of SymmTensor<Type>, and print output
+template<class Type>
+void test_global_funcs(Type)
+{
+    const SymmTensor<Type> sT
+    (
+        Type(1), Type(2), Type(-3),
+                 Type(5), Type(-6),
+                          Type(-9)
+    );
+
+    Info<< "# Operand: " << nl
+        << "  SymmTensor = " << sT << nl << endl;
+
+
+    cmp("  Trace = ", tr(sT), Type(-3));
+    cmp("  Spherical part = ", sph(sT), SphericalTensor<Type>(tr(sT)/Type(3)));
+    cmp("  Symmetric part = ", symm(sT), sT);
+    cmp("  Twice the symmetric part = ", twoSymm(sT), 2*sT);
+    cmp
+    (
+        "  Deviatoric part = ",
+        dev(sT),
+        SymmTensor<Type>
+        (
+            Type(2), Type(2), Type(-3),
+                     Type(6), Type(-6),
+                              Type(-8)
+        )
+    );
+    cmp("  Two-third deviatoric part = ", dev2(sT), sT - 2*sph(sT));
+    cmp("  Determinant = ", det(sT), Type(-17.999999999999996));
+    cmp
+    (
+        "  Cofactor tensor = ",
+        cof(sT),
+        SymmTensor<Type>
+        (
+            Type(-81), Type(36),  Type(3),
+                       Type(-18), Type(0),
+                                  Type(1)
+        )
+    );
+    cmp
+    (
+        "  Inverse = ",
+        inv(sT, det(sT)),
+        SymmTensor<Type>
+        (
+            Type(4.5), Type(-2), Type(-0.16666667),
+                       Type(1),  Type(0),
+                                 Type(-0.05555556)
+        ),
+        1e-8,
+        1e-8
+    );
+    cmp
+    (
+        "  Inverse (another) = ",
+        inv(sT),
+        SymmTensor<Type>
+        (
+            Type(4.5), Type(-2), Type(-0.16666667),
+                       Type(1),  Type(0),
+                                 Type(-0.05555556)
+        ),
+        1e-8,
+        1e-8
+    );
+    cmp("  First invariant = ", invariantI(sT), Type(-3));
+    cmp("  Second invariant = ", invariantII(sT), Type(-98));
+    cmp("  Third invariant = ", invariantIII(sT), Type(-17.999999999999996));
+    cmp
+    (
+        "  Inner-product with self = ",
+        innerSqr(sT),
+        SymmTensor<Type>
+        (
+            Type(14), Type(30), Type(12),
+                      Type(65), Type(18),
+                                Type(126)
+        )
+    );
+    cmp("  Square of Frobenius norm = ", magSqr(sT), Type(205));
+}
+
+
+// Execute each global operator of SymmTensor<Type>, and print output
+template<class Type>
+void test_global_opers(Type)
+{
+    const Tensor<Type> T
+    (
+        Type(1), Type(2), Type(-3),
+        Type(4), Type(5), Type(-6),
+        Type(7), Type(8), Type(-9)
+    );
+    const SymmTensor<Type> sT
+    (
+        Type(1), Type(2), Type(-3),
+                 Type(5), Type(-6),
+                          Type(-9)
+    );
+    const SphericalTensor<Type> spT(Type(1));
+    const Vector<Type> v(Type(3), Type(2), Type(1));
+    const Type x(4);
+
+    Info<< "# Operands:" << nl
+        << "  Tensor = " << T << nl
+        << "  SymmTensor = " << sT << nl
+        << "  SphericalTensor = " << spT << nl
+        << "  Vector = " << v << nl
+        << "  Type = " << x << endl;
+
+
+    cmp
+    (
+        "  Sum of SpTensor-SymmTensor = ",
+        (spT + sT),
+        SymmTensor<Type>
+        (
+            Type(2), Type(2), Type(-3),
+                     Type(6), Type(-6),
+                              Type(-8)
+        )
+    );
+    cmp
+    (
+        "  Sum of SymmTensor-SpTensor = ",
+        (sT + spT),
+        SymmTensor<Type>
+        (
+            Type(2), Type(2), Type(-3),
+                     Type(6), Type(-6),
+                              Type(-8)
+        )
+    );
+    cmp
+    (
+        "  Subtract SymmTensor from SpTensor = ",
+        (spT - sT),
+        SymmTensor<Type>
+        (
+            Type(0), Type(-2), Type(3),
+                     Type(-4), Type(6),
+                               Type(10)
+        )
+    );
+    cmp
+    (
+        "  Subtract SpTensor from SymmTensor = ",
+        (sT - spT),
+        SymmTensor<Type>
+        (
+            Type(0), Type(2), Type(-3),
+                     Type(4), Type(-6),
+                              Type(-10)
+        )
+    );
+    cmp
+    (
+        "  Hodge dual of a SymmTensor",
+        *sT,
+        Vector<Type>(Type(-6), Type(3), Type(2))
+    );
+    cmp
+    (
+        "  Division of a SymmTensor by a Type",
+        sT/x,
+        SymmTensor<Type>
+        (
+            Type(0.25), Type(0.5),  Type(-0.75),
+                        Type(1.25), Type(-1.5),
+                                    Type(-2.25)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SymmTensor-SymmTensor = ",
+        (sT & sT),
+        Tensor<Type>
+        (
+            Type(14), Type(30), Type(12),
+            Type(30), Type(65), Type(18),
+            Type(12), Type(18), Type(126)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SpTensor-SymmTensor = ",
+        (spT & sT),
+        SymmTensor<Type>
+        (
+            Type(1), Type(2), Type(-3),
+                     Type(5), Type(-6),
+                              Type(-9)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SymmTensor-SpTensor = ",
+        (sT & spT),
+        SymmTensor<Type>
+        (
+            Type(1), Type(2), Type(-3),
+                     Type(5), Type(-6),
+                              Type(-9)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SymmTensor-Vector = ",
+        (sT & v),
+        Vector<Type>(Type(4), Type(10), Type(-30)) // Column-vector
+    );
+    cmp
+    (
+        "  Inner-product of Vector-SymmTensor = ",
+        (v & sT),
+        Vector<Type>(Type(4), Type(10), Type(-30)) // Row-vector
+    );
+    cmp("  D-inner-product of SymmTensor-SymmTensor = ", (sT && sT), Type(205));
+    cmp("  D-inner-product of SymmTensor-SpTensor = ", (sT && spT), Type(-3));
+    cmp("  D-inner-product of SpTensor-SymmTensor = ", (spT && sT), Type(-3));
+}
+
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
+    test_member_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global functions: "<< typeID[I] << " ##" << nl;
+    test_global_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global operators: "<< typeID[I] <<" ##" << nl;
+    test_global_opers(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main()
+{
+    const std::tuple<floatScalar, doubleScalar, complex> types
+    (
+        std::make_tuple(Zero, Zero, Zero)
+    );
+
+    const List<word> typeID
+    ({
+        "SymmTensor<floatScalar>",
+        "SymmTensor<doubleScalar>",
+        "SymmTensor<complex>"
+    });
+
+    run_tests(types, typeID);
+
+
+    if (nFail_)
+    {
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
+    }
+
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
+    return 0;
+
+    return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/test/SymmTensor2D/Make/files

@@ -0,0 +1,3 @@
+Test-SymmTensor2D.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-SymmTensor2D

applications/test/SymmTensor2D/Make/options

@@ -0,0 +1,2 @@
+/* EXE_INC = */
+/* EXE_LIBS = */

applications/test/SymmTensor2D/Test-SymmTensor2D.C

@@ -0,0 +1,538 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    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.
+
+    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
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-SymmTensor2D
+
+Description
+    Tests for \c SymmTensor2D constructors, member functions and operators
+    using \c floatScalar, \c doubleScalar, and \c complex base types.
+
+    Eigen decomposition tests for \c symmTensor2D, i.e. SymmTensor2D<scalar>.
+
+    Cross-checks were obtained from 'NumPy 1.15.1' and 'SciPy 1.1.0' if no
+    theoretical cross-check exists (like eigendecomposition relations), and
+    were hard-coded for elementwise comparisons.
+
+    For \c complex base type, the cross-checks do only involve zero imag part.
+
+\*---------------------------------------------------------------------------*/
+
+#include "symmTensor2D.H"
+#include "transform.H"
+#include "Random.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Total number of unit tests
+unsigned nTest_ = 0;
+
+
+// Total number of failed unit tests
+unsigned nFail_ = 0;
+
+
+// Create a random symmTensor2D
+symmTensor2D makeRandomContainer(Random& rnd)
+{
+    symmTensor2D A(Zero);
+    std::generate(A.begin(), A.end(), [&]{ return rnd.GaussNormal<scalar>(); });
+    return A;
+}
+
+
+// Compare two floating point types, and print output.
+// Do ++nFail_ if values of two objects are not equal within a given tolerance.
+// The function is converted from PEP-485.
+template<class Type>
+typename std::enable_if
+<
+    std::is_same<floatScalar, Type>::value ||
+    std::is_same<doubleScalar, Type>::value ||
+    std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,  //<! are values the same within 8 decimals
+    const scalar absTol = 0      //<! useful for cmps near zero
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
+    {
+        ++nFail;
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Compare two containers elementwise, and print output.
+// Do ++nFail_ if two components are not equal within a given tolerance.
+// The function is converted from PEP-485
+template<class Type>
+typename std::enable_if
+<
+    !std::is_same<floatScalar, Type>::value &&
+    !std::is_same<doubleScalar, Type>::value &&
+    !std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,
+    const scalar absTol = 0
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    for (label i = 0; i < pTraits<Type>::nComponents; ++i)
+    {
+        if (max(absTol, relTol*max(mag(x[i]), mag(y[i]))) < mag(x[i] - y[i]))
+        {
+            ++nFail;
+        }
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Create each constructor of SymmTensor2D<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+    {
+        Info<< "# Construct initialized to zero:" << nl;
+        const SymmTensor2D<Type> sT(Zero);
+        Info<< sT << endl;
+    }
+    {
+        Info<< "# Construct given VectorSpace:" << nl;
+        const VectorSpace<SymmTensor2D<Type>, Type, 3> V(Zero);
+        const SymmTensor2D<Type> sT(V);
+        Info<< sT << endl;
+    }
+    {
+        Info<< "# Construct given SphericalTensor2D:" << nl;
+        const SphericalTensor2D<Type> Sp(Type(5));
+        const SymmTensor2D<Type> sT(Sp);
+        Info<< sT << endl;
+    }
+    {
+        Info<< "# Construct given the three components:" << nl;
+        const SymmTensor2D<Type> sT
+        (
+            Type(1), Type(2),
+                     Type(4)
+        );
+        Info<< sT << endl;
+    }
+    {
+        Info<< "# Copy construct:" << nl;
+        const SymmTensor2D<Type> S(Type(1), Type(2), Type(3));
+        const SymmTensor2D<Type> sT(S);
+        Info<< sT << endl;
+    }
+}
+
+
+// Execute each member function of SymmTensor2D<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+    SymmTensor2D<Type> sT(Type(1), Type(2), Type(-3));
+    const SymmTensor2D<Type> csT(Type(-3), Type(2), Type(1));
+
+
+    {
+        Info<< "# Component access:" << nl;
+
+        SymmTensor2D<Type> cpsT
+        (
+            sT.xx(), sT.xy(),
+                     sT.yy()
+        );
+        cmp("  'SymmTensor2D' access:", sT, cpsT);
+        cmp("  xy()=yx():", sT.xy(), sT.yx());
+
+        const SymmTensor2D<Type> cpcsT
+        (
+            csT.xx(), csT.xy(),
+                      csT.yy()
+        );
+        cmp("  'const SymmTensor2D' access:", csT, cpcsT);
+        cmp("  xy()=yx():", sT.xy(), sT.yx());
+    }
+    {
+        Info<< "# Diagonal access:" << nl;
+
+        cmp
+        (
+            "  'SymmTensor2D'.diag():",
+            sT.diag(),
+            Vector2D<Type>(Type(1), Type(-3))
+        );
+        cmp
+        (
+            "  'const SymmTensor2D'.diag():",
+            csT.diag(),
+            Vector2D<Type>(Type(-3), Type(1))
+        );
+
+
+        Info<< "# Diagonal manipulation:" << nl;
+
+        sT.diag(Vector2D<Type>(Type(-10), Type(-15)));
+        cmp
+        (
+            "  'SymmTensor2D'.diag('Vector2D'):",
+            sT.diag(),
+            Vector2D<Type>(Type(-10), Type(-15))
+        );
+    }
+    {
+        Info<< "# Tensor operations:" << nl;
+
+        Info<< "  Transpose:" << nl;
+        cmp("  'SymmTensor2D'.T():", sT.T(), sT);
+    }
+    {
+        Info<< "# Member operators:" << nl;
+
+        sT = SphericalTensor2D<Type>(Type(5));
+        cmp
+        (
+            "  Assign to a SphericalTensor2D:",
+            sT,
+            SymmTensor2D<Type>
+            (
+                Type(5),  Zero,
+                          Type(5)
+            )
+        );
+    }
+}
+
+
+// Execute each global function of SymmTensor2D<Type>, and print output
+template<class Type>
+void test_global_funcs(Type)
+{
+    const SymmTensor2D<Type> sT(Type(1), Type(2), Type(-3));
+    const Vector2D<Type> v(Type(-3), Type(2));
+
+    Info<< "# Operands: " << nl
+        << "  SymmTensor2D<Type> = " << sT << nl
+        << "  Vector2D<Type> = " << v << endl;
+
+
+    cmp("  Trace = ", tr(sT), Type(-2));
+    cmp("  Spherical part = ", sph(sT), SphericalTensor2D<Type>(tr(sT)/Type(2)));
+    cmp("  Symmetric part = ", symm(sT), sT);
+    cmp("  Twice the symmetric part = ", twoSymm(sT), 2*sT);
+    cmp
+    (
+        "  Deviatoric part = ",
+        dev(sT),
+        SymmTensor2D<Type>
+        (
+            Type(2), Type(2),
+                     Type(-2)
+        )
+    );
+    cmp("  Two-third deviatoric part = ", dev2(sT), sT - 2*sph(sT));
+    cmp("  Determinant = ", det(sT), Type(-7.000000000000001));
+    cmp
+    (
+        "  Cofactor tensor = ",
+        cof(sT),
+        SymmTensor2D<Type>
+        (
+            Type(-3), Type(-2),
+                      Type(1)
+        )
+    );
+    cmp
+    (
+        "  Inverse = ",
+        inv(sT, det(sT)),
+        SymmTensor2D<Type>
+        (
+            Type(0.42857143), Type(0.28571429),
+                              Type(-0.14285714)
+        ),
+        1e-8,
+        1e-6
+    );
+    cmp
+    (
+        "  Inverse (another) = ",
+        inv(sT),
+        SymmTensor2D<Type>
+        (
+            Type(0.42857143), Type(0.28571429),
+                              Type(-0.14285714)
+        ),
+        1e-8,
+        1e-6
+    );
+    cmp("  First invariant = ", invariantI(sT), Type(-2));
+    cmp("  Second invariant = ", invariantII(sT), Type(-7));
+    cmp
+    (
+        "  Inner-product with self = ",
+        innerSqr(sT),
+        SymmTensor2D<Type>
+        (
+            Type(5), Type(-4),
+                     Type(13)
+        )
+    );
+    cmp("  Square of Frobenius norm = ", magSqr(sT), Type(17.999999999999996));
+    cmp
+    (
+        "  Outer-product of a Vector2D with itself = ",
+        sqr(v),
+        SymmTensor2D<Type>
+        (
+            Type(9), Type(-6),
+                     Type(4)
+        )
+    );
+}
+
+
+// Execute each global operator of SymmTensor2D<Type>, and print output
+template<class Type>
+void test_global_opers(Type)
+{
+    const Tensor2D<Type> T
+    (
+        Type(1), Type(-2),
+        Type(3), Type(-4)
+    );
+    const SymmTensor2D<Type> sT
+    (
+        Type(1), Type(2),
+                 Type(-4)
+    );
+    const SphericalTensor2D<Type> spT(Type(-1));
+    const Vector2D<Type> v(Type(3), Type(-2));
+    const Type x(-4);
+
+    Info<< "# Operands:" << nl
+        << "  Tensor2D = " << T << nl
+        << "  SymmTensor2D = " << sT << nl
+        << "  SphericalTensor2D = " << spT << nl
+        << "  Vector2D = " << v << nl
+        << "  Type = " << x << endl;
+
+
+    cmp
+    (
+        "  Sum of SpTensor2D-SymmTensor2D = ",
+        (spT + sT),
+        SymmTensor2D<Type>
+        (
+            Type(0), Type(2),
+                     Type(-5)
+        )
+    );
+    cmp
+    (
+        "  Sum of SymmTensor2D-SpTensor2D = ",
+        (sT + spT),
+        SymmTensor2D<Type>
+        (
+            Type(0), Type(2),
+                     Type(-5)
+        )
+    );
+    cmp
+    (
+        "  Subtract SymmTensor2D from SpTensor2D = ",
+        (spT - sT),
+        SymmTensor2D<Type>
+        (
+            Type(-2), Type(-2),
+                      Type(3)
+        )
+    );
+    cmp
+    (
+        "  Subtract SpTensor2D from SymmTensor2D = ",
+        (sT - spT),
+        SymmTensor2D<Type>
+        (
+            Type(2), Type(2),
+                     Type(-3)
+        )
+    );
+    cmp
+    (
+        "  Division of a SymmTensor2D by a Type",
+        sT/x,
+        SymmTensor2D<Type>
+        (
+            Type(-0.25), Type(-0.5),
+                         Type(1)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SymmTensor2D-SymmTensor2D = ",
+        (sT & sT),
+        Tensor2D<Type>
+        (
+            Type(5),  Type(-6),
+            Type(-6), Type(20)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SpTensor2D-SymmTensor2D = ",
+        (spT & sT),
+        SymmTensor2D<Type>
+        (
+            Type(-1), Type(-2),
+                      Type(4)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SymmTensor2D-SpTensor2D = ",
+        (sT & spT),
+        SymmTensor2D<Type>
+        (
+            Type(-1), Type(-2),
+                      Type(4)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SymmTensor2D-Vector2D = ",
+        (sT & v),
+        Vector2D<Type>(Type(-1), Type(14)) // Column-vector
+    );
+    cmp
+    (
+        "  Inner-product of Vector2D-SymmTensor2D = ",
+        (v & sT),
+        Vector2D<Type>(Type(-1), Type(14)) // Row-vector
+    );
+    cmp("  D-inner-prod of SymmTensor2D-SymmTensor2D = ", (sT && sT), Type(25));
+    cmp("  D-inner-prod of SymmTensor2D-SpTensor2D = ", (sT && spT), Type(3));
+    cmp("  D-inner-prod of SpTensor2D-SymmTensor2D = ", (spT && sT), Type(3));
+}
+
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
+    test_member_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global functions: "<< typeID[I] << " ##" << nl;
+    test_global_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global operators: "<< typeID[I] <<" ##" << nl;
+    test_global_opers(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    const std::tuple<floatScalar, doubleScalar, complex> types
+    (
+        std::make_tuple(Zero, Zero, Zero)
+    );
+
+    const List<word> typeID
+    ({
+        "SymmTensor2D<floatScalar>",
+        "SymmTensor2D<doubleScalar>",
+        "SymmTensor2D<complex>"
+    });
+
+    run_tests(types, typeID);
+
+
+    if (nFail_)
+    {
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
+    }
+
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/Tensor/Make/files

@@ -0,0 +1,3 @@
+Test-Tensor.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-Tensor

applications/test/Tensor/Make/options

@@ -0,0 +1,2 @@
+/* EXE_INC = */
+/* EXE_LIBS = */

applications/test/Tensor/Test-Tensor.C

@@ -0,0 +1,964 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2018 OpenFOAM Foundation
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    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.
+
+    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
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-Tensor
+
+Description
+    Tests for \c Tensor constructors, member functions and operators
+    using \c floatScalar, \c doubleScalar, and \c complex base types.
+
+    Eigen decomposition tests for \c tensor, i.e. Tensor<scalar>.
+
+    Cross-checks were obtained from 'NumPy 1.15.1' and 'SciPy 1.1.0' if no
+    theoretical cross-check exists (like eigendecomposition relations), and
+    were hard-coded for elementwise comparisons.
+
+    For \c complex base type, the cross-checks do only involve zero imag part.
+
+\*---------------------------------------------------------------------------*/
+
+#include "tensor.H"
+#include "transform.H"
+#include "Random.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Total number of unit tests
+unsigned nTest_ = 0;
+
+
+// Total number of failed unit tests
+unsigned nFail_ = 0;
+
+
+// Create a random tensor
+tensor makeRandomContainer(Random& rnd)
+{
+    tensor A(Zero);
+    std::generate(A.begin(), A.end(), [&]{ return rnd.GaussNormal<scalar>(); });
+    return A;
+}
+
+
+// Compare two floating point types, and print output.
+// Do ++nFail_ if values of two objects are not equal within a given tolerance.
+// The function is converted from PEP-485.
+template<class Type>
+typename std::enable_if
+<
+    std::is_same<floatScalar, Type>::value ||
+    std::is_same<doubleScalar, Type>::value ||
+    std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,  //<! are values the same within 8 decimals
+    const scalar absTol = 0      //<! useful for cmps near zero
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
+    {
+        ++nFail;
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Compare two containers elementwise, and print output.
+// Do ++nFail_ if two components are not equal within a given tolerance.
+// The function is converted from PEP-485
+template<class Type>
+typename std::enable_if
+<
+    !std::is_same<floatScalar, Type>::value &&
+    !std::is_same<doubleScalar, Type>::value &&
+    !std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,
+    const scalar absTol = 0
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    for (label i = 0; i < pTraits<Type>::nComponents; ++i)
+    {
+        if (max(absTol, relTol*max(mag(x[i]), mag(y[i]))) < mag(x[i] - y[i]))
+        {
+            ++nFail;
+        }
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Create each constructor of Tensor<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+    {
+        Info<< "# Construct initialized to zero:" << nl;
+        const Tensor<Type> T(Zero);
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Construct given MatrixSpace of the same rank:" << nl;
+        const MatrixSpace<Tensor<Type>, Type, 3, 3> M(Zero);
+        const Tensor<Type> T(M);
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Construct given VectorSpace of the same rank:" << nl;
+        const VectorSpace<Tensor<Type>, Type, 9> V(Zero);
+        const Tensor<Type> T(V);
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Construct given SphericalTensor<Type>:" << nl;
+        const SphericalTensor<Type> Sp(Type(5));
+        const Tensor<Type> T(Sp);
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Construct given SymmTensor<Type>:" << nl;
+        const SymmTensor<Type> S
+        (
+            Type(1), Type(2), Type(3),
+                     Type(5), Type(6),
+                              Type(9)
+        );
+        const Tensor<Type> T(S);
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Construct given triad of row vectors,"
+            << " optionally treated as transposed (ie, column vectors)" << nl;
+        const Vector<Vector<Type>> vecs
+        (
+            Vector<Type>(Type(1), Type(2), Type(3)),
+            Vector<Type>(Type(4), Type(5), Type(6)),
+            Vector<Type>(Type(7), Type(8), Type(9))
+        );
+        const Tensor<Type> T(vecs);
+        Info<< T << nl;
+
+        const Tensor<Type> transposedT(vecs, true);
+        Info<< transposedT << endl;
+    }
+    {
+        Info<< "# Construct given the three row vectors,"
+            << " optionally treated as transposed (ie, column vectors)" << nl;
+        const Vector<Type> a(Type(1), Type(2), Type(3));
+        const Vector<Type> b(Type(4), Type(5), Type(6));
+        const Vector<Type> c(Type(7), Type(8), Type(9));
+        const Tensor<Type> T(a, b, c);
+        Info<< T << nl;
+
+        const Tensor<Type> transposedT(a, b, c, true);
+        Info<< transposedT << endl;
+    }
+    {
+        Info<< "# Construct given the nine components:" << nl;
+        const Tensor<Type> T
+        (
+            Type(1), Type(2), Type(-3),
+            Type(4), Type(5), Type(-6),
+            Type(7), Type(8), Type(-9)
+        );
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Copy construct:" << nl;
+        const Tensor<Type> T(Zero);
+        const Tensor<Type> copyT(T);
+        Info<< T << tab << copyT << endl;
+    }
+}
+
+
+// Execute each member function of Tensor<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+    Tensor<Type> T
+    (
+        Type(1), Type(2), Type(-3),
+        Type(4), Type(5), Type(-6),
+        Type(7), Type(8), Type(-9)
+    );
+    Tensor<Type> Tbak = T;
+    const Tensor<Type> cT
+    (
+        Type(-9), Type(8), Type(7),
+        Type(-6), Type(5), Type(4),
+        Type(-3), Type(2), Type(1)
+    );
+
+    Info<< "# Operand: " << nl
+        << "  Tensor = " << T << endl;
+
+
+    {
+        Info<< "# Component access:" << nl;
+
+        Tensor<Type> cpT
+        (
+            T.xx(), T.xy(), T.xz(),
+            T.yx(), T.yy(), T.yz(),
+            T.zx(), T.zy(), T.zz()
+        );
+        cmp("  'Tensor' access:", T, cpT);
+
+        const Tensor<Type> cpcT
+        (
+            cT.xx(), cT.xy(), cT.xz(),
+            cT.yx(), cT.yy(), cT.yz(),
+            cT.zx(), cT.zy(), cT.zz()
+        );
+        cmp("  'const Tensor' access:", cT, cpcT);
+    }
+    {
+        Info<< "# Column-vector access:" << nl;
+        cmp("  cx():", T.cx(), Vector<Type>(Type(1), Type(4), Type(7)));
+        cmp("  cy():", T.cy(), Vector<Type>(Type(2), Type(5), Type(8)));
+        cmp("  cz():", T.cz(), Vector<Type>(Type(-3), Type(-6), Type(-9)));
+        cmp("  col(0):", T.col(0), Vector<Type>(Type(1), Type(4), Type(7)));
+        cmp("  col(1):", T.col(1), Vector<Type>(Type(2), Type(5), Type(8)));
+        cmp("  col(2):", T.col(2), Vector<Type>(Type(-3), Type(-6), Type(-9)));
+        cmp
+        (
+            "  col<0>:",
+            T.template col<0>(),
+            Vector<Type>(Type(1), Type(4), Type(7))
+        );
+        cmp
+        (
+            "  col<1>:",
+            T.template col<1>(),
+            Vector<Type>(Type(2), Type(5), Type(8))
+        );
+        cmp
+        (
+            "  col<2>:",
+            T.template col<2>(),
+            Vector<Type>(Type(-3), Type(-6), Type(-9))
+        );
+        // Compilation error:  Info << "  col<3> = " << T.col<3>() << nl;
+
+
+        Info<< "# Column-vector manipulation:" << nl;
+        T.col(1, Vector<Type>(Type(0), Type(1), Type(99)));
+        cmp
+        (
+            "  col(1, Vector):",
+            T.col(1),
+            Vector<Type>(Type(0), Type(1), Type(99))
+        );
+
+        T.cols
+        (
+            Vector<Type>(Type(1), Type(1), Type(1)),
+            Vector<Type>(Type(-1), Type(1), Type(2)),
+            Vector<Type>(Type(1), Type(1), Type(3))
+        );
+        cmp
+        (
+            "  cols(Vectors):",
+            T,
+            Tensor<Type>
+            (
+                Type(1), Type(-1), Type(1),
+                Type(1), Type(1),  Type(1),
+                Type(1), Type(2),  Type(3)
+            )
+        );
+    }
+    {
+        Info<< "# Row-vector access:" << nl;
+        T = Tbak;
+        cmp("  x():", T.x(), Vector<Type>(Type(1), Type(2), Type(-3)));
+        cmp("  y():", T.y(), Vector<Type>(Type(4), Type(5), Type(-6)));
+        cmp("  z():", T.z(), Vector<Type>(Type(7), Type(8), Type(-9)));
+        cmp("  row(0):", T.row(0), Vector<Type>(Type(1), Type(2), Type(-3)));
+        cmp("  row(1):", T.row(1), Vector<Type>(Type(4), Type(5), Type(-6)));
+        cmp("  row(2):", T.row(2), Vector<Type>(Type(7), Type(8), Type(-9)));
+        cmp
+        (
+            "  row<0>:",
+            T.template row<0>(),
+            Vector<Type>(Type(1), Type(2), Type(-3))
+        );
+        cmp
+        (
+            "  row<1>:",
+            T.template row<1>(),
+            Vector<Type>(Type(4), Type(5), Type(-6))
+        );
+        cmp
+        (
+            "  row<2>:",
+            T.template row<2>(),
+            Vector<Type>(Type(7), Type(8), Type(-9))
+        );
+        // Compilation error:  Info << "  row<3> = " << T.row<3>() << nl;
+
+
+        Info<< "# Row-vector manipulation:" << nl;
+        T.row(1, Vector<Type>(Type(0), Type(1), Type(99)));
+        cmp
+        (
+            "  row(1, Vector):",
+            T.row(1),
+            Vector<Type>(Type(0), Type(1), Type(99))
+        );
+
+        T.rows
+        (
+            Vector<Type>(Type(1), Type(1), Type(1)),
+            Vector<Type>(Type(-1), Type(1), Type(2)),
+            Vector<Type>(Type(1), Type(1), Type(3))
+        );
+        cmp
+        (
+            "  rows(Vectors):",
+            T,
+            Tensor<Type>
+            (
+                Type(1), Type(1), Type(1),
+                Type(-1), Type(1), Type(2),
+                Type(1), Type(1), Type(3)
+            )
+        );
+    }
+    {
+        Info<< "# Diagonal access:" << nl;
+
+        T = Tbak;
+        cmp
+        (
+            "  'Tensor'.diag():",
+            T.diag(),
+            Vector<Type>(Type(1), Type(5), Type(-9))
+        );
+        cmp
+        (
+            "  'const Tensor'.diag():",
+            cT.diag(),
+            Vector<Type>(Type(-9), Type(5), Type(1))
+        );
+
+
+        Info<< "# Diagonal manipulation:" << nl;
+
+        T.diag(Vector<Type>(Type(-10), Type(-15), Type(-20)));
+        cmp
+        (
+            "  'Tensor'.diag('Vector'):",
+            T.diag(),
+            Vector<Type>(Type(-10), Type(-15), Type(-20))
+        );
+    }
+    {
+        Info<< "# Tensor operations:" << nl;
+
+        T = Tbak;
+        cmp("  Transpose:", T, (T.T()).T());
+        cmp    // Singular matrix
+        (
+            "  Inverse:",
+            T.inv(),
+            Tensor<Type>
+            (
+            Type(-4.50359963e+15), Type(9.00719925e+15), Type(-4.50359963e+15),
+            Type(9.00719925e+15), Type(-1.80143985e+16), Type(9.00719925e+15),
+            Type(4.50359963e+15), Type(-9.00719925e+15), Type(4.50359963e+15)
+            )
+        );
+        cmp
+        (
+            "  Inner-product:",
+            T.inner(T),
+            Tensor<Type>
+            (
+                Type(-12), Type(-12), Type(12),
+                Type(-18), Type(-15), Type(12),
+                Type(-24), Type(-18), Type(12)
+            )
+        );
+        cmp
+        (
+            "  Schur-product:",
+            T.schur(T),
+            Tensor<Type>
+            (
+                Type(1),  Type(4),  Type(9),
+                Type(16), Type(25), Type(36),
+                Type(49), Type(64), Type(81)
+            )
+        );
+    }
+    {
+        Info<< "# Member operators:" << nl;
+
+        T = Tbak;
+        T &= T;
+        cmp
+        (
+            "  Assign inner-product of this with another Tensor:",
+            T,
+            (Tbak & Tbak)
+        );
+
+        T = VectorSpace<Tensor<Type>, Type, 9>(Zero);
+        cmp
+        (
+            "  Assign to an equivalent vector space:",
+            T,
+            Tensor<Type>(Zero)
+        );
+
+        T = SphericalTensor<Type>(Type(5));
+        cmp
+        (
+            "  Assign to a SphericalTensor:",
+            T,
+            Tensor<Type>
+            (
+                Type(5), Zero,    Zero,
+                Zero,    Type(5), Zero,
+                Zero,    Zero,    Type(5)
+            )
+        );
+
+        T = SymmTensor<Type>
+        (
+            Type(1), Type(2), Type(-3),
+                     Type(5), Type(-6),
+                              Type(-9)
+        );
+        cmp
+        (
+            "  Assign to a SymmTensor:",
+            T,
+            Tensor<Type>
+            (
+                Type(1),  Type(2),  Type(-3),
+                Type(2),  Type(5),  Type(-6),
+                Type(-3), Type(-6), Type(-9)
+            )
+        );
+
+        T = Vector<Vector<Type>>
+            (
+                Vector<Type>(Type(-1), Type(2), Type(3)),
+                Vector<Type>(Type(-4), Type(5), Type(6)),
+                Vector<Type>(Type(4), Type(-5), Type(6))
+            );
+        cmp
+        (
+            "  Assign to a triad of row vectors:",
+            T,
+            Tensor<Type>
+            (
+                Type(-1), Type(2), Type(3),
+                Type(-4), Type(5), Type(6),
+                Type(4), Type(-5), Type(6)
+            )
+        );
+    }
+}
+
+
+// Execute each global function of Tensor<Type>, and print output
+template<class Type>
+void test_global_funcs(Type)
+{
+    const Tensor<Type> T
+    (
+        Type(-1), Type(2), Type(-3),
+        Type(4), Type(5), Type(-6),
+        Type(7), Type(8), Type(-9)
+    );
+    const SymmTensor<Type> sT
+    (
+        Type(-1), Type(2), Type(-3),
+                  Type(5), Type(-6),
+                           Type(-9)
+    );
+
+    Info<< "# Operands: " << nl
+        << "  Tensor = " << T << nl
+        << "  SymmTensor = " << sT << endl;
+
+
+    cmp("  Trace = ", tr(T), Type(-5));
+    cmp("  Spherical part = ", sph(T), SphericalTensor<Type>(tr(T)/Type(3)));
+    cmp
+    (
+        "  Symmetric part = ",
+        symm(T),
+        SymmTensor<Type>
+        (
+            Type(-1), Type(3), Type(2),
+                      Type(5), Type(1),
+                               Type(-9)
+        )
+    );
+    cmp
+    (
+        "  Twice the symmetric part = ",
+        twoSymm(T),
+        SymmTensor<Type>
+        (
+            Type(-2), Type(6),  Type(4),
+                      Type(10), Type(2),
+                                Type(-18)
+        )
+    );
+    cmp
+    (
+        "  Skew-symmetric part = ",
+        skew(T),
+        Tensor<Type>
+        (
+            Type(0),  Type(-1), Type(-5),
+            Type(1),  Type(0),  Type(-7),
+            Type(5),  Type(7),  Type(0)
+        )
+    );
+    /*
+    // Complex-type is not supported for this function.
+    cmp
+    (
+        "  Skew-symmetric part of a SymmTensor = ",
+        skew(sT),
+        Tensor<Type>(Zero)
+    );
+    */
+    cmp
+    (
+        "  Deviatoric part = ",
+        dev(T),
+        Tensor<Type>
+        (
+            Type(0.66666667), Type(2),          Type(-3),
+            Type(4),          Type(6.66666667), Type(-6),
+            Type(7),          Type(8),          Type(-7.33333333)
+        ),
+        1e-6,
+        1e-6
+    );
+    cmp("  Two-third deviatoric part = ", dev2(T), T - 2*sph(T));
+    cmp("  Determinant = ", det(T), Type(-6.000000000000005));
+    cmp
+    (
+        "  Cofactor tensor = ",
+        cof(T),
+        Tensor<Type>
+        (
+            Type(3),  Type(-6),  Type(-3),
+            Type(-6), Type(30),  Type(22),
+            Type(3),  Type(-18), Type(-13)
+        )
+    );
+    cmp
+    (
+        "  Inverse = ",
+        inv(T, det(T)),
+        Tensor<Type>
+        (
+            Type(-0.5),  Type(1),            Type(-0.5),
+            Type(1),     Type(-5),           Type(3),
+            Type(0.5),   Type(-3.66666667),  Type(2.16666667)
+        ),
+        1e-8,
+        1e-8
+    );
+    cmp
+    (
+        "  Inverse (another) = ",
+        inv(T),
+        Tensor<Type>
+        (
+            Type(-0.5),  Type(1),            Type(-0.5),
+            Type(1),     Type(-5),           Type(3),
+            Type(0.5),   Type(-3.66666667),  Type(2.16666667)
+        ),
+        1e-8,
+        1e-8
+    );
+    cmp
+    (
+        "  Inverse (another) = ",
+        T.inv(),
+        Tensor<Type>
+        (
+            Type(-0.5),  Type(1),            Type(-0.5),
+            Type(1),     Type(-5),           Type(3),
+            Type(0.5),   Type(-3.66666667),  Type(2.16666667)
+        ),
+        1e-8,
+        1e-8
+    );
+    cmp("  First invariant = ", invariantI(T), Type(-5));
+    cmp("  Second invariant = ", invariantII(T), Type(20));
+    cmp("  Third invariant = ", invariantIII(T), Type(-6.000000000000005));
+}
+
+
+// Execute each global operator of Tensor<Type>, and print output
+template<class Type>
+void test_global_opers(Type)
+{
+    const Tensor<Type> T
+    (
+        Type(-1), Type(2), Type(-3),
+        Type(4),  Type(5),  Type(-6),
+        Type(7),  Type(8),  Type(-9)
+    );
+    const SymmTensor<Type> sT
+    (
+        Type(-1), Type(2), Type(-3),
+                 Type(5), Type(-6),
+                          Type(-9)
+    );
+    const SphericalTensor<Type> spT(Type(1));
+    const Vector<Type> v(Type(3), Type(2), Type(1));
+    const Type x(4);
+
+    Info<< "# Operands:" << nl
+        << "  Tensor = " << T << nl
+        << "  SymmTensor = " << sT << nl
+        << "  SphericalTensor = " << spT << nl
+        << "  Vector = " << v << nl
+        << "  Type = " << x << endl;
+
+
+    cmp
+    (
+        "  Sum of SpTensor-Tensor = ",
+        (spT + T),
+        Tensor<Type>
+        (
+            Type(0), Type(2), Type(-3),
+            Type(4), Type(6), Type(-6),
+            Type(7), Type(8), Type(-8)
+        )
+    );
+    cmp
+    (
+        "  Sum of Tensor-SpTensor = ",
+        (T + spT),
+        Tensor<Type>
+        (
+            Type(0), Type(2), Type(-3),
+            Type(4), Type(6), Type(-6),
+            Type(7), Type(8), Type(-8)
+        )
+    );
+    cmp
+    (
+        "  Sum of SymmTensor-Tensor = ",
+        (sT + T),
+        Tensor<Type>
+        (
+            Type(-2), Type(4),  Type(-6),
+            Type(6),  Type(10), Type(-12),
+            Type(4),  Type(2),  Type(-18)
+        )
+    );
+    cmp
+    (
+        "  Sum of Tensor-SymmTensor = ",
+        (T + sT),
+        Tensor<Type>
+        (
+            Type(-2), Type(4),  Type(-6),
+            Type(6),  Type(10), Type(-12),
+            Type(4),  Type(2),  Type(-18)
+        )
+    );
+    cmp
+    (
+        "  Subtract Tensor from SpTensor = ",
+        (spT - T),
+        Tensor<Type>
+        (
+            Type(2),  Type(-2), Type(3),
+            Type(-4), Type(-4), Type(6),
+            Type(-7), Type(-8), Type(10)
+        )
+    );
+    cmp
+    (
+        "  Subtract SpTensor from Tensor = ",
+        (T - spT),
+        Tensor<Type>
+        (
+            Type(-2), Type(2), Type(-3),
+            Type(4),  Type(4), Type(-6),
+            Type(7),  Type(8), Type(-10)
+        )
+    );
+    cmp
+    (
+        "  Subtract Tensor from SymmTensor = ",
+        (sT - T),
+        Tensor<Type>
+        (
+            Type(0),   Type(0),   Type(0),
+            Type(-2),  Type(0),   Type(0),
+            Type(-10), Type(-14), Type(0)
+        )
+    );
+    cmp
+    (
+        "  Subtract SymmTensor from Tensor = ",
+        (T - sT),
+        Tensor<Type>
+        (
+            Type(0),  Type(0),  Type(0),
+            Type(2),  Type(0),  Type(0),
+            Type(10), Type(14), Type(0)
+        )
+    );
+    cmp
+    (
+        "  Hodge dual of a Tensor = ",
+        *T,
+        Vector<Type>(T.yz(), -T.xz(), T.xy())
+    );
+    cmp
+    (
+        "  Hodge dual of a Vector = ",
+        *v,
+        Tensor<Type>
+        (
+             Zero,  -v.z(),  v.y(),
+             v.z(),  Zero,  -v.x(),
+            -v.y(),  v.x(),  Zero
+        )
+    );
+    /*cmp
+    (
+        "  Division of Vector by Tensor = ",
+        (v/T),
+        Tensor<Type>
+        (
+            Type(-3),          Type(1),    Type(-0.33333333),
+            Type(0.75),        Type(0.4),  Type(-0.16666667),
+            Type(0.42857143),  Type(0.25), Type(-0.11111111)
+        )
+    );*/
+    cmp
+    (
+        "  Division of Tensor by Type = ",
+        (T/x),
+        Tensor<Type>
+        (
+            Type(-0.25), Type(0.5),  Type(-0.75),
+            Type(1),     Type(1.25), Type(-1.5),
+            Type(1.75),  Type(2),    Type(-2.25)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of Tensor-Tensor = ",
+        (T & T),
+        Tensor<Type>
+        (
+            Type(-12), Type(-16), Type(18),
+            Type(-26), Type(-15), Type(12),
+            Type(-38), Type(-18), Type(12)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SpTensor-Tensor = ",
+        (spT & T),
+        Tensor<Type>
+        (
+            Type(-1), Type(2),  Type(-3),
+            Type(4),  Type(5),  Type(-6),
+            Type(7),  Type(8),  Type(-9)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of Tensor-SpTensor = ",
+        (T & spT),
+        Tensor<Type>
+        (
+            Type(-1), Type(2), Type(-3),
+            Type(4),  Type(5), Type(-6),
+            Type(7),  Type(8), Type(-9)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SymmTensor-Tensor = ",
+        (sT & T),
+        Tensor<Type>
+        (
+            Type(-12), Type(-16),  Type(18),
+            Type(-24), Type(-19),  Type(18),
+            Type(-84), Type(-108), Type(126)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of Tensor-SymmTensor = ",
+        (T & sT),
+        Tensor<Type>
+        (
+            Type(14), Type(26),  Type(18),
+            Type(24), Type(69),  Type(12),
+            Type(36), Type(108), Type(12)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of Tensor-Vector = ",
+        (T & v),
+        Vector<Type>(Type(-2), Type(16),  Type(28)) // Column-vector
+    );
+    cmp
+    (
+        "  Inner-product of Vector-Tensor = ",
+        (v & T),
+        Vector<Type>(Type(12), Type(24),  Type(-30)) // Row-vector
+    );
+    cmp("  D-inner-product of SpTensor-Tensor = ", (spT && T), Type(-5));
+    cmp("  D-inner-product of Tensor-SpTensor = ", (T && spT), Type(-5));
+    cmp("  D-inner-product of SymmTensor-Tensor = ", (sT && T), Type(95));
+    cmp("  D-inner-product of Tensor-SymmTensor = ", (T && sT), Type(95));
+    cmp
+    (
+        "  Outer-product of Vector-Vector = ",
+        (v*v),
+        Tensor<Type>
+        (
+            Type(9), Type(6), Type(3),
+            Type(6), Type(4), Type(2),
+            Type(3), Type(2), Type(1)
+        )
+    );
+}
+
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
+    test_member_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global functions: "<< typeID[I] << " ##" << nl;
+    test_global_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global operators: "<< typeID[I] <<" ##" << nl;
+    test_global_opers(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main()
+{
+    const std::tuple<floatScalar, doubleScalar, complex> types
+    (
+        std::make_tuple(Zero, Zero, Zero)
+    );
+
+    const List<word> typeID
+    ({
+        "Tensor<floatScalar>",
+        "Tensor<doubleScalar>",
+        "Tensor<complex>"
+    });
+
+    run_tests(types, typeID);
+
+
+    if (nFail_)
+    {
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
+    }
+
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
+    return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/test/Tensor2D/Make/files

@@ -0,0 +1,3 @@
+Test-Tensor2D.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-Tensor2D

applications/test/Tensor2D/Make/options

@@ -0,0 +1,2 @@
+/* EXE_INC = */
+/* EXE_LIBS = */

applications/test/Tensor2D/Test-Tensor2D.C

@@ -0,0 +1,841 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2014 OpenFOAM Foundation
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    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.
+
+    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
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-Tensor2D
+
+Description
+    Tests for \c Tensor2D constructors, member functions and operators
+    using \c floatScalar, \c doubleScalar, and \c complex base types.
+
+    Eigen decomposition tests for \c tensor2D, i.e. Tensor2D<scalar>.
+
+    Cross-checks were obtained from 'NumPy 1.15.1' and 'SciPy 1.1.0' if no
+    theoretical cross-check exists (like eigendecomposition relations), and
+    were hard-coded for elementwise comparisons.
+
+    For \c complex base type, the cross-checks do only involve zero imag part.
+
+\*---------------------------------------------------------------------------*/
+
+#include "vector2DField.H"
+#include "tensor2D.H"
+#include "symmTensor2D.H"
+#include "transform.H"
+#include "Random.H"
+#include "floatScalar.H"
+#include "doubleScalar.H"
+#include "complex.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Total number of unit tests
+unsigned nTest_ = 0;
+
+
+// Total number of failed unit tests
+unsigned nFail_ = 0;
+
+
+// Create a random tensor2D
+tensor2D makeRandomContainer(Random& rnd)
+{
+    tensor2D A(Zero);
+    std::generate(A.begin(), A.end(), [&]{ return rnd.GaussNormal<scalar>(); });
+    return A;
+}
+
+
+// Compare two floating point types, and print output.
+// Do ++nFail_ if values of two objects are not equal within a given tolerance.
+// The function is converted from PEP-485.
+template<class Type>
+typename std::enable_if
+<
+    std::is_same<floatScalar, Type>::value ||
+    std::is_same<doubleScalar, Type>::value ||
+    std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,  //<! are values the same within 8 decimals
+    const scalar absTol = 0      //<! useful for cmps near zero
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
+    {
+        ++nFail;
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Compare two containers elementwise, and print output.
+// Do ++nFail_ if two components are not equal within a given tolerance.
+// The function is converted from PEP-485
+template<class Type>
+typename std::enable_if
+<
+    !std::is_same<floatScalar, Type>::value &&
+    !std::is_same<doubleScalar, Type>::value &&
+    !std::is_same<complex, Type>::value,
+    void
+>::type cmp
+(
+    const word& msg,
+    const Type& x,
+    const Type& y,
+    const scalar relTol = 1e-8,
+    const scalar absTol = 0
+)
+{
+    Info<< msg << x << endl;
+
+    unsigned nFail = 0;
+
+    for (label i = 0; i < pTraits<Type>::nComponents; ++i)
+    {
+        if (max(absTol, relTol*max(mag(x[i]), mag(y[i]))) < mag(x[i] - y[i]))
+        {
+            ++nFail;
+        }
+    }
+
+    if (nFail)
+    {
+        Info<< nl
+            << "        #### Fail in " << nFail << " comps ####" << nl << endl;
+        ++nFail_;
+    }
+    ++nTest_;
+}
+
+
+// Create each constructor of Tensor2D<Type>, and print output
+template<class Type>
+void test_constructors(Type)
+{
+    {
+        Info<< "# Construct initialized to zero:" << nl;
+        const Tensor2D<Type> T(Zero);
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Construct given VectorSpace:" << nl;
+        const VectorSpace<Tensor2D<Type>, Type, 4> V(Zero);
+        const Tensor2D<Type> T(V);
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Construct given SymmTensor2D:" << nl;
+        const SymmTensor2D<Type> S
+        (
+            Type(1), Type(2),
+                     Type(3)
+        );
+        const Tensor2D<Type> T(S);
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Construct given SphericalTensor2D:" << nl;
+        const SphericalTensor2D<Type> Sp(Type(5));
+        const Tensor2D<Type> T(Sp);
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Construct given the two row vectors:" << nl;
+        const Vector2D<Type> x(Type(1), Type(2));
+        const Vector2D<Type> y(Type(3), Type(4));
+        const Tensor2D<Type> T(x, y);
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Construct given the four components:" << nl;
+        const Tensor2D<Type> T
+        (
+            Type(1), Type(2),
+            Type(3), Type(4)
+        );
+        Info<< T << endl;
+    }
+    {
+        Info<< "# Copy construct:" << nl;
+        const Tensor2D<Type> T(Zero);
+        const Tensor2D<Type> Tcopy(T);
+        Info<< T << endl;
+    }
+}
+
+
+// Execute each member function of Tensor2D<Type>, and print output
+template<class Type>
+void test_member_funcs(Type)
+{
+    Tensor2D<Type> T
+    (
+        Type(1), Type(2),
+        Type(4), Type(5)
+    );
+    Tensor2D<Type> Tbak = T;
+    const Tensor2D<Type> cT
+    (
+        Type(-9), Type(8),
+        Type(-6), Type(5)
+    );
+
+    Info<< "# Operand: " << nl
+        << "  Tensor2D = " << T << endl;
+
+
+    {
+        Info<< "# Component access:" << nl;
+
+        Tensor2D<Type> cpT
+        (
+            T.xx(), T.xy(),
+            T.yx(), T.yy()
+        );
+        cmp("  'Tensor2D' access:", T, cpT);
+
+        const Tensor2D<Type> cpcT
+        (
+            cT.xx(), cT.xy(),
+            cT.yx(), cT.yy()
+        );
+        cmp("  'const Tensor2D' access:", cT, cpcT);
+    }
+    {
+        Info<< "# Column-vector access:" << nl;
+        cmp("  cx():", T.cx(), Vector2D<Type>(Type(1), Type(4)));
+        cmp("  cy():", T.cy(), Vector2D<Type>(Type(2), Type(5)));
+        cmp("  col(0):", T.col(0), Vector2D<Type>(Type(1), Type(4)));
+        cmp("  col(1):", T.col(1), Vector2D<Type>(Type(2), Type(5)));
+        cmp
+        (
+            "  col<0>:",
+            T.template col<0>(),
+            Vector2D<Type>(Type(1), Type(4))
+        );
+        cmp
+        (
+            "  col<1>:",
+            T.template col<1>(),
+            Vector2D<Type>(Type(2), Type(5))
+        );
+        // Compilation error:  Info << "  col<2> = " << T.col<2>() << nl;
+
+
+        Info<< "# Column-vector manipulation:" << nl;
+        T.col(1, Vector2D<Type>(Type(0), Type(1)));
+        cmp
+        (
+            "  col(1, Vector):",
+            T.col(1),
+            Vector2D<Type>(Type(0), Type(1))
+        );
+
+        T.cols
+        (
+            Vector2D<Type>(Type(1), Type(1)),
+            Vector2D<Type>(Type(-1), Type(1))
+        );
+        cmp
+        (
+            "  cols(Vectors):",
+            T,
+            Tensor2D<Type>
+            (
+                Type(1), Type(-1),
+                Type(1), Type(1)
+            )
+        );
+    }
+    {
+        Info<< "# Row-vector access:" << nl;
+        T = Tbak;
+        cmp("  x():", T.x(), Vector2D<Type>(Type(1), Type(2)));
+        cmp("  y():", T.y(), Vector2D<Type>(Type(4), Type(5)));
+        cmp("  row(0):", T.row(0), Vector2D<Type>(Type(1), Type(2)));
+        cmp("  row(1):", T.row(1), Vector2D<Type>(Type(4), Type(5)));
+        cmp
+        (
+            "  row<0>:",
+            T.template row<0>(),
+            Vector2D<Type>(Type(1), Type(2))
+        );
+        cmp
+        (
+            "  row<1>:",
+            T.template row<1>(),
+            Vector2D<Type>(Type(4), Type(5))
+        );
+        // Compilation error:  Info << "  row<2> = " << T.row<2>() << nl;
+
+
+        Info<< "# Row-vector manipulation:" << nl;
+        T.row(1, Vector2D<Type>(Type(0), Type(1)));
+        cmp
+        (
+            "  row(1, Vector):",
+            T.row(1),
+            Vector2D<Type>(Type(0), Type(1))
+        );
+
+        T.rows
+        (
+            Vector2D<Type>(Type(1), Type(1)),
+            Vector2D<Type>(Type(-1), Type(1))
+        );
+        cmp
+        (
+            "  rows(Vectors):",
+            T,
+            Tensor2D<Type>
+            (
+                Type(1), Type(1),
+                Type(-1), Type(1)
+            )
+        );
+    }
+    {
+        Info<< "# Diagonal access:" << nl;
+
+        T = Tbak;
+        cmp
+        (
+            "  'Tensor2D'.diag():",
+            T.diag(),
+            Vector2D<Type>(Type(1), Type(5))
+        );
+        cmp
+        (
+            "  'const Tensor2D'.diag():",
+            cT.diag(),
+            Vector2D<Type>(Type(-9), Type(5))
+        );
+
+
+        Info<< "# Diagonal manipulation:" << nl;
+
+        T.diag(Vector2D<Type>(Type(-10), Type(-15)));
+        cmp
+        (
+            "  'Tensor2D'.diag('Vector'):",
+            T.diag(),
+            Vector2D<Type>(Type(-10), Type(-15))
+        );
+    }
+    {
+        Info<< "# Tensor operations:" << nl;
+
+        T = Tbak;
+        cmp("  Transpose:", T, (T.T()).T());
+        cmp
+        (
+            "  Inner-product:",
+            T.inner(T),
+            Tensor2D<Type>
+            (
+                Type(9),  Type(12),
+                Type(24), Type(33)
+            )
+        );
+        cmp
+        (
+            "  Schur-product:",
+            T.schur(T),
+            Tensor2D<Type>
+            (
+                Type(1),  Type(4),
+                Type(16), Type(25)
+            )
+        );
+    }
+    {
+        Info<< "# Member operators:" << nl;
+
+        T = SphericalTensor2D<Type>(Type(5));
+        cmp
+        (
+            "  Assign to a SphericalTensor2D:",
+            T,
+            Tensor2D<Type>
+            (
+                Type(5), Zero,
+                Zero,    Type(5)
+            )
+        );
+
+        T = SymmTensor2D<Type>
+        (
+            Type(1), Type(2),
+                     Type(5)
+        );
+        cmp
+        (
+            "  Assign to a SymmTensor2D:",
+            T,
+            Tensor2D<Type>
+            (
+                Type(1),  Type(2),
+                Type(2),  Type(5)
+            )
+        );
+    }
+}
+
+
+// Execute each global function of Tensor2D<Type>, and print output
+template<class Type>
+void test_global_funcs(Type)
+{
+    const Tensor2D<Type> T
+    (
+        Type(-1), Type(2),
+        Type(4), Type(5)
+    );
+
+    Info<< "# Operand: " << nl
+        << "  Tensor2D = " << T << endl;
+
+
+    cmp("  Trace = ", tr(T), Type(4));
+    cmp("  Spherical part = ", sph(T), SphericalTensor2D<Type>(tr(T)/Type(2)));
+    cmp
+    (
+        "  Symmetric part = ",
+        symm(T),
+        SymmTensor2D<Type>
+        (
+            Type(-1), Type(3),
+                      Type(5)
+        )
+    );
+    cmp
+    (
+        "  Twice the symmetric part = ",
+        twoSymm(T),
+        SymmTensor2D<Type>
+        (
+            Type(-2), Type(6),
+                      Type(10)
+        )
+    );
+    cmp
+    (
+        "  Skew-symmetric part = ",
+        skew(T),
+        Tensor2D<Type>
+        (
+            Type(0),  Type(-1),
+            Type(1),  Type(0)
+        )
+    );
+    cmp
+    (
+        "  Deviatoric part = ",
+        dev(T),
+        Tensor2D<Type>
+        (
+            Type(-3), Type(2),
+            Type(4),  Type(3)
+        )
+    );
+    cmp("  Two-third deviatoric part = ", dev2(T), T - 2*sph(T));
+    cmp("  Determinant = ", det(T), Type(-13));
+    cmp
+    (
+        "  Cofactor tensor2D = ",
+        cof(T),
+        Tensor2D<Type>
+        (
+            Type(5),  Type(-4),
+            Type(-2), Type(-1)
+        )
+    );
+    cmp
+    (
+        "  Inverse = ",
+        inv(T, det(T)),
+        Tensor2D<Type>
+        (
+            Type(-0.38461538),  Type(0.15384615),
+            Type(0.30769231),   Type(0.07692308)
+        ),
+        1e-6,
+        1e-6
+    );
+    cmp
+    (
+        "  Inverse (another) = ",
+        inv(T),
+        Tensor2D<Type>
+        (
+            Type(-0.38461538),  Type(0.15384615),
+            Type(0.30769231),   Type(0.07692308)
+        ),
+        1e-6,
+        1e-6
+    );
+    cmp("  First invariant = ", invariantI(T), Type(4));
+    cmp("  Second invariant = ", invariantII(T), Type(-13));
+}
+
+
+// Execute each global operator of Tensor2D<Type>, and print output
+template<class Type>
+void test_global_opers(Type)
+{
+    const Tensor2D<Type> T
+    (
+        Type(-1), Type(2),
+        Type(4),  Type(5)
+    );
+    const SymmTensor2D<Type> sT
+    (
+        Type(1), Type(2),
+                 Type(5)
+    );
+    const SphericalTensor2D<Type> spT(Type(1));
+    const Vector2D<Type> v(Type(3), Type(2));
+    const Type x(4);
+
+    Info<< "# Operands:" << nl
+        << "  Tensor2D = " << T << nl
+        << "  SymmTensor2D = " << sT << nl
+        << "  SphericalTensor2D = " << spT << nl
+        << "  Vector2D = " << v << nl
+        << "  Type = " << x << endl;
+
+
+    cmp
+    (
+        "  Sum of SpTensor2D-Tensor2D = ",
+        (spT + T),
+        Tensor2D<Type>
+        (
+            Type(0), Type(2),
+            Type(4), Type(6)
+        )
+    );
+    cmp
+    (
+        "  Sum of Tensor2D-SpTensor2D = ",
+        (T + spT),
+        Tensor2D<Type>
+        (
+            Type(0), Type(2),
+            Type(4), Type(6)
+        )
+    );
+    cmp
+    (
+        "  Sum of SymmTensor2D-Tensor2D = ",
+        (sT + T),
+        Tensor2D<Type>
+        (
+            Type(0), Type(4),
+            Type(6),  Type(10)
+        )
+    );
+    cmp
+    (
+        "  Sum of Tensor2D-SymmTensor2D = ",
+        (T + sT),
+        Tensor2D<Type>
+        (
+            Type(0), Type(4),
+            Type(6),  Type(10)
+        )
+    );
+    cmp
+    (
+        "  Subtract Tensor2D from SpTensor2D = ",
+        (spT - T),
+        Tensor2D<Type>
+        (
+            Type(2),   Type(-2),
+            Type(-4),  Type(-4)
+        )
+    );
+    cmp
+    (
+        "  Subtract SpTensor2D from Tensor2D = ",
+        (T - spT),
+        Tensor2D<Type>
+        (
+            Type(-2), Type(2),
+            Type(4),  Type(4)
+        )
+    );
+    cmp
+    (
+        "  Subtract Tensor2D from SymmTensor2D = ",
+        (sT - T),
+        Tensor2D<Type>
+        (
+            Type(2),  Type(0),
+            Type(-2), Type(0)
+        )
+    );
+    cmp
+    (
+        "  Subtract SymmTensor2D from Tensor2D = ",
+        (T - sT),
+        Tensor2D<Type>
+        (
+            Type(-2), Type(0),
+            Type(2), Type(0)
+        )
+    );
+    cmp
+    (
+        "  Division of Tensor2D by Type = ",
+        (T/x),
+        Tensor2D<Type>
+        (
+            Type(-0.25), Type(0.5),
+            Type(1),     Type(1.25)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of Tensor2D-Tensor2D = ",
+        (T & T),
+        Tensor2D<Type>
+        (
+            Type(9),  Type(8),
+            Type(16), Type(33)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SpTensor2D-Tensor2D = ",
+        (spT & T),
+        Tensor2D<Type>
+        (
+            Type(-1), Type(2),
+            Type(4),  Type(5)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of Tensor2D-SpTensor2D = ",
+        (T & spT),
+        Tensor2D<Type>
+        (
+            Type(-1), Type(2),
+            Type(4),  Type(5)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of SymmTensor2D-Tensor2D = ",
+        (sT & T),
+        Tensor2D<Type>
+        (
+            Type(7),  Type(12),
+            Type(18), Type(29)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of Tensor2D-SymmTensor2D = ",
+        (T & sT),
+        Tensor2D<Type>
+        (
+            Type(3),  Type(8),
+            Type(14), Type(33)
+        )
+    );
+    cmp
+    (
+        "  Inner-product of Tensor2D-Vector2D = ",
+        (T & v),
+        Vector2D<Type>(Type(1), Type(22)) // Column-vector
+    );
+    cmp
+    (
+        "  Inner-product of Vector2D-Tensor2D = ",
+        (v & T),
+        Vector2D<Type>(Type(5), Type(16)) // Row-vector
+    );
+    cmp("  D-inner-product of SpTensor2D-Tensor2D = ", (spT && T), Type(4));
+    cmp("  D-inner-product of Tensor2D-SpTensor2D = ", (T && spT), Type(4));
+    cmp("  D-inner-product of SymmTensor2D-Tensor2D = ", (sT && T), Type(36));
+    cmp("  D-inner-product of Tensor2D-SymmTensor2D = ", (T && sT), Type(36));
+    cmp
+    (
+        "  Outer-product of Vector2D-Vector2D = ",
+        (v*v),
+        Tensor2D<Type>
+        (
+            Type(9), Type(6),
+            Type(6), Type(4)
+        )
+    );
+}
+
+
+// Do compile-time recursion over the given types
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I == sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
+
+
+template<std::size_t I = 0, typename... Tp>
+inline typename std::enable_if<I < sizeof...(Tp), void>::type
+run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
+{
+    Info<< nl << "    ## Test constructors: "<< typeID[I] <<" ##" << nl;
+    test_constructors(std::get<I>(types));
+
+    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
+    test_member_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global functions: "<< typeID[I] << " ##" << nl;
+    test_global_funcs(std::get<I>(types));
+
+    Info<< nl << "    ## Test global operators: "<< typeID[I] <<" ##" << nl;
+    test_global_opers(std::get<I>(types));
+
+    run_tests<I + 1, Tp...>(types, typeID);
+}
+
+
+// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+    const std::tuple<floatScalar, doubleScalar, complex> types
+    (
+        std::make_tuple(Zero, Zero, Zero)
+    );
+
+    const List<word> typeID
+    ({
+        "Tensor2D<floatScalar>",
+        "Tensor2D<doubleScalar>",
+        "Tensor2D<complex>"
+    });
+
+    run_tests(types, typeID);
+
+
+    {
+        Info<< "# Pre-v2006 tests:" << nl;
+        vector2D v1(1, 2), v2(3, 4);
+        tensor2D t3 = v1*v2;
+
+        Info<< v1 << "*" << v2 << " = " << t3 << endl;
+
+        {
+            Info<< "rows:" << nl;
+            for (direction i = 0; i < 2; ++i)
+            {
+                Info<< "  (" << i << ") = " << t3.row(i) << nl;
+            }
+        }
+
+        {
+            Info<< "cols:" << nl;
+            for (direction i = 0; i < 2; ++i)
+            {
+                Info<< "  (" << i << ") = " << t3.col(i) << nl;
+            }
+            Info<< "col<0> = " << t3.col<0>() << nl;
+            Info<< "col<1> = " << t3.col<1>() << nl;
+            // Compilation error:  Info << "col<3> = " << t3.col<3>() << nl;
+
+            t3.col<0>({0, 2});
+            Info<< "replaced col<0> = " << t3.col<0>() << nl;
+            Info<< "tensor " << t3 << nl;
+
+            t3.row<1>(Zero);
+            Info<< "replaced row<1> = " << t3.row<1>() << nl;
+            Info<< "tensor " << t3 << nl;
+        }
+
+
+        {
+            vector2DField vfld1(8, Zero);
+
+            forAll(vfld1, i)
+            {
+                vfld1[i] = (i + 1) * ((i % 2) ? v1 : v2);
+            }
+
+            Info<< "vector: " << flatOutput(vfld1) << nl;
+
+            scalarField xvals(8);
+            scalarField yvals(8);
+            unzip(vfld1, xvals, yvals);
+
+            Info<< "unzip" << nl
+                << " x => " << flatOutput(xvals) << nl
+                << " y => " << flatOutput(yvals) << nl;
+
+            reverse(xvals);
+            zip(vfld1, xvals, yvals);
+
+            Info<< "rezip (with reversed x)" << nl
+                << "   => " << flatOutput(vfld1) << nl;
+        }
+    }
+
+
+    if (nFail_)
+    {
+        Info<< nl << "        #### "
+            << "Failed in " << nFail_ << " tests "
+            << "out of total " << nTest_ << " tests "
+            << "####\n" << endl;
+        return 1;
+    }
+
+    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/tensor/Make/files

@@ -1,3 +0,0 @@
-Test-tensor.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-tensor

applications/test/tensor/Test-tensor.C

@@ -1,235 +0,0 @@
-#include "tensor.H"
-#include "triad.H"
-#include "symmTensor.H"
-#include "transform.H"
-#include "stringList.H"
-#include "IOstreams.H"
-
-using namespace Foam;
-
-int main()
-{
-    tensor t1(1, 2, 3, 4, 5, 6, 7, 8, 9);
-    tensor t2(1, 2, 3, 1, 2, 3, 1, 2, 3);
-
-    tensor t3 = t1 + t2;
-
-    Info<< "tensor " << t1 << " + " << t2 << nl
-        << " = " << t3 << nl << nl;
-
-    {
-        Info<< "rows:" << nl;
-        for (direction i=0; i < 3; ++i)
-        {
-            Info<< "  (" << i << ") = " << t3.row(i) << nl;
-        }
-    }
-    {
-        Info<< "cols:" << nl;
-        for (direction i=0; i < 3; ++i)
-        {
-            Info<< "  (" << i << ") = " << t3.col(i) << nl;
-        }
-        Info<< "col<0> = " << t3.col<0>() << nl;
-        Info<< "col<1> = " << t3.col<1>() << nl;
-        Info<< "col<2> = " << t3.col<2>() << nl;
-        // Compilation error:  Info << "col<3> = " << t3.col<3>() << nl;
-
-        t3.col<1>({0, 2, 4});
-        Info<< "replaced col<1> = " << t3.col<1>() << nl;
-        Info<< "tensor " << t3 << nl;
-
-        t3.row<2>(Zero);
-        Info<< "replaced row<2> = " << t3.row<2>() << nl;
-        Info<< "tensor " << t3 << nl;
-
-        triad tr3(t3);
-        Info<< "triad " << tr3 << " :: T() " << tr3.T() << nl;
-    }
-    Info<< nl;
-
-
-    tensor t4(3,-2,1,-2,2,0,1, 0, 4);
-
-    Info<< inv(t4) << endl;
-    Info<< (inv(t4) & t4) << endl;
-
-    Info<< t1.x() << t1.y() << t1.z() << endl;
-
-    tensor t6(1,0,-4,0,5,4,-4,4,3);
-    //tensor t6(1,2,0,2,5,0,0,0,0);
-
-    Info<< "tensor " << t6 << endl;
-    vector e = eigenValues(t6);
-    Info<< "eigenvalues " << e << endl;
-
-    tensor ev = eigenVectors(t6);
-    Info<< "eigenvectors " << ev << endl;
-
-    Info<< "Check determinant " << e.x()*e.y()*e.z() << " " << det(t6) << endl;
-
-    Info<< "Check eigenvectors "
-        << (eigenVectors(t6, e) & t6) << " "
-        << (e.x() * eigenVectors(t6, e)).x()
-        << (e.x() * eigenVectors(t6, e)).y()
-        << (e.z() * eigenVectors(t6, e)).z()
-        << endl;
-
-    Info<< "Check eigenvalues for symmTensor "
-        << eigenValues(symm(t6)) - eigenValues(tensor(symm(t6))) << endl;
-
-    Info<< "Check eigenvectors for symmTensor "
-        << eigenVectors(symm(t6)) - eigenVectors(tensor(symm(t6))) << endl;
-
-    tensor t7(1, 2, 3, 2, 4, 5, 3, 5, 6);
-
-    Info<< "Check transformation "
-        << (t1 & t7 & t1.T()) << " " << transform(t1, t7) << endl;
-
-    symmTensor st1(1, 2, 3, 4, 5, 6);
-    symmTensor st2(7, 8, 9, 10, 11, 12);
-
-    Info<< "Check symmetric transformation "
-        << transform(t1, st1) << endl;
-
-    Info<< "Check dot product of symmetric tensors "
-        << (st1 & st2) << endl;
-
-    Info<< "Check inner sqr of a symmetric tensor "
-        << innerSqr(st1) << " " << innerSqr(st1) - (st1 & st1) << endl;
-
-    Info<< "Check symmetric part of dot product of symmetric tensors "
-        << twoSymm(st1&st2) - ((st1&st2) + (st2&st1)) << endl;
-
-    tensor sk1 = skew(t6);
-    Info<< "Check dot product of symmetric and skew tensors "
-        << twoSymm(st1&sk1) - ((st1&sk1) - (sk1&st1)) << endl;
-
-    vector v1(1, 2, 3);
-
-    Info<< sqr(v1) << endl;
-    Info<< symm(t7) << endl;
-    Info<< twoSymm(t7) << endl;
-    Info<< magSqr(st1) << endl;
-    Info<< mag(st1) << endl;
-
-    Info<< (symm(t7) && t7) - (0.5*(t7 + t7.T()) && t7) << endl;
-    Info<< (t7 && symm(t7)) - (t7 && 0.5*(t7 + t7.T())) << endl;
-
-
-    /*
-    // Lots of awkward eigenvector tests ...
-
-    tensor T_rand_real
-    (
-        0.9999996423721313, 0.3330855667591095, 0.6646450161933899,
-        0.9745196104049683, 0.0369445420801640, 0.0846728682518005,
-        0.6474838852882385, 0.1617118716239929, 0.2041363865137100
-    );
-    DebugVar(T_rand_real);
-    vector L_rand_real(eigenValues(T_rand_real));
-    DebugVar(L_rand_real);
-    tensor U_rand_real(eigenVectors(T_rand_real));
-    DebugVar(U_rand_real);
-
-    Info << endl << endl;
-
-    tensor T_rand_imag
-    (
-        0.8668024539947510, 0.1664607226848602, 0.8925783634185791,
-        0.9126510620117188, 0.7408077120780945, 0.1499115079641342,
-        0.0936608463525772, 0.7615650296211243, 0.8953040242195129
-    );
-    DebugVar(T_rand_imag);
-    vector L_rand_imag(eigenValues(T_rand_imag));
-    DebugVar(L_rand_imag);
-    tensor U_rand_imag(eigenVectors(T_rand_imag));
-    DebugVar(U_rand_imag);
-
-    Info << endl << endl;
-
-    tensor T_rand_symm
-    (
-        1.9999992847442627, 1.3076051771640778, 1.3121289014816284,
-        1.3076051771640778, 0.0738890841603279, 0.2463847398757935,
-        1.3121289014816284, 0.2463847398757935, 0.4082727730274200
-    );
-    DebugVar(T_rand_symm);
-    vector L_rand_symm(eigenValues(T_rand_symm));
-    DebugVar(L_rand_symm);
-    tensor U_rand_symm(eigenVectors(T_rand_symm));
-    DebugVar(U_rand_symm);
-
-    Info << endl << endl;
-
-    symmTensor T_rand_Symm
-    (
-        1.9999992847442627, 1.3076051771640778, 1.3121289014816284,
-                            0.0738890841603279, 0.2463847398757935,
-                                                0.4082727730274200
-    );
-    DebugVar(T_rand_Symm);
-    vector L_rand_Symm(eigenValues(T_rand_Symm));
-    DebugVar(L_rand_Symm);
-    tensor U_rand_Symm(eigenVectors(T_rand_Symm));
-    DebugVar(U_rand_Symm);
-
-    Info << endl << endl;
-
-    tensor T_rand_diag
-    (
-        0.8668024539947510, 0, 0,
-        0, 0.7408077120780945, 0,
-        0, 0, 0.8953040242195129
-    );
-    DebugVar(T_rand_diag);
-    vector L_rand_diag(eigenValues(T_rand_diag));
-    DebugVar(L_rand_diag);
-    tensor U_rand_diag(eigenVectors(T_rand_diag));
-    DebugVar(U_rand_diag);
-
-    Info << endl << endl;
-
-    tensor T_repeated
-    (
-        0, 1, 1,
-        1, 0, 1,
-        1, 1, 0
-    );
-    DebugVar(T_repeated);
-    vector L_repeated(eigenValues(T_repeated));
-    DebugVar(L_repeated);
-    tensor U_repeated(eigenVectors(T_repeated));
-    DebugVar(U_repeated);
-
-    Info << endl << endl;
-
-    tensor T_repeated_zero
-    (
-        1, 1, 1,
-        1, 1, 1,
-        1, 1, 1
-    );
-    DebugVar(T_repeated_zero);
-    vector L_repeated_zero(eigenValues(T_repeated_zero));
-    DebugVar(L_repeated_zero);
-    tensor U_repeated_zero(eigenVectors(T_repeated_zero));
-    DebugVar(U_repeated_zero);
-
-    Info << endl << endl;
-
-    tensor T_triple
-    (
-        2, 0, 0,
-        0, 2, 0,
-        0, 0, 2
-    );
-    DebugVar(T_triple);
-    vector L_triple(eigenValues(T_triple));
-    DebugVar(L_triple);
-    tensor U_triple(eigenVectors(T_triple));
-    DebugVar(U_triple);
-    */
-
-    return 0;
-}

src/OpenFOAM/primitives/DiagTensor/DiagTensor.H

@@ -28,11 +28,11 @@ Class
     Foam::DiagTensor
 
 Description
-    Templated 3D DiagTensor derived from VectorSpace.
+    A templated (3 x 3) diagonal tensor of objects of \<T\>, effectively
+    containing 3 elements, derived from VectorSpace.
 
-    Adding construction from 3 components, element access using xx(), yy()
-    and zz() member functions and the inner-product (dot-product) and
-    outer-product operators.
+See also
+    Test-DiagTensor.C
 
 SourceFiles
     DiagTensorI.H

src/OpenFOAM/primitives/DiagTensor/DiagTensorI.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -116,6 +117,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
+//- Sum of a DiagTensor and a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator+(const DiagTensor<Cmpt>& dt1, const Tensor<Cmpt>& t2)
@@ -129,6 +131,7 @@ operator+(const DiagTensor<Cmpt>& dt1, const Tensor<Cmpt>& t2)
 }
 
 
+//- Sum of a Tensor and a DiagTensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator+(const Tensor<Cmpt>& t1, const DiagTensor<Cmpt>& dt2)
@@ -142,6 +145,7 @@ operator+(const Tensor<Cmpt>& t1, const DiagTensor<Cmpt>& dt2)
 }
 
 
+//- Subtract a Tensor from a DiagTensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator-(const DiagTensor<Cmpt>& dt1, const Tensor<Cmpt>& t2)
@@ -155,6 +159,7 @@ operator-(const DiagTensor<Cmpt>& dt1, const Tensor<Cmpt>& t2)
 }
 
 
+//- Subtract a DiagTensor from a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator-(const Tensor<Cmpt>& t1, const DiagTensor<Cmpt>& dt2)
@@ -168,7 +173,7 @@ operator-(const Tensor<Cmpt>& t1, const DiagTensor<Cmpt>& dt2)
 }
 
 
-//- Inner-product between two diagonal tensors
+//- Inner-product of a DiagTensor and a DiagTensor
 template<class Cmpt>
 inline DiagTensor<Cmpt>
 operator&(const DiagTensor<Cmpt>& dt1, const DiagTensor<Cmpt>& dt2)
@@ -182,7 +187,7 @@ operator&(const DiagTensor<Cmpt>& dt1, const DiagTensor<Cmpt>& dt2)
 }
 
 
-//- Inner-product between a diagonal tensor and a tensor
+//- Inner-product of a DiagTensor and a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator&(const DiagTensor<Cmpt>& dt1, const Tensor<Cmpt>& t2)
@@ -204,7 +209,7 @@ operator&(const DiagTensor<Cmpt>& dt1, const Tensor<Cmpt>& t2)
 }
 
 
-//- Inner-product between a tensor and a diagonal tensor
+//- Inner-product of a Tensor and a DiagTensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator&(const Tensor<Cmpt>& t1, const DiagTensor<Cmpt>& dt2)
@@ -226,7 +231,7 @@ operator&(const Tensor<Cmpt>& t1, const DiagTensor<Cmpt>& dt2)
 }
 
 
-//- Inner-product between a diagonal tensor and a vector
+//- Inner-product of a DiagTensor and a Vector
 template<class Cmpt>
 inline Vector<Cmpt>
 operator&(const DiagTensor<Cmpt>& dt, const Vector<Cmpt>& v)
@@ -240,7 +245,7 @@ operator&(const DiagTensor<Cmpt>& dt, const Vector<Cmpt>& v)
 }
 
 
-//- Inner-product between a vector and a diagonal tensor
+//- Inner-product of a Vector and a DiagTensor
 template<class Cmpt>
 inline Vector<Cmpt>
 operator&(const Vector<Cmpt>& v, const DiagTensor<Cmpt>& dt)
@@ -254,25 +259,67 @@ operator&(const Vector<Cmpt>& v, const DiagTensor<Cmpt>& dt)
 }
 
 
-//- Division of a scalar by a diagonalTensor
+//- Division of a Cmpt by a DiagTensor
 template<class Cmpt>
 inline DiagTensor<Cmpt>
-operator/(const scalar s, const DiagTensor<Cmpt>& dt)
+operator/(const Cmpt s, const DiagTensor<Cmpt>& dt)
 {
+    #ifdef FULLDEBUG
+    if (mag(det(dt)) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "Cmpt = " << s
+            << " is not divisible by the DiagTensor due to a zero element:"
+            << "DiagTensor = " << dt
+            << abort(FatalError);
+    }
+    #endif
+
     return DiagTensor<Cmpt>(s/dt.xx(), s/dt.yy(), s/dt.zz());
 }
 
 
-//- Division of a vector by a diagonalTensor
+//- Division of a DiagTensor by a Cmpt
+template<class Cmpt>
+inline DiagTensor<Cmpt>
+operator/(const DiagTensor<Cmpt>& dt, const Cmpt s)
+{
+    #ifdef FULLDEBUG
+    if (mag(s) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "DiagTensor = " << dt
+            << " is not divisible due to a zero value in Cmpt:"
+            << "Cmpt = " << s
+            << abort(FatalError);
+    }
+    #endif
+
+    return DiagTensor<Cmpt>(dt.xx()/s, dt.yy()/s, dt.zz()/s);
+}
+
+
+//- Division of a Vector by a DiagTensor
 template<class Cmpt>
 inline Vector<Cmpt>
 operator/(const Vector<Cmpt> v, const DiagTensor<Cmpt>& dt)
 {
+    #ifdef FULLDEBUG
+    if (mag(det(dt)) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "Vector = " << v
+            << " is not divisible by the DiagTensor due to a zero element:"
+            << "DiagTensor = " << dt
+            << abort(FatalError);
+    }
+    #endif
+
     return Vector<Cmpt>(v.x()/dt.xx(), v.y()/dt.yy(), v.z()/dt.zz());
 }
 
 
-//- Return the trace of a diagonal tensor
+//- Return the trace of a DiagTensor
 template<class Cmpt>
 inline Cmpt tr(const DiagTensor<Cmpt>& dt)
 {
@@ -280,42 +327,44 @@ inline Cmpt tr(const DiagTensor<Cmpt>& dt)
 }
 
 
-//- Return the spherical part of a diagonal tensor
+//- Return the spherical part of a DiagTensor as a SphericalTensor
 template<class Cmpt>
 inline SphericalTensor<Cmpt> sph(const DiagTensor<Cmpt>& dt)
 {
     return SphericalTensor<Cmpt>
     (
-        1.0/3.0*tr(dt)
+        (1.0/3.0)*tr(dt)
     );
 }
 
 
-//- Return the determinant of a diagonal tensor
+//- Return the determinant of a DiagTensor
 template<class Cmpt>
-inline Cmpt det(const DiagTensor<Cmpt>& t)
+inline Cmpt det(const DiagTensor<Cmpt>& dt)
 {
-    return t.xx()*t.yy()*t.zz();
+    return dt.xx()*dt.yy()*dt.zz();
 }
 
 
-//- Return the inverse of a diagonal tensor
+//- Return the inverse of a DiagTensor as a DiagTensor
 template<class Cmpt>
 inline DiagTensor<Cmpt> inv(const DiagTensor<Cmpt>& dt)
 {
-    return DiagTensor<Cmpt>(1.0/dt.xx(), 1.0/dt.yy(), 1.0/dt.zz());
+    #ifdef FULLDEBUG
+    if (mag(det(dt)) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "DiagTensor is not invertible due to the zero determinant:"
+            << "det(DiagTensor) = " << det(dt)
+            << abort(FatalError);
+    }
+    #endif
+
+    return DiagTensor<Cmpt>(1/dt.xx(), 1/dt.yy(), 1/dt.zz());
 }
 
 
-//- Return the diagonal of a symmetric tensor as a diagonal tensor
-template<class Cmpt>
-inline DiagTensor<Cmpt> diag(const SymmTensor<Cmpt>& t)
-{
-    return DiagTensor<Cmpt>(t.xx(), t.yy(), t.zz());
-}
-
-
-//- Return the diagonal of a tensor as a diagonal tensor
+//- Return the diagonal of a Tensor as a DiagTensor
 template<class Cmpt>
 inline DiagTensor<Cmpt> diag(const Tensor<Cmpt>& t)
 {
@@ -323,6 +372,14 @@ inline DiagTensor<Cmpt> diag(const Tensor<Cmpt>& t)
 }
 
 
+//- Return the diagonal of a SymmTensor as a DiagTensor
+template<class Cmpt>
+inline DiagTensor<Cmpt> diag(const SymmTensor<Cmpt>& st)
+{
+    return DiagTensor<Cmpt>(st.xx(), st.yy(), st.zz());
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/primitives/DiagTensor/diagTensor/diagTensor.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,10 @@ Typedef
     Foam::diagTensor
 
 Description
-    A scalar version of the templated DiagTensor
+    DiagTensor of scalars, i.e. DiagTensor<scalar>.
+
+See also
+    Test-DiagTensor.C
 
 SourceFiles
     diagTensor.C

src/OpenFOAM/primitives/SphericalTensor/SphericalTensor.H

@@ -28,9 +28,11 @@ Class
     Foam::SphericalTensor
 
 Description
-    Templated 3D SphericalTensor derived from VectorSpace adding construction
-    from 1 component, element access using th ii() member function and the
-    inner-product (dot-product) and outer-product operators.
+    A templated (3 x 3) diagonal tensor of objects of \<T\>, effectively
+    containing 1 element, derived from VectorSpace.
+
+See also
+    Test-SphericalTensor.C
 
 SourceFiles
     SphericalTensorI.H
@@ -121,7 +123,7 @@ public:
             inline Cmpt& ii();
 
 
-        //- Transpose (no-op)
+        //- Return non-Hermitian transpose (no-op)
         inline const SphericalTensor<Cmpt>& T() const;
 };
 

src/OpenFOAM/primitives/SphericalTensor/SphericalTensorI.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -92,7 +93,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
-//- Inner-product between two spherical tensors
+//- Inner-product of a SphericalTensor and a SphericalTensor
 template<class Cmpt>
 inline SphericalTensor<Cmpt>
 operator&(const SphericalTensor<Cmpt>& st1, const SphericalTensor<Cmpt>& st2)
@@ -101,7 +102,7 @@ operator&(const SphericalTensor<Cmpt>& st1, const SphericalTensor<Cmpt>& st2)
 }
 
 
-//- Inner-product between a spherical tensor and a vector
+//- Inner-product of a SphericalTensor and a Vector
 template<class Cmpt>
 inline Vector<Cmpt>
 operator&(const SphericalTensor<Cmpt>& st, const Vector<Cmpt>& v)
@@ -115,7 +116,7 @@ operator&(const SphericalTensor<Cmpt>& st, const Vector<Cmpt>& v)
 }
 
 
-//- Inner-product between a vector and a spherical tensor
+//- Inner-product of a Vector and a SphericalTensor
 template<class Cmpt>
 inline Vector<Cmpt>
 operator&(const Vector<Cmpt>& v, const SphericalTensor<Cmpt>& st)
@@ -129,7 +130,7 @@ operator&(const Vector<Cmpt>& v, const SphericalTensor<Cmpt>& st)
 }
 
 
-//- Double-dot-product between a spherical tensor and a spherical tensor
+//- Double-inner-product of a SphericalTensor and a SphericalTensor
 template<class Cmpt>
 inline Cmpt
 operator&&(const SphericalTensor<Cmpt>& st1, const SphericalTensor<Cmpt>& st2)
@@ -138,22 +139,55 @@ operator&&(const SphericalTensor<Cmpt>& st1, const SphericalTensor<Cmpt>& st2)
 }
 
 
-//- Division of a scalar by a sphericalTensor
+//- Division of a Cmpt by a SphericalTensor
 template<class Cmpt>
 inline SphericalTensor<Cmpt>
-operator/(const scalar s, const SphericalTensor<Cmpt>& st)
+operator/(const Cmpt s, const SphericalTensor<Cmpt>& st)
 {
+    #ifdef FULLDEBUG
+    if (mag(st.ii()) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "Cmpt = " << s
+            << " is not divisible due to a zero element in SphericalTensor:"
+            << "SphericalTensor = " << st
+            << abort(FatalError);
+    }
+    #endif
+
     return SphericalTensor<Cmpt>(s/st.ii());
 }
 
 
+//- Division of a SphericalTensor by a Cmpt
 template<class Cmpt>
-inline Cmpt magSqr(const SphericalTensor<Cmpt>& st)
+inline SphericalTensor<Cmpt>
+operator/(const SphericalTensor<Cmpt>& st, const Cmpt s)
 {
-    return 3*magSqr(st.ii());
+    #ifdef FULLDEBUG
+    if (mag(s) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "SphericalTensor = " << st
+            << " is not divisible due to a zero value in Cmpt:"
+            << "Cmpt = " << s
+            << abort(FatalError);
+    }
+    #endif
+
+    return SphericalTensor<Cmpt>(st.ii()/s);
 }
 
 
+//- Return the square of Frobenius norm of a SphericalTensor as a Cmpt
+template<class Cmpt>
+inline Cmpt magSqr(const SphericalTensor<Cmpt>& st)
+{
+    return Cmpt(3*mag(st.ii()*st.ii()));
+}
+
+
+//- Return the max component of a SphericalTensor
 template<class Cmpt>
 inline Cmpt cmptMax(const SphericalTensor<Cmpt>& st)
 {
@@ -161,6 +195,7 @@ inline Cmpt cmptMax(const SphericalTensor<Cmpt>& st)
 }
 
 
+//- Return the min component of a SphericalTensor
 template<class Cmpt>
 inline Cmpt cmptMin(const SphericalTensor<Cmpt>& st)
 {
@@ -168,6 +203,7 @@ inline Cmpt cmptMin(const SphericalTensor<Cmpt>& st)
 }
 
 
+//- Return the sum of components of a SphericalTensor
 template<class Cmpt>
 inline Cmpt cmptSum(const SphericalTensor<Cmpt>& st)
 {
@@ -175,6 +211,7 @@ inline Cmpt cmptSum(const SphericalTensor<Cmpt>& st)
 }
 
 
+//- Return the arithmetic average of components of a SphericalTensor
 template<class Cmpt>
 inline Cmpt cmptAv(const SphericalTensor<Cmpt>& st)
 {
@@ -182,7 +219,7 @@ inline Cmpt cmptAv(const SphericalTensor<Cmpt>& st)
 }
 
 
-//- Return the trace of a spherical tensor
+//- Return the trace of a SphericalTensor
 template<class Cmpt>
 inline Cmpt tr(const SphericalTensor<Cmpt>& st)
 {
@@ -190,7 +227,7 @@ inline Cmpt tr(const SphericalTensor<Cmpt>& st)
 }
 
 
-//- Return the spherical part of a spherical tensor, i.e. itself
+//- Return the spherical part of a SphericalTensor, i.e. itself
 template<class Cmpt>
 inline SphericalTensor<Cmpt> sph(const SphericalTensor<Cmpt>& st)
 {
@@ -198,7 +235,7 @@ inline SphericalTensor<Cmpt> sph(const SphericalTensor<Cmpt>& st)
 }
 
 
-//- Return the determinant of a spherical tensor
+//- Return the determinant of a SphericalTensor
 template<class Cmpt>
 inline Cmpt det(const SphericalTensor<Cmpt>& st)
 {
@@ -206,14 +243,26 @@ inline Cmpt det(const SphericalTensor<Cmpt>& st)
 }
 
 
-//- Return the inverse of a spherical tensor
+//- Return the inverse of a SphericalTensor
 template<class Cmpt>
 inline SphericalTensor<Cmpt> inv(const SphericalTensor<Cmpt>& st)
 {
-    return SphericalTensor<Cmpt>(1.0/st.ii());
+    #ifdef FULLDEBUG
+    if (mag(st.ii()) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "SphericalTensor is not invertible due to the zero determinant:"
+            << "det(SphericalTensor) = " << det(st)
+            << abort(FatalError);
+    }
+    #endif
+
+    return SphericalTensor<Cmpt>(1/st.ii());
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 template<class Cmpt>
 class outerProduct<SphericalTensor<Cmpt>, Cmpt>
 {

src/OpenFOAM/primitives/SphericalTensor/labelSphericalTensor/labelSphericalTensor.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,10 @@ Typedef
     Foam::labelSphericalTensor
 
 Description
-    SphericalTensor of labels.
+    SphericalTensor of labels, i.e. SphericalTensor<label>.
+
+See also
+    Test-SphericalTensor.C
 
 SourceFiles
     labelSphericalTensor.C

src/OpenFOAM/primitives/SphericalTensor/sphericalTensor/sphericalTensor.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,10 @@ Typedef
     Foam::sphericalTensor
 
 Description
-    SphericalTensor of scalars.
+    SphericalTensor of scalars, i.e. SphericalTensor<scalar>.
+
+See also
+    Test-SphericalTensor.C
 
 SourceFiles
     sphericalTensor.C

src/OpenFOAM/primitives/SphericalTensor2D/SphericalTensor2D.H

@@ -28,9 +28,11 @@ Class
     Foam::SphericalTensor2D
 
 Description
-    Templated 2D sphericalTensor derived from VectorSpace adding construction
-    from 1 component, element access using ii() member function and the
-    inner-product (dot-product) and outer-product operators.
+    A templated (2 x 2) diagonal tensor of objects of \<T\>, effectively
+    containing 1 element, derived from VectorSpace.
+
+See also
+    Test-SphericalTensor2D.C
 
 SourceFiles
     SphericalTensor2DI.H

src/OpenFOAM/primitives/SphericalTensor2D/SphericalTensor2DI.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,7 +84,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
-//- Inner-product between two spherical tensors
+//- Inner-product of a SphericalTensor2D and a SphericalTensor2D
 template<class Cmpt>
 inline SphericalTensor2D<Cmpt>
 operator&
@@ -96,7 +97,7 @@ operator&
 }
 
 
-//- Inner-product between a spherical tensor and a vector
+//- Inner-product of a SphericalTensor2D and a Vector2D
 template<class Cmpt>
 inline Vector2D<Cmpt>
 operator&(const SphericalTensor2D<Cmpt>& st, const Vector2D<Cmpt>& v)
@@ -109,7 +110,7 @@ operator&(const SphericalTensor2D<Cmpt>& st, const Vector2D<Cmpt>& v)
 }
 
 
-//- Inner-product between a vector and a spherical tensor
+//- Inner-product of a Vector2D and a SphericalTensor2D
 template<class Cmpt>
 inline Vector2D<Cmpt>
 operator&(const Vector2D<Cmpt>& v, const SphericalTensor2D<Cmpt>& st)
@@ -122,16 +123,47 @@ operator&(const Vector2D<Cmpt>& v, const SphericalTensor2D<Cmpt>& st)
 }
 
 
-//- Division of a scalar by a sphericalTensor2D
+//- Division of a Cmpt by a SphericalTensor2D
 template<class Cmpt>
 inline SphericalTensor2D<Cmpt>
-operator/(const scalar s, const SphericalTensor2D<Cmpt>& st)
+operator/(const Cmpt s, const SphericalTensor2D<Cmpt>& st)
 {
+    #ifdef FULLDEBUG
+    if (mag(st.ii()) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "Cmpt = " << s
+            << " is not divisible due to a zero element in SphericalTensor2D:"
+            << "SphericalTensor2D = " << st
+            << abort(FatalError);
+    }
+    #endif
+
     return SphericalTensor2D<Cmpt>(s/st.ii());
 }
 
 
-//- Return the trace of a spherical tensor
+//- Division of a SphericalTensor2D by a Cmpt
+template<class Cmpt>
+inline SphericalTensor2D<Cmpt>
+operator/(const SphericalTensor2D<Cmpt>& st, const Cmpt s)
+{
+    #ifdef FULLDEBUG
+    if (mag(s) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "SphericalTensor2D = " << st
+            << " is not divisible due to a zero value in Cmpt:"
+            << "Cmpt = " << s
+            << abort(FatalError);
+    }
+    #endif
+
+    return SphericalTensor2D<Cmpt>(st.ii()/s);
+}
+
+
+//- Return the trace of a SphericalTensor2D
 template<class Cmpt>
 inline Cmpt tr(const SphericalTensor2D<Cmpt>& st)
 {
@@ -139,7 +171,7 @@ inline Cmpt tr(const SphericalTensor2D<Cmpt>& st)
 }
 
 
-//- Return the spherical part of a spherical tensor, i.e. itself
+//- Return the spherical part of a SphericalTensor2D, i.e. itself
 template<class Cmpt>
 inline SphericalTensor2D<Cmpt> sph(const SphericalTensor2D<Cmpt>& st)
 {
@@ -147,7 +179,7 @@ inline SphericalTensor2D<Cmpt> sph(const SphericalTensor2D<Cmpt>& st)
 }
 
 
-//- Return the determinant of a spherical tensor
+//- Return the determinant of a SphericalTensor2D
 template<class Cmpt>
 inline Cmpt det(const SphericalTensor2D<Cmpt>& st)
 {
@@ -155,14 +187,26 @@ inline Cmpt det(const SphericalTensor2D<Cmpt>& st)
 }
 
 
-//- Return the inverse of a symmetric tensor
+//- Return the inverse of a SphericalTensor2D
 template<class Cmpt>
 inline SphericalTensor2D<Cmpt> inv(const SphericalTensor2D<Cmpt>& st)
 {
-    return SphericalTensor2D<Cmpt>(1.0/st.ii());
+    #ifdef FULLDEBUG
+    if (mag(st.ii()) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "SphericalTensor2D is not invertible due to zero determinant:"
+            << "det(SphericalTensor2D) = " << det(st)
+            << abort(FatalError);
+    }
+    #endif
+
+    return SphericalTensor2D<Cmpt>(1/st.ii());
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 template<class Cmpt>
 class outerProduct<SphericalTensor2D<Cmpt>, Cmpt>
 {

src/OpenFOAM/primitives/SphericalTensor2D/sphericalTensor2D/sphericalTensor2D.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,10 @@ Typedef
     Foam::sphericalTensor2D
 
 Description
-    SphericalTensor2D of scalars.
+    SphericalTensor2D of scalars, i.e. SphericalTensor2D<scalar>.
+
+See also
+    Test-SphericalTensor2D.C
 
 SourceFiles
     sphericalTensor2D.C

src/OpenFOAM/primitives/SymmTensor/SymmTensor.H

@@ -28,10 +28,11 @@ Class
     Foam::SymmTensor
 
 Description
-    Templated 3D symmetric tensor derived from VectorSpace adding construction
-    from 6 components, element access using xx(), xy() etc. member functions
-    and the inner-product (dot-product) and outer-product of two Vectors
-    (tensor-product) operators.
+    A templated (3 x 3) symmetric tensor of objects of \<T\>, effectively
+    containing 6 elements, derived from VectorSpace.
+
+See also
+    Test-SymmTensor.C
 
 SourceFiles
     SymmTensorI.H
@@ -125,19 +126,25 @@ public:
             inline const Cmpt& xx() const;
             inline const Cmpt& xy() const;
             inline const Cmpt& xz() const;
+            inline const Cmpt& yx() const;
             inline const Cmpt& yy() const;
             inline const Cmpt& yz() const;
+            inline const Cmpt& zx() const;
+            inline const Cmpt& zy() const;
             inline const Cmpt& zz() const;
 
             inline Cmpt& xx();
             inline Cmpt& xy();
             inline Cmpt& xz();
+            inline Cmpt& yx();
             inline Cmpt& yy();
             inline Cmpt& yz();
+            inline Cmpt& zx();
+            inline Cmpt& zy();
             inline Cmpt& zz();
 
 
-        // Diagonal access.
+        // Diagonal access and manipulation
 
             //- Extract the diagonal as a vector
             inline Vector<Cmpt> diag() const;
@@ -148,7 +155,7 @@ public:
 
     // Tensor Operations
 
-        //- Transpose
+        //- Return non-Hermitian transpose
         inline const SymmTensor<Cmpt>& T() const;
 
 

src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H

@@ -51,8 +51,8 @@ inline Foam::SymmTensor<Cmpt>::SymmTensor
 template<class Cmpt>
 inline Foam::SymmTensor<Cmpt>::SymmTensor(const SphericalTensor<Cmpt>& st)
 {
-    this->v_[XX] = st.ii(); this->v_[XY] = 0; this->v_[XZ] = 0;
-    this->v_[YY] = st.ii(); this->v_[YZ] = 0;
+    this->v_[XX] = st.ii(); this->v_[XY] = Zero; this->v_[XZ] = Zero;
+    this->v_[YY] = st.ii(); this->v_[YZ] = Zero;
     this->v_[ZZ] = st.ii();
 }
 
@@ -98,6 +98,12 @@ inline const Cmpt& Foam::SymmTensor<Cmpt>::xz() const
     return this->v_[XZ];
 }
 
+template<class Cmpt>
+inline const Cmpt& Foam::SymmTensor<Cmpt>::yx() const
+{
+    return this->v_[XY];
+}
+
 template<class Cmpt>
 inline const Cmpt& Foam::SymmTensor<Cmpt>::yy() const
 {
@@ -110,6 +116,18 @@ inline const Cmpt& Foam::SymmTensor<Cmpt>::yz() const
     return this->v_[YZ];
 }
 
+template<class Cmpt>
+inline const Cmpt& Foam::SymmTensor<Cmpt>::zx() const
+{
+    return this->v_[XZ];
+}
+
+template<class Cmpt>
+inline const Cmpt& Foam::SymmTensor<Cmpt>::zy() const
+{
+    return this->v_[YZ];
+}
+
 template<class Cmpt>
 inline const Cmpt& Foam::SymmTensor<Cmpt>::zz() const
 {
@@ -135,6 +153,12 @@ inline Cmpt& Foam::SymmTensor<Cmpt>::xz()
     return this->v_[XZ];
 }
 
+template<class Cmpt>
+inline Cmpt& Foam::SymmTensor<Cmpt>::yx()
+{
+    return this->v_[XY];
+}
+
 template<class Cmpt>
 inline Cmpt& Foam::SymmTensor<Cmpt>::yy()
 {
@@ -147,6 +171,18 @@ inline Cmpt& Foam::SymmTensor<Cmpt>::yz()
     return this->v_[YZ];
 }
 
+template<class Cmpt>
+inline Cmpt& Foam::SymmTensor<Cmpt>::zx()
+{
+    return this->v_[XZ];
+}
+
+template<class Cmpt>
+inline Cmpt& Foam::SymmTensor<Cmpt>::zy()
+{
+    return this->v_[YZ];
+}
+
 template<class Cmpt>
 inline Cmpt& Foam::SymmTensor<Cmpt>::zz()
 {
@@ -180,8 +216,8 @@ inline const Foam::SymmTensor<Cmpt>& Foam::SymmTensor<Cmpt>::T() const
 template<class Cmpt>
 inline void Foam::SymmTensor<Cmpt>::operator=(const SphericalTensor<Cmpt>& st)
 {
-    this->v_[XX] = st.ii(); this->v_[XY] = 0;       this->v_[XZ] = 0;
-                            this->v_[YY] = st.ii(); this->v_[YZ] = 0;
+    this->v_[XX] = st.ii(); this->v_[XY] = Zero;    this->v_[XZ] = Zero;
+                            this->v_[YY] = st.ii(); this->v_[YZ] = Zero;
                                                     this->v_[ZZ] = st.ii();
 }
 
@@ -193,7 +229,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
-//- Hodge Dual operator (tensor -> vector)
+//- Return the Hodge dual of a SymmTensor as a Vector
 template<class Cmpt>
 inline Vector<Cmpt> operator*(const SymmTensor<Cmpt>& st)
 {
@@ -201,7 +237,7 @@ inline Vector<Cmpt> operator*(const SymmTensor<Cmpt>& st)
 }
 
 
-//- Inner-product between two symmetric tensors
+//- Inner-product of a SymmTensor and a SymmTensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator&(const SymmTensor<Cmpt>& st1, const SymmTensor<Cmpt>& st2)
@@ -223,7 +259,7 @@ operator&(const SymmTensor<Cmpt>& st1, const SymmTensor<Cmpt>& st2)
 }
 
 
-//- Double-dot-product between a symmetric tensor and a symmetric tensor
+//- Double-inner-product of a SymmTensor and a SymmTensor
 template<class Cmpt>
 inline Cmpt
 operator&&(const SymmTensor<Cmpt>& st1, const SymmTensor<Cmpt>& st2)
@@ -237,7 +273,7 @@ operator&&(const SymmTensor<Cmpt>& st1, const SymmTensor<Cmpt>& st2)
 }
 
 
-//- Inner-product between a symmetric tensor and a vector
+//- Inner-product of a SymmTensor and a Vector
 template<class Cmpt>
 inline Vector<Cmpt>
 operator&(const SymmTensor<Cmpt>& st, const Vector<Cmpt>& v)
@@ -251,7 +287,7 @@ operator&(const SymmTensor<Cmpt>& st, const Vector<Cmpt>& v)
 }
 
 
-//- Inner-product between a vector and a symmetric tensor
+//- Inner-product of a Vector and a SymmTensor
 template<class Cmpt>
 inline Vector<Cmpt>
 operator&(const Vector<Cmpt>& v, const SymmTensor<Cmpt>& st)
@@ -265,7 +301,7 @@ operator&(const Vector<Cmpt>& v, const SymmTensor<Cmpt>& st)
 }
 
 
-//- Inner-sqr of a symmetric tensor
+//- Return the inner-product of a SymmTensor with itself
 template<class Cmpt>
 inline SymmTensor<Cmpt>
 innerSqr(const SymmTensor<Cmpt>& st)
@@ -284,19 +320,20 @@ innerSqr(const SymmTensor<Cmpt>& st)
 }
 
 
+//- Return the square of Frobenius norm of a SymmTensor as a Cmpt
 template<class Cmpt>
 inline Cmpt magSqr(const SymmTensor<Cmpt>& st)
 {
-    return
+    return Cmpt
     (
-        magSqr(st.xx()) + 2*magSqr(st.xy()) + 2*magSqr(st.xz())
-                        +   magSqr(st.yy()) + 2*magSqr(st.yz())
-                                            +   magSqr(st.zz())
+        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())
     );
 }
 
 
-//- Return the trace of a symmetric tensor
+//- Return the trace of a SymmTensor
 template<class Cmpt>
 inline Cmpt tr(const SymmTensor<Cmpt>& st)
 {
@@ -304,7 +341,7 @@ inline Cmpt tr(const SymmTensor<Cmpt>& st)
 }
 
 
-//- Return the spherical part of a symmetric tensor
+//- Return the spherical part of a SymmTensor
 template<class Cmpt>
 inline SphericalTensor<Cmpt> sph(const SymmTensor<Cmpt>& st)
 {
@@ -315,7 +352,7 @@ inline SphericalTensor<Cmpt> sph(const SymmTensor<Cmpt>& st)
 }
 
 
-//- Return the symmetric part of a symmetric tensor, i.e. itself
+//- Return the symmetric part of a SymmTensor, i.e. itself
 template<class Cmpt>
 inline const SymmTensor<Cmpt>& symm(const SymmTensor<Cmpt>& st)
 {
@@ -323,7 +360,7 @@ inline const SymmTensor<Cmpt>& symm(const SymmTensor<Cmpt>& st)
 }
 
 
-//- Return twice the symmetric part of a symmetric tensor
+//- Return twice the symmetric part of a SymmTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt> twoSymm(const SymmTensor<Cmpt>& st)
 {
@@ -331,23 +368,23 @@ inline SymmTensor<Cmpt> twoSymm(const SymmTensor<Cmpt>& st)
 }
 
 
-//- Return the deviatoric part of a symmetric tensor
+//- Return the deviatoric part of a SymmTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt> dev(const SymmTensor<Cmpt>& st)
 {
-    return st - SphericalTensor<Cmpt>::oneThirdI*tr(st);
+    return st - sph(st);
 }
 
 
-//- Return the deviatoric part of a symmetric tensor
+//- Return the two-third deviatoric part of a SymmTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt> dev2(const SymmTensor<Cmpt>& st)
 {
-    return st - SphericalTensor<Cmpt>::twoThirdsI*tr(st);
+    return st - 2*sph(st);
 }
 
 
-//- Return the determinant of a symmetric tensor
+//- Return the determinant of a SymmTensor
 template<class Cmpt>
 inline Cmpt det(const SymmTensor<Cmpt>& st)
 {
@@ -360,7 +397,7 @@ inline Cmpt det(const SymmTensor<Cmpt>& st)
 }
 
 
-//- Return the cofactor symmetric tensor of a symmetric tensor
+//- Return the cofactor of a SymmTensor as a SymmTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt> cof(const SymmTensor<Cmpt>& st)
 {
@@ -378,25 +415,25 @@ inline SymmTensor<Cmpt> cof(const SymmTensor<Cmpt>& st)
 }
 
 
-//- Return the inverse of a symmetric tensor give the determinant
+//- Return the inverse of a SymmTensor by using the given determinant
 template<class Cmpt>
 inline SymmTensor<Cmpt> inv(const SymmTensor<Cmpt>& st, const Cmpt detst)
 {
-    return SymmTensor<Cmpt>
-    (
-        st.yy()*st.zz() - st.yz()*st.yz(),
-        st.xz()*st.yz() - st.xy()*st.zz(),
-        st.xy()*st.yz() - st.xz()*st.yy(),
+    #ifdef FULLDEBUG
+    if (mag(detst) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "SymmTensor is not invertible due to the zero determinant:"
+            << "det(symmTensor) = " << mag(detst)
+            << abort(FatalError);
+    }
+    #endif
 
-        st.xx()*st.zz() - st.xz()*st.xz(),
-        st.xy()*st.xz() - st.xx()*st.yz(),
-
-        st.xx()*st.yy() - st.xy()*st.xy()
-    )/detst;
+    return cof(st).T()/detst;
 }
 
 
-//- Return the inverse of a symmetric tensor
+//- Return the inverse of a SymmTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt> inv(const SymmTensor<Cmpt>& st)
 {
@@ -404,7 +441,7 @@ inline SymmTensor<Cmpt> inv(const SymmTensor<Cmpt>& st)
 }
 
 
-//- Return the 1st invariant of a symmetric tensor
+//- Return the 1st invariant of a SymmTensor
 template<class Cmpt>
 inline Cmpt invariantI(const SymmTensor<Cmpt>& st)
 {
@@ -412,19 +449,19 @@ inline Cmpt invariantI(const SymmTensor<Cmpt>& st)
 }
 
 
-//- Return the 2nd invariant of a symmetric tensor
+//- Return the 2nd invariant of a SymmTensor
 template<class Cmpt>
 inline Cmpt invariantII(const SymmTensor<Cmpt>& st)
 {
     return
     (
         st.xx()*st.yy() + st.yy()*st.zz() + st.xx()*st.zz()
-      - sqr(st.xy()) - sqr(st.yz()) - sqr(st.xz())
+      - st.xy()*st.xy() - st.yz()*st.yz() - st.xz()*st.xz()
     );
 }
 
 
-//- Return the 3rd invariant of a symmetric tensor
+//- Return the 3rd invariant of a SymmTensor
 template<class Cmpt>
 inline Cmpt invariantIII(const SymmTensor<Cmpt>& st)
 {
@@ -432,6 +469,7 @@ inline Cmpt invariantIII(const SymmTensor<Cmpt>& st)
 }
 
 
+//- Sum of a SphericalTensor and a SymmTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt>
 operator+(const SphericalTensor<Cmpt>& spt1, const SymmTensor<Cmpt>& st2)
@@ -445,6 +483,7 @@ operator+(const SphericalTensor<Cmpt>& spt1, const SymmTensor<Cmpt>& st2)
 }
 
 
+//- Sum of a SymmTensor and a SphericalTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt>
 operator+(const SymmTensor<Cmpt>& st1, const SphericalTensor<Cmpt>& spt2)
@@ -458,6 +497,7 @@ operator+(const SymmTensor<Cmpt>& st1, const SphericalTensor<Cmpt>& spt2)
 }
 
 
+//- Subtract a SymmTensor from a SphericalTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt>
 operator-(const SphericalTensor<Cmpt>& spt1, const SymmTensor<Cmpt>& st2)
@@ -471,6 +511,7 @@ operator-(const SphericalTensor<Cmpt>& spt1, const SymmTensor<Cmpt>& st2)
 }
 
 
+//- Subtract a SphericalTensor from a SymmTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt>
 operator-(const SymmTensor<Cmpt>& st1, const SphericalTensor<Cmpt>& spt2)
@@ -484,7 +525,7 @@ operator-(const SymmTensor<Cmpt>& st1, const SphericalTensor<Cmpt>& spt2)
 }
 
 
-//- Inner-product between a spherical symmetric tensor and a symmetric tensor
+//- Inner-product of a SphericalTensor and a SymmTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt>
 operator&(const SphericalTensor<Cmpt>& spt1, const SymmTensor<Cmpt>& st2)
@@ -498,7 +539,7 @@ operator&(const SphericalTensor<Cmpt>& spt1, const SymmTensor<Cmpt>& st2)
 }
 
 
-//- Inner-product between a tensor and a spherical tensor
+//- Inner-product of a SymmTensor and a SphericalTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt>
 operator&(const SymmTensor<Cmpt>& st1, const SphericalTensor<Cmpt>& spt2)
@@ -512,24 +553,25 @@ operator&(const SymmTensor<Cmpt>& st1, const SphericalTensor<Cmpt>& spt2)
 }
 
 
-//- Double-dot-product between a spherical tensor and a symmetric tensor
+//- Double-inner-product of a SphericalTensor and a SymmTensor
 template<class Cmpt>
 inline Cmpt
 operator&&(const SphericalTensor<Cmpt>& spt1, const SymmTensor<Cmpt>& st2)
 {
-    return(spt1.ii()*st2.xx() + spt1.ii()*st2.yy() + spt1.ii()*st2.zz());
+    return (spt1.ii()*st2.xx() + spt1.ii()*st2.yy() + spt1.ii()*st2.zz());
 }
 
 
-//- Double-dot-product between a tensor and a spherical tensor
+//- Double-inner-product of a SymmTensor and a SphericalTensor
 template<class Cmpt>
 inline Cmpt
 operator&&(const SymmTensor<Cmpt>& st1, const SphericalTensor<Cmpt>& spt2)
 {
-    return(st1.xx()*spt2.ii() + st1.yy()*spt2.ii() + st1.zz()*spt2.ii());
+    return (st1.xx()*spt2.ii() + st1.yy()*spt2.ii() + st1.zz()*spt2.ii());
 }
 
 
+//- Return the square of a Vector as a SymmTensor
 template<class Cmpt>
 inline SymmTensor<Cmpt> sqr(const Vector<Cmpt>& v)
 {
@@ -542,6 +584,22 @@ inline SymmTensor<Cmpt> sqr(const Vector<Cmpt>& v)
 }
 
 
+//- Division of a SymmTensor by a Cmpt
+template<class Cmpt>
+inline SymmTensor<Cmpt>
+operator/(const SymmTensor<Cmpt>& st, const Cmpt s)
+{
+    return SymmTensor<Cmpt>
+    (
+        st.xx()/s, st.xy()/s, st.xz()/s,
+                   st.yy()/s, st.yz()/s,
+                              st.zz()/s
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 template<class Cmpt>
 class outerProduct<SymmTensor<Cmpt>, Cmpt>
 {

src/OpenFOAM/primitives/SymmTensor/labelSymmTensor/labelSymmTensor.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,10 @@ Typedef
     Foam::labelSymmTensor
 
 Description
-    SymmTensor of labels.
+    SymmTensor of labels, i.e. SymmTensor<label>.
+
+See also
+    Test-SymmTensor.C
 
 SourceFiles
     labelSymmTensor.C

src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2D.H

@@ -28,10 +28,11 @@ Class
     Foam::SymmTensor2D
 
 Description
-    Templated 2D symmetric tensor derived from VectorSpace adding construction
-    from 4 components, element access using xx(), xy() etc. member functions
-    and the inner-product (dot-product) and outer-product of two Vectors
-    (tensor-product) operators.
+    A templated (2 x 2) symmetric tensor of objects of \<T\>, effectively
+    containing 3 elements, derived from VectorSpace.
+
+See also
+    Test-SymmTensor2D.C
 
 SourceFiles
     SymmTensor2DI.H
@@ -122,14 +123,16 @@ public:
 
             inline const Cmpt& xx() const;
             inline const Cmpt& xy() const;
+            inline const Cmpt& yx() const;
             inline const Cmpt& yy() const;
 
             inline Cmpt& xx();
             inline Cmpt& xy();
+            inline Cmpt& yx();
             inline Cmpt& yy();
 
 
-        // Diagonal access.
+        // Diagonal access and manipulation
 
             //- Extract the diagonal as a vector
             inline Vector2D<Cmpt> diag() const;
@@ -138,7 +141,7 @@ public:
             inline void diag(const Vector2D<Cmpt>& v);
 
 
-        //- Transpose
+        //- Return non-Hermitian transpose
         inline const SymmTensor2D<Cmpt>& T() const;
 
 

src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2DI.H

@@ -50,7 +50,7 @@ inline Foam::SymmTensor2D<Cmpt>::SymmTensor2D
 template<class Cmpt>
 inline Foam::SymmTensor2D<Cmpt>::SymmTensor2D(const SphericalTensor2D<Cmpt>& st)
 {
-    this->v_[XX] = st.ii(); this->v_[XY] = 0;
+    this->v_[XX] = st.ii(); this->v_[XY] = Zero;
     this->v_[YY] = st.ii();
 }
 
@@ -88,6 +88,12 @@ inline const Cmpt& Foam::SymmTensor2D<Cmpt>::xy() const
     return this->v_[XY];
 }
 
+template<class Cmpt>
+inline const Cmpt& Foam::SymmTensor2D<Cmpt>::yx() const
+{
+    return this->v_[XY];
+}
+
 template<class Cmpt>
 inline const Cmpt& Foam::SymmTensor2D<Cmpt>::yy() const
 {
@@ -107,6 +113,12 @@ inline Cmpt& Foam::SymmTensor2D<Cmpt>::xy()
     return this->v_[XY];
 }
 
+template<class Cmpt>
+inline Cmpt& Foam::SymmTensor2D<Cmpt>::yx()
+{
+    return this->v_[XY];
+}
+
 template<class Cmpt>
 inline Cmpt& Foam::SymmTensor2D<Cmpt>::yy()
 {
@@ -143,7 +155,7 @@ inline void Foam::SymmTensor2D<Cmpt>::operator=
     const SphericalTensor2D<Cmpt>& st
 )
 {
-    this->v_[XX] = st.ii(); this->v_[XY] = 0;
+    this->v_[XX] = st.ii(); this->v_[XY] = Zero;
                             this->v_[YY] = st.ii();
 }
 
@@ -155,7 +167,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
-//- Inner-product between two symmetric tensors
+//- Inner-product of a SymmTensor2D and a SymmTensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator&(const SymmTensor2D<Cmpt>& st1, const SymmTensor2D<Cmpt>& st2)
@@ -171,7 +183,7 @@ operator&(const SymmTensor2D<Cmpt>& st1, const SymmTensor2D<Cmpt>& st2)
 }
 
 
-//- Double-dot-product between a symmetric tensor and a symmetric tensor
+//- Double-inner-product of a SymmTensor2D and a SymmTensor2D
 template<class Cmpt>
 inline Cmpt
 operator&&(const SymmTensor2D<Cmpt>& st1, const SymmTensor2D<Cmpt>& st2)
@@ -184,7 +196,7 @@ operator&&(const SymmTensor2D<Cmpt>& st1, const SymmTensor2D<Cmpt>& st2)
 }
 
 
-//- Inner-product between a symmetric tensor and a vector
+//- Inner-product of a SymmTensor2D and a Vector2D
 template<class Cmpt>
 inline Vector2D<Cmpt>
 operator&(const SymmTensor2D<Cmpt>& st, const Vector2D<Cmpt>& v)
@@ -197,7 +209,7 @@ operator&(const SymmTensor2D<Cmpt>& st, const Vector2D<Cmpt>& v)
 }
 
 
-//- Inner-product between a vector and a symmetric tensor
+//- Inner-product of a Vector2D and a SymmTensor2D
 template<class Cmpt>
 inline Vector2D<Cmpt>
 operator&(const Vector2D<Cmpt>& v, const SymmTensor2D<Cmpt>& st)
@@ -210,7 +222,7 @@ operator&(const Vector2D<Cmpt>& v, const SymmTensor2D<Cmpt>& st)
 }
 
 
-//- Inner-sqr of a symmetric tensor
+//- Return the inner-product of a SymmTensor2D with itself
 template<class Cmpt>
 inline SymmTensor2D<Cmpt>
 innerSqr(const SymmTensor2D<Cmpt>& st)
@@ -224,18 +236,19 @@ innerSqr(const SymmTensor2D<Cmpt>& st)
 }
 
 
+//- Return the square of Frobenius norm of a SymmTensor2D as a Cmpt
 template<class Cmpt>
 inline Cmpt magSqr(const SymmTensor2D<Cmpt>& st)
 {
-    return
+    return Cmpt
     (
-        magSqr(st.xx()) + 2*magSqr(st.xy())
-                        +   magSqr(st.yy())
+        mag(st.xx()*st.xx()) + 2*mag(st.xy()*st.xy())
+                             +   mag(st.yy()*st.yy())
     );
 }
 
 
-//- Return the trace of a symmetric tensor
+//- Return the trace of a SymmTensor2D
 template<class Cmpt>
 inline Cmpt tr(const SymmTensor2D<Cmpt>& st)
 {
@@ -243,7 +256,7 @@ inline Cmpt tr(const SymmTensor2D<Cmpt>& st)
 }
 
 
-//- Return the spherical part of a symmetric tensor
+//- Return the spherical part of a SymmTensor2D
 template<class Cmpt>
 inline SphericalTensor2D<Cmpt> sph(const SymmTensor2D<Cmpt>& st)
 {
@@ -254,7 +267,7 @@ inline SphericalTensor2D<Cmpt> sph(const SymmTensor2D<Cmpt>& st)
 }
 
 
-//- Return the symmetric part of a symmetric tensor, i.e. itself
+//- Return the symmetric part of a SymmTensor2D, i.e. itself
 template<class Cmpt>
 inline const SymmTensor2D<Cmpt>& symm(const SymmTensor2D<Cmpt>& st)
 {
@@ -262,7 +275,7 @@ inline const SymmTensor2D<Cmpt>& symm(const SymmTensor2D<Cmpt>& st)
 }
 
 
-//- Return twice the symmetric part of a symmetric tensor
+//- Return twice the symmetric part of a SymmTensor2D, i.e. twice itself
 template<class Cmpt>
 inline SymmTensor2D<Cmpt> twoSymm(const SymmTensor2D<Cmpt>& st)
 {
@@ -270,54 +283,61 @@ inline SymmTensor2D<Cmpt> twoSymm(const SymmTensor2D<Cmpt>& st)
 }
 
 
-//- Return the deviatoric part of a symmetric tensor
+//- Return the deviatoric part of a SymmTensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt> dev(const SymmTensor2D<Cmpt>& st)
 {
-    return st - SphericalTensor2D<Cmpt>::oneThirdI*tr(st);
+    return st - sph(st);
 }
 
 
-//- Return the deviatoric part of a symmetric tensor
+//- Return the two-third deviatoric part of a SymmTensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt> dev2(const SymmTensor2D<Cmpt>& st)
 {
-    return st - SphericalTensor2D<Cmpt>::twoThirdsI*tr(st);
+    return st - 2*sph(st);
 }
 
 
-//- Return the determinant of a symmetric tensor
+//- Return the determinant of a SymmTensor2D
 template<class Cmpt>
 inline Cmpt det(const SymmTensor2D<Cmpt>& st)
 {
-    return
-    (
-        st.xx()*st.yy() - st.xy()*st.xy()
-    );
+    return (st.xx()*st.yy() - st.xy()*st.xy());
 }
 
 
-//- Return the cofactor symmetric tensor of a symmetric tensor
+//- Return the cofactor SymmTensor2D of a SymmTensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt> cof(const SymmTensor2D<Cmpt>& st)
 {
     return SymmTensor2D<Cmpt>
     (
-        st.yy(), -st.xy(),
+        st.yy(), -st.yx(),
                   st.xx()
     );
 }
 
 
-//- Return the inverse of a symmetric tensor give the determinant
+//- Return the inverse of a SymmTensor2D using a given determinant
 template<class Cmpt>
 inline SymmTensor2D<Cmpt> inv(const SymmTensor2D<Cmpt>& st, const Cmpt detst)
 {
+    #ifdef FULLDEBUG
+    if (mag(detst) < SMALL)
+    {
+        FatalErrorInFunction
+            << "SymmTensor2D is not invertible due to the zero determinant:"
+            << "det(SymmTensor2D) = " << mag(detst)
+            << abort(FatalError);
+    }
+    #endif
+    // cofactor and adjugate matrices of SymmTensor are the same, hence no .T()
     return cof(st)/detst;
 }
 
 
-//- Return the inverse of a symmetric tensor
+//- Return the inverse of a SymmTensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt> inv(const SymmTensor2D<Cmpt>& st)
 {
@@ -325,7 +345,7 @@ inline SymmTensor2D<Cmpt> inv(const SymmTensor2D<Cmpt>& st)
 }
 
 
-//- Return the 1st invariant of a symmetric tensor
+//- Return the 1st invariant of a SymmTensor2D
 template<class Cmpt>
 inline Cmpt invariantI(const SymmTensor2D<Cmpt>& st)
 {
@@ -333,30 +353,15 @@ inline Cmpt invariantI(const SymmTensor2D<Cmpt>& st)
 }
 
 
-//- Return the 2nd invariant of a symmetric tensor
+//- Return the 2nd invariant of a SymmTensor2D
 template<class Cmpt>
 inline Cmpt invariantII(const SymmTensor2D<Cmpt>& st)
-{
-    return
-    (
-        0.5*sqr(tr(st))
-      - 0.5*
-        (
-           st.xx()*st.xx() + st.xy()*st.xy()
-         + st.xy()*st.xy() + st.yy()*st.yy()
-        )
-    );
-}
-
-
-//- Return the 3rd invariant of a symmetric tensor
-template<class Cmpt>
-inline Cmpt invariantIII(const SymmTensor2D<Cmpt>& st)
 {
     return det(st);
 }
 
 
+//- Sum of a SphericalTensor2D and a SymmTensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt>
 operator+(const SphericalTensor2D<Cmpt>& spt1, const SymmTensor2D<Cmpt>& st2)
@@ -369,6 +374,7 @@ operator+(const SphericalTensor2D<Cmpt>& spt1, const SymmTensor2D<Cmpt>& st2)
 }
 
 
+//- Sum of a SymmTensor2D and a SphericalTensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt>
 operator+(const SymmTensor2D<Cmpt>& st1, const SphericalTensor2D<Cmpt>& spt2)
@@ -381,6 +387,7 @@ operator+(const SymmTensor2D<Cmpt>& st1, const SphericalTensor2D<Cmpt>& spt2)
 }
 
 
+//- Subtract a SymmTensor2D from a SphericalTensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt>
 operator-(const SphericalTensor2D<Cmpt>& spt1, const SymmTensor2D<Cmpt>& st2)
@@ -393,6 +400,7 @@ operator-(const SphericalTensor2D<Cmpt>& spt1, const SymmTensor2D<Cmpt>& st2)
 }
 
 
+//- Subtract a SphericalTensor2D from a SymmTensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt>
 operator-(const SymmTensor2D<Cmpt>& st1, const SphericalTensor2D<Cmpt>& spt2)
@@ -405,7 +413,7 @@ operator-(const SymmTensor2D<Cmpt>& st1, const SphericalTensor2D<Cmpt>& spt2)
 }
 
 
-//- Inner-product between a spherical symmetric tensor and a symmetric tensor
+//- Inner-product of a SphericalTensor2D and a SymmTensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt>
 operator&(const SphericalTensor2D<Cmpt>& spt1, const SymmTensor2D<Cmpt>& st2)
@@ -418,7 +426,7 @@ operator&(const SphericalTensor2D<Cmpt>& spt1, const SymmTensor2D<Cmpt>& st2)
 }
 
 
-//- Inner-product between a tensor and a spherical tensor
+//- Inner-product of a SymmTensor2D and a SphericalTensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt>
 operator&(const SymmTensor2D<Cmpt>& st1, const SphericalTensor2D<Cmpt>& spt2)
@@ -431,24 +439,25 @@ operator&(const SymmTensor2D<Cmpt>& st1, const SphericalTensor2D<Cmpt>& spt2)
 }
 
 
-//- Double-dot-product between a spherical tensor and a symmetric tensor
+//- Double-inner-product of a SphericalTensor2D and a SymmTensor2D
 template<class Cmpt>
 inline Cmpt
 operator&&(const SphericalTensor2D<Cmpt>& spt1, const SymmTensor2D<Cmpt>& st2)
 {
-    return(spt1.ii()*st2.xx() + spt1.ii()*st2.yy());
+    return (spt1.ii()*st2.xx() + spt1.ii()*st2.yy());
 }
 
 
-//- Double-dot-product between a tensor and a spherical tensor
+//- Double-inner-product of a SymmTensor2D and a SphericalTensor2D
 template<class Cmpt>
 inline Cmpt
 operator&&(const SymmTensor2D<Cmpt>& st1, const SphericalTensor2D<Cmpt>& spt2)
 {
-    return(st1.xx()*spt2.ii() + st1.yy()*spt2.ii());
+    return (st1.xx()*spt2.ii() + st1.yy()*spt2.ii());
 }
 
 
+//- Outer-product of a Vector2D with itself
 template<class Cmpt>
 inline SymmTensor2D<Cmpt> sqr(const Vector2D<Cmpt>& v)
 {
@@ -460,6 +469,21 @@ inline SymmTensor2D<Cmpt> sqr(const Vector2D<Cmpt>& v)
 }
 
 
+//- Division of a SymmTensor2D by a Cmpt
+template<class Cmpt>
+inline SymmTensor2D<Cmpt>
+operator/(const SymmTensor2D<Cmpt>& st, const Cmpt s)
+{
+    return SymmTensor2D<Cmpt>
+    (
+        st.xx()/s, st.xy()/s,
+                   st.yy()/s
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 template<class Cmpt>
 class outerProduct<SymmTensor2D<Cmpt>, Cmpt>
 {

src/OpenFOAM/primitives/Tensor/Tensor.H

@@ -28,10 +28,10 @@ Class
     Foam::Tensor
 
 Description
-    Templated 3D tensor derived from MatrixSpace adding construction from
-    9 components, element access using xx(), xy() etc. member functions and
-    the inner-product (dot-product) and outer-product of two Vectors
-    (tensor-product) operators.
+    A templated (3 x 3) tensor of objects of \<T\> derived from MatrixSpace.
+
+See also
+    Test-Tensor.C
 
 SourceFiles
     TensorI.H
@@ -264,7 +264,7 @@ public:
             );
 
 
-        // Diagonal access.
+        // Diagonal access and manipulation
 
             //- Extract the diagonal as a vector
             inline Vector<Cmpt> diag() const;
@@ -275,7 +275,7 @@ public:
 
     // Tensor Operations
 
-        //- Return transpose
+        //- Return non-Hermitian transpose
         inline Tensor<Cmpt> T() const;
 
         //- Return inverse

src/OpenFOAM/primitives/Tensor/TensorI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,6 +25,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 \*---------------------------------------------------------------------------*/
+#include <type_traits>
 
 #include "SymmTensor.H"
 
@@ -62,9 +63,9 @@ inline Foam::Tensor<Cmpt>::Tensor
 template<class Cmpt>
 inline Foam::Tensor<Cmpt>::Tensor(const SphericalTensor<Cmpt>& st)
 {
-    this->v_[XX] = st.ii(); this->v_[XY] = 0;       this->v_[XZ] = 0;
-    this->v_[YX] = 0;       this->v_[YY] = st.ii(); this->v_[YZ] = 0;
-    this->v_[ZX] = 0;       this->v_[ZY] = 0;       this->v_[ZZ] = st.ii();
+    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();
 }
 
 
@@ -572,9 +573,9 @@ inline void Foam::Tensor<Cmpt>::operator=
 template<class Cmpt>
 inline void Foam::Tensor<Cmpt>::operator=(const SphericalTensor<Cmpt>& st)
 {
-    this->v_[XX] = st.ii(); this->v_[XY] = 0; this->v_[XZ] = 0;
-    this->v_[YX] = 0; this->v_[YY] = st.ii(); this->v_[YZ] = 0;
-    this->v_[ZX] = 0; this->v_[ZY] = 0; this->v_[ZZ] = st.ii();
+    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();
 }
 
 
@@ -611,6 +612,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
+//- Return the Hodge dual of a Tensor as a Vector
 template<class Cmpt>
 inline Vector<Cmpt> operator*(const Tensor<Cmpt>& t)
 {
@@ -618,18 +620,20 @@ inline Vector<Cmpt> operator*(const Tensor<Cmpt>& t)
 }
 
 
+//- Return the Hodge dual of a Vector as a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt> operator*(const Vector<Cmpt>& v)
 {
     return Tensor<Cmpt>
     (
-             0, -v.z(),   v.y(),
-         v.z(),      0,  -v.x(),
-        -v.y(),  v.x(),       0
+         Zero,  -v.z(),  v.y(),
+         v.z(),  Zero,  -v.x(),
+        -v.y(),  v.x(),  Zero
     );
 }
 
 
+//- Inner-product of a Tensor and a Tensor
 template<class Cmpt>
 inline typename innerProduct<Tensor<Cmpt>, Tensor<Cmpt>>::type
 operator&(const Tensor<Cmpt>& t1, const Tensor<Cmpt>& t2)
@@ -638,6 +642,7 @@ operator&(const Tensor<Cmpt>& t1, const Tensor<Cmpt>& t2)
 }
 
 
+//- 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)
@@ -651,6 +656,7 @@ operator&(const Tensor<Cmpt>& t, const Vector<Cmpt>& v)
 }
 
 
+//- 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)
@@ -664,6 +670,7 @@ operator&(const Vector<Cmpt>& v, const Tensor<Cmpt>& t)
 }
 
 
+//- Outer-product of a Vector and a Vector
 template<class Cmpt>
 inline typename outerProduct<Vector<Cmpt>, Vector<Cmpt>>::type
 operator*(const Vector<Cmpt>& v1, const Vector<Cmpt>& v2)
@@ -677,6 +684,7 @@ operator*(const Vector<Cmpt>& v1, const Vector<Cmpt>& v2)
 }
 
 
+//- Division of a Vector by a Tensor
 template<class Cmpt>
 inline typename innerProduct<Vector<Cmpt>, Tensor<Cmpt>>::type
 operator/(const Vector<Cmpt>& v, const Tensor<Cmpt>& t)
@@ -685,9 +693,34 @@ operator/(const Vector<Cmpt>& v, const Tensor<Cmpt>& t)
 }
 
 
+//- Division of a Tensor by a Cmpt
+template<class Cmpt>
+inline Tensor<Cmpt>
+operator/(const Tensor<Cmpt>& t, const Cmpt s)
+{
+    #ifdef FULLDEBUG
+    if (mag(s) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "Tensor = " << t
+            << " is not divisible due to a zero value in Cmpt:"
+            << "Cmpt = " << s
+            << abort(FatalError);
+    }
+    #endif
+
+    return Tensor<Cmpt>
+    (
+        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
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 
-//- Return the trace of a tensor
+//- Return the trace of a Tensor
 template<class Cmpt>
 inline Cmpt tr(const Tensor<Cmpt>& t)
 {
@@ -695,7 +728,7 @@ inline Cmpt tr(const Tensor<Cmpt>& t)
 }
 
 
-//- Return the spherical part of a tensor
+//- Return the spherical part of a Tensor
 template<class Cmpt>
 inline SphericalTensor<Cmpt> sph(const Tensor<Cmpt>& t)
 {
@@ -706,7 +739,7 @@ inline SphericalTensor<Cmpt> sph(const Tensor<Cmpt>& t)
 }
 
 
-//- Return the symmetric part of a tensor
+//- Return the symmetric part of a Tensor
 template<class Cmpt>
 inline SymmTensor<Cmpt> symm(const Tensor<Cmpt>& t)
 {
@@ -719,7 +752,7 @@ inline SymmTensor<Cmpt> symm(const Tensor<Cmpt>& t)
 }
 
 
-//- Return twice the symmetric part of a tensor
+//- Return twice the symmetric part of a Tensor
 template<class Cmpt>
 inline SymmTensor<Cmpt> twoSymm(const Tensor<Cmpt>& t)
 {
@@ -732,20 +765,20 @@ inline SymmTensor<Cmpt> twoSymm(const Tensor<Cmpt>& t)
 }
 
 
-//- Return the skew-symmetric part of a tensor
+//- Return the skew-symmetric part of a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt> skew(const Tensor<Cmpt>& t)
 {
     return Tensor<Cmpt>
     (
-        0.0, 0.5*(t.xy() - t.yx()), 0.5*(t.xz() - t.zx()),
-        0.5*(t.yx() - t.xy()), 0.0, 0.5*(t.yz() - t.zy()),
-        0.5*(t.zx() - t.xz()), 0.5*(t.zy() - t.yz()), 0.0
+        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
     );
 }
 
 
-//- Return the skew-symmetric part of a symmetric tensor
+//- Return the skew-symmetric part of a SymmTensor as a Tensor
 template<class Cmpt>
 inline const Tensor<Cmpt>& skew(const SymmTensor<Cmpt>& st)
 {
@@ -753,23 +786,23 @@ inline const Tensor<Cmpt>& skew(const SymmTensor<Cmpt>& st)
 }
 
 
-//- Return the deviatoric part of a tensor
+//- Return the deviatoric part of a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt> dev(const Tensor<Cmpt>& t)
 {
-    return t - SphericalTensor<Cmpt>::oneThirdI*tr(t);
+    return t - sph(t);
 }
 
 
-//- Return the deviatoric part of a tensor
+//- Return the two-third deviatoric part of a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt> dev2(const Tensor<Cmpt>& t)
 {
-    return t - SphericalTensor<Cmpt>::twoThirdsI*tr(t);
+    return t - 2*sph(t);
 }
 
 
-//- Return the determinant of a tensor
+//- Return the determinant of a Tensor
 template<class Cmpt>
 inline Cmpt det(const Tensor<Cmpt>& t)
 {
@@ -782,7 +815,7 @@ inline Cmpt det(const Tensor<Cmpt>& t)
 }
 
 
-//- Return the cofactor tensor of a tensor
+//- Return the cofactor Tensor of a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt> cof(const Tensor<Cmpt>& t)
 {
@@ -803,28 +836,25 @@ inline Tensor<Cmpt> cof(const Tensor<Cmpt>& t)
 }
 
 
-//- Return the inverse of a tensor given the determinant
+//- Return the inverse of a Tensor by using the given determinant
 template<class Cmpt>
 inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t, const Cmpt dett)
 {
-    return Tensor<Cmpt>
-    (
-        t.yy()*t.zz() - t.zy()*t.yz(),
-        t.xz()*t.zy() - t.xy()*t.zz(),
-        t.xy()*t.yz() - t.xz()*t.yy(),
+    #ifdef FULLDEBUG
+    if (mag(dett) < SMALL)
+    {
+        FatalErrorInFunction
+            << "Tensor is not invertible due to the zero determinant:"
+            << "det(Tensor) = " << mag(dett)
+            << abort(FatalError);
+    }
+    #endif
 
-        t.zx()*t.yz() - t.yx()*t.zz(),
-        t.xx()*t.zz() - t.xz()*t.zx(),
-        t.yx()*t.xz() - t.xx()*t.yz(),
-
-        t.yx()*t.zy() - t.yy()*t.zx(),
-        t.xy()*t.zx() - t.xx()*t.zy(),
-        t.xx()*t.yy() - t.yx()*t.xy()
-    )/dett;
+    return cof(t).T()/dett;
 }
 
 
-//- Return the inverse of a tensor
+//- Return the inverse of a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t)
 {
@@ -832,6 +862,7 @@ inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t)
 }
 
 
+//- Return the inverse of this Tensor
 template<class Cmpt>
 inline Tensor<Cmpt> Tensor<Cmpt>::inv() const
 {
@@ -839,7 +870,7 @@ inline Tensor<Cmpt> Tensor<Cmpt>::inv() const
 }
 
 
-//- Return the 1st invariant of a tensor
+//- Return the 1st invariant of a Tensor
 template<class Cmpt>
 inline Cmpt invariantI(const Tensor<Cmpt>& t)
 {
@@ -847,7 +878,7 @@ inline Cmpt invariantI(const Tensor<Cmpt>& t)
 }
 
 
-//- Return the 2nd invariant of a tensor
+//- Return the 2nd invariant of a Tensor
 template<class Cmpt>
 inline Cmpt invariantII(const Tensor<Cmpt>& t)
 {
@@ -859,7 +890,7 @@ inline Cmpt invariantII(const Tensor<Cmpt>& t)
 }
 
 
-//- Return the 3rd invariant of a tensor
+//- Return the 3rd invariant of a Tensor
 template<class Cmpt>
 inline Cmpt invariantIII(const Tensor<Cmpt>& t)
 {
@@ -869,6 +900,7 @@ inline Cmpt invariantIII(const Tensor<Cmpt>& t)
 
 // * * * * * * * * * Mixed Tensor SphericalTensor Operators  * * * * * * * * //
 
+//- Sum of a SphericalTensor and a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator+(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
@@ -882,6 +914,7 @@ operator+(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 }
 
 
+//- Sum of a Tensor and a SphericalTensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator+(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
@@ -895,6 +928,7 @@ operator+(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
 }
 
 
+//- Subtract a Tensor from a SphericalTensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator-(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
@@ -908,6 +942,7 @@ operator-(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 }
 
 
+//- Subtract a SphericalTensor from a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator-(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
@@ -921,7 +956,7 @@ operator-(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
 }
 
 
-//- Inner-product between a spherical tensor and a tensor
+//- Inner-product of a SphericalTensor and a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator&(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
@@ -935,7 +970,7 @@ operator&(const SphericalTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 }
 
 
-//- Inner-product between a tensor and a spherical tensor
+//- Inner-product of a Tensor and a SphericalTensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator&(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
@@ -949,21 +984,21 @@ operator&(const Tensor<Cmpt>& t1, const SphericalTensor<Cmpt>& st2)
 }
 
 
-//- Double-dot-product between a spherical tensor and a tensor
+//- 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());
+    return (st1.ii()*t2.xx() + st1.ii()*t2.yy() + st1.ii()*t2.zz());
 }
 
 
-//- Double-dot-product between a tensor and a spherical tensor
+//- 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());
+    return (t1.xx()*st2.ii() + t1.yy()*st2.ii() + t1.zz()*st2.ii());
 }
 
 
@@ -1005,6 +1040,7 @@ public:
 
 // * * * * * * * * * * Mixed Tensor SymmTensor Operators * * * * * * * * * * //
 
+//- Sum of a SymmTensor and a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator+(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
@@ -1018,6 +1054,7 @@ operator+(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 }
 
 
+//- Sum of a Tensor and a SymmTensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator+(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
@@ -1031,6 +1068,7 @@ operator+(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
 }
 
 
+//- Subtract a Tensor from a SymmTensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator-(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
@@ -1044,6 +1082,7 @@ operator-(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 }
 
 
+//- Subtract a SymmTensor from a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator-(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
@@ -1057,7 +1096,7 @@ operator-(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
 }
 
 
-//- Inner-product between a symmetric tensor and a tensor
+//- Inner-product of a SymmTensor and a Tensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator&(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
@@ -1079,7 +1118,7 @@ operator&(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 }
 
 
-//- Inner-product between a tensor and a symmetric tensor
+//- Inner-product of a Tensor and a SymmTensor
 template<class Cmpt>
 inline Tensor<Cmpt>
 operator&(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
@@ -1101,7 +1140,7 @@ operator&(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)
 }
 
 
-//- Double-dot-product between a symmetric tensor and a tensor
+//- Double-inner-product of a SymmTensor and a Tensor
 template<class Cmpt>
 inline Cmpt
 operator&&(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
@@ -1115,7 +1154,7 @@ operator&&(const SymmTensor<Cmpt>& st1, const Tensor<Cmpt>& t2)
 }
 
 
-//- Double-dot-product between a tensor and a symmetric tensor
+//- Double-inner-product of a Tensor and a SymmTensor
 template<class Cmpt>
 inline Cmpt
 operator&&(const Tensor<Cmpt>& t1, const SymmTensor<Cmpt>& st2)

src/OpenFOAM/primitives/Tensor2D/Tensor2D.H

@@ -28,10 +28,10 @@ Class
     Foam::Tensor2D
 
 Description
-    Templated 2D tensor derived from VectorSpace adding construction from
-    4 components, element access using xx(), xy(), yx() and yy() member
-    functions and the iner-product (dot-product) and outer-product of two
-    Vector2Ds (tensor-product) operators.
+    A templated (2 x 2) tensor of objects of \<T\> derived from VectorSpace.
+
+See also
+    Test-Tensor2D.C
 
 SourceFiles
     Tensor2DI.H
@@ -203,7 +203,7 @@ public:
             inline void rows(const Vector2D<Cmpt>& x, const Vector2D<Cmpt>& y);
 
 
-        // Diagonal access.
+        // Diagonal access and manipulation
 
             //- Extract the diagonal as a vector
             inline Vector2D<Cmpt> diag() const;
@@ -214,7 +214,7 @@ public:
 
     // Tensor Operations
 
-        //- Return transpose
+        //- Return non-Hermitian transpose
         inline Tensor2D<Cmpt> T() const;
 
         //- Inner-product of this with another Tensor2D.

src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,8 +56,8 @@ inline Foam::Tensor2D<Cmpt>::Tensor2D(const SymmTensor2D<Cmpt>& st)
 template<class Cmpt>
 inline Foam::Tensor2D<Cmpt>::Tensor2D(const SphericalTensor2D<Cmpt>& st)
 {
-    this->v_[XX] = st.ii(); this->v_[XY] = 0;
-    this->v_[YX] = 0;       this->v_[YY] = st.ii();
+    this->v_[XX] = st.ii(); this->v_[XY] = Zero;
+    this->v_[YX] = Zero;    this->v_[YY] = st.ii();
 }
 
 
@@ -395,8 +395,8 @@ inline void Foam::Tensor2D<Cmpt>::operator=(const SymmTensor2D<Cmpt>& st)
 template<class Cmpt>
 inline void Foam::Tensor2D<Cmpt>::operator=(const SphericalTensor2D<Cmpt>& st)
 {
-    this->v_[XX] = st.ii(); this->v_[XY] = 0;
-    this->v_[YX] = 0; this->v_[YY] = st.ii();
+    this->v_[XX] = st.ii(); this->v_[XY] = Zero;
+    this->v_[YX] = Zero;    this->v_[YY] = st.ii();
 }
 
 
@@ -407,7 +407,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
-//- Inner-product between two tensors
+//- Inner-product of a Tensor2D and a Tensor2D
 template<class Cmpt>
 inline typename innerProduct<Tensor2D<Cmpt>, Tensor2D<Cmpt>>::type
 operator&(const Tensor2D<Cmpt>& t1, const Tensor2D<Cmpt>& t2)
@@ -416,7 +416,7 @@ operator&(const Tensor2D<Cmpt>& t1, const Tensor2D<Cmpt>& t2)
 }
 
 
-//- Inner-product between a tensor and a vector
+//- Inner-product of a Tensor2D and a Vector2D
 template<class Cmpt>
 inline typename innerProduct<Tensor2D<Cmpt>, Vector2D<Cmpt>>::type
 operator&(const Tensor2D<Cmpt>& t, const Vector2D<Cmpt>& v)
@@ -428,7 +428,7 @@ operator&(const Tensor2D<Cmpt>& t, const Vector2D<Cmpt>& v)
     );
 }
 
-//- Inner-product between a vector and a tensor
+//- Inner-product of a Vector2D and a Tensor2D
 template<class Cmpt>
 inline typename innerProduct<Vector2D<Cmpt>, Tensor2D<Cmpt>>::type
 operator&(const Vector2D<Cmpt>& v, const Tensor2D<Cmpt>& t)
@@ -440,7 +440,7 @@ operator&(const Vector2D<Cmpt>& v, const Tensor2D<Cmpt>& t)
     );
 }
 
-//- Outer-product between two vectors
+//- Outer-product of a Vector2D and a Vector2D
 template<class Cmpt>
 inline typename outerProduct<Vector2D<Cmpt>, Vector2D<Cmpt>>::type
 operator*(const Vector2D<Cmpt>& v1, const Vector2D<Cmpt>& v2)
@@ -453,7 +453,7 @@ operator*(const Vector2D<Cmpt>& v1, const Vector2D<Cmpt>& v2)
 }
 
 
-//- Return the trace of a tensor
+//- Return the trace of a Tensor2D
 template<class Cmpt>
 inline Cmpt tr(const Tensor2D<Cmpt>& t)
 {
@@ -461,7 +461,7 @@ inline Cmpt tr(const Tensor2D<Cmpt>& t)
 }
 
 
-//- Return the spherical part of a tensor
+//- Return the spherical part of a Tensor2D
 template<class Cmpt>
 inline SphericalTensor2D<Cmpt> sph(const Tensor2D<Cmpt>& t)
 {
@@ -472,7 +472,7 @@ inline SphericalTensor2D<Cmpt> sph(const Tensor2D<Cmpt>& t)
 }
 
 
-//- Return the symmetric part of a tensor
+//- Return the symmetric part of a Tensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt> symm(const Tensor2D<Cmpt>& t)
 {
@@ -483,7 +483,8 @@ inline SymmTensor2D<Cmpt> symm(const Tensor2D<Cmpt>& t)
     );
 }
 
-//- Return the twice the symmetric part of a tensor
+
+//- Return the twice the symmetric part of a Tensor2D
 template<class Cmpt>
 inline SymmTensor2D<Cmpt> twoSymm(const Tensor2D<Cmpt>& t)
 {
@@ -494,63 +495,74 @@ inline SymmTensor2D<Cmpt> twoSymm(const Tensor2D<Cmpt>& t)
     );
 }
 
-//- Return the skew-symmetric part of a tensor
+
+//- Return the skew-symmetric part of a Tensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt> skew(const Tensor2D<Cmpt>& t)
 {
     return Tensor2D<Cmpt>
     (
-        0.0, 0.5*(t.xy() - t.yx()),
-        0.5*(t.yx() - t.xy()), 0.0
+        Zero,                  0.5*(t.xy() - t.yx()),
+        0.5*(t.yx() - t.xy()), Zero
     );
 }
 
 
-//- Return the deviatoric part of a tensor
+//- Return the deviatoric part of a Tensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt> dev(const Tensor2D<Cmpt>& t)
 {
-    return t - SphericalTensor2D<Cmpt>::oneThirdI*tr(t);
+    return t - sph(t);
 }
 
 
-//- Return the deviatoric part of a tensor
+//- Return the two-third deviatoric part of a Tensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt> dev2(const Tensor2D<Cmpt>& t)
 {
-    return t - SphericalTensor2D<Cmpt>::twoThirdsI*tr(t);
+    return t - 2*sph(t);
 }
 
 
-//- Return the determinant of a tensor
+//- Return the determinant of a Tensor2D
 template<class Cmpt>
 inline Cmpt det(const Tensor2D<Cmpt>& t)
 {
-    return(t.xx()*t.yy() - t.xy()*t.yx());
+    return t.xx()*t.yy() - t.xy()*t.yx();
 }
 
 
-//- Return the cofactor tensor of a tensor
+//- Return the cofactor Tensor2D of a Tensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt> cof(const Tensor2D<Cmpt>& t)
 {
     return Tensor2D<Cmpt>
     (
-        t.yy(), -t.xy(),
-       -t.yx(),  t.xx()
+        t.yy(), -t.yx(),
+       -t.xy(),  t.xx()
     );
 }
 
 
-//- Return the inverse of a tensor given the determinant
+//- Return the inverse of a Tensor2D using a given determinant
 template<class Cmpt>
 inline Tensor2D<Cmpt> inv(const Tensor2D<Cmpt>& t, const Cmpt dett)
 {
-    return cof(t)/dett;
+    #ifdef FULLDEBUG
+    if (mag(dett) < SMALL)
+    {
+        FatalErrorInFunction
+            << "Tensor2D is not invertible due to the zero determinant:"
+            << "det(Tensor2D) = " << mag(dett)
+            << abort(FatalError);
+    }
+    #endif
+
+    return cof(t).T()/dett;
 }
 
 
-//- Return the inverse of a tensor
+//- Return the inverse of a Tensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt> inv(const Tensor2D<Cmpt>& t)
 {
@@ -558,7 +570,7 @@ inline Tensor2D<Cmpt> inv(const Tensor2D<Cmpt>& t)
 }
 
 
-//- Return the 1st invariant of a tensor
+//- Return the 1st invariant of a Tensor2D
 template<class Cmpt>
 inline Cmpt invariantI(const Tensor2D<Cmpt>& t)
 {
@@ -566,25 +578,9 @@ inline Cmpt invariantI(const Tensor2D<Cmpt>& t)
 }
 
 
-//- Return the 2nd invariant of a tensor
+//- Return the 2nd invariant of a Tensor2D
 template<class Cmpt>
 inline Cmpt invariantII(const Tensor2D<Cmpt>& t)
-{
-    return
-    (
-        0.5*sqr(tr(t))
-      - 0.5*
-        (
-           t.xx()*t.xx() + t.xy()*t.xy()
-         + t.yx()*t.yx() + t.yy()*t.yy()
-        )
-    );
-}
-
-
-//- Return the 3rd invariant of a tensor
-template<class Cmpt>
-inline Cmpt invariantIII(const Tensor2D<Cmpt>& t)
 {
     return det(t);
 }
@@ -592,6 +588,7 @@ inline Cmpt invariantIII(const Tensor2D<Cmpt>& t)
 
 // * * * * * * * * * Mixed Tensor SphericalTensor Operators  * * * * * * * * //
 
+//- Sum of a SphericalTensor2D and a Tensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator+(const SphericalTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
@@ -604,6 +601,7 @@ operator+(const SphericalTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
 }
 
 
+//- Sum of a Tensor2D and a SphericalTensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator+(const Tensor2D<Cmpt>& t1, const SphericalTensor2D<Cmpt>& st2)
@@ -616,6 +614,7 @@ operator+(const Tensor2D<Cmpt>& t1, const SphericalTensor2D<Cmpt>& st2)
 }
 
 
+//- Subtract a Tensor2D from a SphericalTensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator-(const SphericalTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
@@ -628,6 +627,7 @@ operator-(const SphericalTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
 }
 
 
+//- Subtract a SphericalTensor2D from a Tensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator-(const Tensor2D<Cmpt>& t1, const SphericalTensor2D<Cmpt>& st2)
@@ -640,7 +640,7 @@ operator-(const Tensor2D<Cmpt>& t1, const SphericalTensor2D<Cmpt>& st2)
 }
 
 
-//- Inner-product between a spherical tensor and a tensor
+//- Inner-product of a SphericalTensor2D and Tensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator&(const SphericalTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
@@ -655,7 +655,7 @@ operator&(const SphericalTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
 }
 
 
-//- Inner-product between a tensor and a spherical tensor
+//- Inner-product of a Tensor2D and a SphericalTensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator&(const Tensor2D<Cmpt>& t1, const SphericalTensor2D<Cmpt>& st2)
@@ -671,26 +671,27 @@ operator&(const Tensor2D<Cmpt>& t1, const SphericalTensor2D<Cmpt>& st2)
 }
 
 
-//- Double-dot-product between a spherical tensor and a tensor
+//- Double-inner-product of a SphericalTensor2D and a Tensor2D
 template<class Cmpt>
 inline Cmpt
 operator&&(const SphericalTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
 {
-    return(st1.ii()*t2.xx() + st1.ii()*t2.yy());
+    return (st1.ii()*t2.xx() + st1.ii()*t2.yy());
 }
 
 
-//- Double-dot-product between a tensor and a spherical tensor
+//- Double-inner-product of a Tensor2D and a SphericalTensor2D
 template<class Cmpt>
 inline Cmpt
 operator&&(const Tensor2D<Cmpt>& t1, const SphericalTensor2D<Cmpt>& st2)
 {
-    return(t1.xx()*st2.ii() + t1.yy()*st2.ii());
+    return (t1.xx()*st2.ii() + t1.yy()*st2.ii());
 }
 
 
 // * * * * * * * * * * Mixed Tensor SymmTensor Operators * * * * * * * * * * //
 
+//- Sum of a SymmTensor2D and a Tensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator+(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
@@ -703,6 +704,7 @@ operator+(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
 }
 
 
+//- Sum of a Tensor2D and a SymmTensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator+(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
@@ -715,6 +717,7 @@ operator+(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
 }
 
 
+//- Subtract a Tensor2D from a SymmTensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator-(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
@@ -727,6 +730,7 @@ operator-(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
 }
 
 
+//- Subtract a SymmTensor2D from a Tensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator-(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
@@ -739,7 +743,7 @@ operator-(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
 }
 
 
-//- Inner-product between a spherical tensor and a tensor
+//- Inner-product of a SymmTensor2D and Tensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator&(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
@@ -755,7 +759,7 @@ operator&(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
 }
 
 
-//- Inner-product between a tensor and a spherical tensor
+//- Inner-product of a Tensor2D and a SymmTensor2D
 template<class Cmpt>
 inline Tensor2D<Cmpt>
 operator&(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
@@ -771,7 +775,7 @@ operator&(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
 }
 
 
-//- Double-dot-product between a spherical tensor and a tensor
+//- Double-inner-product of a SymmTensor2D and a Tensor2D
 template<class Cmpt>
 inline Cmpt
 operator&&(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
@@ -784,7 +788,7 @@ operator&&(const SymmTensor2D<Cmpt>& st1, const Tensor2D<Cmpt>& t2)
 }
 
 
-//- Double-dot-product between a tensor and a spherical tensor
+//- Double-inner-product of a Tensor2D and a SymmTensor2D
 template<class Cmpt>
 inline Cmpt
 operator&&(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
@@ -797,6 +801,32 @@ operator&&(const Tensor2D<Cmpt>& t1, const SymmTensor2D<Cmpt>& st2)
 }
 
 
+//- Division of a Tensor2D by a Cmpt
+template<class Cmpt>
+inline Tensor2D<Cmpt>
+operator/(const Tensor2D<Cmpt>& t, const Cmpt s)
+{
+    #ifdef FULLDEBUG
+    if (mag(s) < VSMALL)
+    {
+        FatalErrorInFunction
+            << "Tensor2D = " << t
+            << " is not divisible due to a zero value in Cmpt:"
+            << "Cmpt = " << s
+            << abort(FatalError);
+    }
+    #endif
+
+    return Tensor2D<Cmpt>
+    (
+        t.xx()/s, t.xy()/s,
+        t.yx()/s, t.yy()/s
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 template<class Cmpt>
 class typeOfSum<SphericalTensor2D<Cmpt>, Tensor2D<Cmpt>>
 {

src/OpenFOAM/primitives/Vector/VectorI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -176,6 +176,20 @@ inline Vector<Cmpt> operator^(const Vector<Cmpt>& v1, const Vector<Cmpt>& v2)
 }
 
 
+template<class Cmpt>
+inline Vector<Cmpt> operator*(const Cmpt& s, const Vector<Cmpt>& v)
+{
+    return Vector<Cmpt>(s*v.x(), s*v.y(), s*v.z());
+}
+
+
+template<class Cmpt>
+inline Vector<Cmpt> operator*(const Vector<Cmpt>& v, const Cmpt& s)
+{
+    return s*v;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/primitives/Vector2D/Vector2DI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -130,6 +130,20 @@ inline scalar Vector2D<Cmpt>::perp(const Vector2D<Cmpt>& b) const
 }
 
 
+template<class Cmpt>
+inline Vector2D<Cmpt> operator*(const Cmpt& s, const Vector2D<Cmpt>& v)
+{
+    return Vector2D<Cmpt>(s*v.x(), s*v.y());
+}
+
+
+template<class Cmpt>
+inline Vector2D<Cmpt> operator*(const Vector2D<Cmpt>& v, const Cmpt& s)
+{
+    return s*v;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam