diff --git a/applications/test/SymmTensor/Test-SymmTensor.C b/applications/test/SymmTensor/Test-SymmTensor.C index 5e48ef7c34..2d1aaa970a 100644 --- a/applications/test/SymmTensor/Test-SymmTensor.C +++ b/applications/test/SymmTensor/Test-SymmTensor.C @@ -512,6 +512,78 @@ void test_global_opers(Type) } +// Return false if given eigenvalues fail to satisy eigenvalue relations +// Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307) +void test_eigenvalues(const symmTensor& sT, const vector& EVals) +{ + { + const scalar determinant = det(sT); + const scalar EValsProd = EVals.x()*EVals.y()*EVals.z(); + cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant); + } + { + const scalar trace = tr(sT); + scalar EValsSum = 0.0; + for (const auto& val : EVals) + { + EValsSum += val; + } + cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace); + } +} + + +// Return false if a given eigenvalue-eigenvector pair +// fails to satisfy the characteristic equation +void test_characteristic_equation +( + const symmTensor& sT, + const vector& EVals, + const tensor& EVecs +) +{ + Info<< "# Characteristic equation:" << nl; + for (direction dir = 0; dir < pTraits::nComponents; ++dir) + { + Info<< "EVal = " << EVals[dir] << nl + << "EVec = " << EVecs.row(dir) << endl; + const vector leftSide(sT & EVecs.row(dir)); + const vector rightSide(EVals[dir]*EVecs.row(dir)); + const vector X(leftSide - rightSide); + + for (const auto x : X) + { + cmp(" (sT & EVec - EVal*EVec) = 0:", x, 0.0, 1e-8, 1e-6); + } + } +} + + +// Return false if the eigen functions fail to satisfy relations +void test_eigen_funcs(const symmTensor& sT) +{ + Info<< "# Operand:" << nl + << " symmTensor = " << sT << nl; + + + Info<< "# Return eigenvalues of a given symmTensor:" << nl; + const vector EVals(eigenValues(sT)); + Info<< EVals << endl; + test_eigenvalues(sT, EVals); + + Info<< "# Return eigenvectors of a given symmTensor corresponding to" + << " given eigenvalues:" << nl; + const tensor EVecs0(eigenVectors(sT, EVals)); + Info<< EVecs0 << endl; + test_characteristic_equation(sT, EVals, EVecs0); + + Info<< "# Return eigenvectors of a given symmTensor by computing" + << " the eigenvalues of the symmTensor in the background:" << nl; + const tensor EVecs1(eigenVectors(sT)); + Info<< EVecs1 << endl; +} + + // Do compile-time recursion over the given types template inline typename std::enable_if::type @@ -557,6 +629,74 @@ int main() run_tests(types, typeID); + Info<< nl << " ## Test SymmTensor eigen functions: ##" << nl; + const label numberOfTests = 10000; + Random rndGen(1234); + + for (label i = 0; i < numberOfTests; ++i) + { + const symmTensor T(makeRandomContainer(rndGen)); + test_eigen_funcs(T); + } + + { + Info<< nl << " ## Test eigen functions by a zero tensor: ##"<< nl; + const symmTensor zeroT(Zero); + test_eigen_funcs(zeroT); + } + { + Info<< nl + << " ## Test eigen functions by a tensor consisting of" + << " 2 repeated eigenvalues: ##" + << nl; + const symmTensor T + ( + 1.0, 0.0, Foam::sqrt(2.0), + 2.0, 0.0, + 0.0 + ); + test_eigen_funcs(T); + } + { + Info<< nl + << " ## Test eigen functions by a tensor consisting of" + << " 3 repeated eigenvalues: ##" + << nl; + const symmTensor T + ( + 0.023215, -5.0739e-09, -7.0012e-09, + 0.023215, -8.148e-10, + 0.023215 + ); + test_eigen_funcs(T); + } + { + Info<< nl + << " ## Test eigen functions by a tensor consisting of" + << " SMALL off-diagonal elements: ##" + << nl; + + for (label i = 0; i < numberOfTests; ++i) + { + symmTensor T(makeRandomContainer(rndGen)); + T.xy() = SMALL*rndGen.GaussNormal(); + T.xz() = SMALL*rndGen.GaussNormal(); + T.yz() = SMALL*rndGen.GaussNormal(); + test_eigen_funcs(T); + } + } + { + Info<< nl << " ## Test eigen functions by a stiff tensor: ##" << nl; + const symmTensor stiff + ( + pow(10.0, 10), pow(10.0, 8), pow(10.0, -8), + pow(10.0, -8), pow(10.0, 8), + pow(10.0, 7) + ); + test_eigen_funcs(stiff); + } + + if (nFail_) { Info<< nl << " #### " @@ -568,8 +708,6 @@ int main() Info<< nl << " #### Passed all " << nTest_ <<" tests ####\n" << endl; return 0; - - return 0; } diff --git a/applications/test/SymmTensor2D/Test-SymmTensor2D.C b/applications/test/SymmTensor2D/Test-SymmTensor2D.C index 236516fb04..82bbc77745 100644 --- a/applications/test/SymmTensor2D/Test-SymmTensor2D.C +++ b/applications/test/SymmTensor2D/Test-SymmTensor2D.C @@ -476,6 +476,80 @@ void test_global_opers(Type) } +// Return false if given eigenvalues fail to satisy eigenvalue relations +// Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307) +void test_eigenvalues(const symmTensor2D& T, const vector2D& EVals) +{ + { + Info<< "# Product of eigenvalues = det(T):" << nl; + const scalar determinant = det(T); + const scalar EValsProd = EVals.x()*EVals.y(); + cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant); + } + + { + const scalar trace = tr(T); + scalar EValsSum = 0.0; + for (const auto& val : EVals) + { + EValsSum += val; + } + cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace); + } +} + + +// Return false if a given eigenvalue-eigenvector pair +// fails to satisfy the characteristic equation +void test_characteristic_equation +( + const symmTensor2D& T, + const vector2D& EVals, + const tensor2D& EVecs +) +{ + Info<< "# Characteristic equation:" << nl; + for (direction dir = 0; dir < pTraits::nComponents; ++dir) + { + Info<< "EVal = " << EVals[dir] << nl + << "EVec = " << EVecs.row(dir) << nl; + const vector2D leftSide(T & EVecs.row(dir)); + const vector2D rightSide(EVals[dir]*EVecs.row(dir)); + const vector2D X(leftSide - rightSide); + + for (const auto x : X) + { + cmp(" (sT & EVec - EVal*EVec) = 0:", x, 0.0, 1e-8, 1e-6); + } + } +} + + +// Return false if the eigen functions fail to satisfy relations +void test_eigen_funcs(const symmTensor2D& T) +{ + Info<< "# Operand:" << nl + << " symmTensor2D = " << T << nl; + + + Info<< "# Return eigenvalues of a given symmTensor2D:" << nl; + const vector2D EVals(eigenValues(T)); + Info<< EVals << endl; + test_eigenvalues(T, EVals); + + Info<< "# Return eigenvectors of a given symmTensor2D corresponding to" + << " given eigenvalues:" << nl; + const tensor2D EVecs0(eigenVectors(T, EVals)); + Info<< EVecs0 << endl; + test_characteristic_equation(T, EVals, EVecs0); + + Info<< "# Return eigenvectors of a given symmTensor2D by computing" + << " the eigenvalues of the symmTensor2D in the background:" << nl; + const tensor2D EVecs1(eigenVectors(T)); + Info<< EVecs1 << endl; +} + + // Do compile-time recursion over the given types template inline typename std::enable_if::type @@ -521,6 +595,52 @@ int main(int argc, char *argv[]) run_tests(types, typeID); + Info<< nl << " ## Test symmTensor2D eigen functions: ##" << nl; + const label numberOfTests = 10000; + Random rndGen(1234); + + for (label i = 0; i < numberOfTests; ++i) + { + const symmTensor2D T(makeRandomContainer(rndGen)); + test_eigen_funcs(T); + } + + Info<< nl << " ## Test symmTensor2D eigen functions" + << " with T.xy = VSMALL: ##" << nl; + for (label i = 0; i < numberOfTests; ++i) + { + symmTensor2D T(makeRandomContainer(rndGen)); + T.xy() = VSMALL; + test_eigen_funcs(T); + } + + Info<< nl << " ## Test symmTensor2D eigen functions" + << " with T.xx = T.yy: ##" << nl; + for (label i = 0; i < numberOfTests; ++i) + { + symmTensor2D T(makeRandomContainer(rndGen)); + T.xx() = T.yy(); + test_eigen_funcs(T); + } + + { + Info<< nl + << " ## Test eigen functions by a zero symmTensor2D: ##"<< nl; + const symmTensor2D zeroT(Zero); + test_eigen_funcs(zeroT); + } + + { + Info<< nl << " ## Test eigen functions by a stiff tensor2D: ##"<< nl; + const symmTensor2D stiff + ( + pow(10.0, 10), pow(10.0, -8), + pow(10.0, 9) + ); + test_eigen_funcs(stiff); + } + + if (nFail_) { Info<< nl << " #### " diff --git a/applications/test/Tensor/Test-Tensor.C b/applications/test/Tensor/Test-Tensor.C index 9bc1066726..18f8d2b51c 100644 --- a/applications/test/Tensor/Test-Tensor.C +++ b/applications/test/Tensor/Test-Tensor.C @@ -902,6 +902,104 @@ void test_global_opers(Type) } +// Return false if given eigenvalues fail to satisy eigenvalue relations +// Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307) +void test_eigenvalues +( + const tensor& T, + const Vector& EVals, + const bool prod = true +) +{ + if (prod) + { + const scalar determinant = det(T); + // In case of complex EVals, the production is effectively scalar + // due to the (complex*complex conjugate) results in zero imag part + const scalar EValsProd = ((EVals.x()*EVals.y()*EVals.z()).real()); + cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant); + } + + { + const scalar trace = tr(T); + scalar EValsSum = 0.0; + // In case of complex EVals, the summation is effectively scalar + // due to the (complex+complex conjugate) results in zero imag part + for (const auto& val : EVals) + { + EValsSum += val.real(); + } + cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace); + } +} + + +// Return false if a given eigenvalue-eigenvector pair +// fails to satisfy the characteristic equation +void test_characteristic_equation +( + const tensor& T, + const Vector& EVals, + const Tensor& EVecs +) +{ + Info<< "# Characteristic equation:" << nl; + + Tensor Tc(Zero); + forAll(T, i) + { + Tc[i] = complex(T[i], 0); + } + + for (direction dir = 0; dir < pTraits::nComponents; ++dir) + { + const Vector leftSide(Tc & EVecs.row(dir)); + const Vector rightSide(EVals[dir]*EVecs.row(dir)); + const Vector X(leftSide - rightSide); + + for (const auto x : X) + { + cmp(" (sT & EVec - EVal*EVec) = 0:", mag(x), 0.0, 1e-8, 1e-6); + } + } +} + + +// Return false if the eigen functions fail to satisfy relations +void test_eigen_funcs(const tensor& T, const bool prod = true) +{ + Info<< "# Operand:" << nl + << " tensor = " << T << nl; + + + Info<< "# Return eigenvalues of a given tensor:" << nl; + const Vector EVals(eigenValues(T)); + Info<< EVals << endl; + test_eigenvalues(T, EVals, prod); + + Info<< "# Return an eigenvector of a given tensor in a given direction" + << " corresponding to a given eigenvalue:" << nl; + const Vector standardBasis1(Zero, pTraits::one, Zero); + const Vector standardBasis2(Zero, Zero, pTraits::one); + const Vector EVec + ( + eigenVector(T, EVals.x(), standardBasis1, standardBasis2) + ); + Info<< EVec << endl; + + Info<< "# Return eigenvectors of a given tensor corresponding to" + << " given eigenvalues:" << nl; + const Tensor EVecs0(eigenVectors(T, EVals)); + Info<< EVecs0 << endl; + test_characteristic_equation(T, EVals, EVecs0); + + Info<< "# Return eigenvectors of a given tensor by computing" + << " the eigenvalues of the tensor in the background:" << nl; + const Tensor EVecs1(eigenVectors(T)); + Info<< EVecs1 << endl; +} + + // Do compile-time recursion over the given types template inline typename std::enable_if::type @@ -947,6 +1045,53 @@ int main() run_tests(types, typeID); + Info<< nl << " ## Test tensor eigen functions: ##" << nl; + const label numberOfTests = 10000; + Random rndGen(1234); + + for (label i = 0; i < numberOfTests; ++i) + { + const tensor T(makeRandomContainer(rndGen)); + test_eigen_funcs(T); + } + + { + Info<< nl << " ## Test eigen functions by a zero tensor: ##"<< nl; + const tensor zeroT(Zero); + test_eigen_funcs(zeroT); + } + { + Info<< nl + << " ## Test eigen functions by a skew-symmetric tensor" + << " consisting of no-real eigenvalues: ##" + << nl; + const tensor T + ( + 0, 1, 1, + -1, 0, 1, + -1, -1, 0 + ); + test_eigen_funcs(T); + } + { + Info<< nl + << " ## Test eigen functions by a stiff tensor: ##" + << nl; + const tensor stiff + ( + pow(10.0, 10), pow(10.0, 8), pow(10.0, 7), + pow(10.0, -8), pow(10.0, 9), pow(10.0, -8), + pow(10.0, 10), pow(10.0, 8), pow(10.0, 7) + ); + // Although eigendecomposition is successful for the stiff tensors, + // cross-check between prod(eigenvalues) ?= det(stiff) is inherently + // problematic; therefore, eigenvalues of the stiff tensors are + // cross-checked by only sum(eigenvalues) ?= trace(stiff) + const bool testProd = false; + test_eigen_funcs(stiff, testProd); + } + + if (nFail_) { Info<< nl << " #### " diff --git a/applications/test/Tensor2D/Test-Tensor2D.C b/applications/test/Tensor2D/Test-Tensor2D.C index 4b8950e611..b504e502fe 100644 --- a/applications/test/Tensor2D/Test-Tensor2D.C +++ b/applications/test/Tensor2D/Test-Tensor2D.C @@ -716,6 +716,94 @@ void test_global_opers(Type) } +// Return false if given eigenvalues fail to satisy eigenvalue relations +// Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307) +void test_eigenvalues(const tensor2D& T, const Vector2D& EVals) +{ + { + const scalar determinant = det(T); + // In case of complex EVals, the production is effectively scalar + // due to the (complex*complex conjugate) results in zero imag part + const scalar EValsProd = ((EVals.x()*EVals.y()).real()); + cmp("# Product of eigenvalues = det(sT):", EValsProd, determinant); + } + + { + const scalar trace = tr(T); + scalar EValsSum = 0.0; + // In case of complex EVals, the summation is effectively scalar + // due to the (complex+complex conjugate) results in zero imag part + for (const auto& val : EVals) + { + EValsSum += val.real(); + } + cmp("# Sum of eigenvalues = trace(sT):", EValsSum, trace, 1e-8, 1e-8); + } +} + + +// Return false if a given eigenvalue-eigenvector pair +// fails to satisfy the characteristic equation +void test_characteristic_equation +( + const tensor2D& T, + const Vector2D& EVals, + const Tensor2D& EVecs +) +{ + Info<< "# Characteristic equation:" << nl; + + Tensor2D Tc(Zero); + forAll(T, i) + { + Tc[i] = complex(T[i], 0); + } + + for (direction dir = 0; dir < pTraits::nComponents; ++dir) + { + const Vector2D leftSide(Tc & EVecs.row(dir)); + const Vector2D rightSide(EVals[dir]*EVecs.row(dir)); + const Vector2D X(leftSide - rightSide); + + for (const auto x : X) + { + cmp(" (Tc & EVec - EVal*EVec) = 0:", mag(x), 0.0, 1e-8, 1e-6); + } + } +} + + +// Return false if the eigen functions fail to satisfy relations +void test_eigen_funcs(const tensor2D& T) +{ + Info<< "# Operand:" << nl + << " tensor2D = " << T << nl; + + + Info<< "# Return eigenvalues of a given tensor2D:" << nl; + const Vector2D EVals(eigenValues(T)); + Info<< EVals << endl; + test_eigenvalues(T, EVals); + + Info<< "# Return an eigenvector of a given tensor2D in a given direction" + << " corresponding to a given eigenvalue:" << nl; + const Vector2D standardBasis(pTraits::one, Zero); + const Vector2D EVec(eigenVector(T, EVals.x(), standardBasis)); + Info<< EVec << endl; + + Info<< "# Return eigenvectors of a given tensor2D corresponding to" + << " given eigenvalues:" << nl; + const Tensor2D EVecs0(eigenVectors(T, EVals)); + Info<< EVecs0 << endl; + test_characteristic_equation(T, EVals, EVecs0); + + Info<< "# Return eigenvectors of a given tensor2D by computing" + << " the eigenvalues of the tensor2D in the background:" << nl; + const Tensor2D EVecs1(eigenVectors(T)); + Info<< EVecs1 << endl; +} + + // Do compile-time recursion over the given types template inline typename std::enable_if::type @@ -762,6 +850,58 @@ int main(int argc, char *argv[]) run_tests(types, typeID); + Info<< nl << " ## Test tensor2D eigen functions: ##" << nl; + const label numberOfTests = 10000; + Random rndGen(1234); + + for (label i = 0; i < numberOfTests; ++i) + { + const tensor2D T(makeRandomContainer(rndGen)); + test_eigen_funcs(T); + } + + { + Info<< nl << " ## Test eigen functions by a zero tensor2D: ##"<< nl; + const tensor2D zeroT(Zero); + test_eigen_funcs(zeroT); + } + { + Info<< nl + << " ## Test eigen functions by a tensor2D consisting of" + << " repeated eigenvalues: ##" + << nl; + const tensor2D T + ( + -1.0, 2.0, + 0.0, -1.0 + ); + test_eigen_funcs(T); + } + { + Info<< nl + << " ## Test eigen functions by a skew-symmetric tensor2D" + << " consisting of no-real eigenvalues: ##" + << nl; + const tensor2D T + ( + 0.0, 1.0, + -1.0, 0.0 + ); + test_eigen_funcs(T); + } + { + Info<< nl + << " ## Test eigen functions by a stiff tensor2D: ##" + << nl; + const tensor2D stiff + ( + pow(10.0, 10), pow(10.0, 8), + pow(10.0, -8), pow(10.0, 9) + ); + test_eigen_funcs(stiff); + } + + { Info<< "# Pre-v2006 tests:" << nl; vector2D v1(1, 2), v2(3, 4); @@ -823,7 +963,6 @@ int main(int argc, char *argv[]) } } - if (nFail_) { Info<< nl << " #### " diff --git a/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.C b/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.C index 3d4fe827af..3be61102ab 100644 --- a/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.C +++ b/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.C @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2013 OpenFOAM Foundation + Copyright (C) 2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -145,28 +146,6 @@ dimensionedTensor skew(const dimensionedTensor& dt) } -dimensionedVector eigenValues(const dimensionedTensor& dt) -{ - return dimensionedVector - ( - "eigenValues("+dt.name()+')', - dt.dimensions(), - eigenValues(dt.value()) - ); -} - - -dimensionedTensor eigenVectors(const dimensionedTensor& dt) -{ - return dimensionedTensor - ( - "eigenVectors("+dt.name()+')', - dimless, - eigenVectors(dt.value()) - ); -} - - dimensionedVector eigenValues(const dimensionedSymmTensor& dt) { return dimensionedVector diff --git a/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.H b/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.H index ec9341de60..ad51af2fd0 100644 --- a/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.H +++ b/src/OpenFOAM/dimensionedTypes/dimensionedTensor/dimensionedTensor.H @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation + Copyright (C) 2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -63,9 +64,6 @@ dimensionedSymmTensor symm(const dimensionedTensor&); dimensionedSymmTensor twoSymm(const dimensionedTensor&); dimensionedTensor skew(const dimensionedTensor&); -dimensionedVector eigenValues(const dimensionedTensor&); -dimensionedTensor eigenVectors(const dimensionedTensor&); - dimensionedVector eigenValues(const dimensionedSymmTensor&); dimensionedTensor eigenVectors(const dimensionedSymmTensor&); diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedTensorField/DimensionedTensorField.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedTensorField/DimensionedTensorField.C index f50a4ee0df..612916f800 100644 --- a/src/OpenFOAM/fields/DimensionedFields/DimensionedTensorField/DimensionedTensorField.C +++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedTensorField/DimensionedTensorField.C @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2013 OpenFOAM Foundation + Copyright (C) 2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -48,8 +49,6 @@ UNARY_FUNCTION(tensor, tensor, dev2, transform) UNARY_FUNCTION(scalar, tensor, det, pow3) UNARY_FUNCTION(tensor, tensor, cof, pow2) UNARY_FUNCTION(tensor, tensor, inv, inv) -UNARY_FUNCTION(vector, tensor, eigenValues, transform) -UNARY_FUNCTION(tensor, tensor, eigenVectors, sign) UNARY_FUNCTION(vector, symmTensor, eigenValues, transform) UNARY_FUNCTION(symmTensor, symmTensor, eigenVectors, sign) diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedTensorField/DimensionedTensorField.H b/src/OpenFOAM/fields/DimensionedFields/DimensionedTensorField/DimensionedTensorField.H index 28b1e965a6..e05f243543 100644 --- a/src/OpenFOAM/fields/DimensionedFields/DimensionedTensorField/DimensionedTensorField.H +++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedTensorField/DimensionedTensorField.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. @@ -61,8 +62,6 @@ UNARY_FUNCTION(tensor, tensor, dev2, transform) UNARY_FUNCTION(scalar, tensor, det, transform) UNARY_FUNCTION(tensor, tensor, cof, cof) UNARY_FUNCTION(tensor, tensor, inv, inv) -UNARY_FUNCTION(vector, tensor, eigenValues, transform) -UNARY_FUNCTION(tensor, tensor, eigenVectors, sign) UNARY_FUNCTION(vector, symmTensor, eigenValues, transform) UNARY_FUNCTION(symmTensor, symmTensor, eigenVectors, sign) diff --git a/src/OpenFOAM/fields/FieldFields/tensorFieldField/tensorFieldField.C b/src/OpenFOAM/fields/FieldFields/tensorFieldField/tensorFieldField.C index 62747be172..b0e46f32b0 100644 --- a/src/OpenFOAM/fields/FieldFields/tensorFieldField/tensorFieldField.C +++ b/src/OpenFOAM/fields/FieldFields/tensorFieldField/tensorFieldField.C @@ -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. @@ -217,8 +217,6 @@ UNARY_FUNCTION(tensor, tensor, dev2) UNARY_FUNCTION(scalar, tensor, det) UNARY_FUNCTION(tensor, tensor, cof) UNARY_FUNCTION(tensor, tensor, inv) -UNARY_FUNCTION(vector, tensor, eigenValues) -UNARY_FUNCTION(tensor, tensor, eigenVectors) UNARY_FUNCTION(vector, symmTensor, eigenValues) UNARY_FUNCTION(symmTensor, symmTensor, eigenVectors) diff --git a/src/OpenFOAM/fields/FieldFields/tensorFieldField/tensorFieldField.H b/src/OpenFOAM/fields/FieldFields/tensorFieldField/tensorFieldField.H index af817dddd1..34e9966545 100644 --- a/src/OpenFOAM/fields/FieldFields/tensorFieldField/tensorFieldField.H +++ b/src/OpenFOAM/fields/FieldFields/tensorFieldField/tensorFieldField.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. @@ -166,8 +166,6 @@ UNARY_FUNCTION(tensor, tensor, dev2) UNARY_FUNCTION(scalar, tensor, det) UNARY_FUNCTION(tensor, tensor, cof) UNARY_FUNCTION(tensor, tensor, inv) -UNARY_FUNCTION(vector, tensor, eigenValues) -UNARY_FUNCTION(tensor, tensor, eigenVectors) UNARY_FUNCTION(vector, symmTensor, eigenValues) UNARY_FUNCTION(symmTensor, symmTensor, eigenVectors) diff --git a/src/OpenFOAM/fields/Fields/tensorField/tensorField.C b/src/OpenFOAM/fields/Fields/tensorField/tensorField.C index 10e0112a17..37df235b32 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 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -122,9 +122,6 @@ tmp inv(const tmp& tf) return tres; } -UNARY_FUNCTION(vector, tensor, eigenValues) -UNARY_FUNCTION(tensor, tensor, eigenVectors) - UNARY_FUNCTION(vector, symmTensor, eigenValues) UNARY_FUNCTION(tensor, symmTensor, eigenVectors) diff --git a/src/OpenFOAM/fields/Fields/tensorField/tensorField.H b/src/OpenFOAM/fields/Fields/tensorField/tensorField.H index 8ad14710f8..386d8ce648 100644 --- a/src/OpenFOAM/fields/Fields/tensorField/tensorField.H +++ b/src/OpenFOAM/fields/Fields/tensorField/tensorField.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. @@ -190,8 +190,6 @@ UNARY_FUNCTION(tensor, tensor, dev2) UNARY_FUNCTION(scalar, tensor, det) UNARY_FUNCTION(tensor, tensor, cof) UNARY_FUNCTION(tensor, tensor, inv) -UNARY_FUNCTION(vector, tensor, eigenValues) -UNARY_FUNCTION(tensor, tensor, eigenVectors) UNARY_FUNCTION(vector, symmTensor, eigenValues) UNARY_FUNCTION(tensor, symmTensor, eigenVectors) diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricTensorField/GeometricTensorField.C b/src/OpenFOAM/fields/GeometricFields/GeometricTensorField/GeometricTensorField.C index 8c1e0eada8..5e470bbb61 100644 --- a/src/OpenFOAM/fields/GeometricFields/GeometricTensorField/GeometricTensorField.C +++ b/src/OpenFOAM/fields/GeometricFields/GeometricTensorField/GeometricTensorField.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2013 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -267,8 +267,6 @@ UNARY_FUNCTION(tensor, tensor, dev2, transform) UNARY_FUNCTION(scalar, tensor, det, pow3) UNARY_FUNCTION(tensor, tensor, cof, pow2) UNARY_FUNCTION(tensor, tensor, inv, inv) -UNARY_FUNCTION(vector, tensor, eigenValues, transform) -UNARY_FUNCTION(tensor, tensor, eigenVectors, sign) UNARY_FUNCTION(vector, symmTensor, eigenValues, transform) UNARY_FUNCTION(symmTensor, symmTensor, eigenVectors, sign) diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricTensorField/GeometricTensorField.H b/src/OpenFOAM/fields/GeometricFields/GeometricTensorField/GeometricTensorField.H index f9808b90ca..45626d2aea 100644 --- a/src/OpenFOAM/fields/GeometricFields/GeometricTensorField/GeometricTensorField.H +++ b/src/OpenFOAM/fields/GeometricFields/GeometricTensorField/GeometricTensorField.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. @@ -167,8 +167,6 @@ UNARY_FUNCTION(tensor, tensor, dev2, transform) UNARY_FUNCTION(scalar, tensor, det, transform) UNARY_FUNCTION(tensor, tensor, cof, cof) UNARY_FUNCTION(tensor, tensor, inv, inv) -UNARY_FUNCTION(vector, tensor, eigenValues, transform) -UNARY_FUNCTION(tensor, tensor, eigenVectors, sign) UNARY_FUNCTION(vector, symmTensor, eigenValues, transform) UNARY_FUNCTION(symmTensor, symmTensor, eigenVectors, sign) diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldFunctions.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldFunctions.H index 8818963dbc..6624ba47a4 100644 --- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldFunctions.H +++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchFieldFunctions.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. @@ -346,20 +346,6 @@ inline void skew ) {} -inline void eigenValues -( - pointPatchField&, - const pointPatchField& -) -{} - -inline void eigenVectors -( - pointPatchField&, - const pointPatchField& -) -{} - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C index 30582a25ea..b7a33e05e8 100644 --- a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C +++ b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation + Copyright (C) 2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -26,6 +27,11 @@ License \*---------------------------------------------------------------------------*/ #include "symmTensor.H" +#include "cubicEqn.H" +#include "mathematicalConstants.H" + +using namespace Foam::constant::mathematical; + // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -85,4 +91,248 @@ const Foam::symmTensor Foam::symmTensor::I ); +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::vector Foam::eigenValues(const symmTensor& t) +{ + // Coefficients of the characteristic cubic polynomial (a = 1) + const scalar b = + - t.xx() - t.yy() - t.zz(); + const scalar c = + t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz() + - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz(); + const scalar d = + - t.xx()*t.yy()*t.zz() + - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx() + + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx(); + + // Determine the roots of the characteristic cubic polynomial + Roots<3> roots(cubicEqn(1, b, c, d).roots()); + + vector eVals(Zero); + + // Check the root types + forAll(roots, i) + { + switch (roots.type(i)) + { + case roots::real: + eVals[i] = roots[i]; + break; + case roots::complex: + case roots::posInf: + case roots::negInf: + case roots::nan: + FatalErrorInFunction + << "Eigenvalue computation fails for symmTensor: " << t + << "due to the non-real root = " << roots[i] + << exit(FatalError); + } + } + + // Sort the eigenvalues into ascending order + if (eVals.x() > eVals.y()) + { + Swap(eVals.x(), eVals.y()); + } + if (eVals.y() > eVals.z()) + { + Swap(eVals.y(), eVals.z()); + } + if (eVals.x() > eVals.y()) + { + Swap(eVals.x(), eVals.y()); + } + + return eVals; +} + + +Foam::vector Foam::eigenVector +( + const symmTensor& T, + const scalar eVal, + const vector& standardBasis1, + const vector& standardBasis2 +) +{ + // Construct the characteristic equation system for this eigenvalue + const tensor A(T - eVal*I); + + // Determinants of the 2x2 sub-matrices used to find the eigenvectors + // Sub-determinants for a unique eigenvenvalue + scalar sd0 = A.yy()*A.zz() - A.yz()*A.zy(); + scalar sd1 = A.zz()*A.xx() - A.zx()*A.xz(); + scalar sd2 = A.xx()*A.yy() - A.xy()*A.yx(); + scalar magSd0 = mag(sd0); + scalar magSd1 = mag(sd1); + scalar magSd2 = mag(sd2); + + // Evaluate the eigenvector using the largest sub-determinant + if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL) + { + const vector eVec + ( + 1, + (A.yz()*A.zx() - A.zz()*A.yx())/sd0, + (A.zy()*A.yx() - A.yy()*A.zx())/sd0 + ); + + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); + } + else if (magSd1 >= magSd2 && magSd1 > SMALL) + { + const vector eVec + ( + (A.xz()*A.zy() - A.zz()*A.xy())/sd1, + 1, + (A.zx()*A.xy() - A.xx()*A.zy())/sd1 + ); + + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); + } + else if (magSd2 > SMALL) + { + const vector eVec + ( + (A.xy()*A.yz() - A.yy()*A.xz())/sd2, + (A.yx()*A.xz() - A.xx()*A.yz())/sd2, + 1 + ); + + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); + } + + // Sub-determinants for a repeated eigenvalue + sd0 = A.yy()*standardBasis1.z() - A.yz()*standardBasis1.y(); + sd1 = A.zz()*standardBasis1.x() - A.zx()*standardBasis1.z(); + sd2 = A.xx()*standardBasis1.y() - A.xy()*standardBasis1.x(); + magSd0 = mag(sd0); + magSd1 = mag(sd1); + magSd2 = mag(sd2); + + // Evaluate the eigenvector using the largest sub-determinant + if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL) + { + const vector eVec + ( + 1, + (A.yz()*standardBasis1.x() - standardBasis1.z()*A.yx())/sd0, + (standardBasis1.y()*A.yx() - A.yy()*standardBasis1.x())/sd0 + ); + + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); + } + else if (magSd1 >= magSd2 && magSd1 > SMALL) + { + const vector eVec + ( + (standardBasis1.z()*A.zy() - A.zz()*standardBasis1.y())/sd1, + 1, + (A.zx()*standardBasis1.y() - standardBasis1.x()*A.zy())/sd1 + ); + + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); + } + else if (magSd2 > SMALL) + { + const vector eVec + ( + (A.xy()*standardBasis1.z() - standardBasis1.y()*A.xz())/sd2, + (standardBasis1.x()*A.xz() - A.xx()*standardBasis1.z())/sd2, + 1 + ); + + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); + } + + // Triple eigenvalue + return standardBasis1^standardBasis2; +} + + +Foam::tensor Foam::eigenVectors +( + const symmTensor& T, + const vector& eVals +) +{ + vector Ux(1, 0, 0), Uy(0, 1, 0), Uz(0, 0, 1); + + Ux = eigenVector(T, eVals.x(), Uy, Uz); + Uy = eigenVector(T, eVals.y(), Uz, Ux); + Uz = eigenVector(T, eVals.z(), Ux, Uy); + + return tensor(Ux, Uy, Uz); +} + + +Foam::tensor Foam::eigenVectors(const symmTensor& T) +{ + const vector eVals(eigenValues(T)); + + return eigenVectors(T, eVals); +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.H b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.H index b37e59a566..9625b332ca 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 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -28,7 +28,13 @@ Typedef Foam::symmTensor Description - SymmTensor of scalars. + SymmTensor of scalars, i.e. SymmTensor. + + Analytical functions for the computation of real eigenvalues and + real eigenvectors from a given symmTensor. + +See also + Test-SymmTensor.C SourceFiles symmTensor.C @@ -39,6 +45,9 @@ SourceFiles #define symmTensor_H #include "SymmTensor.H" +#include "vector.H" +#include "sphericalTensor.H" +#include "tensor.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -49,6 +58,55 @@ namespace Foam typedef SymmTensor symmTensor; +typedef Tensor tensor; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//- Return real ascending-order eigenvalues of a given symmTensor +// \param T symmTensor +// +// \return vector real eigenvalues +vector eigenValues(const symmTensor& T); + + +//- Return a real eigenvector corresponding to +//- a given real eigenvalue of a given symmTensor +// \param T symmTensor +// \param eVal real eigenvalue +// \param standardBasis1 symmTensor orthogonal component 1 +// \param standardBasis2 symmTensor orthogonal component 2 +// +// \return vector real eigenvector +vector eigenVector +( + const symmTensor& T, + const scalar eVal, + const vector& standardBasis1, + const vector& standardBasis2 +); + + +//- Return real eigenvectors corresponding to +//- given real eigenvalues of a given symmTensor +// \param T symmTensor +// \param eVals real eigenvalues +// +// \return tensor real eigenvectors, each row is an eigenvector +tensor eigenVectors +( + const symmTensor& T, + const vector& eVals +); + + +//- Return real eigenvectors of a given symmTensor by computing +//- the real eigenvalues of the tensor in the background +// \param T symmTensor +// +// \return tensor real eigenvectors, each row is an eigenvector +tensor eigenVectors(const symmTensor& T); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.C b/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.C index fc993e2386..7ff6ade25e 100644 --- a/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.C +++ b/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.C @@ -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,4 +84,101 @@ const Foam::symmTensor2D Foam::symmTensor2D::I ); +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::vector2D Foam::eigenValues(const symmTensor2D& T) +{ + //(K:Eqs. 3.2-3.3) + const scalar skewTrace = T.xx() - T.yy(); + const scalar trace = tr(T); + const scalar gap = sign(skewTrace)*hypot(skewTrace, 2.0*T.xy()); + + return vector2D (0.5*(trace + gap), 0.5*(trace - gap)); +} + + +Foam::vector2D Foam::eigenVector +( + const symmTensor2D& T, + const scalar eVal, + const vector2D& standardBasis +) +{ + // Construct the characteristic equation system for this eigenvalue + const tensor2D A(T - eVal*tensor2D::I); + + // Evaluate the eigenvector using the largest divisor + if (mag(A.yy()) > mag(A.xx()) && mag(A.yy()) > SMALL) + { + const vector2D eVec(1, -A.yx()/A.yy()); + + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); + } + else if (mag(A.xx()) > SMALL) + { + const vector2D eVec(-A.xy()/A.xx(), 1); + + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); + } + + // Repeated eigenvalue + return vector2D(-standardBasis.y(), standardBasis.x()); +} + + +Foam::tensor2D Foam::eigenVectors +( + const symmTensor2D& T, + const vector2D& eVals +) +{ + // (K:Eq. 3.5) + const scalar skewTrace = T.xx() - T.yy(); + + if (mag(skewTrace) > SMALL) + { + const scalar phi = 0.5*atan(2.0*T.xy()/skewTrace); + const scalar cphi = cos(phi); + const scalar sphi = sin(phi); + return tensor2D(cphi, sphi, -sphi, cphi); + } + else if (mag(T.xy()) > SMALL) + { + const scalar a = 0.70710678; // phi ~ 45deg + return tensor2D(a, sign(T.xy())*a, -1*sign(T.xy())*a, a); + } + + // (K:p. 3) + return tensor2D(1, 0, 0, 1); +} + + +Foam::tensor2D Foam::eigenVectors(const symmTensor2D& T) +{ + const vector2D eVals(eigenValues(T)); + + return eigenVectors(T, eVals); +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.H b/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.H index a772673b6b..30d0fc935b 100644 --- a/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.H +++ b/src/OpenFOAM/primitives/SymmTensor2D/symmTensor2D/symmTensor2D.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2013 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -28,7 +28,22 @@ Typedef Foam::symmTensor2D Description - SymmTensor2D of scalars. + SymmTensor2D of scalars, i.e. SymmTensor2D. + + Analytical functions for the computation of real eigenvalues + and real eigenvectors from a given symmTensor2D. + + Reference: + \verbatim + 2-by-2 eigenvalue algorithm (tag:K): + Kronenburg, M. J. (2015). + A method for fast diagonalization of + a 2x2 or 3x3 real symmetric matrix. + arXiv:1306.6291. + \endverbatim + +See also + Test-SymmTensor2D.C SourceFiles symmTensor2D.C @@ -39,6 +54,7 @@ SourceFiles #define symmTensor2D_H #include "SymmTensor2D.H" +#include "vector2D.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -49,8 +65,49 @@ namespace Foam typedef SymmTensor2D symmTensor2D; +typedef Tensor2D tensor2D; + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +//- Return real arbitrary-order eigenvalues of a given symmTensor2D +// \param T symmTensor2D +// +// \return vector2D real eigenvalues +vector2D eigenValues(const symmTensor2D& T); + + +//- Return a real eigenvector corresponding to +//- a given real eigenvalue of a given symmTensor2D +// \param T symmTensor2D +// \param eVal real eigenvalue +// \param standardBasis symmTensor2D orthogonal component, e.g. vector2D(1, 0) +// +// \return vector2D real eigenvector +vector2D eigenVector +( + const symmTensor2D& T, + const scalar eVal, + const vector2D& standardBasis +); + + +//- Return real eigenvectors corresponding to +//- given real eigenvalues of a given symmTensor2D +// \param T symmTensor2D +// \param eVals real eigenvalues +// +// \return tensor2D real eigenvectors, each row is an eigenvector +tensor2D eigenVectors(const symmTensor2D& T, const vector2D& eVals); + + +//- Return real eigenvectors of a given symmTensor2D by computing +//- the real eigenvalues of the symmTensor2D in the background +// \param T symmTensor2D +// +// \return tensor2D real eigenvectors, each row is an eigenvector +tensor2D eigenVectors(const symmTensor2D& T); + + } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C index 4aef370e09..69471de43f 100644 --- a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C +++ b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation + Copyright (C) 2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -73,7 +74,7 @@ const Foam::tensor Foam::tensor::I // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -Foam::vector Foam::eigenValues(const tensor& t) +Foam::Vector Foam::eigenValues(const tensor& t) { // Coefficients of the characteristic cubic polynomial (a = 1) const scalar b = @@ -86,117 +87,147 @@ Foam::vector Foam::eigenValues(const tensor& t) - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx() + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx(); - // Solve - Roots<3> roots = cubicEqn(1, b, c, d).roots(); + // Determine the roots of the characteristic cubic polynomial + const Roots<3> roots(cubicEqn(1, b, c, d).roots()); // Check the root types - vector lambda = vector::zero; + bool isComplex = false; forAll(roots, i) { switch (roots.type(i)) { - case roots::real: - lambda[i] = roots[i]; - break; case roots::complex: - WarningInFunction - << "Complex eigenvalues detected for tensor: " << t - << endl; - lambda[i] = 0; + isComplex = true; break; case roots::posInf: - lambda[i] = VGREAT; - break; case roots::negInf: - lambda[i] = - VGREAT; - break; case roots::nan: FatalErrorInFunction - << "Eigenvalue calculation failed for tensor: " << t + << "Eigenvalue computation fails for tensor: " << t + << "due to the not-a-number root = " << roots[i] << exit(FatalError); + case roots::real: + break; } } - // Sort the eigenvalues into ascending order - if (lambda.x() > lambda.y()) + if (isComplex) { - Swap(lambda.x(), lambda.y()); - } - if (lambda.y() > lambda.z()) - { - Swap(lambda.y(), lambda.z()); - } - if (lambda.x() > lambda.y()) - { - Swap(lambda.x(), lambda.y()); + return + Vector + ( + complex(roots[0], 0), + complex(roots[1], roots[2]), + complex(roots[1], -roots[2]) + ); } - return lambda; + return + Vector + ( + complex(roots[0], 0), + complex(roots[1], 0), + complex(roots[2], 0) + ); } -Foam::vector Foam::eigenVector +Foam::Vector Foam::eigenVector ( const tensor& T, - const scalar lambda, - const vector& direction1, - const vector& direction2 + const complex& eVal, + const Vector& standardBasis1, + const Vector& standardBasis2 ) { - // Construct the linear system for this eigenvalue - tensor A(T - lambda*I); + // Construct the characteristic equation system for this eigenvalue + Tensor A(Zero); + forAll(A, i) + { + A[i] = complex(T[i], 0); + } + A.xx() -= eVal; + A.yy() -= eVal; + A.zz() -= eVal; // Determinants of the 2x2 sub-matrices used to find the eigenvectors - scalar sd0, sd1, sd2; - scalar magSd0, magSd1, magSd2; - - // Sub-determinants for a unique eivenvalue - sd0 = A.yy()*A.zz() - A.yz()*A.zy(); - sd1 = A.zz()*A.xx() - A.zx()*A.xz(); - sd2 = A.xx()*A.yy() - A.xy()*A.yx(); - magSd0 = mag(sd0); - magSd1 = mag(sd1); - magSd2 = mag(sd2); + // Sub-determinants for a unique eigenvenvalue + complex sd0 = A.yy()*A.zz() - A.yz()*A.zy(); + complex sd1 = A.zz()*A.xx() - A.zx()*A.xz(); + complex sd2 = A.xx()*A.yy() - A.xy()*A.yx(); + scalar magSd0 = mag(sd0); + scalar magSd1 = mag(sd1); + scalar magSd2 = mag(sd2); // Evaluate the eigenvector using the largest sub-determinant if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL) { - vector ev + const Vector eVec ( - 1, + complex(1, 0), (A.yz()*A.zx() - A.zz()*A.yx())/sd0, (A.zy()*A.yx() - A.yy()*A.zx())/sd0 ); - return ev/mag(ev); + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); } else if (magSd1 >= magSd2 && magSd1 > SMALL) { - vector ev + const Vector eVec ( (A.xz()*A.zy() - A.zz()*A.xy())/sd1, - 1, + complex(1, 0), (A.zx()*A.xy() - A.xx()*A.zy())/sd1 ); - return ev/mag(ev); + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); } else if (magSd2 > SMALL) { - vector ev + const Vector eVec ( (A.xy()*A.yz() - A.yy()*A.xz())/sd2, (A.yx()*A.xz() - A.xx()*A.yz())/sd2, - 1 + complex(1, 0) ); - return ev/mag(ev); + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); } // Sub-determinants for a repeated eigenvalue - sd0 = A.yy()*direction1.z() - A.yz()*direction1.y(); - sd1 = A.zz()*direction1.x() - A.zx()*direction1.z(); - sd2 = A.xx()*direction1.y() - A.xy()*direction1.x(); + sd0 = A.yy()*standardBasis1.z() - A.yz()*standardBasis1.y(); + sd1 = A.zz()*standardBasis1.x() - A.zx()*standardBasis1.z(); + sd2 = A.xx()*standardBasis1.y() - A.xy()*standardBasis1.x(); magSd0 = mag(sd0); magSd1 = mag(sd1); magSd2 = mag(sd2); @@ -204,90 +235,96 @@ Foam::vector Foam::eigenVector // Evaluate the eigenvector using the largest sub-determinant if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL) { - vector ev + const Vector eVec ( - 1, - (A.yz()*direction1.x() - direction1.z()*A.yx())/sd0, - (direction1.y()*A.yx() - A.yy()*direction1.x())/sd0 + complex(1, 0), + (A.yz()*standardBasis1.x() - standardBasis1.z()*A.yx())/sd0, + (standardBasis1.y()*A.yx() - A.yy()*standardBasis1.x())/sd0 ); - return ev/mag(ev); + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); } else if (magSd1 >= magSd2 && magSd1 > SMALL) { - vector ev + const Vector eVec ( - (direction1.z()*A.zy() - A.zz()*direction1.y())/sd1, - 1, - (A.zx()*direction1.y() - direction1.x()*A.zy())/sd1 + (standardBasis1.z()*A.zy() - A.zz()*standardBasis1.y())/sd1, + complex(1, 0), + (A.zx()*standardBasis1.y() - standardBasis1.x()*A.zy())/sd1 ); - return ev/mag(ev); + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); } else if (magSd2 > SMALL) { - vector ev + const Vector eVec ( - (A.xy()*direction1.z() - direction1.y()*A.xz())/sd2, - (direction1.x()*A.xz() - A.xx()*direction1.z())/sd2, - 1 + (A.xy()*standardBasis1.z() - standardBasis1.y()*A.xz())/sd2, + (standardBasis1.x()*A.xz() - A.xx()*standardBasis1.z())/sd2, + complex(1, 0) ); - return ev/mag(ev); + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); } // Triple eigenvalue - return direction1^direction2; + return standardBasis1^standardBasis2; } -Foam::tensor Foam::eigenVectors(const tensor& T, const vector& lambdas) -{ - vector Ux(1, 0, 0), Uy(0, 1, 0), Uz(0, 0, 1); - - Ux = eigenVector(T, lambdas.x(), Uy, Uz); - Uy = eigenVector(T, lambdas.y(), Uz, Ux); - Uz = eigenVector(T, lambdas.z(), Ux, Uy); - - return tensor(Ux, Uy, Uz); -} - - -Foam::tensor Foam::eigenVectors(const tensor& T) -{ - const vector lambdas(eigenValues(T)); - - return eigenVectors(T, lambdas); -} - - -Foam::vector Foam::eigenValues(const symmTensor& T) -{ - return eigenValues(tensor(T)); -} - - -Foam::vector Foam::eigenVector +Foam::Tensor Foam::eigenVectors ( - const symmTensor& T, - const scalar lambda, - const vector& direction1, - const vector& direction2 + const tensor& T, + const Vector& eVals ) { - return eigenVector(tensor(T), lambda, direction1, direction2); + Vector Ux(complex(1, 0), Zero, Zero); + Vector Uy(Zero, complex(1, 0), Zero); + Vector Uz(Zero, Zero, complex(1, 0)); + + Ux = eigenVector(T, eVals.x(), Uy, Uz); + Uy = eigenVector(T, eVals.y(), Uz, Ux); + Uz = eigenVector(T, eVals.z(), Ux, Uy); + + return Tensor(Ux, Uy, Uz); } -Foam::tensor Foam::eigenVectors(const symmTensor& T, const vector& lambdas) +Foam::Tensor Foam::eigenVectors(const tensor& T) { - return eigenVectors(tensor(T), lambdas); -} + const Vector eVals(eigenValues(T)); - -Foam::tensor Foam::eigenVectors(const symmTensor& T) -{ - return eigenVectors(tensor(T)); + return eigenVectors(T, eVals); } diff --git a/src/OpenFOAM/primitives/Tensor/tensor/tensor.H b/src/OpenFOAM/primitives/Tensor/tensor/tensor.H index c729200cf4..95c7560a3e 100644 --- a/src/OpenFOAM/primitives/Tensor/tensor/tensor.H +++ b/src/OpenFOAM/primitives/Tensor/tensor/tensor.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2014 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -28,7 +28,13 @@ Typedef Foam::tensor Description - Tensor of scalars. + Tensor of scalars, i.e. Tensor. + + Analytical functions for the computation of complex eigenvalues and + complex eigenvectors from a given tensor. + +See also + Test-Tensor.C SourceFiles tensor.C @@ -42,6 +48,7 @@ SourceFiles #include "vector.H" #include "sphericalTensor.H" #include "symmTensor.H" +#include "complex.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -52,27 +59,51 @@ namespace Foam typedef Tensor tensor; -vector eigenValues(const tensor& T); -vector eigenVector +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//- Return complex eigenvalues of a given tensor +// \param T tensor +// +// \return Vector eigenvalues +Vector eigenValues(const tensor& T); + + +//- Return a complex eigenvector corresponding to +//- a given complex eigenvalue of a given tensor +// \param T tensor +// \param eVal complex eigenvalue +// \param standardBasis1 tensor orthogonal component 1 +// \param standardBasis2 tensor orthogonal component 2 +// +// \return Vector eigenvector +Vector eigenVector ( const tensor& T, - const scalar lambda, - const vector& direction1, - const vector& direction2 + const complex& eVal, + const Vector& standardBasis1, + const Vector& standardBasis2 ); -tensor eigenVectors(const tensor& T, const vector& lambdas); -tensor eigenVectors(const tensor& T); -vector eigenValues(const symmTensor& T); -vector eigenVector + +//- Return complex eigenvectors corresponding to +//- given complex eigenvalues of a given tensor +// \param T tensor +// \param eVals complex eigenvalues +// +// \return Tensor eigenvectors, each row is an eigenvector +Tensor eigenVectors ( - const symmTensor& T, - const scalar lambda, - const vector& direction1, - const vector& direction2 + const tensor& T, + const Vector& eVals ); -tensor eigenVectors(const symmTensor& T, const vector& lambdas); -tensor eigenVectors(const symmTensor& T); + + +//- Return complex eigenvectors of a given tensor by computing +//- the complex eigenvalues of the tensor in the background +// \param T tensor +// +// \return Tensor complex eigenvectors, each row is an eigenvector +Tensor eigenVectors(const tensor& T); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C b/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C index 689837ca99..4cfca5d440 100644 --- a/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C +++ b/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation + Copyright (C) 2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -86,98 +87,133 @@ const Foam::tensor2D Foam::tensor2D::I // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -Foam::vector2D Foam::eigenValues(const tensor2D& t) +Foam::Vector2D Foam::eigenValues(const tensor2D& T) { - // Coefficients of the characteristic quadratic polynomial (a = 1) - const scalar b = - t.xx() - t.yy(); - const scalar c = t.xx()*t.yy() - t.xy()*t.yx(); + const scalar a = T.xx(); + const scalar b = T.xy(); + const scalar c = T.yx(); + const scalar d = T.yy(); - // Solve - Roots<2> roots = quadraticEqn(1, b, c).roots(); + const scalar trace = a + d; - // Check the root types - vector2D lambda = vector2D::zero; - forAll(roots, i) + // (JLM:p. 2246) + scalar w = b*c; + scalar e = std::fma(-b, c, w); + scalar f = std::fma(a, d, -w); + const scalar determinant = f + e; + + // Square-distance between two eigenvalues + const scalar gapSqr = std::fma(-4.0, determinant, sqr(trace)); + + // (F:Sec. 8.4.2.) + // Eigenvalues are effectively real + if (0 <= gapSqr) { - switch (roots.type(i)) + scalar firstRoot = 0.5*(trace - sign(-trace)*Foam::sqrt(gapSqr)); + + if (mag(firstRoot) < SMALL) { - case roots::real: - lambda[i] = roots[i]; - break; - case roots::complex: - WarningInFunction - << "Complex eigenvalues detected for tensor: " << t - << endl; - lambda[i] = 0; - break; - case roots::posInf: - lambda[i] = VGREAT; - break; - case roots::negInf: - lambda[i] = - VGREAT; - break; - case roots::nan: - FatalErrorInFunction - << "Eigenvalue calculation failed for tensor: " << t - << exit(FatalError); + WarningInFunction + << "Zero-valued root is found. Adding SMALL to the root " + << "to avoid floating-point exception." << nl; + firstRoot = SMALL; } - } - // Sort the eigenvalues into ascending order - if (lambda.x() > lambda.y()) + Vector2D eVals + ( + complex(firstRoot, 0), + complex(determinant/firstRoot, 0) + ); + + // Sort the eigenvalues into ascending order + if (mag(eVals.x()) > mag(eVals.y())) + { + Swap(eVals.x(), eVals.y()); + } + + return eVals; + } + // Eigenvalues are complex + else { - Swap(lambda.x(), lambda.y()); - } + const complex eVal(0.5*trace, 0.5*Foam::sqrt(mag(gapSqr))); - return lambda; + return Vector2D + ( + eVal, + eVal.conjugate() + ); + } } -Foam::vector2D Foam::eigenVector +Foam::Vector2D Foam::eigenVector ( const tensor2D& T, - const scalar lambda, - const vector2D& direction1 + const complex& eVal, + const Vector2D& standardBasis ) { - // Construct the linear system for this eigenvalue - tensor2D A(T - lambda*tensor2D::I); - + // (K:p. 47-48) // Evaluate the eigenvector using the largest divisor - if (mag(A.yy()) > mag(A.xx()) && mag(A.yy()) > SMALL) + if (mag(T.yx()) > mag(T.xy()) && mag(T.yx()) > VSMALL) { - vector2D ev(1, - A.yx()/A.yy()); + const Vector2D eVec(eVal - T.yy(), complex(T.yx())); - return ev/mag(ev); + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); } - else if (mag(A.xx()) > SMALL) + else if (mag(T.xy()) > VSMALL) { - vector2D ev(- A.xy()/A.xx(), 1); + const Vector2D eVec(complex(T.xy()), eVal - T.xx()); - return ev/mag(ev); + #ifdef FULLDEBUG + if (mag(eVec) < SMALL) + { + FatalErrorInFunction + << "Eigenvector magnitude should be non-zero:" + << "mag(eigenvector) = " << mag(eVec) + << abort(FatalError); + } + #endif + + return eVec/mag(eVec); } - // Repeated eigenvalue - return vector2D(- direction1.y(), direction1.x()); + return standardBasis; } -Foam::tensor2D Foam::eigenVectors(const tensor2D& T, const vector2D& lambdas) +Foam::Tensor2D Foam::eigenVectors +( + const tensor2D& T, + const Vector2D& eVals +) { - vector2D Ux(1, 0), Uy(0, 1); + Vector2D Ux(pTraits::one, Zero); + Vector2D Uy(Zero, pTraits::one); - Ux = eigenVector(T, lambdas.x(), Uy); - Uy = eigenVector(T, lambdas.y(), Ux); + Ux = eigenVector(T, eVals.x(), Uy); + Uy = eigenVector(T, eVals.y(), Ux); - return tensor2D(Ux, Uy); + return Tensor2D(Ux, Uy); } -Foam::tensor2D Foam::eigenVectors(const tensor2D& T) +Foam::Tensor2D Foam::eigenVectors(const tensor2D& T) { - const vector2D lambdas(eigenValues(T)); + const Vector2D eVals(eigenValues(T)); - return eigenVectors(T, lambdas); + return eigenVectors(T, eVals); } diff --git a/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.H b/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.H index 8b7edcef31..dcc7ac172a 100644 --- a/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.H +++ b/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -28,7 +28,39 @@ Typedef Foam::tensor2D Description - Tensor2D of scalars. + Tensor2D of scalars, i.e. Tensor2D. + + Analytical functions for the computation of complex eigenvalues and + complex eigenvectors from a given tensor2D. + + Reference: + \verbatim + 2-by-2 eigenvalue algorithm (tags:F; B): + Ford, W. (2014). + Numerical linear algebra with applications: Using MATLAB. + London: Elsevier/Academic Press. + DOI:10.1016/C2011-0-07533-6 + + Blinn, J. (1996). + Consider the lowly 2 x 2 matrix. + IEEE Computer Graphics and Applications, 16(2), 82-88. + DOI:10.1109/38.486688 + + 2-by-2 eigenvector algorithm (tag:K): + Knill, O. (2004). + Mathematics Math21b Fall 2004. + bit.ly/2kjPVlX (Retrieved:06-09-19) + + Kahan summation algorithm for 2-by-2 matrix determinants (tag:JLM): + Jeannerod, C.-P., Louvet, N., & Muller, J.-M., (2013). + Further analysis of Kahan's algorithm for the accurate computation + of 2x2 determinants. + Math. Comp. 82 (2013), 2245-2264. + DOI:10.1090/S0025-5718-2013-02679-8 + \endverbatim + +See also + Test-Tensor2D.C SourceFiles tensor2D.C @@ -39,7 +71,7 @@ SourceFiles #define tensor2D_H #include "Tensor2D.H" -#include "vector2D.H" +#include "complex.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -50,15 +82,49 @@ namespace Foam typedef Tensor2D tensor2D; -vector2D eigenValues(const tensor2D& t); -vector2D eigenVector +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//- Return complex eigenvalues of a given tensor2D +// \param T tensor2D +// +// \return Vector2D eigenvalues +Vector2D eigenValues(const tensor2D& T); + + +//- Return a complex eigenvector corresponding to +//- a given complex eigenvalue of a given tensor2D +// \param T tensor2D +// \param eVal complex eigenvalue +// \param standardBasis tensor2D orthogonal component, e.g. vector2D(1, 0) +// +// \return Vector2D eigenvector +Vector2D eigenVector ( - const tensor2D& t, - const scalar lambda, - const vector2D& direction1 + const tensor2D& T, + const complex& eVal, + const Vector2D& standardBasis ); -tensor2D eigenVectors(const tensor2D& t, const vector2D& lambdas); -tensor2D eigenVectors(const tensor2D& t); + + +//- Return complex eigenvectors corresponding to +//- given complex eigenvalues of a given tensor2D +// \param T tensor2D +// \param eVals eigenvalues +// +// \return Tensor2D eigenvectors, each row is an eigenvector +Tensor2D eigenVectors +( + const tensor2D& T, + const Vector2D& eVals +); + + +//- Return complex eigenvectors of a given tensor2D by computing +//- the complex eigenvalues of the tensor2D in the background +// \param T tensor2D +// +// \return Tensor2D eigenvectors, each row is an eigenvector +Tensor2D eigenVectors(const tensor2D& T); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C index f5ce6889e5..2e63f1f461 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -386,7 +386,7 @@ void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio { const pointField& pts = mesh_.points(); - tensor J = f.inertia(pts, fC); + symmTensor J(symm(f.inertia(pts, fC))); // Find the dominant collapse direction by finding the eigenvector // that corresponds to the normal direction, discarding it. The @@ -423,7 +423,7 @@ void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio } else { - vector eVals = eigenValues(J); + vector eVals(eigenValues(J)); if (mag(eVals.y() - eVals.x()) < 100*SMALL) { diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.C index 045069926e..a72e883cbb 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2015 OpenFOAM Foundation - Copyright (C) 2016 OpenCFD Ltd. + Copyright (C) 2016-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -145,7 +145,7 @@ Foam::eddy::eddy dir1_(0) { // Principal stresses - eigenvalues returned in ascending order - vector lambda = eigenValues(R); + vector lambda(eigenValues(R)); // Eddy rotation from principal-to-global axes // - given by the 3 eigenvectors of the Reynold stress tensor as rows in diff --git a/src/functionObjects/field/Lambda2/Lambda2.C b/src/functionObjects/field/Lambda2/Lambda2.C index de01ad8dde..c1aa6ec28f 100644 --- a/src/functionObjects/field/Lambda2/Lambda2.C +++ b/src/functionObjects/field/Lambda2/Lambda2.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2013-2016 OpenFOAM Foundation - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -58,10 +58,13 @@ bool Foam::functionObjects::Lambda2::calc() const tmp tgradU(fvc::grad(U)); const volTensorField& gradU = tgradU(); - const volTensorField SSplusWW + const volSymmTensorField SSplusWW ( - (symm(gradU) & symm(gradU)) - + (skew(gradU) & skew(gradU)) + symm + ( + (symm(gradU) & symm(gradU)) + + (skew(gradU) & skew(gradU)) + ) ); return store diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H index 919e775c9b..c168ea69af 100644 --- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H +++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H @@ -6,6 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation + Copyright (C) 2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -139,17 +140,17 @@ inline Foam::molecule::constantProperties::constantProperties // Calculate the inertia tensor in the current orientation - tensor momOfInertia(Zero); + symmTensor momOfInertia(Zero); forAll(siteReferencePositions_, i) { const vector& p(siteReferencePositions_[i]); - momOfInertia += siteMasses_[i]*tensor + momOfInertia += siteMasses_[i]*symmTensor ( p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(), - -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(), - -p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y() + p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(), + p.x()*p.x() + p.y()*p.y() ); } @@ -166,7 +167,7 @@ inline Foam::molecule::constantProperties::constantProperties momOfInertia /= eigenValues(momOfInertia).x(); - tensor e = eigenVectors(momOfInertia); + tensor e(eigenVectors(momOfInertia)); // Calculate the transformation between the principle axes to the // global axes @@ -203,11 +204,11 @@ inline Foam::molecule::constantProperties::constantProperties { const vector& p(siteReferencePositions_[i]); - momOfInertia += siteMasses_[i]*tensor + momOfInertia += siteMasses_[i]*symmTensor ( p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(), - -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(), - -p.z()*p.x(), -p.z()*p.y(), p.x()*p.x() + p.y()*p.y() + p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(), + p.x()*p.x() + p.y()*p.y() ); } diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C index a447be4101..e4f1d5d5d4 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2017 OpenCFD Ltd. + Copyright (C) 2017-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -309,7 +309,7 @@ Foam::triSurfaceTools::curvatures { pointFundamentalTensors[pI] /= (accumulatedWeights[pI] + SMALL); - vector2D principalCurvatures = eigenValues(pointFundamentalTensors[pI]); + vector2D principalCurvatures(eigenValues(pointFundamentalTensors[pI])); //scalar curvature = // (principalCurvatures[0] + principalCurvatures[1])/2;