cubicEqn.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2017 OpenFOAM Foundation
9 Copyright (C) 2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "linearEqn.H"
30#include "quadraticEqn.H"
31#include "cubicEqn.H"
32
33// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
34
36{
37 const scalar a = this->a();
38 const scalar b = this->b();
39 const scalar c = this->c();
40 const scalar d = this->d();
41
42 #ifdef FULLDEBUG
43 Info<< "#DEBUG#" << nl
44 << "Coefficients of the characteristic cubic polynomial:" << nl
45 << "a = " << a << nl
46 << "b = " << b << nl
47 << "c = " << c << nl
48 << "d = " << d << nl
49 << "#######" << endl;
50 #endif
51
52 // Check the leading term in the cubic eqn exists
53 if (mag(a) < VSMALL)
54 {
55 return Roots<3>(quadraticEqn(b, c, d).roots(), roots::nan, 0);
56 }
57
58 // (JLM:p. 2246) [p = a*c - b*b/3]
59 const scalar w = a*c;
60 const scalar p = -(fma(-a, c, w) + fma(b, b/3.0, -w));
61 const scalar q = b*b*b*scalar(2)/27 - b*c*a/3 + d*a*a;
62 const scalar numDiscr = p*p*p/27 + q*q/4;
63 const scalar discr = (mag(numDiscr) > VSMALL) ? numDiscr : 0;
64
65 // Determine the number and types of the roots
66 const bool threeReal = discr < 0;
67 const bool oneRealTwoComplex = discr > 0;
68 const bool twoReal = //p != 0; && discr == 0;
69 (mag(p) > sqrt(SMALL)) && !(threeReal || oneRealTwoComplex);
70 // const bool oneReal = p == 0 && discr == 0;
71
72 #ifdef FULLDEBUG
73 Info<< "#DEBUG#" << nl
74 << "Numerical discriminant:" << tab << numDiscr << nl
75 << "Adjusted discriminant:" << tab << discr << nl
76 << "Number and types of the roots:" << nl
77 << "threeReal = " << threeReal << nl
78 << "oneRealTwoComplex = " << oneRealTwoComplex << nl
79 << "twoReal = " << twoReal << nl
80 << "oneReal = " << !(threeReal || oneRealTwoComplex || twoReal) << nl
81 << "#######" << endl;
82 #endif
83
84 static const scalar sqrt3 = sqrt(3.0);
85
86 scalar x = 0;
87
88 if (threeReal)
89 {
90 const scalar wCbRe = -q/2;
91 const scalar wCbIm = sqrt(-discr);
92 const scalar wAbs = cbrt(hypot(wCbRe, wCbIm));
93 const scalar wArg = atan2(wCbIm, wCbRe)/3;
94 const scalar wRe = wAbs*cos(wArg);
95 const scalar wIm = wAbs*sin(wArg);
96
97 if (b > 0)
98 {
99 x = -wRe - mag(wIm)*sqrt3 - b/3;
100 }
101 else
102 {
103 x = 2*wRe - b/3;
104 }
105 }
106 else if (oneRealTwoComplex)
107 {
108 const scalar wCb = -q/2 - sign(q)*sqrt(discr);
109 const scalar w = cbrt(wCb);
110 const scalar t = w - p/(3*w);
111
112 if (p + t*b < 0)
113 {
114 x = t - b/3;
115 }
116 else
117 {
118 const scalar xRe = -t/2 - b/3;
119 const scalar xIm = sqrt3/2*(w + p/3/w);
120
121 return
123 (
124 Roots<1>(roots::real, -a*d/(xRe*xRe + xIm*xIm)),
126 (
129 )
130 );
131 }
132 }
133 else if (twoReal)
134 {
135 if (q*b > 0)
136 {
137 x = -2*cbrt(q/2) - b/3;
138 }
139 else
140 {
141 x = cbrt(q/2) - b/3;
142 const Roots<1> r(linearEqn(-a, x).roots());
143 return Roots<3>(Roots<2>(r, r), linearEqn(x*x, a*d).roots());
144 }
145 }
146 else // (oneReal)
147 {
148 const Roots<1> r(linearEqn(a, b/3).roots());
149 return Roots<3>(r.type(0), r[0]);
150 }
151
152 #if FULLDEBUG
153 Info<< "#DEBUG#" << nl
154 << "x = " << x << nl
155 << "#######" << endl;
156 #endif
157
158 return
160 (
161 linearEqn(-a, x).roots(),
162 quadraticEqn(-x*x, c*x + a*d, d*x).roots()
163 );
164}
165
166
167// ************************************************************************* //
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
Definition: Roots.H:73
void type(const direction i, const roots::type t)
Set the type of the i-th root.
Definition: RootsI.H:122
scalar d() const
Definition: cubicEqnI.H:75
scalar c() const
Definition: cubicEqnI.H:69
Roots< 3 > roots() const
Return the roots of the cubic equation with no particular order.
Definition: cubicEqn.C:35
scalar a() const
Definition: cubicEqnI.H:57
scalar b() const
Definition: cubicEqnI.H:63
Container to encapsulate various operations for linear equation of the forms with real coefficients:
Definition: linearEqn.H:63
Container to encapsulate various operations for quadratic equation of the forms with real coefficient...
Definition: quadraticEqn.H:91
volScalarField & p
@ complex
Definition: Roots.H:57
@ real
Definition: Roots.H:56
dimensionedScalar sign(const dimensionedScalar &ds)
dimensionedScalar sin(const dimensionedScalar &ds)
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensionedScalar atan2(const dimensionedScalar &x, const dimensionedScalar &y)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
dimensionedScalar hypot(const dimensionedScalar &x, const dimensionedScalar &y)
dimensionedScalar cbrt(const dimensionedScalar &ds)
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
constexpr char tab
The tab '\t' character(0x09)
Definition: Ostream.H:52