Polynomial.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "Polynomial.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<int PolySize>
32 :
33  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
34  logActive_(false),
35  logCoeff_(0.0)
36 {
37  for (int i = 0; i < PolySize; ++i)
38  {
39  this->v_[i] = 0.0;
40  }
41 }
42 
43 
44 template<int PolySize>
45 Foam::Polynomial<PolySize>::Polynomial(const scalar coeffs[PolySize])
46 :
47  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
48  logActive_(false),
49  logCoeff_(0.0)
50 {
51  for (int i=0; i<PolySize; i++)
52  {
53  this->v_[i] = coeffs[i];
54  }
55 }
56 
57 
58 template<int PolySize>
60 :
61  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
62  logActive_(false),
63  logCoeff_(0.0)
64 {
65  if (coeffs.size() != PolySize)
66  {
68  << "Size mismatch: Needed " << PolySize
69  << " but given " << coeffs.size()
70  << nl << exit(FatalError);
71  }
72 
73  for (int i = 0; i < PolySize; ++i)
74  {
75  this->v_[i] = coeffs[i];
76  }
77 }
78 
79 
80 template<int PolySize>
82 :
83  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(is),
84  logActive_(false),
85  logCoeff_(0.0)
86 {}
87 
88 
89 template<int PolySize>
91 :
92  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
93  logActive_(false),
94  logCoeff_(0.0)
95 {
96  word isName(is);
97 
98  if (isName != name)
99  {
101  << "Expected polynomial name " << name << " but read " << isName
102  << nl << exit(FatalError);
103  }
104 
105  VectorSpace<Polynomial<PolySize>, scalar, PolySize>::
106  operator=(VectorSpace<Polynomial<PolySize>, scalar, PolySize>(is));
107 
108  if (this->size() == 0)
109  {
111  << "Polynomial coefficients for entry " << isName
112  << " are invalid (empty)" << nl << exit(FatalError);
113  }
114 }
115 
116 
117 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
118 
119 template<int PolySize>
121 {
122  return logActive_;
123 }
124 
125 
126 template<int PolySize>
128 {
129  return logCoeff_;
130 }
131 
132 
133 template<int PolySize>
134 Foam::scalar Foam::Polynomial<PolySize>::value(const scalar x) const
135 {
136  scalar val = this->v_[0];
137 
138  scalar powX = 1;
139  for (label i=1; i<PolySize; ++i)
140  {
141  powX *= x;
142  val += this->v_[i]*powX;
143  }
144 
145  if (logActive_)
146  {
147  val += logCoeff_*log(x);
148  }
149 
150  return val;
151 }
152 
153 
154 template<int PolySize>
155 Foam::scalar Foam::Polynomial<PolySize>::derivative(const scalar x) const
156 {
157  scalar deriv = 0;
158 
159  if (PolySize > 1)
160  {
161  deriv += this->v_[1];
162 
163  scalar powX = 1;
164  for (label i=2; i<PolySize; ++i)
165  {
166  powX *= x;
167  deriv += i*this->v_[i]*powX;
168  }
169  }
170 
171  if (logActive_)
172  {
173  deriv += logCoeff_/x;
174  }
175 
176  return deriv;
177 }
178 
179 
180 template<int PolySize>
182 (
183  const scalar x1,
184  const scalar x2
185 ) const
186 {
187  scalar powX1 = x1;
188  scalar powX2 = x2;
189 
190  scalar integ = this->v_[0]*(powX2 - powX1);
191  for (label i=1; i<PolySize; ++i)
192  {
193  powX1 *= x1;
194  powX2 *= x2;
195  integ += this->v_[i]/(i + 1)*(powX2 - powX1);
196  }
197 
198  if (logActive_)
199  {
200  integ += logCoeff_*((x2*log(x2) - x2) - (x1*log(x1) - x1));
201  }
202 
203  return integ;
204 }
205 
206 
207 template<int PolySize>
209 Foam::Polynomial<PolySize>::integral(const scalar intConstant) const
210 {
211  intPolyType newCoeffs;
212 
213  newCoeffs[0] = intConstant;
214  forAll(*this, i)
215  {
216  newCoeffs[i+1] = this->v_[i]/(i + 1);
217  }
218 
219  return newCoeffs;
220 }
221 
222 
223 template<int PolySize>
225 Foam::Polynomial<PolySize>::integralMinus1(const scalar intConstant) const
226 {
227  polyType newCoeffs;
228 
229  if (this->v_[0] > vSmall)
230  {
231  newCoeffs.logActive_ = true;
232  newCoeffs.logCoeff_ = this->v_[0];
233  }
234 
235  newCoeffs[0] = intConstant;
236  for (label i=1; i<PolySize; ++i)
237  {
238  newCoeffs[i] = this->v_[i]/i;
239  }
240 
241  return newCoeffs;
242 }
243 
244 
245 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
bool logActive() const
Return true if the log term is active.
Definition: Polynomial.C:120
dimensionedScalar log(const dimensionedScalar &ds)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
Templated vector space.
Definition: VectorSpace.H:53
scalar integral(const scalar x1, const scalar x2) const
Return integral between two values.
Definition: Polynomial.C:182
Polynomial()
Construct null, with all coefficients = 0.0.
Definition: Polynomial.C:31
A class for handling words, derived from string.
Definition: word.H:59
scalar value(const scalar x) const
Return polynomial value.
Definition: Polynomial.C:134
static const char nl
Definition: Ostream.H:260
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Polynomial templated on size (order):
Definition: Polynomial.H:65
scalar logCoeff() const
Return the log coefficient.
Definition: Polynomial.C:127
polyType integralMinus1(const scalar intConstant=0.0) const
Return integral coefficients when lowest order is -1.
Definition: Polynomial.C:225
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
scalar derivative(const scalar x) const
Return derivative of the polynomial at the given x.
Definition: Polynomial.C:155