Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "Polynomial.H"
00027
00028
00029
00030 template<int PolySize>
00031 Foam::Polynomial<PolySize>::Polynomial()
00032 :
00033 VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
00034 logActive_(false),
00035 logCoeff_(0.0)
00036 {}
00037
00038
00039 template<int PolySize>
00040 Foam::Polynomial<PolySize>::Polynomial(const word& name, Istream& is)
00041 :
00042 VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
00043 logActive_(false),
00044 logCoeff_(0.0)
00045 {
00046 word isName(is);
00047
00048 if (isName != name)
00049 {
00050 FatalErrorIn
00051 (
00052 "Polynomial<PolySize>::Polynomial(const word&, Istream&)"
00053 ) << "Expected polynomial name " << name << " but read " << isName
00054 << nl << exit(FatalError);
00055 }
00056
00057 VectorSpace<Polynomial<PolySize>, scalar, PolySize>::
00058 operator=(VectorSpace<Polynomial<PolySize>, scalar, PolySize>(is));
00059
00060 if (this->size() == 0)
00061 {
00062 FatalErrorIn
00063 (
00064 "Polynomial<PolySize>::Polynomial(const word&, Istream&)"
00065 ) << "Polynomial coefficients for entry " << isName
00066 << " are invalid (empty)" << nl << exit(FatalError);
00067 }
00068 }
00069
00070
00071 template<int PolySize>
00072 Foam::Polynomial<PolySize>::Polynomial
00073 (
00074 const Polynomial<PolySize>& poly
00075 )
00076 :
00077 VectorSpace<Polynomial<PolySize>, scalar, PolySize>(poly),
00078 logActive_(poly.logActive_),
00079 logCoeff_(poly.logCoeff_)
00080 {}
00081
00082
00083
00084
00085 template<int PolySize>
00086 bool& Foam::Polynomial<PolySize>::logActive()
00087 {
00088 return logActive_;
00089 }
00090
00091
00092 template<int PolySize>
00093 Foam::scalar& Foam::Polynomial<PolySize>::logCoeff()
00094 {
00095 return logCoeff_;
00096 }
00097
00098
00099 template<int PolySize>
00100 Foam::scalar Foam::Polynomial<PolySize>::evaluate(const scalar x) const
00101 {
00102 scalar y = this->v_[0];
00103
00104 for (label i=1; i<PolySize; i++)
00105 {
00106 y += this->v_[i]*pow(x, i);
00107 }
00108
00109 if (logActive_)
00110 {
00111 y += logCoeff_*log(x);
00112 }
00113
00114 return y;
00115 }
00116
00117
00118 template<int PolySize>
00119 Foam::scalar Foam::Polynomial<PolySize>::integrateLimits
00120 (
00121 const scalar x1,
00122 const scalar x2
00123 ) const
00124 {
00125 if (logActive_)
00126 {
00127 FatalErrorIn
00128 (
00129 "scalar Polynomial<PolySize>::integrateLimits"
00130 "("
00131 "const scalar, "
00132 "const scalar"
00133 ") const"
00134 ) << "Cannot integrate polynomial with logarithmic coefficients"
00135 << nl << abort(FatalError);
00136 }
00137
00138 intPolyType poly = this->integrate();
00139
00140 return poly.evaluate(x2) - poly.evaluate(x1);
00141 }
00142
00143
00144 template<int PolySize>
00145 typename Foam::Polynomial<PolySize>::intPolyType
00146 Foam::Polynomial<PolySize>::integrate(const scalar intConstant)
00147 {
00148 intPolyType newCoeffs;
00149
00150 newCoeffs[0] = intConstant;
00151 forAll(*this, i)
00152 {
00153 newCoeffs[i + 1] = this->v_[i]/(i + 1);
00154 }
00155
00156 return newCoeffs;
00157 }
00158
00159
00160 template<int PolySize>
00161 typename Foam::Polynomial<PolySize>::polyType
00162 Foam::Polynomial<PolySize>::integrateMinus1(const scalar intConstant)
00163 {
00164 polyType newCoeffs;
00165
00166 if (this->v_[0] > VSMALL)
00167 {
00168 newCoeffs.logActive() = true;
00169 newCoeffs.logCoeff() = this->v_[0];
00170 }
00171
00172 newCoeffs[0] = intConstant;
00173
00174 if (PolySize > 0)
00175 {
00176 for (label i=1; i<PolySize; i++)
00177 {
00178 newCoeffs[i] = this->v_[i]/i;
00179 }
00180 }
00181
00182 return newCoeffs;
00183 }
00184
00185
00186