FreeFOAM The Cross-Platform CFD Toolkit
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) 2008-2010 OpenCFD Ltd.
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 
38 
39 template<int PolySize>
41 :
42  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
43  logActive_(false),
44  logCoeff_(0.0)
45 {
46  word isName(is);
47 
48  if (isName != name)
49  {
51  (
52  "Polynomial<PolySize>::Polynomial(const word&, Istream&)"
53  ) << "Expected polynomial name " << name << " but read " << isName
54  << nl << exit(FatalError);
55  }
56 
57  VectorSpace<Polynomial<PolySize>, scalar, PolySize>::
58  operator=(VectorSpace<Polynomial<PolySize>, scalar, PolySize>(is));
59 
60  if (this->size() == 0)
61  {
63  (
64  "Polynomial<PolySize>::Polynomial(const word&, Istream&)"
65  ) << "Polynomial coefficients for entry " << isName
66  << " are invalid (empty)" << nl << exit(FatalError);
67  }
68 }
69 
70 
71 template<int PolySize>
73 (
74  const Polynomial<PolySize>& poly
75 )
76 :
77  VectorSpace<Polynomial<PolySize>, scalar, PolySize>(poly),
78  logActive_(poly.logActive_),
79  logCoeff_(poly.logCoeff_)
80 {}
81 
82 
83 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
84 
85 template<int PolySize>
87 {
88  return logActive_;
89 }
90 
91 
92 template<int PolySize>
94 {
95  return logCoeff_;
96 }
97 
98 
99 template<int PolySize>
100 Foam::scalar Foam::Polynomial<PolySize>::evaluate(const scalar x) const
101 {
102  scalar y = this->v_[0];
103 
104  for (label i=1; i<PolySize; i++)
105  {
106  y += this->v_[i]*pow(x, i);
107  }
108 
109  if (logActive_)
110  {
111  y += logCoeff_*log(x);
112  }
113 
114  return y;
115 }
116 
117 
118 template<int PolySize>
120 (
121  const scalar x1,
122  const scalar x2
123 ) const
124 {
125  if (logActive_)
126  {
128  (
129  "scalar Polynomial<PolySize>::integrateLimits"
130  "("
131  "const scalar, "
132  "const scalar"
133  ") const"
134  ) << "Cannot integrate polynomial with logarithmic coefficients"
135  << nl << abort(FatalError);
136  }
137 
138  intPolyType poly = this->integrate();
139 
140  return poly.evaluate(x2) - poly.evaluate(x1);
141 }
142 
143 
144 template<int PolySize>
146 Foam::Polynomial<PolySize>::integrate(const scalar intConstant)
147 {
148  intPolyType newCoeffs;
149 
150  newCoeffs[0] = intConstant;
151  forAll(*this, i)
152  {
153  newCoeffs[i + 1] = this->v_[i]/(i + 1);
154  }
155 
156  return newCoeffs;
157 }
158 
159 
160 template<int PolySize>
163 {
164  polyType newCoeffs;
165 
166  if (this->v_[0] > VSMALL)
167  {
168  newCoeffs.logActive() = true;
169  newCoeffs.logCoeff() = this->v_[0];
170  }
171 
172  newCoeffs[0] = intConstant;
173 
174  if (PolySize > 0)
175  {
176  for (label i=1; i<PolySize; i++)
177  {
178  newCoeffs[i] = this->v_[i]/i;
179  }
180  }
181 
182  return newCoeffs;
183 }
184 
185 
186 // ************************ vim: set sw=4 sts=4 et: ************************ //