ENH: add StaticAssert to Polynomial - positive number of terms only

ENH: allow construct from UList of coefficients, from C-arrays
     avoid uninitialized values for null constructor

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;

applications/test/Polynomial/polyTestInput

@@ -1,11 +0,0 @@
-    testPoly
-    (
-        0.11
-        0.45
-       -0.94
-        1.58
-       -2.58
-        0.08
-        3.15
-       -4.78
-    )

src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C

@@ -33,7 +33,50 @@ Foam::Polynomial<PolySize>::Polynomial()
     VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
     logActive_(false),
     logCoeff_(0.0)
-{}
+{
+    for (int i = 0; i < PolySize; ++i)
+    {
+        this->v_[i] = 0.0;
+    }
+}
+
+
+template<int PolySize>
+Foam::Polynomial<PolySize>::Polynomial(const scalar coeffs[PolySize])
+:
+    VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
+    logActive_(false),
+    logCoeff_(0.0)
+{
+    for (int i=0; i<PolySize; i++)
+    {
+        this->v_[i] = coeffs[i];
+    }
+}
+
+
+template<int PolySize>
+Foam::Polynomial<PolySize>::Polynomial(const UList<scalar>& coeffs)
+:
+    VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
+    logActive_(false),
+    logCoeff_(0.0)
+{
+    if (coeffs.size() != PolySize)
+    {
+        FatalErrorIn
+        (
+            "Polynomial<PolySize>::Polynomial(const UList<scalar>&)"
+        )   << "Size mismatch: Needed " << PolySize
+            << " but got " << coeffs.size()
+            << nl << exit(FatalError);
+    }
+
+    for (int i = 0; i < PolySize; ++i)
+    {
+        this->v_[i] = coeffs[i];
+    }
+}
 
 
 template<int PolySize>
@@ -173,13 +216,9 @@ Foam::Polynomial<PolySize>::integrateMinus1(const scalar intConstant)
     }
 
     newCoeffs[0] = intConstant;
-
-    if (PolySize > 0)
+    for (label i=1; i<PolySize; ++i)
     {
-        for (label i=1; i<PolySize; i++)
-        {
-            newCoeffs[i] = this->v_[i]/i;
-        }
+        newCoeffs[i] = this->v_[i]/i;
     }
 
     return newCoeffs;

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<Polynomial<PolySize>, 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<scalar>& coeffs);
+
         //- Construct from name and Istream
         Polynomial(const word& name, Istream& is);