songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ENH: robuster handling of inv() of fields with singular tensors (#2724)

Browse Source 
- replace the "one-size-fits-all" approach of tensor field inv()
  with individual 'failsafe' inverts.

  The inv() field function historically just checked the first entry
  to detect 2D cases and adjusted/readjusted *all* tensors accordingly
  (to avoid singularity tensors and/or noisy inversions).

  This seems to have worked reasonably well with 3D volume meshes, but
  breaks down for 2D area meshes, which can be axis-aligned
  differently on different sections of the mesh.
This commit is contained in:
Mark Olesen 2023-04-13 11:25:03 +02:00
parent 3947f3c441
commit a96dcf706b
21 changed files with 446 additions and 326 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
applications/test/invTensor/Make/files Normal file
View File

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

2
applications/test/invTensor/Make/options Normal file
View File

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

96
applications/test/invTensor/Test-invTensor.C Normal file
View File

@ -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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

4
src/OpenFOAM/dimensionedTypes/dimensionedSymmTensor/dimensionedSymmTensor.H
View File

@ -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&);

5
src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.H
View File

@ -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&);

51
src/OpenFOAM/fields/Fields/symmTensorField/symmTensorField.C
View File

@ -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)

51
src/OpenFOAM/fields/Fields/tensorField/tensorField.C
View File

@ -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)

17
src/OpenFOAM/primitives/SymmTensor/SymmTensor.H
View File

@ -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, with optional (ad hoc) failsafe handling
//- of 2D tensors
inline SymmTensor<Cmpt> inv(const bool failsafe=false) const;
//- Return inverse
inline SymmTensor<Cmpt> inv() 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;

162
src/OpenFOAM/primitives/SymmTensor/SymmTensorI.H
View File

@ -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(),
xx()*zz() - xz()*xz(),
xy()*xz() - xx()*yz(),
xx()*yy() - xy()*xy()
yy()*zz() - yz()*yz(), xz()*yz() - xy()*zz(), xy()*yz() - xz()*yy(),
xx()*zz() - xz()*xz(), xy()*xz() - xx()*yz(),
xx()*yy() - xy()*xy()
);
}
template<class Cmpt>
inline Foam::SymmTensor<Cmpt> Foam::SymmTensor<Cmpt>::cof2D
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>::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();
}

6
src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
View File

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

5
src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2D.H
View File

@ -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

45
src/OpenFOAM/primitives/SymmTensor2D/SymmTensor2DI.H
View File

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

17
src/OpenFOAM/primitives/Tensor/Tensor.H
View File

@ -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, with optional (ad hoc) failsafe handling
//- of 2D tensors
inline Tensor<Cmpt> inv(const bool failsafe=false) const;
//- Return inverse
inline Tensor<Cmpt> inv() 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;

177
src/OpenFOAM/primitives/Tensor/TensorI.H
View File

@ -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(),
xz()*zy() - xy()*zz(),
xx()*zz() - xz()*zx(),
xy()*zx() - xx()*zy(),
xy()*yz() - xz()*yy(),
yx()*xz() - xx()*yz(),
xx()*yy() - yx()*xy()
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()
);
}
template<class Cmpt>
inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::cof2D
inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::cof() const
{
return this->adjunct().T();
}
template<class Cmpt>
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();
}

6
src/OpenFOAM/primitives/Tensor/floats/tensor.C
View File

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

5
src/OpenFOAM/primitives/Tensor2D/Tensor2D.H
View File

@ -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

47
src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H
View File

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

22
src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C
View File

@ -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

47
src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/invDistLeastSquaresVectors.C
View File

@ -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];
const scalar magSqrDist = d.magSqr();
patchLsP[patchFacei] = (invDd[faceCells[patchFacei]] & d)/magSqr(d);
if (magSqrDist > ROOTVSMALL)
{
patchLsP[patchFacei] =
(invDd[faceCells[patchFacei]] & d)/magSqrDist;
}
else
{
patchLsP[patchFacei] = Zero;
}
}
}

2
src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/leastSquaresVectors.C
View File

@ -151,7 +151,7 @@ void Foam::leastSquaresVectors::calcLeastSquaresVectors()
}
// Invert the dd tensor
// Invert the dd tensors - including failsafe checks
const symmTensorField invDd(inv(dd));

2
src/finiteVolume/finiteVolume/gradSchemes/leastSquaresGrad/unweightedLeastSquaresVectors.C
View File

@ -129,7 +129,7 @@ void Foam::leastSquaresVectors::calcLeastSquaresVectors()
}
// Invert the dd tensor
// Invert the dd tensors - including failsafe checks
const symmTensorField invDd(inv(dd));