diff --git a/applications/test/Polynomial/PolynomialTest.C b/applications/test/Polynomial/PolynomialTest.C index d3fa03e595..fbfbcc2000 100644 --- a/applications/test/Polynomial/PolynomialTest.C +++ b/applications/test/Polynomial/PolynomialTest.C @@ -29,29 +29,30 @@ Description \*---------------------------------------------------------------------------*/ -#include "IFstream.H" +#include "IStringStream.H" #include "Polynomial.H" #include "Random.H" #include "cpuTime.H" using namespace Foam; -const int nCoeffs = 8; -const scalar coeff[] = { 0.11, 0.45, -0.94, 1.58, 2.58, 0.08, 3.15, -4.78 }; +const int PolySize = 8; +const scalar coeff[] = { 0.11, 0.45, -0.94, 1.58, -2.58, 0.08, 3.15, -4.78 }; +const char* polyDef = "testPoly (0.11 0.45 -0.94 1.58 -2.58 0.08 3.15 -4.78)"; scalar polyValue(const scalar x) { // Hard-coded polynomial 8 coeff (7th order) return - 0.11 - + 0.45*x - - 0.94*sqr(x) - + 1.58*pow3(x) - - 2.58*pow4(x) - + 0.08*pow5(x) - + 3.15*pow6(x) - - 4.78*x*pow6(x); + coeff[0] + + coeff[1]*x + + coeff[2]*sqr(x) + + coeff[3]*pow3(x) + + coeff[4]*pow4(x) + + coeff[5]*pow5(x) + + coeff[6]*pow6(x) + + coeff[7]*x*pow6(x); } @@ -59,14 +60,14 @@ scalar intPolyValue(const scalar x) { // Hard-coded integrated form of above polynomial return - 0.11*x - + 0.45/2.0*sqr(x) - - 0.94/3.0*pow3(x) - + 1.58/4.0*pow4(x) - - 2.58/5.0*pow5(x) - + 0.08/6.0*pow6(x) - + 3.15/7.0*x*pow6(x) - - 4.78/8.0*x*x*pow6(x); + coeff[0]*x + + coeff[1]/2.0*sqr(x) + + coeff[2]/3.0*pow3(x) + + coeff[3]/4.0*pow4(x) + + coeff[4]/5.0*pow5(x) + + coeff[5]/6.0*pow6(x) + + coeff[6]/7.0*x*pow6(x) + + coeff[7]/8.0*x*x*pow6(x); } @@ -75,7 +76,7 @@ scalar polyValue1(const scalar x) // "normal" evaluation using pow() scalar value = coeff[0]; - for (int i=1; i < nCoeffs; ++i) + for (int i=1; i < PolySize; ++i) { value += coeff[i]*pow(x, i); } @@ -84,15 +85,16 @@ scalar polyValue1(const scalar x) } -// calculation avoiding pow() scalar polyValue2(const scalar x) { + // calculation avoiding pow() scalar value = coeff[0]; scalar powX = x; - for (int i=1; i < nCoeffs; ++i, powX *= x) + for (int i=1; i < PolySize; ++i) { value += coeff[i] * powX; + powX *= x; } return value; @@ -107,9 +109,10 @@ int main(int argc, char *argv[]) const label nIters = 1000; scalar sum = 0.0; - IFstream is("polyTestInput"); + Info<< "null poly = " << (Polynomial<8>()) << nl << endl; - Polynomial<8> poly("testPoly", is); + // Polynomial<8> poly("testPoly", IStringStream(polyDef)()); + Polynomial<8> poly(coeff); Polynomial<9> intPoly(poly.integrate(0.0)); Info<< "poly = " << poly << endl; diff --git a/applications/test/Polynomial/polyTestInput b/applications/test/Polynomial/polyTestInput deleted file mode 100644 index 4ba8f65e51..0000000000 --- a/applications/test/Polynomial/polyTestInput +++ /dev/null @@ -1,11 +0,0 @@ - testPoly - ( - 0.11 - 0.45 - -0.94 - 1.58 - -2.58 - 0.08 - 3.15 - -4.78 - ) diff --git a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C index 903bdffe2b..cff3744644 100644 --- a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C +++ b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C @@ -33,7 +33,50 @@ Foam::Polynomial::Polynomial() VectorSpace, scalar, PolySize>(), logActive_(false), logCoeff_(0.0) -{} +{ + for (int i = 0; i < PolySize; ++i) + { + this->v_[i] = 0.0; + } +} + + +template +Foam::Polynomial::Polynomial(const scalar coeffs[PolySize]) +: + VectorSpace, scalar, PolySize>(), + logActive_(false), + logCoeff_(0.0) +{ + for (int i=0; iv_[i] = coeffs[i]; + } +} + + +template +Foam::Polynomial::Polynomial(const UList& coeffs) +: + VectorSpace, scalar, PolySize>(), + logActive_(false), + logCoeff_(0.0) +{ + if (coeffs.size() != PolySize) + { + FatalErrorIn + ( + "Polynomial::Polynomial(const UList&)" + ) << "Size mismatch: Needed " << PolySize + << " but got " << coeffs.size() + << nl << exit(FatalError); + } + + for (int i = 0; i < PolySize; ++i) + { + this->v_[i] = coeffs[i]; + } +} template @@ -173,13 +216,9 @@ Foam::Polynomial::integrateMinus1(const scalar intConstant) } newCoeffs[0] = intConstant; - - if (PolySize > 0) + for (label i=1; iv_[i]/i; - } + newCoeffs[i] = this->v_[i]/i; } return newCoeffs; diff --git a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H index cbfe045803..f2449d1835 100644 --- a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H +++ b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H @@ -29,7 +29,7 @@ Description poly = logCoeff*log(x) + sum(coeff_[i]*x^i) - where 0 <= i <= n + where 0 \<= i \<= n - integer powers, starting at zero - evaluate(x) to evaluate the poly for a given value @@ -51,6 +51,7 @@ SourceFiles #include "scalar.H" #include "Ostream.H" #include "VectorSpace.H" +#include "StaticAssert.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -79,6 +80,9 @@ class Polynomial : public VectorSpace, scalar, PolySize> { + //- Size must be positive (non-zero) + StaticAssert(PolySize > 0); + // Private data //- Include the log term? - only activated using integrateMinus1() @@ -97,9 +101,15 @@ public: // Constructors - //- Construct null + //- Construct null, with all coefficients = 0.0 Polynomial(); + //- Construct from C-array of coefficients + explicit Polynomial(const scalar coeffs[PolySize]); + + //- Construct from a list of coefficients + explicit Polynomial(const UList& coeffs); + //- Construct from name and Istream Polynomial(const word& name, Istream& is);