diff --git a/applications/test/invTensor/Make/files b/applications/test/invTensor/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..c99507752a8dfbb9a22d31479d899f68d3d40dd7
--- /dev/null
+++ b/applications/test/invTensor/Make/files
@@ -0,0 +1,3 @@
+Test-invTensor.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-invTensor
diff --git a/applications/test/invTensor/Make/options b/applications/test/invTensor/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..18e6fe47afacb902cddccf82632772447704fd88
--- /dev/null
+++ b/applications/test/invTensor/Make/options
@@ -0,0 +1,2 @@
+/* EXE_INC = */
+/* EXE_LIBS = */
diff --git a/applications/test/invTensor/Test-invTensor.C b/applications/test/invTensor/Test-invTensor.C
new file mode 100644
index 0000000000000000000000000000000000000000..b16b7d36662a76393f36367183f91cba96f11291
--- /dev/null
+++ b/applications/test/invTensor/Test-invTensor.C
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2023 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-invTensor
+
+Description
+ Tests for regular and corner cases of tensor inversion.
+
+\*---------------------------------------------------------------------------*/
+
+#include "tensor.H"
+#include "symmTensor.H"
+#include "transform.H"
+#include "unitConversion.H"
+#include "Random.H"
+#include "scalar.H"
+#include "complex.H"
+#include "sigFpe.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class TensorType>
+void test_invert(const TensorType& tt)
+{
+ Info<< pTraits<TensorType>::typeName << nl
+ << "ten : " << tt << nl;
+
+ // Non-failsafe. try/catch does not work here
+ if (tt.det() < ROOTVSMALL)
+ {
+ Info<< "inv : " << "nan/inf" << nl;
+ }
+ else
+ {
+ Info<< "inv : " << tt.inv() << nl;
+ }
+
+ // Failsafe
+ Info<< "inv : " << tt.safeInv() << nl;
+
+ Info<< nl;
+}
+
+
+// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+ // Run with FOAM_SIGFPE=true to ensure we see divide by zero
+ Foam::sigFpe::set(true);
+
+ {
+ symmTensor st(1e-3, 0, 0, 1e-3, 0, 1e-6);
+ test_invert(st);
+ }
+
+ {
+ symmTensor st(1e-3, 0, 0, 1e-3, 0, 1e-30);
+ test_invert(st);
+ }
+
+ {
+ symmTensor st(1e-3, 0, 0, 1e-3, 0, 0);
+ test_invert(st);
+ }
+
+ return 0;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/dimensionedTypes/dimensionedSymmTensor/dimensionedSymmTensor.H b/src/OpenFOAM/dimensionedTypes/dimensionedSymmTensor/dimensionedSymmTensor.H
index 2d7cbf60163a9f4231feb24adf117a92f9f37bab..a43534155124d42092a3c9e1be91a5f95618eb52 100644
--- a/src/OpenFOAM/dimensionedTypes/dimensionedSymmTensor/dimensionedSymmTensor.H
+++ b/src/OpenFOAM/dimensionedTypes/dimensionedSymmTensor/dimensionedSymmTensor.H
@@ -50,7 +50,7 @@ namespace Foam
typedef dimensioned<symmTensor> dimensionedSymmTensor;
-// global functions
+// Global Functions
dimensionedSymmTensor sqr(const dimensionedVector&);
dimensionedSymmTensor innerSqr(const dimensionedSymmTensor&);
@@ -65,7 +65,7 @@ dimensionedSymmTensor cof(const dimensionedSymmTensor&);
dimensionedSymmTensor inv(const dimensionedSymmTensor&);
-// global operators
+// Global Operators
//- Hodge Dual operator (tensor -> vector)
dimensionedVector operator*(const dimensionedSymmTensor&);
diff --git a/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.H b/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.H
index 8de0683efe3a05a224c92a483113818910865e83..a66aea3c07f71026e9068eee4a90f7fc737a24bf 100644
--- a/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.H
+++ b/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.H
@@ -51,8 +51,7 @@ namespace Foam
typedef dimensioned<tensor> dimensionedTensor;
-
-// global functions
+// Global Functions
dimensionedScalar tr(const dimensionedTensor&);
dimensionedTensor dev(const dimensionedTensor&);
@@ -68,7 +67,7 @@ dimensionedVector eigenValues(const dimensionedSymmTensor&);
dimensionedTensor eigenVectors(const dimensionedSymmTensor&);
-// global operators
+// Global Operators
//- Hodge Dual operator (tensor -> vector)
dimensionedVector operator*(const dimensionedTensor&);
diff --git a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
index 1be5183d4c650e0a929e00a8f445c02ce5abbfcd..afdf86c186e6f2d6641cef97b2f55b855ded80e6 100644
--- a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
+++ b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
@@ -51,55 +51,10 @@ UNARY_FUNCTION(symmTensor, symmTensor, dev2)
UNARY_FUNCTION(scalar, symmTensor, det)
UNARY_FUNCTION(symmTensor, symmTensor, cof)
-void inv(Field<symmTensor>& result, const UList<symmTensor>& tf1)
+void inv(Field<symmTensor>& result, const UList<symmTensor>& f1)
{
- if (result.empty() || tf1.empty())
- {
- return;
- }
-
- // Attempting to identify 2-D cases
- const scalar minThreshold = SMALL * magSqr(tf1[0]);
-
- const bool small_xx = (magSqr(tf1[0].xx()) < minThreshold);
- const bool small_yy = (magSqr(tf1[0].yy()) < minThreshold);
- const bool small_zz = (magSqr(tf1[0].zz()) < minThreshold);
-
- if (small_xx || small_yy || small_zz)
- {
- const vector adjust
- (
- (small_xx ? 1 : 0),
- (small_yy ? 1 : 0),
- (small_zz ? 1 : 0)
- );
-
- // Cannot use TFOR_ALL_F_OP_FUNC_F (additional operations)
-
- const label loopLen = (result).size();
-
- /* pragmas... */
- for (label i = 0; i < loopLen; ++i)
- {
- symmTensor work(tf1[i]);
- work.addDiag(adjust);
-
- result[i] = Foam::inv(work);
- result[i].subtractDiag(adjust);
- }
- }
- else
- {
- // Same as TFOR_ALL_F_OP_FUNC_F
-
- const label loopLen = (result).size();
-
- /* pragmas... */
- for (label i = 0; i < loopLen; ++i)
- {
- result[i] = Foam::inv(tf1[i]);
- }
- }
+ // With 'failsafe' invert
+ TFOR_ALL_F_OP_F_FUNC(symmTensor, result, =, symmTensor, f1, safeInv)
}
tmp<symmTensorField> inv(const UList<symmTensor>& tf)
diff --git a/src/OpenFOAM/fields/Fields/tensorField/tensorField.C b/src/OpenFOAM/fields/Fields/tensorField/tensorField.C
index c4f76fc9788b76a95aa3918a325c7167cf8876c4..9b343fa5a0389597f15439b1f715ef369f58d283 100644
--- a/src/OpenFOAM/fields/Fields/tensorField/tensorField.C
+++ b/src/OpenFOAM/fields/Fields/tensorField/tensorField.C
@@ -50,55 +50,10 @@ UNARY_FUNCTION(tensor, tensor, dev2)
UNARY_FUNCTION(scalar, tensor, det)
UNARY_FUNCTION(tensor, tensor, cof)
-void inv(Field<tensor>& result, const UList<tensor>& tf1)
+void inv(Field<tensor>& result, const UList<tensor>& f1)
{
- if (result.empty() || tf1.empty())
- {
- return;
- }
-
- // Attempting to identify 2-D cases
- const scalar minThreshold = SMALL * magSqr(tf1[0]);
-
- const bool small_xx = (magSqr(tf1[0].xx()) < minThreshold);
- const bool small_yy = (magSqr(tf1[0].yy()) < minThreshold);
- const bool small_zz = (magSqr(tf1[0].zz()) < minThreshold);
-
- if (small_xx || small_yy || small_zz)
- {
- const vector adjust
- (
- (small_xx ? 1 : 0),
- (small_yy ? 1 : 0),
- (small_zz ? 1 : 0)
- );
-
- // Cannot use TFOR_ALL_F_OP_FUNC_F (additional operations)
-
- const label loopLen = (result).size();
-
- /* pragmas... */
- for (label i = 0; i < loopLen; ++i)
- {
- tensor work(tf1[i]);
- work.addDiag(adjust);
-
- result[i] = Foam::inv(work);
- result[i].subtractDiag(adjust);
- }
- }
- else
- {
- // Same as TFOR_ALL_F_OP_FUNC_F
-
- const label loopLen = (result).size();
-
- /* pragmas... */
- for (label i = 0; i < loopLen; ++i)
- {
- result[i] = Foam::inv(tf1[i]);
- }
- }
+ // With 'failsafe' invert
+ TFOR_ALL_F_OP_F_FUNC(tensor, result, =, tensor, f1, safeInv)
}
tmp<tensorField> inv(const UList<tensor>& tf)
diff --git a/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H b/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H
index e3650810231eea727c5477625309d3d1b1e5b50e..2f38d3e481aed30348dc7c40ece150976ed451fb 100644
--- a/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H
+++ b/src/OpenFOAM/primitives/SymmTensor/SymmTensor.H
@@ -245,15 +245,20 @@ public:
//- The 2D determinant by excluding given direction
inline Cmpt det2D(const direction excludeCmpt) const;
- //- Return cofactor matrix
+ //- Return adjunct matrix (transpose of cofactor matrix)
+ inline SymmTensor<Cmpt> adjunct() const;
+
+ //- Return cofactor matrix (transpose of adjunct matrix)
inline SymmTensor<Cmpt> cof() const;
- //- Return 2D cofactor matrix by excluding given direction
- inline SymmTensor<Cmpt> cof2D(const direction excludeCmpt) const;
+ //- Return 2D adjunct matrix by excluding given direction
+ inline SymmTensor<Cmpt> adjunct2D(const direction excludeCmpt) const;
+
+ //- Return inverse
+ inline SymmTensor<Cmpt> inv() const;
- //- Return inverse, with optional (ad hoc) failsafe handling
- //- of 2D tensors
- inline SymmTensor<Cmpt> inv(const bool failsafe=false) const;
+ //- Return inverse, with (ad hoc) failsafe handling of 2D tensors
+ inline SymmTensor<Cmpt> safeInv() const;
//- Return inverse of 2D tensor (by excluding given direction)
inline SymmTensor<Cmpt> inv2D(const direction excludeCmpt) const;
diff --git a/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H b/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H
index 47c10fa13b8fe17c7d1b336f2f2a6238e6a0b3d2..5c28c6c2bc6ff3f1fd3cc6279f57780f9453fe7f 100644
--- a/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H
+++ b/src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H
@@ -272,24 +272,28 @@ inline Cmpt Foam::SymmTensor<Cmpt>::det2D(const direction excludeCmpt) const
template<class Cmpt>
-inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::cof() const
+inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::adjunct() const
{
+ // symmetric: cof() == adjunct()
return SymmTensor<Cmpt>
(
- yy()*zz() - yz()*yz(),
- xz()*yz() - xy()*zz(),
- xy()*yz() - xz()*yy(),
+ yy()*zz() - yz()*yz(), xz()*yz() - xy()*zz(), xy()*yz() - xz()*yy(),
+ xx()*zz() - xz()*xz(), xy()*xz() - xx()*yz(),
+ xx()*yy() - xy()*xy()
+ );
+}
- xx()*zz() - xz()*xz(),
- xy()*xz() - xx()*yz(),
- xx()*yy() - xy()*xy()
- );
+template<class Cmpt>
+inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::cof() const
+{
+ // symmetric: cof() == adjunct()
+ return this->adjunct();
}
template<class Cmpt>
-inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::cof2D
+inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::adjunct2D
(
const direction excludeCmpt
) const
@@ -335,7 +339,84 @@ inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::inv2D
{
const Cmpt detval = this->det2D(excludeCmpt);
- return this->cof2D(excludeCmpt).T()/detval;
+ return this->adjunct2D(excludeCmpt)/detval;
+}
+
+
+// Invert without much error handling
+template<class Cmpt>
+inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::inv() const
+{
+ const Cmpt detval = this->det();
+
+ #ifdef FULLDEBUG
+ if (mag(detval) < VSMALL)
+ {
+ FatalErrorInFunction
+ << "Tensor not properly invertible, determinant:"
+ << detval << " tensor:" << *this << nl
+ << abort(FatalError);
+ }
+ #endif
+
+ return this->adjunct()/detval;
+}
+
+
+// Invert with some error handling
+template<class Cmpt>
+inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::safeInv() const
+{
+ {
+ // Attempt to identify and handle 2-D cases.
+ // - use diagSqr instead of magSqr for fewer operations
+ const scalar magSqr_xx = Foam::magSqr(xx());
+ const scalar magSqr_yy = Foam::magSqr(yy());
+ const scalar magSqr_zz = Foam::magSqr(zz());
+
+ // SMALL: 1e-15 (double), 1e-6 (float), but 1e-6 may be adequate
+
+ const scalar threshold = SMALL * (magSqr_xx + magSqr_yy + magSqr_zz);
+
+ const bool small_xx = (magSqr_xx < threshold);
+ const bool small_yy = (magSqr_yy < threshold);
+ const bool small_zz = (magSqr_zz < threshold);
+
+ if (small_xx || small_yy || small_zz)
+ {
+ SymmTensor<Cmpt> work(*this);
+
+ if (small_xx) { work.xx() += pTraits<Cmpt>::one; }
+ if (small_yy) { work.yy() += pTraits<Cmpt>::one; }
+ if (small_zz) { work.zz() += pTraits<Cmpt>::one; }
+
+ const Cmpt detval = work.det();
+
+ if (mag(detval) < ROOTVSMALL)
+ {
+ // Appears to be nearly zero - leave untouched?
+ return SymmTensor<Cmpt>(Zero);
+ }
+
+ work = work.adjunct()/detval;
+
+ if (small_xx) { work.xx() -= pTraits<Cmpt>::one; }
+ if (small_yy) { work.yy() -= pTraits<Cmpt>::one; }
+ if (small_zz) { work.zz() -= pTraits<Cmpt>::one; }
+
+ return work;
+ }
+ }
+
+ const Cmpt detval = this->det();
+
+ if (mag(detval) < ROOTVSMALL)
+ {
+ // Appears to be nearly zero - leave untouched?
+ return SymmTensor<Cmpt>(Zero);
+ }
+
+ return this->adjunct()/detval;
}
@@ -438,7 +519,7 @@ inline SymmTensor<Cmpt> inv(const SymmTensor<Cmpt>& st, const Cmpt detval)
}
#endif
- return st.cof().T()/detval;
+ return st.adjunct()/detval;
}
@@ -446,62 +527,7 @@ inline SymmTensor<Cmpt> inv(const SymmTensor<Cmpt>& st, const Cmpt detval)
template<class Cmpt>
inline SymmTensor<Cmpt> inv(const SymmTensor<Cmpt>& st)
{
- return inv(st, st.det());
-}
-
-
-template<class Cmpt>
-inline SymmTensor<Cmpt> SymmTensor<Cmpt>::inv(const bool failsafe) const
-{
- const Cmpt detval = this->det();
-
- if (failsafe)
- {
- // Attempt to identify and handle 2-D cases.
- // - use diagSqr instead of magSqr for fewer operations
- const scalar magSqr_xx = Foam::magSqr(xx());
- const scalar magSqr_yy = Foam::magSqr(yy());
- const scalar magSqr_zz = Foam::magSqr(zz());
-
- const scalar threshold = SMALL * (magSqr_xx + magSqr_yy + magSqr_zz);
-
- const bool small_xx = (magSqr_xx < threshold);
- const bool small_yy = (magSqr_yy < threshold);
- const bool small_zz = (magSqr_zz < threshold);
-
- if (small_xx || small_yy || small_zz)
- {
- SymmTensor<Cmpt> work(*this);
-
- if (small_xx) { work.xx() += pTraits<Cmpt>::one; }
- if (small_yy) { work.yy() += pTraits<Cmpt>::one; }
- if (small_zz) { work.zz() += pTraits<Cmpt>::one; }
-
- work = Foam::inv(work, work.det());
-
- if (small_xx) { work.xx() -= pTraits<Cmpt>::one; }
- if (small_yy) { work.yy() -= pTraits<Cmpt>::one; }
- if (small_zz) { work.zz() -= pTraits<Cmpt>::one; }
-
- return work;
- }
- else if (mag(detval) < VSMALL)
- {
- // Really appears to be zero - leave untouched?
- return SymmTensor<Cmpt>(Zero);
- }
- }
- #ifdef FULLDEBUG
- else if (mag(detval) < VSMALL)
- {
- FatalErrorInFunction
- << "Tensor not properly invertible, determinant:"
- << detval << " tensor:" << *this << nl
- << abort(FatalError);
- }
- #endif
-
- return this->cof().T()/detval;
+ return st.inv();
}
diff --git a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
index 8869e9fb56a4b0a55e0bf974beabfb65335ccc3a..fb41001b1f8009a25c3a12572172d697bd4b794c 100644
--- a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
+++ b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
@@ -340,9 +340,9 @@ Foam::tensor Foam::eigenVectors(const symmTensor& T)
Foam::symmTensor Foam::pinv(const symmTensor& st)
{
- const scalar dt = det(st);
+ const scalar detval = st.det();
- if (dt < ROOTVSMALL)
+ if (detval < ROOTVSMALL)
{
// Fall back to pseudo inverse
scalarRectangularMatrix mat(3, 3);
@@ -361,7 +361,7 @@ Foam::symmTensor Foam::pinv(const symmTensor& st)
);
}
- return inv(st, dt);
+ return Foam::inv(st, detval);
}
diff --git a/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2D.H b/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2D.H
index 1308b015e9d851e74a768f2dac237c0680d988b9..e70a951c9b9c1b2419fb5dd67e33d8051162e8a9 100644
--- a/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2D.H
+++ b/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2D.H
@@ -152,7 +152,10 @@ public:
//- The determinate
inline Cmpt det() const;
- //- Return cofactor matrix
+ //- Return adjunct matrix (transpose of cofactor matrix)
+ inline SymmTensor2D<Cmpt> adjunct() const;
+
+ //- Return cofactor matrix (transpose of adjunct matrix)
inline SymmTensor2D<Cmpt> cof() const;
//- Return inverse
diff --git a/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2DI.H b/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2DI.H
index b583b669c5ef7a3c37fe289c9c0ec07f70c0f81f..4eba6e6cc3bcc5fddee58091313f7137fa209867 100644
--- a/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2DI.H
+++ b/src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2DI.H
@@ -109,16 +109,45 @@ inline Cmpt Foam::SymmTensor2D<Cmpt>::det() const
template<class Cmpt>
-inline Foam::SymmTensor2D<Cmpt> Foam::SymmTensor2D<Cmpt>::cof() const
+inline Foam::SymmTensor2D<Cmpt> Foam::SymmTensor2D<Cmpt>::adjunct() const
{
+ // symmetric: cof() == adjunct()
return SymmTensor2D<Cmpt>
(
- yy(), -yx(),
+ yy(), -xy(),
xx()
);
}
+template<class Cmpt>
+inline Foam::SymmTensor2D<Cmpt> Foam::SymmTensor2D<Cmpt>::cof() const
+{
+ // symmetric: cof() == adjunct()
+ return this->adjunct();
+}
+
+
+// Invert without much error handling
+template<class Cmpt>
+inline Foam::SymmTensor2D<Cmpt> Foam::SymmTensor2D<Cmpt>::inv() const
+{
+ const Cmpt detval = this->det();
+
+ #ifdef FULLDEBUG
+ if (mag(detval) < SMALL)
+ {
+ FatalErrorInFunction
+ << "SymmTensor2D not properly invertible, determinant:"
+ << detval << " tensor:" << *this << nl
+ << abort(FatalError);
+ }
+ #endif
+
+ return this->adjunct()/detval;
+}
+
+
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Cmpt>
@@ -220,7 +249,7 @@ inline SymmTensor2D<Cmpt> inv(const SymmTensor2D<Cmpt>& st, const Cmpt detval)
}
#endif
- return st.cof().T()/detval;
+ return st.adjunct()/detval;
}
@@ -228,15 +257,7 @@ inline SymmTensor2D<Cmpt> inv(const SymmTensor2D<Cmpt>& st, const Cmpt detval)
template<class Cmpt>
inline SymmTensor2D<Cmpt> inv(const SymmTensor2D<Cmpt>& st)
{
- return inv(st, st.det());
-}
-
-
-// The inverse of this SymmTensor2D
-template<class Cmpt>
-inline SymmTensor2D<Cmpt> SymmTensor2D<Cmpt>::inv() const
-{
- return Foam::inv(*this, this->det());
+ return st.inv();
}
diff --git a/src/OpenFOAM/primitives/Tensor/Tensor.H b/src/OpenFOAM/primitives/Tensor/Tensor.H
index 7391d06865bd8d02d4a58f5dc19058adbd549e7b..601e00441674f65891255b40e3bd2bc542d2bcaa 100644
--- a/src/OpenFOAM/primitives/Tensor/Tensor.H
+++ b/src/OpenFOAM/primitives/Tensor/Tensor.H
@@ -287,15 +287,20 @@ public:
//- The 2D determinant by excluding given direction
inline Cmpt det2D(const direction excludeCmpt) const;
- //- Return cofactor matrix
+ //- Return adjunct matrix (transpose of cofactor matrix)
+ inline Tensor<Cmpt> adjunct() const;
+
+ //- Return cofactor matrix (transpose of adjunct matrix)
inline Tensor<Cmpt> cof() const;
- //- Return 2D cofactor matrix by excluding given direction
- inline Tensor<Cmpt> cof2D(const direction excludeCmpt) const;
+ //- Return 2D adjunct matrix by excluding given direction
+ inline Tensor<Cmpt> adjunct2D(const direction excludeCmpt) const;
+
+ //- Return inverse
+ inline Tensor<Cmpt> inv() const;
- //- Return inverse, with optional (ad hoc) failsafe handling
- //- of 2D tensors
- inline Tensor<Cmpt> inv(const bool failsafe=false) const;
+ //- Return inverse, with (ad hoc) failsafe handling of 2D tensors
+ inline Tensor<Cmpt> safeInv() const;
//- Return inverse of 2D tensor (by excluding given direction)
inline Tensor<Cmpt> inv2D(const direction excludeCmpt) const;
diff --git a/src/OpenFOAM/primitives/Tensor/TensorI.H b/src/OpenFOAM/primitives/Tensor/TensorI.H
index 6e24c0f506cb408dacf6a927fb7d868c485cce0f..b389f12ee13596c278cec9ad40f6c3b6d795ae72 100644
--- a/src/OpenFOAM/primitives/Tensor/TensorI.H
+++ b/src/OpenFOAM/primitives/Tensor/TensorI.H
@@ -468,27 +468,26 @@ inline Cmpt Foam::Tensor<Cmpt>::det2D(const direction excludeCmpt) const
template<class Cmpt>
-inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::cof() const
+inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::adjunct() const
{
return Tensor<Cmpt>
(
- yy()*zz() - zy()*yz(),
- zx()*yz() - yx()*zz(),
- yx()*zy() - yy()*zx(),
+ yy()*zz() - zy()*yz(), xz()*zy() - xy()*zz(), xy()*yz() - xz()*yy(),
+ zx()*yz() - yx()*zz(), xx()*zz() - xz()*zx(), yx()*xz() - xx()*yz(),
+ yx()*zy() - yy()*zx(), xy()*zx() - xx()*zy(), xx()*yy() - yx()*xy()
+ );
+}
- xz()*zy() - xy()*zz(),
- xx()*zz() - xz()*zx(),
- xy()*zx() - xx()*zy(),
- xy()*yz() - xz()*yy(),
- yx()*xz() - xx()*yz(),
- xx()*yy() - yx()*xy()
- );
+template<class Cmpt>
+inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::cof() const
+{
+ return this->adjunct().T();
}
template<class Cmpt>
-inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::cof2D
+inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::adjunct2D
(
const direction excludeCmpt
) const
@@ -500,8 +499,8 @@ inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::cof2D
return Tensor<Cmpt>
(
Zero, Zero, Zero,
- Zero, zz(), -zy(),
- Zero, -yz(), yy()
+ Zero, zz(), -yz(),
+ Zero, -zy(), yy()
);
}
@@ -509,9 +508,9 @@ inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::cof2D
{
return Tensor<Cmpt>
(
- zz(), Zero, -zx(),
+ zz(), Zero, -xz(),
Zero, Zero, Zero,
- -xz(), Zero, xx()
+ -zx(), Zero, xx()
);
}
}
@@ -519,8 +518,8 @@ inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::cof2D
// Fall-through: Eliminate z
return Tensor<Cmpt>
(
- yy(), -yx(), Zero,
- -xy(), xx(), Zero,
+ yy(), -xy(), Zero,
+ -yx(), xx(), Zero,
Zero, Zero, Zero
);
}
@@ -534,7 +533,7 @@ inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::inv2D
{
const Cmpt detval = this->det2D(excludeCmpt);
- return this->cof2D(excludeCmpt).T()/detval;
+ return this->adjunct2D(excludeCmpt)/detval;
}
@@ -576,6 +575,84 @@ Foam::Tensor<Cmpt>::schur(const Tensor<Cmpt>& t2) const
}
+// Invert without much error handling
+template<class Cmpt>
+inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::inv() const
+{
+ const Cmpt detval = this->det();
+
+ #ifdef FULLDEBUG
+ if (mag(detval) < VSMALL)
+ {
+ FatalErrorInFunction
+ << "Tensor not properly invertible, determinant:"
+ << detval << " tensor:" << *this << nl
+ << abort(FatalError);
+ }
+ #endif
+
+ return this->adjunct()/detval;
+}
+
+
+// Invert with some error handling
+template<class Cmpt>
+inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::safeInv() const
+{
+ {
+ // Attempt to identify and handle 2-D cases
+ // - use diagSqr instead of magSqr for fewer operations
+
+ const scalar magSqr_xx = Foam::magSqr(xx());
+ const scalar magSqr_yy = Foam::magSqr(yy());
+ const scalar magSqr_zz = Foam::magSqr(zz());
+
+ // SMALL: 1e-15 (double), 1e-6 (float), but 1e-6 may be adequate
+
+ const scalar threshold = SMALL * (magSqr_xx + magSqr_yy + magSqr_zz);
+
+ const bool small_xx = (magSqr_xx < threshold);
+ const bool small_yy = (magSqr_yy < threshold);
+ const bool small_zz = (magSqr_zz < threshold);
+
+ if (small_xx || small_yy || small_zz)
+ {
+ Tensor<Cmpt> work(*this);
+
+ if (small_xx) { work.xx() += pTraits<Cmpt>::one; }
+ if (small_yy) { work.yy() += pTraits<Cmpt>::one; }
+ if (small_zz) { work.zz() += pTraits<Cmpt>::one; }
+
+ const Cmpt detval = work.det();
+
+ if (mag(detval) < ROOTVSMALL)
+ {
+ // Appears to be nearly zero - leave untouched?
+ return Tensor<Cmpt>(Zero);
+ }
+
+ work = work.adjunct()/detval;
+
+ if (small_xx) { work.xx() -= pTraits<Cmpt>::one; }
+ if (small_yy) { work.yy() -= pTraits<Cmpt>::one; }
+ if (small_zz) { work.zz() -= pTraits<Cmpt>::one; }
+
+ return work;
+ }
+ }
+
+ const Cmpt detval = this->det();
+
+ if (mag(detval) < ROOTVSMALL)
+ {
+ // Appears to be nearly zero - leave untouched?
+ return Tensor<Cmpt>(Zero);
+ }
+
+ return this->adjunct()/detval;
+}
+
+
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Cmpt>
@@ -750,7 +827,7 @@ inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t, const Cmpt detval)
}
#endif
- return t.cof().T()/detval;
+ return t.adjunct()/detval;
}
@@ -758,63 +835,7 @@ inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t, const Cmpt detval)
template<class Cmpt>
inline Tensor<Cmpt> inv(const Tensor<Cmpt>& t)
{
- return inv(t, t.det());
-}
-
-
-template<class Cmpt>
-inline Tensor<Cmpt> Tensor<Cmpt>::inv(const bool failsafe) const
-{
- const Cmpt detval = this->det();
-
- if (failsafe)
- {
- // Attempt to identify and handle 2-D cases
- // - use diagSqr instead of magSqr for fewer operations
-
- const scalar magSqr_xx = Foam::magSqr(xx());
- const scalar magSqr_yy = Foam::magSqr(yy());
- const scalar magSqr_zz = Foam::magSqr(zz());
-
- const scalar threshold = SMALL * (magSqr_xx + magSqr_yy + magSqr_zz);
-
- const bool small_xx = (magSqr_xx < threshold);
- const bool small_yy = (magSqr_yy < threshold);
- const bool small_zz = (magSqr_zz < threshold);
-
- if (small_xx || small_yy || small_zz)
- {
- Tensor<Cmpt> work(*this);
-
- if (small_xx) { work.xx() += pTraits<Cmpt>::one; }
- if (small_yy) { work.yy() += pTraits<Cmpt>::one; }
- if (small_zz) { work.zz() += pTraits<Cmpt>::one; }
-
- work = Foam::inv(work, work.det());
-
- if (small_xx) { work.xx() -= pTraits<Cmpt>::one; }
- if (small_yy) { work.yy() -= pTraits<Cmpt>::one; }
- if (small_zz) { work.zz() -= pTraits<Cmpt>::one; }
-
- return work;
- }
- else if (mag(detval) < VSMALL)
- {
- // Really appears to be zero - leave untouched?
- return Tensor<Cmpt>(Zero);
- }
- }
- #ifdef FULLDEBUG
- else if (mag(detval) < VSMALL)
- {
- FatalErrorInFunction
- << "Tensor not properly invertible, determinant:"
- << detval << " tensor:" << *this << nl
- << abort(FatalError);
- }
- #endif
-
- return this->cof().T()/detval;
+ return t.inv();
}
diff --git a/src/OpenFOAM/primitives/Tensor/floats/tensor.C b/src/OpenFOAM/primitives/Tensor/floats/tensor.C
index dfe5851bd0aea0797ebb56448d492795934ac975..d26aad40a4147af6c05f97316a986534f555d6cf 100644
--- a/src/OpenFOAM/primitives/Tensor/floats/tensor.C
+++ b/src/OpenFOAM/primitives/Tensor/floats/tensor.C
@@ -304,9 +304,9 @@ Foam::Tensor<Foam::complex> Foam::eigenVectors(const tensor& T)
Foam::tensor Foam::pinv(const tensor& t)
{
- const scalar dt = det(t);
+ const scalar detval = t.det();
- if (dt < ROOTVSMALL)
+ if (detval < ROOTVSMALL)
{
// Fall back to pseudo inverse
scalarRectangularMatrix mat(3, 3);
@@ -325,7 +325,7 @@ Foam::tensor Foam::pinv(const tensor& t)
);
}
- return inv(t, dt);
+ return Foam::inv(t, detval);
}
diff --git a/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H b/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H
index f12b68104f9ff89b632d28e31fdc25fe655918f1..69f8efd09111b440d4cd91bc300602a88540690f 100644
--- a/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H
+++ b/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H
@@ -212,7 +212,10 @@ public:
//- The determinate
inline Cmpt det() const;
- //- Return cofactor matrix
+ //- Return adjunct matrix (transpose of cofactor matrix)
+ inline Tensor2D<Cmpt> adjunct() const;
+
+ //- Return cofactor matrix (transpose of adjunct matrix)
inline Tensor2D<Cmpt> cof() const;
//- Return inverse
diff --git a/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H b/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H
index 9367f4c82bf07a11945b9edc244fa3ce02beb22c..1acbb876a8c78aab927ad61f5d34ff329a6d0578 100644
--- a/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H
+++ b/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H
@@ -318,16 +318,23 @@ inline Cmpt Foam::Tensor2D<Cmpt>::det() const
template<class Cmpt>
-inline Foam::Tensor2D<Cmpt> Foam::Tensor2D<Cmpt>::cof() const
+inline Foam::Tensor2D<Cmpt> Foam::Tensor2D<Cmpt>::adjunct() const
{
return Tensor2D<Cmpt>
(
- yy(), -yx(),
- -xy(), xx()
+ yy(), -xy(),
+ -yx(), xx()
);
}
+template<class Cmpt>
+inline Foam::Tensor2D<Cmpt> Foam::Tensor2D<Cmpt>::cof() const
+{
+ return this->adjunct().T();
+}
+
+
template<class Cmpt>
inline Foam::Tensor2D<Cmpt>
Foam::Tensor2D<Cmpt>::inner(const Tensor2D<Cmpt>& t2) const
@@ -359,6 +366,26 @@ Foam::Tensor2D<Cmpt>::schur(const Tensor2D<Cmpt>& t2) const
}
+// Invert without much error handling
+template<class Cmpt>
+inline Foam::Tensor2D<Cmpt> Foam::Tensor2D<Cmpt>::inv() const
+{
+ const Cmpt detval = this->det();
+
+ #ifdef FULLDEBUG
+ if (mag(detval) < SMALL)
+ {
+ FatalErrorInFunction
+ << "SymmTensor2D not properly invertible, determinant:"
+ << detval << " tensor:" << *this << nl
+ << abort(FatalError);
+ }
+ #endif
+
+ return this->adjunct()/detval;
+}
+
+
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Cmpt>
@@ -463,7 +490,7 @@ inline Cmpt det(const Tensor2D<Cmpt>& t)
}
-//- Return the cofactor Tensor2D of a Tensor2D
+//- Return the cofactor of a Tensor2D
template<class Cmpt>
inline Tensor2D<Cmpt> cof(const Tensor2D<Cmpt>& t)
{
@@ -485,7 +512,7 @@ inline Tensor2D<Cmpt> inv(const Tensor2D<Cmpt>& t, const Cmpt detval)
}
#endif
- return t.cof().T()/detval;
+ return t.adjunct()/detval;
}
@@ -493,15 +520,7 @@ inline Tensor2D<Cmpt> inv(const Tensor2D<Cmpt>& t, const Cmpt detval)
template<class Cmpt>
inline Tensor2D<Cmpt> inv(const Tensor2D<Cmpt>& t)
{
- return inv(t, t.det());
-}
-
-
-// Return the inverse of this Tensor2D
-template<class Cmpt>
-inline Tensor2D<Cmpt> Tensor2D<Cmpt>::inv() const
-{
- return Foam::inv(*this, this->det());
+ return t.inv();
}
diff --git a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C
index f2bd352bbd066ac922be82dccedd009fc0029b81..e40618b0a079ab0ce61629b10497677c3a601a07 100644
--- a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C
+++ b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C
@@ -148,26 +148,8 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
}
- // Invert the dd tensor.
-
- // Cannot rely on the usual field inv() since that only uses the
- // first element to guess if the remaining tensors are 2D (or
- // singular). We, however, can have a mixture (eg, skipping over
- // zero-length edges can yield a zero). Instead use the
- // 'failsafe' inv() that checks each tensor (#2724)
-
- // Fragile: const symmTensorField invDd(inv(dd));
-
- symmTensorField invDd(dd.size());
- {
- const label loopLen = dd.size();
-
- /* pragmas... */
- for (label i = 0; i < loopLen; ++i)
- {
- invDd[i] = dd[i].inv(true); // With 'failsafe' handling
- }
- }
+ // Invert the dd tensors - including failsafe checks
+ const symmTensorField invDd(inv(dd));
// Revisit all faces and calculate the lsP and lsN vectors
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/invDistLeastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/invDistLeastSquaresVectors.C
index 5bfe9a515fa72299078fa47e07e6b4b34f3d79e2..584968db7d3fcbda94907c84a69e3f2a384eb850 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/invDistLeastSquaresVectors.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/invDistLeastSquaresVectors.C
@@ -104,14 +104,17 @@ void Foam::leastSquaresVectors::calcLeastSquaresVectors()
const label nei = neighbour[facei];
const vector d(C[nei] - C[own]);
- const symmTensor wdd(sqr(d)/magSqr(d));
- dd[own] += wdd;
- dd[nei] += wdd;
- }
+ const scalar magSqrDist = d.magSqr();
+ if (magSqrDist > ROOTVSMALL)
+ {
+ const symmTensor wdd(sqr(d)/magSqrDist);
+ dd[own] += wdd;
+ dd[nei] += wdd;
+ }
+ }
- surfaceVectorField::Boundary& blsP =
- pVectors_.boundaryField();
+ auto& blsP = pVectors_.boundaryField();
forAll(blsP, patchi)
{
@@ -126,13 +129,17 @@ void Foam::leastSquaresVectors::calcLeastSquaresVectors()
forAll(pd, patchFacei)
{
const vector& d = pd[patchFacei];
+ const scalar magSqrDist = d.magSqr();
- dd[faceCells[patchFacei]] += sqr(d)/magSqr(d);
+ if (magSqrDist > ROOTVSMALL)
+ {
+ dd[faceCells[patchFacei]] += sqr(d)/magSqrDist;
+ }
}
}
- // Invert the dd tensor
+ // Invert the dd tensors - including failsafe checks
const symmTensorField invDd(inv(dd));
@@ -143,9 +150,18 @@ void Foam::leastSquaresVectors::calcLeastSquaresVectors()
const label nei = neighbour[facei];
const vector d(C[nei] - C[own]);
+ const scalar magSqrDist = d.magSqr();
- pVectors_[facei] = (invDd[own] & d)/magSqr(d);
- nVectors_[facei] = -(invDd[nei] & d)/magSqr(d);
+ if (magSqrDist > ROOTVSMALL)
+ {
+ pVectors_[facei] = (invDd[own] & d)/magSqrDist;
+ nVectors_[facei] = -(invDd[nei] & d)/magSqrDist;
+ }
+ else
+ {
+ pVectors_[facei] = Zero;
+ nVectors_[facei] = Zero;
+ }
}
forAll(blsP, patchi)
@@ -161,8 +177,17 @@ void Foam::leastSquaresVectors::calcLeastSquaresVectors()
forAll(pd, patchFacei)
{
const vector& d = pd[patchFacei];
-
- patchLsP[patchFacei] = (invDd[faceCells[patchFacei]] & d)/magSqr(d);
+ const scalar magSqrDist = d.magSqr();
+
+ if (magSqrDist > ROOTVSMALL)
+ {
+ patchLsP[patchFacei] =
+ (invDd[faceCells[patchFacei]] & d)/magSqrDist;
+ }
+ else
+ {
+ patchLsP[patchFacei] = Zero;
+ }
}
}
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
index 40dfc42445da7ead4c948bdf416681c703e0f8c9..c8f5b010511ed8372bc579ace2b020645988d932 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
@@ -151,7 +151,7 @@ void Foam::leastSquaresVectors::calcLeastSquaresVectors()
}
- // Invert the dd tensor
+ // Invert the dd tensors - including failsafe checks
const symmTensorField invDd(inv(dd));
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/unweightedLeastSquaresVectors.C b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/unweightedLeastSquaresVectors.C
index 20c8bf54060c120581d5477e44e10a44f4f5990f..c19c1ba86b29734578404c0cbe4e4bd0336acbc8 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/unweightedLeastSquaresVectors.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/unweightedLeastSquaresVectors.C
@@ -129,7 +129,7 @@ void Foam::leastSquaresVectors::calcLeastSquaresVectors()
}
- // Invert the dd tensor
+ // Invert the dd tensors - including failsafe checks
const symmTensorField invDd(inv(dd));