Polynomial.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 | Copyright (C) 2011-2012 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>
46 (
47  const Polynomial<PolySize>& poly
48 )
49 :
50  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(poly),
51  logActive_(poly.logActive_),
52  logCoeff_(poly.logCoeff_)
53 {}
54 
55 
56 template<int PolySize>
57 Foam::Polynomial<PolySize>::Polynomial(const scalar coeffs[PolySize])
58 :
59  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
60  logActive_(false),
61  logCoeff_(0.0)
62 {
63  for (int i=0; i<PolySize; i++)
64  {
65  this->v_[i] = coeffs[i];
66  }
67 }
68 
69 
70 template<int PolySize>
72 :
73  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
74  logActive_(false),
75  logCoeff_(0.0)
76 {
77  if (coeffs.size() != PolySize)
78  {
80  (
81  "Polynomial<PolySize>::Polynomial(const UList<scalar>&)"
82  ) << "Size mismatch: Needed " << PolySize
83  << " but given " << coeffs.size()
84  << nl << exit(FatalError);
85  }
86 
87  for (int i = 0; i < PolySize; ++i)
88  {
89  this->v_[i] = coeffs[i];
90  }
91 }
92 
93 
94 template<int PolySize>
96 :
97  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(is),
98  logActive_(false),
99  logCoeff_(0.0)
100 {}
101 
102 
103 template<int PolySize>
105 :
106  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
107  logActive_(false),
108  logCoeff_(0.0)
109 {
110  word isName(is);
111 
112  if (isName != name)
113  {
115  (
116  "Polynomial<PolySize>::Polynomial(const word&, Istream&)"
117  ) << "Expected polynomial name " << name << " but read " << isName
118  << nl << exit(FatalError);
119  }
120 
121  VectorSpace<Polynomial<PolySize>, scalar, PolySize>::
122  operator=(VectorSpace<Polynomial<PolySize>, scalar, PolySize>(is));
123 
124  if (this->size() == 0)
125  {
127  (
128  "Polynomial<PolySize>::Polynomial(const word&, Istream&)"
129  ) << "Polynomial coefficients for entry " << isName
130  << " are invalid (empty)" << nl << exit(FatalError);
131  }
132 }
133 
134 
135 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
136 
137 template<int PolySize>
139 {
140  return logActive_;
141 }
142 
143 
144 template<int PolySize>
146 {
147  return logCoeff_;
148 }
149 
150 
151 template<int PolySize>
152 Foam::scalar Foam::Polynomial<PolySize>::value(const scalar x) const
153 {
154  scalar val = this->v_[0];
155 
156  // avoid costly pow() in calculation
157  scalar powX = 1;
158  for (label i=1; i<PolySize; ++i)
159  {
160  powX *= x;
161  val += this->v_[i]*powX;
162  }
163 
164  if (logActive_)
165  {
166  val += logCoeff_*log(x);
167  }
168 
169  return val;
170 }
171 
172 
173 template<int PolySize>
174 Foam::scalar Foam::Polynomial<PolySize>::derivative(const scalar x) const
175 {
176  scalar deriv = 0;
177 
178  if (PolySize > 1)
179  {
180  // avoid costly pow() in calculation
181  deriv += this->v_[1];
182 
183  scalar powX = 1;
184  for (label i=2; i<PolySize; ++i)
185  {
186  powX *= x;
187  deriv += i*this->v_[i]*powX;
188  }
189  }
190 
191  if (logActive_)
192  {
193  deriv += logCoeff_/x;
194  }
195 
196  return deriv;
197 }
198 
199 
200 template<int PolySize>
202 (
203  const scalar x1,
204  const scalar x2
205 ) const
206 {
207  // avoid costly pow() in calculation
208  scalar powX1 = x1;
209  scalar powX2 = x2;
210 
211  scalar integ = this->v_[0]*(powX2 - powX1);
212  for (label i=1; i<PolySize; ++i)
213  {
214  powX1 *= x1;
215  powX2 *= x2;
216  integ += this->v_[i]/(i + 1)*(powX2 - powX1);
217  }
218 
219  if (logActive_)
220  {
221  integ += logCoeff_*((x2*log(x2) - x2) - (x1*log(x1) - x1));
222  }
223 
224  return integ;
225 }
226 
227 
228 template<int PolySize>
230 Foam::Polynomial<PolySize>::integral(const scalar intConstant) const
231 {
232  intPolyType newCoeffs;
233 
234  newCoeffs[0] = intConstant;
235  forAll(*this, i)
236  {
237  newCoeffs[i+1] = this->v_[i]/(i + 1);
238  }
239 
240  return newCoeffs;
241 }
242 
243 
244 template<int PolySize>
246 Foam::Polynomial<PolySize>::integralMinus1(const scalar intConstant) const
247 {
248  polyType newCoeffs;
249 
250  if (this->v_[0] > VSMALL)
251  {
252  newCoeffs.logActive_ = true;
253  newCoeffs.logCoeff_ = this->v_[0];
254  }
255 
256  newCoeffs[0] = intConstant;
257  for (label i=1; i<PolySize; ++i)
258  {
259  newCoeffs[i] = this->v_[i]/i;
260  }
261 
262  return newCoeffs;
263 }
264 
265 
266 // ************************************************************************* //
label size() const
Return the number of elements in the VectorSpace = nCmpt.
Polynomial templated on size (order):
Definition: Polynomial.H:63
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
A class for handling words, derived from string.
Definition: word.H:59
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
scalar logCoeff() const
Return the log coefficient.
Definition: Polynomial.C:145
polyType integralMinus1(const scalar intConstant=0.0) const
Return integral coefficients when lowest order is -1.
Definition: Polynomial.C:246
bool logActive() const
Return true if the log term is active.
Definition: Polynomial.C:138
Polynomial()
Construct null, with all coefficients = 0.0.
Definition: Polynomial.C:31
dimensionedScalar log(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:260
scalar v_[nCmpt]
The components of this vector space.
Definition: VectorSpace.H:95
#define forAll(list, i)
Definition: UList.H:421
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
scalar value(const scalar x) const
Return polynomial value.
Definition: Polynomial.C:152
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
scalar integral(const scalar x1, const scalar x2) const
Return integral between two values.
Definition: Polynomial.C:202
error FatalError
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
Templated vector space.
Definition: VectorSpace.H:52
scalar derivative(const scalar x) const
Return derivative of the polynomial at the given x.
Definition: Polynomial.C:174