zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-12-28 03:37:59 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'feature-fa-least-squares' into 'develop'

Browse Source 
ENH: inv: fall back to pseudo-inverse for singular tensors

See merge request Development/openfoam!574
This commit is contained in:
Andrew Heather
2022-12-01 12:09:05 +00:00
parent 0b83d39500 e190df9033
commit 7fe8bdcf99
6 changed files with 80 additions and 151 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

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

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

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

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

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

10
src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.H
View File

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

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

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

8
src/OpenFOAM/primitives/Tensor/floats/tensor.H
View File

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