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 -------------------------------------------------------------------------------
11 License
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
122  Roots<3>
123  (
124  Roots<1>(roots::real, -a*d/(xRe*xRe + xIm*xIm)),
125  Roots<2>
126  (
127  Roots<1>(roots::complex, xRe),
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
159  Roots<3>
160  (
161  linearEqn(-a, x).roots(),
162  quadraticEqn(-x*x, c*x + a*d, d*x).roots()
163  );
164 }
165 
166 
167 // ************************************************************************* //
Foam::Roots::type
void type(const direction i, const roots::type t)
Set the type of the i-th root.
Definition: RootsI.H:122
Foam::roots::complex
Definition: Roots.H:57
p
volScalarField & p
Definition: createFieldRefs.H:8
cubicEqn.H
Foam::roots::nan
Definition: Roots.H:60
Foam::cubicEqn::b
scalar b() const
Definition: cubicEqnI.H:63
Foam::sin
dimensionedScalar sin(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:264
Foam::atan2
dimensionedScalar atan2(const dimensionedScalar &x, const dimensionedScalar &y)
Definition: dimensionedScalar.C:312
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::cubicEqn::roots
Roots< 3 > roots() const
Return the roots of the cubic equation with no particular order.
Definition: cubicEqn.C:35
Foam::sign
dimensionedScalar sign(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:166
Foam::cubicEqn::a
scalar a() const
Definition: cubicEqnI.H:57
Foam::linearEqn
Container to encapsulate various operations for linear equation of the forms with real coefficients:
Definition: linearEqn.H:60
Foam::hypot
dimensionedScalar hypot(const dimensionedScalar &x, const dimensionedScalar &y)
Definition: dimensionedScalar.C:327
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::roots::real
Definition: Roots.H:56
linearEqn.H
Foam::cubicEqn::d
scalar d() const
Definition: cubicEqnI.H:75
Foam::tab
constexpr char tab
Definition: Ostream.H:403
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::cubicEqn::c
scalar c() const
Definition: cubicEqnI.H:69
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
x
x
Definition: LISASMDCalcMethod2.H:52
Foam::Roots
Templated storage for the roots of polynomial equations, plus flags to indicate the nature of the roo...
Definition: Roots.H:70
quadraticEqn.H
Foam::cbrt
dimensionedScalar cbrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:155
Foam::quadraticEqn
Container to encapsulate various operations for quadratic equation of the forms with real coefficient...
Definition: quadraticEqn.H:88
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:265