diff --git a/applications/test/tensor/Test-tensor.C b/applications/test/tensor/Test-tensor.C
index 6da37d33fffad0bd9cbcb01a64bc64b1b384ec89..e406ab7f3bffe1b00b3b9e13014137e211934373 100644
--- a/applications/test/tensor/Test-tensor.C
+++ b/applications/test/tensor/Test-tensor.C
@@ -69,5 +69,121 @@ int main()
Info<< (symm(t7) && t7) - (0.5*(t7 + t7.T()) && t7) << endl;
Info<< (t7 && symm(t7)) - (t7 && 0.5*(t7 + t7.T())) << endl;
+
+ /*
+ // Lots of awkward eigenvector tests ...
+
+ tensor T_rand_real
+ (
+ 0.9999996423721313, 0.3330855667591095, 0.6646450161933899,
+ 0.9745196104049683, 0.0369445420801640, 0.0846728682518005,
+ 0.6474838852882385, 0.1617118716239929, 0.2041363865137100
+ );
+ Debug(T_rand_real);
+ vector L_rand_real(eigenValues(T_rand_real));
+ Debug(L_rand_real);
+ tensor U_rand_real(eigenVectors(T_rand_real));
+ Debug(U_rand_real);
+
+ Info << endl << endl;
+
+ tensor T_rand_imag
+ (
+ 0.8668024539947510, 0.1664607226848602, 0.8925783634185791,
+ 0.9126510620117188, 0.7408077120780945, 0.1499115079641342,
+ 0.0936608463525772, 0.7615650296211243, 0.8953040242195129
+ );
+ Debug(T_rand_imag);
+ vector L_rand_imag(eigenValues(T_rand_imag));
+ Debug(L_rand_imag);
+ tensor U_rand_imag(eigenVectors(T_rand_imag));
+ Debug(U_rand_imag);
+
+ Info << endl << endl;
+
+ tensor T_rand_symm
+ (
+ 1.9999992847442627, 1.3076051771640778, 1.3121289014816284,
+ 1.3076051771640778, 0.0738890841603279, 0.2463847398757935,
+ 1.3121289014816284, 0.2463847398757935, 0.4082727730274200
+ );
+ Debug(T_rand_symm);
+ vector L_rand_symm(eigenValues(T_rand_symm));
+ Debug(L_rand_symm);
+ tensor U_rand_symm(eigenVectors(T_rand_symm));
+ Debug(U_rand_symm);
+
+ Info << endl << endl;
+
+ symmTensor T_rand_Symm
+ (
+ 1.9999992847442627, 1.3076051771640778, 1.3121289014816284,
+ 0.0738890841603279, 0.2463847398757935,
+ 0.4082727730274200
+ );
+ Debug(T_rand_Symm);
+ vector L_rand_Symm(eigenValues(T_rand_Symm));
+ Debug(L_rand_Symm);
+ tensor U_rand_Symm(eigenVectors(T_rand_Symm));
+ Debug(U_rand_Symm);
+
+ Info << endl << endl;
+
+ tensor T_rand_diag
+ (
+ 0.8668024539947510, 0, 0,
+ 0, 0.7408077120780945, 0,
+ 0, 0, 0.8953040242195129
+ );
+ Debug(T_rand_diag);
+ vector L_rand_diag(eigenValues(T_rand_diag));
+ Debug(L_rand_diag);
+ tensor U_rand_diag(eigenVectors(T_rand_diag));
+ Debug(U_rand_diag);
+
+ Info << endl << endl;
+
+ tensor T_repeated
+ (
+ 0, 1, 1,
+ 1, 0, 1,
+ 1, 1, 0
+ );
+ Debug(T_repeated);
+ vector L_repeated(eigenValues(T_repeated));
+ Debug(L_repeated);
+ tensor U_repeated(eigenVectors(T_repeated));
+ Debug(U_repeated);
+
+ Info << endl << endl;
+
+ tensor T_repeated_zero
+ (
+ 1, 1, 1,
+ 1, 1, 1,
+ 1, 1, 1
+ );
+ Debug(T_repeated_zero);
+ vector L_repeated_zero(eigenValues(T_repeated_zero));
+ Debug(L_repeated_zero);
+ tensor U_repeated_zero(eigenVectors(T_repeated_zero));
+ Debug(U_repeated_zero);
+
+ Info << endl << endl;
+
+ tensor T_triple
+ (
+ 2, 0, 0,
+ 0, 2, 0,
+ 0, 0, 2
+ );
+ Debug(T_triple);
+ vector L_triple(eigenValues(T_triple));
+ Debug(L_triple);
+ tensor U_triple(eigenVectors(T_triple));
+ Debug(U_triple);
+ */
+
+
return 0;
}
diff --git a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
index 7ec9702310a3ac0719ba486ea543f5a2d1680ce4..30f5ee246fb55666caeea2558fb30c7cfe3dbfbf 100644
--- a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
+++ b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
@@ -89,99 +89,80 @@ namespace Foam
Foam::vector Foam::eigenValues(const tensor& t)
{
- scalar i = 0;
- scalar ii = 0;
- scalar iii = 0;
+ // The eigenvalues
+ scalar i, ii, iii;
- if
- (
- (
- mag(t.xy()) + mag(t.xz()) + mag(t.yx())
- + mag(t.yz()) + mag(t.zx()) + mag(t.zy())
- )
- < SMALL
- )
+ // Coefficients of the characteristic polynmial
+ // x^3 + a*x^2 + b*x + c = 0
+ scalar a =
+ - 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();
+
+ 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();
+
+ // Auxillary variables
+ scalar aBy3 = a/3;
+
+ scalar P = (a*a - 3*b)/9; // == -p_wikipedia/3
+ scalar PPP = P*P*P;
+
+ scalar Q = (2*a*a*a - 9*a*b + 27*c)/54; // == q_wikipedia/2
+ scalar QQ = Q*Q;
+
+ // Three identical roots
+ if (mag(P) < SMALL && mag(Q) < SMALL)
{
- // diagonal matrix
- i = t.xx();
- ii = t.yy();
- iii = t.zz();
+ return vector(- aBy3, - aBy3, - aBy3);
}
- else
+
+ // Two identical roots and one distinct root
+ else if (mag(PPP/QQ - 1) < SMALL)
+ {
+ scalar sqrtP = sqrt(P);
+ scalar signQ = sign(Q);
+
+ i = ii = signQ*sqrtP - aBy3;
+ iii = - 2*signQ*sqrtP - aBy3;
+ }
+
+ // Three distinct roots
+ else if (PPP > QQ)
{
- scalar a = -t.xx() - t.yy() - t.zz();
+ scalar sqrtP = sqrt(P);
+ scalar value = cos(acos(Q/sqrt(PPP))/3);
+ scalar delta = sqrt(3 - 3*value*value);
- scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
- - t.xy()*t.yx() - t.xz()*t.zx() - t.yz()*t.zy();
+ i = - 2*sqrtP*value - aBy3;
+ ii = sqrtP*(value + delta) - aBy3;
+ iii = sqrtP*(value - delta) - aBy3;
+ }
- scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.zx()
- - t.xz()*t.yx()*t.zy() + t.xz()*t.yy()*t.zx()
- + t.xy()*t.yx()*t.zz() + t.xx()*t.yz()*t.zy();
+ // One real root, two imaginary roots
+ // based on the above logic, PPP must be less than QQ
+ else
+ {
+ WarningIn("eigenValues(const tensor&)")
+ << "complex eigenvalues detected for tensor: " << t
+ << endl;
- // If there is a zero root
- if (mag(c) < 1e-100)
+ if (mag(P) < SMALL)
{
- scalar disc = sqr(a) - 4*b;
-
- if (disc >= -SMALL)
- {
- scalar q = -0.5*sqrt(max(0.0, disc));
-
- i = 0;
- ii = -0.5*a + q;
- iii = -0.5*a - q;
- }
- else
- {
- FatalErrorIn("eigenValues(const tensor&)")
- << "zero and complex eigenvalues in tensor: " << t
- << abort(FatalError);
- }
+ i = cbrt(QQ/2);
}
else
{
- scalar Q = (a*a - 3*b)/9;
- scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
-
- scalar R2 = sqr(R);
- scalar Q3 = pow3(Q);
-
- // Three different real roots
- if (R2 < Q3)
- {
- scalar sqrtQ = sqrt(Q);
- scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ))));
-
- scalar m2SqrtQ = -2*sqrtQ;
- scalar aBy3 = a/3;
-
- i = m2SqrtQ*cos(theta/3) - aBy3;
- ii = m2SqrtQ*cos((theta + twoPi)/3) - aBy3;
- iii = m2SqrtQ*cos((theta - twoPi)/3) - aBy3;
- }
- else
- {
- scalar A = cbrt(R + sqrt(R2 - Q3));
-
- // Three equal real roots
- if (A < SMALL)
- {
- scalar root = -a/3;
- return vector(root, root, root);
- }
- else
- {
- // Complex roots
- WarningIn("eigenValues(const tensor&)")
- << "complex eigenvalues detected for tensor: " << t
- << endl;
-
- return vector::zero;
- }
- }
+ scalar w = cbrt(- Q - sqrt(QQ - PPP));
+ i = w + P/w - aBy3;
}
- }
+ return vector(-VGREAT, i, VGREAT);
+ }
// Sort the eigenvalues into ascending order
if (i > ii)
@@ -203,24 +184,35 @@ Foam::vector Foam::eigenValues(const tensor& t)
}
-Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
+Foam::vector Foam::eigenVector
+(
+ const tensor& t,
+ const scalar lambda
+)
{
- if (mag(lambda) < SMALL)
- {
- return vector::zero;
- }
-
- // Construct the matrix for the eigenvector problem
+ // Constantly rotating direction ensures different eigenvectors are
+ // generated when called sequentially with a multiple eigenvalue
+ static vector direction(0,0,1);
+ 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);
- // Calculate the sub-determinants of the 3 components
- scalar sd0 = A.yy()*A.zz() - A.yz()*A.zy();
- scalar sd1 = A.xx()*A.zz() - A.xz()*A.zx();
- scalar sd2 = A.xx()*A.yy() - A.xy()*A.yx();
+ // Determinants of the 2x2 sub-matrices used to find the eigenvectors
+ scalar sd0, sd1, sd2;
+ scalar magSd0, magSd1, magSd2;
- scalar magSd0 = mag(sd0);
- scalar magSd1 = mag(sd1);
- scalar magSd2 = mag(sd2);
+ // 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);
// Evaluate the eigenvector using the largest sub-determinant
if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
@@ -231,9 +223,8 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
(A.yz()*A.zx() - A.zz()*A.yx())/sd0,
(A.zy()*A.yx() - A.yy()*A.zx())/sd0
);
- ev /= mag(ev);
- return ev;
+ return ev/mag(ev);
}
else if (magSd1 >= magSd2 && magSd1 > SMALL)
{
@@ -243,9 +234,8 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
1,
(A.zx()*A.xy() - A.xx()*A.zy())/sd1
);
- ev /= mag(ev);
- return ev;
+ return ev/mag(ev);
}
else if (magSd2 > SMALL)
{
@@ -255,167 +245,17 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
(A.yx()*A.xz() - A.xx()*A.yz())/sd2,
1
);
- ev /= mag(ev);
-
- return ev;
- }
- else
- {
- return vector::zero;
- }
-}
-
-Foam::tensor Foam::eigenVectors(const tensor& t)
-{
- vector evals(eigenValues(t));
-
- tensor evs
- (
- eigenVector(t, evals.x()),
- eigenVector(t, evals.y()),
- eigenVector(t, evals.z())
- );
-
- return evs;
-}
-
-
-// Return eigenvalues in ascending order of absolute values
-Foam::vector Foam::eigenValues(const symmTensor& t)
-{
- scalar i = 0;
- scalar ii = 0;
- scalar iii = 0;
-
- if
- (
- (
- mag(t.xy()) + mag(t.xz()) + mag(t.xy())
- + mag(t.yz()) + mag(t.xz()) + mag(t.yz())
- )
- < SMALL
- )
- {
- // diagonal matrix
- i = t.xx();
- ii = t.yy();
- iii = t.zz();
+ return ev/mag(ev);
}
- else
- {
- scalar a = -t.xx() - t.yy() - t.zz();
-
- scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
- - t.xy()*t.xy() - t.xz()*t.xz() - t.yz()*t.yz();
- scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.xz()
- - t.xz()*t.xy()*t.yz() + t.xz()*t.yy()*t.xz()
- + t.xy()*t.xy()*t.zz() + t.xx()*t.yz()*t.yz();
-
- // If there is a zero root
- if (mag(c) < 1e-100)
- {
- scalar disc = sqr(a) - 4*b;
-
- if (disc >= -SMALL)
- {
- scalar q = -0.5*sqrt(max(0.0, disc));
-
- i = 0;
- ii = -0.5*a + q;
- iii = -0.5*a - q;
- }
- else
- {
- FatalErrorIn("eigenValues(const tensor&)")
- << "zero and complex eigenvalues in tensor: " << t
- << abort(FatalError);
- }
- }
- else
- {
- scalar Q = (a*a - 3*b)/9;
- scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
-
- scalar R2 = sqr(R);
- scalar Q3 = pow3(Q);
-
- // Three different real roots
- if (R2 < Q3)
- {
- scalar sqrtQ = sqrt(Q);
- scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ))));
-
- scalar m2SqrtQ = -2*sqrtQ;
- scalar aBy3 = a/3;
-
- i = m2SqrtQ*cos(theta/3) - aBy3;
- ii = m2SqrtQ*cos((theta + twoPi)/3) - aBy3;
- iii = m2SqrtQ*cos((theta - twoPi)/3) - aBy3;
- }
- else
- {
- scalar A = cbrt(R + sqrt(R2 - Q3));
-
- // Three equal real roots
- if (A < SMALL)
- {
- scalar root = -a/3;
- return vector(root, root, root);
- }
- else
- {
- // Complex roots
- WarningIn("eigenValues(const symmTensor&)")
- << "complex eigenvalues detected for symmTensor: " << t
- << endl;
-
- return vector::zero;
- }
- }
- }
- }
-
-
- // Sort the eigenvalues into ascending order
- if (i > ii)
- {
- Swap(i, ii);
- }
-
- if (ii > iii)
- {
- Swap(ii, iii);
- }
-
- if (i > ii)
- {
- Swap(i, ii);
- }
-
- return vector(i, ii, iii);
-}
-
-
-Foam::vector Foam::eigenVector(const symmTensor& t, const scalar lambda)
-{
- if (mag(lambda) < SMALL)
- {
- return vector::zero;
- }
-
- // Construct the matrix for the eigenvector problem
- symmTensor A(t - lambda*I);
-
- // Calculate the sub-determinants of the 3 components
- scalar sd0 = A.yy()*A.zz() - A.yz()*A.yz();
- scalar sd1 = A.xx()*A.zz() - A.xz()*A.xz();
- scalar sd2 = A.xx()*A.yy() - A.xy()*A.xy();
-
- scalar magSd0 = mag(sd0);
- scalar magSd1 = mag(sd1);
- scalar magSd2 = mag(sd2);
+ // 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();
+ magSd0 = mag(sd0);
+ magSd1 = mag(sd1);
+ magSd2 = mag(sd2);
// Evaluate the eigenvector using the largest sub-determinant
if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
@@ -423,45 +263,41 @@ Foam::vector Foam::eigenVector(const symmTensor& t, const scalar lambda)
vector ev
(
1,
- (A.yz()*A.xz() - A.zz()*A.xy())/sd0,
- (A.yz()*A.xy() - A.yy()*A.xz())/sd0
+ (A.yz()*direction.x() - direction.z()*A.yx())/sd0,
+ (direction.y()*A.yx() - A.yy()*direction.x())/sd0
);
- ev /= mag(ev);
- return ev;
+ return ev/mag(ev);
}
else if (magSd1 >= magSd2 && magSd1 > SMALL)
{
vector ev
(
- (A.xz()*A.yz() - A.zz()*A.xy())/sd1,
+ (direction.z()*A.zy() - A.zz()*direction.y())/sd1,
1,
- (A.xz()*A.xy() - A.xx()*A.yz())/sd1
+ (A.zx()*direction.y() - direction.x()*A.zy())/sd1
);
- ev /= mag(ev);
- return ev;
+ return ev/mag(ev);
}
else if (magSd2 > SMALL)
{
vector ev
(
- (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
- (A.xy()*A.xz() - A.xx()*A.yz())/sd2,
+ (A.xy()*direction.z() - direction.y()*A.xz())/sd2,
+ (direction.x()*A.xz() - A.xx()*direction.z())/sd2,
1
);
- ev /= mag(ev);
- return ev;
- }
- else
- {
- return vector::zero;
+ return ev/mag(ev);
}
+
+ // Triple eigenvalue
+ return oldDirection;
}
-Foam::tensor Foam::eigenVectors(const symmTensor& t)
+Foam::tensor Foam::eigenVectors(const tensor& t)
{
vector evals(eigenValues(t));
@@ -476,4 +312,22 @@ Foam::tensor Foam::eigenVectors(const symmTensor& t)
}
+Foam::vector Foam::eigenValues(const symmTensor& t)
+{
+ return eigenValues(tensor(t));
+}
+
+
+Foam::vector Foam::eigenVector(const symmTensor& t, const scalar lambda)
+{
+ return eigenVector(tensor(t), lambda);
+}
+
+
+Foam::tensor Foam::eigenVectors(const symmTensor& t)
+{
+ return eigenVectors(tensor(t));
+}
+
+
// ************************************************************************* //