diff --git a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
index f68690da9b344ccd35f70f5abab33c1a78a4e5ab..0699c2e85e1cb290a64d77b4ee8b2b8203da8189 100644
--- a/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
+++ b/src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
- Copyright (C) 2019 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -50,79 +50,7 @@ UNARY_FUNCTION(symmTensor, symmTensor, dev)
UNARY_FUNCTION(symmTensor, symmTensor, dev2)
UNARY_FUNCTION(scalar, symmTensor, det)
UNARY_FUNCTION(symmTensor, symmTensor, cof)
-
-void inv(Field<symmTensor>& tf, const UList<symmTensor>& tf1)
-{
- if (tf.empty())
- {
- return;
- }
-
- scalar scale = magSqr(tf1[0]);
- Vector<bool> removeCmpts
- (
- magSqr(tf1[0].xx())/scale < SMALL,
- magSqr(tf1[0].yy())/scale < SMALL,
- magSqr(tf1[0].zz())/scale < SMALL
- );
-
- if (removeCmpts.x() || removeCmpts.y() || removeCmpts.z())
- {
- symmTensorField tf1Plus(tf1);
-
- if (removeCmpts.x())
- {
- tf1Plus += symmTensor(1,0,0,0,0,0);
- }
-
- if (removeCmpts.y())
- {
- tf1Plus += symmTensor(0,0,0,1,0,0);
- }
-
- if (removeCmpts.z())
- {
- tf1Plus += symmTensor(0,0,0,0,0,1);
- }
-
- TFOR_ALL_F_OP_FUNC_F(symmTensor, tf, =, inv, symmTensor, tf1Plus)
-
- if (removeCmpts.x())
- {
- tf -= symmTensor(1,0,0,0,0,0);
- }
-
- if (removeCmpts.y())
- {
- tf -= symmTensor(0,0,0,1,0,0);
- }
-
- if (removeCmpts.z())
- {
- tf -= symmTensor(0,0,0,0,0,1);
- }
- }
- else
- {
- TFOR_ALL_F_OP_FUNC_F(symmTensor, tf, =, inv, symmTensor, tf1)
- }
-}
-
-tmp<symmTensorField> inv(const UList<symmTensor>& tf)
-{
- auto tresult = tmp<symmTensorField>::New(tf.size());
- inv(tresult.ref(), tf);
- return tresult;
-}
-
-tmp<symmTensorField> inv(const tmp<symmTensorField>& tf)
-{
- tmp<symmTensorField> tresult = New(tf);
- inv(tresult.ref(), tf());
- tf.clear();
- return tresult;
-}
-
+UNARY_FUNCTION(symmTensor, symmTensor, inv)
template<>
tmp<Field<symmTensor>> transformFieldMask<symmTensor>
diff --git a/src/OpenFOAM/fields/Fields/tensorField/tensorField.C b/src/OpenFOAM/fields/Fields/tensorField/tensorField.C
index 37df235b320afe8d94cef77f614f88482cf9bfb3..9522fa6308d9efe9ac545b3948f4a950943e1489 100644
--- a/src/OpenFOAM/fields/Fields/tensorField/tensorField.C
+++ b/src/OpenFOAM/fields/Fields/tensorField/tensorField.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
- Copyright (C) 2019-2020 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -49,79 +49,7 @@ UNARY_FUNCTION(tensor, tensor, dev)
UNARY_FUNCTION(tensor, tensor, dev2)
UNARY_FUNCTION(scalar, tensor, det)
UNARY_FUNCTION(tensor, tensor, cof)
-
-void inv(Field<tensor>& tf, const UList<tensor>& tf1)
-{
- if (tf.empty())
- {
- return;
- }
-
- scalar scale = magSqr(tf1[0]);
- Vector<bool> removeCmpts
- (
- magSqr(tf1[0].xx())/scale < SMALL,
- magSqr(tf1[0].yy())/scale < SMALL,
- magSqr(tf1[0].zz())/scale < SMALL
- );
-
- if (removeCmpts.x() || removeCmpts.y() || removeCmpts.z())
- {
- tensorField tf1Plus(tf1);
-
- if (removeCmpts.x())
- {
- tf1Plus += tensor(1,0,0,0,0,0,0,0,0);
- }
-
- if (removeCmpts.y())
- {
- tf1Plus += tensor(0,0,0,0,1,0,0,0,0);
- }
-
- if (removeCmpts.z())
- {
- tf1Plus += tensor(0,0,0,0,0,0,0,0,1);
- }
-
- TFOR_ALL_F_OP_FUNC_F(tensor, tf, =, inv, tensor, tf1Plus)
-
- if (removeCmpts.x())
- {
- tf -= tensor(1,0,0,0,0,0,0,0,0);
- }
-
- if (removeCmpts.y())
- {
- tf -= tensor(0,0,0,0,1,0,0,0,0);
- }
-
- if (removeCmpts.z())
- {
- tf -= tensor(0,0,0,0,0,0,0,0,1);
- }
- }
- else
- {
- TFOR_ALL_F_OP_FUNC_F(tensor, tf, =, inv, tensor, tf1)
- }
-}
-
-tmp<tensorField> inv(const UList<tensor>& tf)
-{
- auto tres = tmp<tensorField>::New(tf.size());
- inv(tres.ref(), tf);
- return tres;
-}
-
-tmp<tensorField> inv(const tmp<tensorField>& tf)
-{
- auto tres = New(tf);
- inv(tres.ref(), tf());
- tf.clear();
- return tres;
-}
-
+UNARY_FUNCTION(tensor, tensor, inv)
UNARY_FUNCTION(vector, symmTensor, eigenValues)
UNARY_FUNCTION(tensor, symmTensor, eigenVectors)
diff --git a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
index e8212aed3c1e99994e9a895f6db92486f99477b1..2dc1de226b060bf55342e69fb767fed95f1c7ab5 100644
--- a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
+++ b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
- Copyright (C) 2020 OpenCFD Ltd.
+ Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -29,6 +29,7 @@ License
#include "symmTensor.H"
#include "cubicEqn.H"
#include "mathematicalConstants.H"
+#include "SVD.H"
using namespace Foam::constant::mathematical;
@@ -337,4 +338,34 @@ Foam::tensor Foam::eigenVectors(const symmTensor& T)
}
+Foam::symmTensor Foam::inv(const symmTensor& st)
+{
+ const scalar dt = det(st);
+
+ if (dt < ROOTVSMALL)
+ {
+ // Fall back to pseudo inverse
+ scalarRectangularMatrix pinv(3, 3, Zero);
+
+ pinv(0,0) = st.xx();
+ pinv(0,1) = st.xy();
+ pinv(0,2) = st.xz();
+ pinv(1,1) = st.yy();
+ pinv(1,2) = st.yz();
+ pinv(2,2) = st.zz();
+
+ pinv = SVDinv(pinv);
+
+ return symmTensor
+ (
+ pinv(0,0), pinv(0,1), pinv(0,2),
+ pinv(1,1), pinv(1,2),
+ pinv(2,2)
+ );
+ }
+
+ return inv(st, dt);
+}
+
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.H b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.H
index 9625b332ca17fe370efc895393e918ba686ff5eb..338ff2f7d92a0436bd549d87e11a7bc6134ffd35 100644
--- a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.H
+++ b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2012 OpenFOAM Foundation
- Copyright (C) 2019-2020 OpenCFD Ltd.
+ Copyright (C) 2019-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -107,6 +107,14 @@ tensor eigenVectors
tensor eigenVectors(const symmTensor& T);
+//- Return inverse of a given symmTensor, and fall back
+//- to pseudo-inverse if the symmTensor is singular
+// \param st symmTensor
+//
+// \return symmTensor inverse of symmTensor
+symmTensor inv(const symmTensor& st);
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
diff --git a/src/OpenFOAM/primitives/Tensor/floats/tensor.C b/src/OpenFOAM/primitives/Tensor/floats/tensor.C
index 3916def7608d4d817a4d4e91a2617c9bc5a1cdf7..ced500c18d7e6148adc956856cf9910ab1cadcf2 100644
--- a/src/OpenFOAM/primitives/Tensor/floats/tensor.C
+++ b/src/OpenFOAM/primitives/Tensor/floats/tensor.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
- Copyright (C) 2020 OpenCFD Ltd.
+ Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -29,6 +29,7 @@ License
#include "tensor.H"
#include "cubicEqn.H"
#include "mathematicalConstants.H"
+#include "SVD.H"
using namespace Foam::constant::mathematical;
@@ -301,4 +302,29 @@ Foam::Tensor<Foam::complex> Foam::eigenVectors(const tensor& T)
}
+Foam::tensor Foam::inv(const tensor& t)
+{
+ const scalar dt = det(t);
+
+ if (dt < ROOTVSMALL)
+ {
+ // Fall back to pseudo inverse
+ scalarRectangularMatrix pinv(3, 3);
+
+ std::copy(t.cbegin(), t.cend(), pinv.begin());
+
+ pinv = SVDinv(pinv);
+
+ return tensor
+ (
+ pinv(0,0), pinv(0,1), pinv(0,2),
+ pinv(1,0), pinv(1,1), pinv(1,2),
+ pinv(2,0), pinv(2,1), pinv(2,2)
+ );
+ }
+
+ return inv(t, dt);
+}
+
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Tensor/floats/tensor.H b/src/OpenFOAM/primitives/Tensor/floats/tensor.H
index a2724ad568bcd86a7da263455e433a83de14c722..e1c289b53e0a718ecd3dc05d46660b06fa61e886 100644
--- a/src/OpenFOAM/primitives/Tensor/floats/tensor.H
+++ b/src/OpenFOAM/primitives/Tensor/floats/tensor.H
@@ -117,6 +117,14 @@ Tensor<complex> eigenVectors
Tensor<complex> eigenVectors(const tensor& T);
+//- Return inverse of a given tensor, and fall back
+//- into pseudo-inverse if the tensor is singular
+// \param t tensor
+//
+// \return tensor inverse of tensor
+tensor inv(const tensor& t);
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam