From de6313e4ceeaabdbd092538d721a1bc40a0cc1a7 Mon Sep 17 00:00:00 2001
From: andy <andy@shelob.opencfd.co.uk>
Date: Mon, 22 Feb 2016 12:14:53 +0000
Subject: [PATCH] ENH: Migrated eigen vector/value updates from old development
line and updated dependencies
---
.../primitives/Tensor/tensor/tensor.C | 134 ++++++++++--------
.../primitives/Tensor/tensor/tensor.H | 30 ++--
.../polyTopoChange/edgeCollapser.C | 2 +-
3 files changed, 98 insertions(+), 68 deletions(-)
diff --git a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
index cadbef6a9b3..f84d2e53160 100644
--- a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
+++ b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
@@ -70,24 +70,30 @@ const Foam::tensor Foam::tensor::I
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-Foam::vector Foam::eigenValues(const tensor& t)
+Foam::vector Foam::eigenValues(const tensor& T)
{
// The eigenvalues
scalar i, ii, iii;
// diagonal matrix
- if
- (
+ const scalar onDiagMagSum =
(
- mag(t.xy()) + mag(t.xz()) + mag(t.yx())
- + mag(t.yz()) + mag(t.zx()) + mag(t.zy())
- )
- < SMALL
- )
+ mag(T.xx()) + mag(T.yy()) + mag(T.zz())
+ );
+
+ const scalar offDiagMagSum =
+ (
+ mag(T.xy()) + mag(T.xz()) + mag(T.yx())
+ + mag(T.yz()) + mag(T.zx()) + mag(T.zy())
+ );
+
+ const scalar magSum = onDiagMagSum + offDiagMagSum;
+
+ if (offDiagMagSum < max(VSMALL, SMALL*magSum))
{
- i = t.xx();
- ii = t.yy();
- iii = t.zz();
+ i = T.xx();
+ ii = T.yy();
+ iii = T.zz();
}
// non-diagonal matrix
@@ -96,16 +102,16 @@ Foam::vector Foam::eigenValues(const tensor& t)
// Coefficients of the characteristic polynmial
// x^3 + a*x^2 + b*x + c = 0
scalar a =
- - t.xx() - t.yy() - t.zz();
+ - T.xx() - T.yy() - T.zz();
scalar b =
- t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
- - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz();
+ T.xx()*T.yy() + T.xx()*T.zz() + T.yy()*T.zz()
+ - T.xy()*T.yx() - T.yz()*T.zy() - T.zx()*T.xz();
scalar c =
- - 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();
+ - 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();
// Auxillary variables
scalar aBy3 = a/3;
@@ -117,13 +123,13 @@ Foam::vector Foam::eigenValues(const tensor& t)
scalar QQ = Q*Q;
// Three identical roots
- if (mag(P) < SMALL && mag(Q) < SMALL)
+ if (mag(P) < SMALL*sqr(magSum) && mag(Q) < SMALL*pow3(magSum))
{
return vector(- aBy3, - aBy3, - aBy3);
}
// Two identical roots and one distinct root
- else if (mag(PPP/QQ - 1) < SMALL)
+ else if (mag(PPP - QQ) < SMALL*pow6(magSum))
{
scalar sqrtP = sqrt(P);
scalar signQ = sign(Q);
@@ -148,11 +154,11 @@ Foam::vector Foam::eigenValues(const tensor& t)
// based on the above logic, PPP must be less than QQ
else
{
- WarningInFunction
- << "complex eigenvalues detected for tensor: " << t
+ WarningIn("eigenValues(const tensor&)")
+ << "complex eigenvalues detected for tensor: " << T
<< endl;
- if (mag(P) < SMALL)
+ if (mag(P) < SMALL*sqr(magSum))
{
i = cbrt(QQ/2);
}
@@ -188,21 +194,14 @@ Foam::vector Foam::eigenValues(const tensor& t)
Foam::vector Foam::eigenVector
(
- const tensor& t,
- const scalar lambda
+ const tensor& T,
+ const scalar lambda,
+ const vector& direction1,
+ const vector& direction2
)
{
- // Constantly rotating direction ensures different eigenvectors are
- // generated when called sequentially with a multiple eigenvalue
- static vector direction(1,0,0);
- vector oldDirection(direction);
- scalar temp = direction[2];
- direction[2] = direction[1];
- direction[1] = direction[0];
- direction[0] = temp;
-
// Construct the linear system for this eigenvalue
- tensor A(t - lambda*I);
+ tensor A(T - lambda*I);
// Determinants of the 2x2 sub-matrices used to find the eigenvectors
scalar sd0, sd1, sd2;
@@ -252,9 +251,9 @@ Foam::vector Foam::eigenVector
}
// Sub-determinants for a repeated eigenvalue
- sd0 = A.yy()*direction.z() - A.yz()*direction.y();
- sd1 = A.zz()*direction.x() - A.zx()*direction.z();
- sd2 = A.xx()*direction.y() - A.xy()*direction.x();
+ 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();
magSd0 = mag(sd0);
magSd1 = mag(sd1);
magSd2 = mag(sd2);
@@ -265,8 +264,8 @@ Foam::vector Foam::eigenVector
vector ev
(
1,
- (A.yz()*direction.x() - direction.z()*A.yx())/sd0,
- (direction.y()*A.yx() - A.yy()*direction.x())/sd0
+ (A.yz()*direction1.x() - direction1.z()*A.yx())/sd0,
+ (direction1.y()*A.yx() - A.yy()*direction1.x())/sd0
);
return ev/mag(ev);
@@ -275,9 +274,9 @@ Foam::vector Foam::eigenVector
{
vector ev
(
- (direction.z()*A.zy() - A.zz()*direction.y())/sd1,
+ (direction1.z()*A.zy() - A.zz()*direction1.y())/sd1,
1,
- (A.zx()*direction.y() - direction.x()*A.zy())/sd1
+ (A.zx()*direction1.y() - direction1.x()*A.zy())/sd1
);
return ev/mag(ev);
@@ -286,8 +285,8 @@ Foam::vector Foam::eigenVector
{
vector ev
(
- (A.xy()*direction.z() - direction.y()*A.xz())/sd2,
- (direction.x()*A.xz() - A.xx()*direction.z())/sd2,
+ (A.xy()*direction1.z() - direction1.y()*A.xz())/sd2,
+ (direction1.x()*A.xz() - A.xx()*direction1.z())/sd2,
1
);
@@ -295,40 +294,57 @@ Foam::vector Foam::eigenVector
}
// Triple eigenvalue
- return oldDirection;
+ return direction1^direction2;
}
-Foam::tensor Foam::eigenVectors(const tensor& t)
+Foam::tensor Foam::eigenVectors(const tensor& T, const vector& lambdas)
{
- vector evals(eigenValues(t));
+ vector Ux(1, 0, 0), Uy(0, 1, 0), Uz(0, 0, 1);
- tensor evs
- (
- eigenVector(t, evals.x()),
- eigenVector(t, evals.y()),
- eigenVector(t, evals.z())
- );
+ Ux = eigenVector(T, lambdas.x(), Uy, Uz);
+ Uy = eigenVector(T, lambdas.y(), Uz, Ux);
+ Uz = eigenVector(T, lambdas.z(), Ux, Uy);
- return evs;
+ return tensor(Ux, Uy, Uz);
}
-Foam::vector Foam::eigenValues(const symmTensor& t)
+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
+(
+ const symmTensor& T,
+ const scalar lambda,
+ const vector& direction1,
+ const vector& direction2
+)
{
- return eigenValues(tensor(t));
+ return eigenVector(tensor(T), lambda, direction1, direction2);
}
-Foam::vector Foam::eigenVector(const symmTensor& t, const scalar lambda)
+Foam::tensor Foam::eigenVectors(const symmTensor& T, const vector& lambdas)
{
- return eigenVector(tensor(t), lambda);
+ return eigenVectors(tensor(T), lambdas);
}
-Foam::tensor Foam::eigenVectors(const symmTensor& t)
+Foam::tensor Foam::eigenVectors(const symmTensor& T)
{
- return eigenVectors(tensor(t));
+ return eigenVectors(tensor(T));
}
diff --git a/src/OpenFOAM/primitives/Tensor/tensor/tensor.H b/src/OpenFOAM/primitives/Tensor/tensor/tensor.H
index b92cc8e64a1..94d724c9c0b 100644
--- a/src/OpenFOAM/primitives/Tensor/tensor/tensor.H
+++ b/src/OpenFOAM/primitives/Tensor/tensor/tensor.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -50,13 +50,27 @@ namespace Foam
typedef Tensor<scalar> tensor;
-vector eigenValues(const tensor&);
-vector eigenVector(const tensor&, const scalar lambda);
-tensor eigenVectors(const tensor&);
-
-vector eigenValues(const symmTensor&);
-vector eigenVector(const symmTensor&, const scalar lambda);
-tensor eigenVectors(const symmTensor&);
+vector eigenValues(const tensor& T);
+vector eigenVector
+(
+ const tensor& T,
+ const scalar lambda,
+ const vector& direction1,
+ const vector& direction2
+);
+tensor eigenVectors(const tensor& T, const vector& lambdas);
+tensor eigenVectors(const tensor& T);
+
+vector eigenValues(const symmTensor& T);
+vector eigenVector
+(
+ const symmTensor& T,
+ const scalar lambda,
+ const vector& direction1,
+ const vector& direction2
+);
+tensor eigenVectors(const symmTensor& T, const vector& lambdas);
+tensor eigenVectors(const symmTensor& T);
//- Data associated with tensor type are contiguous
template<>
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
index 5136f383762..c64fa8fd3c8 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
@@ -458,7 +458,7 @@ void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio
// normal, as it has the greatest value. The minimum eigenvalue
// is the dominant collapse axis for high aspect ratio faces.
- collapseAxis = eigenVector(J, eVals.x());
+ collapseAxis = eigenVectors(J, eVals).x();
// The inertia calculation describes the mass distribution as a
// function of distance squared to the axis, so the square root of
--
GitLab