From 0be6269add80fd3037e64de38bcd9d5a64e6d1cf Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Wed, 22 Mar 2017 15:11:54 +0000 Subject: [PATCH] Added robust primitive cubic/quadratic/linear equation solutions. Applied to eigen-value calculations. Fixed repeated-eigen-value issues in eigen-vector generation. --- applications/test/cubicEqn/Make/files | 3 + applications/test/cubicEqn/Make/options | 3 + applications/test/cubicEqn/Test-cubicEqn.C | 101 ++++++++ src/OpenFOAM/Make/files | 3 + .../primitives/Tensor/tensor/tensor.C | 229 +++++++----------- .../primitives/Tensor/tensor/tensor.H | 30 ++- .../primitives/Tensor2D/tensor2D/tensor2D.C | 130 +++++----- .../primitives/Tensor2D/tensor2D/tensor2D.H | 10 +- .../primitives/polynomialEqns/Roots.H | 131 ++++++++++ .../primitives/polynomialEqns/RootsI.H | 136 +++++++++++ .../polynomialEqns/cubicEqn/cubicEqn.C | 164 +++++++++++++ .../polynomialEqns/cubicEqn/cubicEqn.H | 115 +++++++++ .../polynomialEqns/cubicEqn/cubicEqnI.H | 115 +++++++++ .../polynomialEqns/linearEqn/linearEqn.H | 104 ++++++++ .../polynomialEqns/linearEqn/linearEqnI.H | 111 +++++++++ .../quadraticEqn/quadraticEqn.C | 88 +++++++ .../quadraticEqn/quadraticEqn.H | 107 ++++++++ .../quadraticEqn/quadraticEqnI.H | 101 ++++++++ .../polyTopoChange/edgeCollapser.C | 4 +- 19 files changed, 1473 insertions(+), 212 deletions(-) create mode 100644 applications/test/cubicEqn/Make/files create mode 100644 applications/test/cubicEqn/Make/options create mode 100644 applications/test/cubicEqn/Test-cubicEqn.C create mode 100644 src/OpenFOAM/primitives/polynomialEqns/Roots.H create mode 100644 src/OpenFOAM/primitives/polynomialEqns/RootsI.H create mode 100644 src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C create mode 100644 src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H create mode 100644 src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H create mode 100644 src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H create mode 100644 src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H create mode 100644 src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C create mode 100644 src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H create mode 100644 src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H diff --git a/applications/test/cubicEqn/Make/files b/applications/test/cubicEqn/Make/files new file mode 100644 index 0000000000..91b9e84dcc --- /dev/null +++ b/applications/test/cubicEqn/Make/files @@ -0,0 +1,3 @@ +Test-cubicEqn.C + +EXE = $(FOAM_USER_APPBIN)/Test-cubicEqn diff --git a/applications/test/cubicEqn/Make/options b/applications/test/cubicEqn/Make/options new file mode 100644 index 0000000000..4c3dd783cb --- /dev/null +++ b/applications/test/cubicEqn/Make/options @@ -0,0 +1,3 @@ +EXE_INC = + +EXE_LIBS = diff --git a/applications/test/cubicEqn/Test-cubicEqn.C b/applications/test/cubicEqn/Test-cubicEqn.C new file mode 100644 index 0000000000..d0ab606666 --- /dev/null +++ b/applications/test/cubicEqn/Test-cubicEqn.C @@ -0,0 +1,101 @@ +#include +#include + +#include "cubicEqn.H" +#include "IOstreams.H" +#include "stringList.H" + +using namespace Foam; + +scalar randomScalar(const scalar min, const scalar max) +{ + static_assert + ( + sizeof(long) == sizeof(scalar), + "Scalar and long are not the same size" + ); + static std::default_random_engine generator(std::time(0)); + static std::uniform_int_distribution + distribution + ( + std::numeric_limits::min(), + std::numeric_limits::max() + ); + scalar x; + do + { + long i = distribution(generator); + x = reinterpret_cast(i); + } + while (min > mag(x) || mag(x) > max || !std::isfinite(x)); + return x; +}; + +template +void test(const Type& polynomialEqn, const scalar tol) +{ + Roots r = polynomialEqn.roots(); + + const scalar nan = std::numeric_limits::quiet_NaN(); + const scalar inf = std::numeric_limits::infinity(); + + FixedList t; + FixedList v(nan); + FixedList e(nan); + bool ok = true; + forAll(r, i) + { + t[i] = r.type(i); + switch (t[i]) + { + case roots::real: + v[i] = polynomialEqn.value(r[i]); + e[i] = polynomialEqn.error(r[i]); + ok = ok && mag(v[i]) < tol*mag(e[i]); + break; + case roots::posInf: + v[i] = + inf; + e[i] = nan; + break; + case roots::negInf: + v[i] = - inf; + e[i] = nan; + break; + default: + v[i] = e[i] = nan; + break; + } + } + + if (!ok) + { + Info<< "Coeffs: " << polynomialEqn << endl + << " Types: " << t << endl + << " Roots: " << r << endl + << "Values: " << v << endl + << "Errors: " << e << endl << endl; + } +} + +int main() +{ + const int nTests = 1000000; + for (int t = 0; t < nTests; ++ t) + { + test + ( + cubicEqn + ( + randomScalar(1e-50, 1e+50), + randomScalar(1e-50, 1e+50), + randomScalar(1e-50, 1e+50), + randomScalar(1e-50, 1e+50) + ), + 100 + ); + } + + Info << nTests << " cubics tested" << endl; + + return 0; +} diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 84f83567d0..44ef405e41 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -125,6 +125,9 @@ $(spatialVectorAlgebra)/SpatialVector/spatialVector/spatialVector.C $(spatialVectorAlgebra)/SpatialTensor/spatialTensor/spatialTensor.C $(spatialVectorAlgebra)/CompactSpatialTensor/compactSpatialTensor/compactSpatialTensor.C +primitives/polynomialEqns/cubicEqn/cubicEqn.C +primitives/polynomialEqns/quadraticEqn/quadraticEqn.C + containers/HashTables/HashTable/HashTableCore.C containers/HashTables/StaticHashTable/StaticHashTableCore.C containers/Lists/SortableList/ParSortableListName.C diff --git a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C index 2d8da5e71c..0914c279a4 100644 --- a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C +++ b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "tensor.H" +#include "cubicEqn.H" #include "mathematicalConstants.H" using namespace Foam::constant::mathematical; @@ -72,137 +73,76 @@ const Foam::tensor Foam::tensor::I Foam::vector Foam::eigenValues(const tensor& t) { - // The eigenvalues - scalar i, ii, iii; - - // diagonal matrix - if - ( - ( - mag(t.xy()) + mag(t.xz()) + mag(t.yx()) - + mag(t.yz()) + mag(t.zx()) + mag(t.zy()) - ) - < SMALL - ) - { - i = t.xx(); - ii = t.yy(); - iii = t.zz(); - } - - // non-diagonal matrix - else + // 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(); + + // Solve + Roots<3> roots = cubicEqn(1, b, c, d).roots(); + + // Check the root types + vector lambda = vector::zero; + forAll(roots, i) { - // 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) + switch (roots.type(i)) { - return vector(- aBy3, - aBy3, - aBy3); - } - - // Two identical roots and one distinct root - else if (mag(QQ) > SMALL && 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 sqrtP = sqrt(P); - scalar value = cos(acos(Q/sqrt(PPP))/3); - scalar delta = sqrt(3 - 3*value*value); - - i = - 2*sqrtP*value - aBy3; - ii = sqrtP*(value + delta) - aBy3; - iii = sqrtP*(value - delta) - aBy3; - } - - // One real root, two imaginary roots - // based on the above logic, PPP must be less than QQ - else - { - WarningInFunction - << "complex eigenvalues detected for tensor: " << t - << endl; - - if (mag(P) < SMALL) - { - i = cbrt(QQ/2); - } - else - { - scalar w = cbrt(- Q - sqrt(QQ - PPP)); - i = w + P/w - aBy3; - } - - return vector(-VGREAT, i, VGREAT); + 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); } } // Sort the eigenvalues into ascending order - if (i > ii) + if (lambda.x() > lambda.y()) { - Swap(i, ii); + Swap(lambda.x(), lambda.y()); } - - if (ii > iii) + if (lambda.y() > lambda.z()) { - Swap(ii, iii); + Swap(lambda.y(), lambda.z()); } - - if (i > ii) + if (lambda.x() > lambda.y()) { - Swap(i, ii); + Swap(lambda.x(), lambda.y()); } - return vector(i, ii, iii); + return lambda; } 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 +192,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 +205,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 +215,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 +226,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 +235,57 @@ Foam::vector Foam::eigenVector } // Triple eigenvalue - return oldDirection; + return direction1^direction2; +} + + +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) +Foam::tensor Foam::eigenVectors(const tensor& T) { - vector evals(eigenValues(t)); + const vector lambdas(eigenValues(T)); + + return eigenVectors(T, lambdas); +} - tensor evs - ( - eigenVector(t, evals.x()), - eigenVector(t, evals.y()), - eigenVector(t, evals.z()) - ); - return evs; +Foam::vector Foam::eigenValues(const symmTensor& T) +{ + return eigenValues(tensor(T)); } -Foam::vector Foam::eigenValues(const symmTensor& 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 b92cc8e64a..89966a9f72 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-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -50,13 +50,27 @@ namespace Foam typedef Tensor 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/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C b/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C index 36f9725255..bf5a12aa7f 100644 --- a/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C +++ b/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "tensor2D.H" +#include "quadraticEqn.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -85,95 +86,96 @@ const Foam::tensor2D Foam::tensor2D::I Foam::vector2D Foam::eigenValues(const tensor2D& t) { - scalar i = 0; - scalar ii = 0; + // 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(); - if (mag(t.xy()) < SMALL && mag(t.yx()) < SMALL) - { - i = t.xx(); - ii = t.yy(); - } - else - { - scalar mb = t.xx() + t.yy(); - scalar c = t.xx()*t.yy() - t.xy()*t.yx(); + // Solve + Roots<2> roots = quadraticEqn(1, b, c).roots(); - // If there is a zero root - if (mag(c) < SMALL) - { - i = 0; - ii = mb; - } - else + // Check the root types + vector2D lambda = vector2D::zero; + forAll(roots, i) + { + switch (roots.type(i)) { - scalar disc = sqr(mb) - 4*c; - - if (disc > 0) - { - scalar q = sqrt(disc); - - i = 0.5*(mb - q); - ii = 0.5*(mb + q); - } - else - { + 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 - << "zero and complex eigenvalues in tensor2D: " << t - << abort(FatalError); - } + << "Eigenvalue calculation failed for tensor: " << t + << exit(FatalError); } } // Sort the eigenvalues into ascending order - if (i > ii) + if (lambda.x() > lambda.y()) { - Swap(i, ii); + Swap(lambda.x(), lambda.y()); } - return vector2D(i, ii); + return lambda; } -Foam::vector2D Foam::eigenVector(const tensor2D& t, const scalar lambda) +Foam::vector2D Foam::eigenVector +( + const tensor2D& T, + const scalar lambda, + const vector2D& direction1 +) { - if (lambda < SMALL) - { - return vector2D::zero; - } + // Construct the linear system for this eigenvalue + tensor2D A(T - lambda*tensor2D::I); - if (mag(t.xy()) < SMALL && mag(t.yx()) < SMALL) + // Evaluate the eigenvector using the largest divisor + if (mag(A.yy()) > mag(A.xx()) && mag(A.yy()) > SMALL) { - if (lambda > min(t.xx(), t.yy())) - { - return vector2D(1, 0); - } - else - { - return vector2D(0, 1); - } - } - else if (mag(t.xy()) < SMALL) - { - return vector2D(lambda - t.yy(), t.yx()); + vector2D ev(1, - A.yx()/A.yy()); + + return ev/mag(ev); } - else + else if (mag(A.xx()) > SMALL) { - return vector2D(t.xy(), lambda - t.yy()); + vector2D ev(- A.xy()/A.xx(), 1); + + return ev/mag(ev); } + + // Repeated eigenvalue + return vector2D(- direction1.y(), direction1.x()); } -Foam::tensor2D Foam::eigenVectors(const tensor2D& t) +Foam::tensor2D Foam::eigenVectors(const tensor2D& T, const vector2D& lambdas) { - vector2D evals(eigenValues(t)); + vector2D Ux(1, 0), Uy(0, 1); + + Ux = eigenVector(T, lambdas.x(), Uy); + Uy = eigenVector(T, lambdas.y(), Ux); + + return tensor2D(Ux, Uy); +} + - tensor2D evs - ( - eigenVector(t, evals.x()), - eigenVector(t, evals.y()) - ); +Foam::tensor2D Foam::eigenVectors(const tensor2D& T) +{ + const vector2D lambdas(eigenValues(T)); - return evs; + return eigenVectors(T, lambdas); } diff --git a/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.H b/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.H index 1b7cf2f600..feb7f0e275 100644 --- a/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.H +++ b/src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.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-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -49,7 +49,13 @@ namespace Foam typedef Tensor2D tensor2D; vector2D eigenValues(const tensor2D& t); -vector2D eigenVector(const tensor2D& t, const scalar lambda); +vector2D eigenVector +( + const tensor2D& t, + const scalar lambda, + const vector2D& direction1 +); +tensor2D eigenVectors(const tensor2D& t, const vector2D& lambdas); tensor2D eigenVectors(const tensor2D& t); //- Data associated with tensor2D type are contiguous diff --git a/src/OpenFOAM/primitives/polynomialEqns/Roots.H b/src/OpenFOAM/primitives/polynomialEqns/Roots.H new file mode 100644 index 0000000000..27829eee46 --- /dev/null +++ b/src/OpenFOAM/primitives/polynomialEqns/Roots.H @@ -0,0 +1,131 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::Roots + +Description + Templated storage for the roots of polynomial equations, plus flags to + indicate the nature of the roots. + +SourceFiles + RootsI.H + Roots.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Roots_H +#define Roots_H + +#include "VectorSpace.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace roots +{ + +//- Types of root +enum type +{ + real = 0, + complex, + posInf, + negInf, + nan +}; + +} + +/*---------------------------------------------------------------------------*\ + Class Roots Declaration +\*---------------------------------------------------------------------------*/ + +template +class Roots +: + public VectorSpace, scalar, N> +{ + // Private data + + //- Root types, encoded into a single integer + label types_; + +public: + + // Constructors + + //- Construct null + inline Roots(); + + //- Construct with a uniform value + inline Roots(const roots::type t, const scalar x); + + //- Construct by concatenation + inline Roots + ( + const roots::type t, + const scalar x, + const Roots& xs + ); + + //- Construct by concatenation + inline Roots + ( + const Roots& xs, + const roots::type t, + const scalar x + ); + + //- Construct by concatenation + template + inline Roots(const Roots& xs, const Roots& ys); + + + // Member Functions + + //- Set the type of the i-th root + inline void type(const direction i, const roots::type t); + + //- Return the type of the i-th root + inline roots::type type(const direction i) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "RootsI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/polynomialEqns/RootsI.H b/src/OpenFOAM/primitives/polynomialEqns/RootsI.H new file mode 100644 index 0000000000..dd1052f0af --- /dev/null +++ b/src/OpenFOAM/primitives/polynomialEqns/RootsI.H @@ -0,0 +1,136 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +inline Foam::Roots::Roots() +: + types_(0) +{ + forAll(*this, i) + { + type(i, roots::nan); + } +} + + +template +inline Foam::Roots::Roots(const roots::type t, const scalar x) +: + types_(0) +{ + forAll(*this, i) + { + this->v_[i] = x; + type(i, t); + } +} + + +template +inline Foam::Roots::Roots +( + const roots::type t, + const scalar x, + const Roots& xs +) +: + types_(0) +{ + this->v_[0] = x; + type(0, t); + forAll(xs, i) + { + this->v_[i+1] = xs[i]; + type(i + 1, xs.type(i)); + } +} + + +template +inline Foam::Roots::Roots +( + const Roots& xs, + const roots::type t, + const scalar x +) +: + types_(0) +{ + forAll(xs, i) + { + this->v_[i] = xs[i]; + type(i, xs.type(i)); + } + this->v_[N-1] = x; + type(N - 1, t); +} + + +template +template +inline Foam::Roots::Roots +( + const Roots& xs, + const Roots& ys +) +: + types_(0) +{ + forAll(xs, i) + { + this->v_[i] = xs[i]; + type(i, xs.type(i)); + } + forAll(ys, i) + { + this->v_[i + M] = ys[i]; + type(i + M, ys.type(i)); + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +inline void Foam::Roots::type +( + const direction i, + const roots::type t +) +{ + types_ += (t - type(i)) << 3*i; +} + + +template +inline Foam::roots::type Foam::Roots::type(const direction i) const +{ + return static_cast((types_ >> 3*i) & 7); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C new file mode 100644 index 0000000000..c46b794c38 --- /dev/null +++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C @@ -0,0 +1,164 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "linearEqn.H" +#include "quadraticEqn.H" +#include "cubicEqn.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::Roots<3> Foam::cubicEqn::roots() const +{ + /* + + This function solves a cubic equation of the following form: + + a*x^3 + b*x^2 + c*x + d = 0 + x^3 + B*x^2 + C*x + D = 0 + + The following two substitutions are used: + + x = t - B/3 + t = w - P/3/w + + This reduces the problem to a quadratic in w^3. + + w^6 + Q*w^3 - P = 0 + + Where Q and P are given in the code below. + + The properties of the cubic can be related to the properties of this + quadratic in w^3. If it has a repeated root a zero, the cubic has a tripl + root. If it has a repeated root not at zero, the cubic has two real roots, + one repeated and one not. If it has two complex roots, the cubic has three + real roots. If it has two real roots, then the cubic has one real root and + two complex roots. + + This is solved for the most numerically accurate value of w^3. See the + quadratic function for details on how to pick a value. This single value of + w^3 can yield up to three cube roots for w, which relate to the three + solutions for x. + + Only a single root, or pair of conjugate roots, is directly evaluated; the + one, or ones with the lowest relative numerical error. Root identities are + then used to recover the remaining roots, possibly utilising a quadratic + and/or linear solution. This seems to be a good way of maintaining the + accuracy of roots at very different magnitudes. + + */ + + const scalar a = this->a(); + const scalar b = this->b(); + const scalar c = this->c(); + const scalar d = this->d(); + + if (a == 0) + { + return Roots<3>(quadraticEqn(b, c, d).roots(), roots::nan, 0); + } + + // This is assumed not to over- or under-flow. If it does, all bets are off. + const scalar p = c*a - b*b/3; + const scalar q = b*b*b*(2.0/27.0) - b*c*a/3 + d*a*a; + const scalar disc = p*p*p/27 + q*q/4; + + // How many roots of what types are available? + const bool oneReal = disc == 0 && p == 0; + const bool twoReal = disc == 0 && p != 0; + const bool threeReal = disc < 0; + //const bool oneRealTwoComplex = disc > 0; + + static const scalar sqrt3 = sqrt(3.0); + + scalar x; + + if (oneReal) + { + const Roots<1> r = linearEqn(- a, b/3).roots(); + return Roots<3>(r.type(0), r[0]); + } + else if (twoReal) + { + if (q*b > 0) + { + x = - 2*cbrt(q/2) - b/3; + } + else + { + x = cbrt(q/2) - b/3; + const Roots<1> r = linearEqn(- a, x).roots(); + return Roots<3>(Roots<2>(r, r), linearEqn(x*x, a*d).roots()); + } + } + else if (threeReal) + { + const scalar wCbRe = - q/2, wCbIm = sqrt(- disc); + const scalar wAbs = cbrt(hypot(wCbRe, wCbIm)); + const scalar wArg = atan2(wCbIm, wCbRe)/3; + const scalar wRe = wAbs*cos(wArg), wIm = wAbs*sin(wArg); + if (b > 0) + { + x = - wRe - mag(wIm)*sqrt3 - b/3; + } + else + { + x = 2*wRe - b/3; + } + } + else // if (oneRealTwoComplex) + { + const scalar wCb = - q/2 - sign(q)*sqrt(disc); + const scalar w = cbrt(wCb); + const scalar t = w - p/(3*w); + if (p + t*b < 0) + { + x = t - b/3; + } + else + { + const scalar xRe = - t/2 - b/3, xIm = sqrt3/2*(w + p/3/w); + x = - a*a*d/(xRe*xRe + xIm*xIm); + + // This form of solving for the remaining roots seems more stable + // for this case. This has been determined by trial and error. + return + Roots<3> + ( + linearEqn(- a, x).roots(), + quadraticEqn(a*x, x*x + b*x, - a*d).roots() + ); + } + } + + return + Roots<3> + ( + linearEqn(- a, x).roots(), + quadraticEqn(- x*x, c*x + a*d, d*x).roots() + ); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H new file mode 100644 index 0000000000..db92376145 --- /dev/null +++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.H @@ -0,0 +1,115 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::cubicEqn + +Description + Cubic equation of the form a*x^3 + b*x^2 + c*x + d = 0 + +SourceFiles + cubicEqnI.H + cubicEqn.C + +\*---------------------------------------------------------------------------*/ + +#ifndef cubicEqn_H +#define cubicEqn_H + +#include "Roots.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class cubicEqn Declaration +\*---------------------------------------------------------------------------*/ + +class cubicEqn +: + public VectorSpace +{ +public: + + //- Component labeling enumeration + enum components { A, B, C, D }; + + + // Constructors + + //- Construct null + inline cubicEqn(); + + //- Construct initialized to zero + inline cubicEqn(const Foam::zero); + + //- Construct from components + inline cubicEqn + ( + const scalar a, + const scalar b, + const scalar c, + const scalar d + ); + + + // Member Functions + + // Access + + inline scalar a() const; + inline scalar b() const; + inline scalar c() const; + inline scalar d() const; + + inline scalar& a(); + inline scalar& b(); + inline scalar& c(); + inline scalar& d(); + + //- Evaluate at x + inline scalar value(const scalar x) const; + + //- Estimate the error of evaluation at x + inline scalar error(const scalar x) const; + + //- Get the roots + Roots<3> roots() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "cubicEqnI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H new file mode 100644 index 0000000000..eb6a3283e4 --- /dev/null +++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H @@ -0,0 +1,115 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +inline Foam::cubicEqn::cubicEqn() +{} + + +inline Foam::cubicEqn::cubicEqn(const Foam::zero) +: + VectorSpace(Foam::zero()) +{} + + +inline Foam::cubicEqn::cubicEqn +( + const scalar a, + const scalar b, + const scalar c, + const scalar d +) +{ + this->v_[A] = a; + this->v_[B] = b; + this->v_[C] = c; + this->v_[D] = d; +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline Foam::scalar Foam::cubicEqn::a() const +{ + return this->v_[A]; +} + + +inline Foam::scalar Foam::cubicEqn::b() const +{ + return this->v_[B]; +} + + +inline Foam::scalar Foam::cubicEqn::c() const +{ + return this->v_[C]; +} + + +inline Foam::scalar Foam::cubicEqn::d() const +{ + return this->v_[D]; +} + + +inline Foam::scalar& Foam::cubicEqn::a() +{ + return this->v_[A]; +} + + +inline Foam::scalar& Foam::cubicEqn::b() +{ + return this->v_[B]; +} + + +inline Foam::scalar& Foam::cubicEqn::c() +{ + return this->v_[C]; +} + + +inline Foam::scalar& Foam::cubicEqn::d() +{ + return this->v_[D]; +} + + +inline Foam::scalar Foam::cubicEqn::value(const scalar x) const +{ + return x*(x*(x*a() + b()) + c()) + d(); +} + + +inline Foam::scalar Foam::cubicEqn::error(const scalar x) const +{ + return mag(SMALL*x*(x*(x*3*a() + 2*b()) + c())); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H new file mode 100644 index 0000000000..debc52cc01 --- /dev/null +++ b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqn.H @@ -0,0 +1,104 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::linearEqn + +Description + Linear equation of the form a*x + b = 0 + +SourceFiles + linearEqnI.H + +\*---------------------------------------------------------------------------*/ + +#ifndef linearEqn_H +#define linearEqn_H + +#include "Roots.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class linearEqn Declaration +\*---------------------------------------------------------------------------*/ + +class linearEqn +: + public VectorSpace +{ +public: + + //- Component labeling enumeration + enum components { A, B }; + + + // Constructors + + //- Construct null + inline linearEqn(); + + //- Construct initialized to zero + inline linearEqn(const Foam::zero); + + //- Construct from components + inline linearEqn(const scalar a, const scalar b); + + + // Member Functions + + // Access + + inline scalar a() const; + inline scalar b() const; + + inline scalar& a(); + inline scalar& b(); + + //- Evaluate at x + inline scalar value(const scalar x) const; + + //- Estimate the error of evaluation at x + inline scalar error(const scalar x) const; + + //- Get the roots + inline Roots<1> roots() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "linearEqnI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H new file mode 100644 index 0000000000..e86ccbdd0a --- /dev/null +++ b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H @@ -0,0 +1,111 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +inline Foam::linearEqn::linearEqn() +{} + + +inline Foam::linearEqn::linearEqn(const Foam::zero) +: + VectorSpace(Foam::zero()) +{} + + +inline Foam::linearEqn::linearEqn(const scalar a, const scalar b) +{ + this->v_[A] = a; + this->v_[B] = b; +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline Foam::scalar Foam::linearEqn::a() const +{ + return this->v_[A]; +} + + +inline Foam::scalar Foam::linearEqn::b() const +{ + return this->v_[B]; +} + + +inline Foam::scalar& Foam::linearEqn::a() +{ + return this->v_[A]; +} + + +inline Foam::scalar& Foam::linearEqn::b() +{ + return this->v_[B]; +} + + +inline Foam::scalar Foam::linearEqn::value(const scalar x) const +{ + return x*a() + b(); +} + + +inline Foam::scalar Foam::linearEqn::error(const scalar x) const +{ + return mag(SMALL*x*a()); +} + + +inline Foam::Roots<1> Foam::linearEqn::roots() const +{ + /* + + This function solves a linear equation of the following form: + + a*x + b = 0 + x + B = 0 + + */ + + const scalar a = this->a(); + const scalar b = this->b(); + + if (a == 0) + { + return Roots<1>(roots::nan, 0); + } + + if (mag(b/VGREAT) >= mag(a)) + { + return Roots<1>(sign(a) == sign(b) ? roots::negInf : roots::posInf, 0); + } + + return Roots<1>(roots::real, - b/a); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C new file mode 100644 index 0000000000..d84c2928c2 --- /dev/null +++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C @@ -0,0 +1,88 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "linearEqn.H" +#include "quadraticEqn.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::Roots<2> Foam::quadraticEqn::roots() const +{ + /* + + This function solves a quadraticEqn equation of the following form: + + a*x^2 + b*x + c = 0 + x^2 + B*x + C = 0 + + The quadraticEqn formula is as follows: + + x = - B/2 +- sqrt(B*B - 4*C)/2 + + If the sqrt generates a complex number, this provides the result. If not + then the real root with the smallest floating point error is calculated. + + x0 = - B/2 - sign(B)*sqrt(B*B - 4*C)/2 + + The other root is the obtained using an identity. + + x1 = C/x0 + + */ + + const scalar a = this->a(); + const scalar b = this->b(); + const scalar c = this->c(); + + if (a == 0) + { + return Roots<2>(linearEqn(b, c).roots(), roots::nan, 0); + } + + // This is assumed not to over- or under-flow. If it does, all bets are off. + const scalar disc = b*b/4 - a*c; + + // How many roots of what types are available? + const bool oneReal = disc == 0; + const bool twoReal = disc > 0; + //const bool twoComplex = disc < 0; + + if (oneReal) + { + const Roots<1> r = linearEqn(- a, b/2).roots(); + return Roots<2>(r, r); + } + else if (twoReal) + { + const scalar x = - b/2 - sign(b)*sqrt(disc); + return Roots<2>(linearEqn(- a, x).roots(), linearEqn(- x, c).roots()); + } + else // if (twoComplex) + { + return Roots<2>(roots::complex, 0); + } +} + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H new file mode 100644 index 0000000000..1327dc15c9 --- /dev/null +++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.H @@ -0,0 +1,107 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::quadraticEqn + +Description + Quadratic equation of the form a*x^2 + b*x + c = 0 + +SourceFiles + quadraticEqnI.H + quadraticEqn.C + +\*---------------------------------------------------------------------------*/ + +#ifndef quadraticEqn_H +#define quadraticEqn_H + +#include "Roots.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class quadraticEqn Declaration +\*---------------------------------------------------------------------------*/ + +class quadraticEqn +: + public VectorSpace +{ +public: + + //- Component labeling enumeration + enum components { A, B, C }; + + + // Constructors + + //- Construct null + inline quadraticEqn(); + + //- Construct initialized to zero + inline quadraticEqn(const Foam::zero); + + //- Construct from components + inline quadraticEqn(const scalar a, const scalar b, const scalar c); + + + // Member Functions + + // Access + + inline scalar a() const; + inline scalar b() const; + inline scalar c() const; + + inline scalar& a(); + inline scalar& b(); + inline scalar& c(); + + //- Evaluate at x + inline scalar value(const scalar x) const; + + //- Estimate the error of evaluation at x + inline scalar error(const scalar x) const; + + //- Get the roots + Roots<2> roots() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "quadraticEqnI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H new file mode 100644 index 0000000000..4b4aacae61 --- /dev/null +++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H @@ -0,0 +1,101 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +inline Foam::quadraticEqn::quadraticEqn() +{} + + +inline Foam::quadraticEqn::quadraticEqn(const Foam::zero) +: + VectorSpace(Foam::zero()) +{} + + +inline Foam::quadraticEqn::quadraticEqn +( + const scalar a, + const scalar b, + const scalar c +) +{ + this->v_[A] = a; + this->v_[B] = b; + this->v_[C] = c; +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline Foam::scalar Foam::quadraticEqn::a() const +{ + return this->v_[A]; +} + + +inline Foam::scalar Foam::quadraticEqn::b() const +{ + return this->v_[B]; +} + + +inline Foam::scalar Foam::quadraticEqn::c() const +{ + return this->v_[C]; +} + + +inline Foam::scalar& Foam::quadraticEqn::a() +{ + return this->v_[A]; +} + + +inline Foam::scalar& Foam::quadraticEqn::b() +{ + return this->v_[B]; +} + + +inline Foam::scalar& Foam::quadraticEqn::c() +{ + return this->v_[C]; +} + + +inline Foam::scalar Foam::quadraticEqn::value(const scalar x) const +{ + return x*(x*a() + b()) + c(); +} + + +inline Foam::scalar Foam::quadraticEqn::error(const scalar x) const +{ + return mag(SMALL*x*(x*2*a() + b())); +} + + +// ************************************************************************* // diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C index c0eb30678c..9933c780d4 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -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