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-2015 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  << "Size mismatch: Needed " << PolySize
81  << " but given " << coeffs.size()
82  << nl << exit(FatalError);
83  }
84 
85  for (int i = 0; i < PolySize; ++i)
86  {
87  this->v_[i] = coeffs[i];
88  }
89 }
90 
91 
92 template<int PolySize>
94 :
95  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(is),
96  logActive_(false),
97  logCoeff_(0.0)
98 {}
99 
100 
101 template<int PolySize>
103 :
104  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
105  logActive_(false),
106  logCoeff_(0.0)
107 {
108  word isName(is);
109 
110  if (isName != name)
111  {
113  << "Expected polynomial name " << name << " but read " << isName
114  << nl << exit(FatalError);
115  }
116 
117  VectorSpace<Polynomial<PolySize>, scalar, PolySize>::
118  operator=(VectorSpace<Polynomial<PolySize>, scalar, PolySize>(is));
119 
120  if (this->size() == 0)
121  {
123  << "Polynomial coefficients for entry " << isName
124  << " are invalid (empty)" << nl << exit(FatalError);
125  }
126 }
127 
128 
129 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
130 
131 template<int PolySize>
133 {
134  return logActive_;
135 }
136 
137 
138 template<int PolySize>
140 {
141  return logCoeff_;
142 }
143 
144 
145 template<int PolySize>
146 Foam::scalar Foam::Polynomial<PolySize>::value(const scalar x) const
147 {
148  scalar val = this->v_[0];
149 
150  // avoid costly pow() in calculation
151  scalar powX = 1;
152  for (label i=1; i<PolySize; ++i)
153  {
154  powX *= x;
155  val += this->v_[i]*powX;
156  }
157 
158  if (logActive_)
159  {
160  val += logCoeff_*log(x);
161  }
162 
163  return val;
164 }
165 
166 
167 template<int PolySize>
168 Foam::scalar Foam::Polynomial<PolySize>::derivative(const scalar x) const
169 {
170  scalar deriv = 0;
171 
172  if (PolySize > 1)
173  {
174  // avoid costly pow() in calculation
175  deriv += this->v_[1];
176 
177  scalar powX = 1;
178  for (label i=2; i<PolySize; ++i)
179  {
180  powX *= x;
181  deriv += i*this->v_[i]*powX;
182  }
183  }
184 
185  if (logActive_)
186  {
187  deriv += logCoeff_/x;
188  }
189 
190  return deriv;
191 }
192 
193 
194 template<int PolySize>
196 (
197  const scalar x1,
198  const scalar x2
199 ) const
200 {
201  // avoid costly pow() in calculation
202  scalar powX1 = x1;
203  scalar powX2 = x2;
204 
205  scalar integ = this->v_[0]*(powX2 - powX1);
206  for (label i=1; i<PolySize; ++i)
207  {
208  powX1 *= x1;
209  powX2 *= x2;
210  integ += this->v_[i]/(i + 1)*(powX2 - powX1);
211  }
212 
213  if (logActive_)
214  {
215  integ += logCoeff_*((x2*log(x2) - x2) - (x1*log(x1) - x1));
216  }
217 
218  return integ;
219 }
220 
221 
222 template<int PolySize>
224 Foam::Polynomial<PolySize>::integral(const scalar intConstant) const
225 {
226  intPolyType newCoeffs;
227 
228  newCoeffs[0] = intConstant;
229  forAll(*this, i)
230  {
231  newCoeffs[i+1] = this->v_[i]/(i + 1);
232  }
233 
234  return newCoeffs;
235 }
236 
237 
238 template<int PolySize>
240 Foam::Polynomial<PolySize>::integralMinus1(const scalar intConstant) const
241 {
242  polyType newCoeffs;
243 
244  if (this->v_[0] > VSMALL)
245  {
246  newCoeffs.logActive_ = true;
247  newCoeffs.logCoeff_ = this->v_[0];
248  }
249 
250  newCoeffs[0] = intConstant;
251  for (label i=1; i<PolySize; ++i)
252  {
253  newCoeffs[i] = this->v_[i]/i;
254  }
255 
256  return newCoeffs;
257 }
258 
259 
260 // ************************************************************************* //
scalar logCoeff() const
Return the log coefficient.
Definition: Polynomial.C:139
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
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
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 derivative(const scalar x) const
Return derivative of the polynomial at the given x.
Definition: Polynomial.C:168
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
bool logActive() const
Return true if the log term is active.
Definition: Polynomial.C:132
static direction size()
Return the number of elements in the VectorSpace = Ncmpts.
static const char nl
Definition: Ostream.H:262
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
Polynomial templated on size (order):
Definition: Polynomial.H:63
friend Ostream & operator(Ostream &, const Polynomial &)
Ostream Operator.
scalar value(const scalar x) const
Return polynomial value.
Definition: Polynomial.C:146
polyType integralMinus1(const scalar intConstant=0.0) const
Return integral coefficients when lowest order is -1.
Definition: Polynomial.C:240
scalar v_[Ncmpts]
The components of this vector space.
Definition: VectorSpace.H:81
scalar integral(const scalar x1, const scalar x2) const
Return integral between two values.
Definition: Polynomial.C:196