diff --git a/applications/test/Polynomial/Test-Polynomial.C b/applications/test/Polynomial/Test-Polynomial.C
index 70fbf4ff57..2ddf858ffe 100644
--- a/applications/test/Polynomial/Test-Polynomial.C
+++ b/applications/test/Polynomial/Test-Polynomial.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,24 +25,29 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    PolynomialTest
+    Test-Polynomial
 
 Description
     Test application for the templated Polynomial class
 
 \*---------------------------------------------------------------------------*/
 
-#include "StringStream.H"
 #include "Polynomial.H"
+#include "FixedList.H"
 #include "polynomialFunction.H"
+#include "ITstream.H"
+#include "OTstream.H"
 #include "Random.H"
 #include "cpuTime.H"
 
 using namespace Foam;
 
-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 = "(0.11 0.45 -0.94 1.58 -2.58 0.08 3.15 -4.78)";
+std::initializer_list<scalar> coeffs
+{
+    0.11, 0.45, -0.94, 1.58, -2.58, 0.08, 3.15, -4.78
+};
+
+const FixedList<scalar, 8> coeff(coeffs);
 
 
 scalar polyValue(const scalar x)
@@ -76,10 +82,10 @@ scalar intPolyValue(const scalar x)
 
 scalar polyValue1(const scalar x)
 {
-    // "normal" evaluation using pow()
+    // Naive evaluation using pow()
     scalar value = coeff[0];
 
-    for (int i=1; i < PolySize; ++i)
+    for (label i=1; i < coeff.size(); ++i)
     {
         value += coeff[i]*pow(x, i);
     }
@@ -90,11 +96,11 @@ scalar polyValue1(const scalar x)
 
 scalar polyValue2(const scalar x)
 {
-    // calculation avoiding pow()
+    // Calculation avoiding pow()
     scalar value = coeff[0];
 
     scalar powX = x;
-    for (int i=1; i < PolySize; ++i)
+    for (label i=1; i < coeff.size(); ++i)
     {
         value += coeff[i] * powX;
         powX *= x;
@@ -108,19 +114,27 @@ scalar polyValue2(const scalar x)
 
 int main(int argc, char *argv[])
 {
-    const label n = 10000;
-    const label nIters = 1000;
-    scalar sum = 0.0;
+    constexpr label n = 10000;
+    constexpr label nIters = 1000;
+    scalar sum = 0;
 
-    Info<< "null poly = " << (Polynomial<8>()) << nl
-        << "null poly = " << (polynomialFunction(8)) << nl
+    Info<< "null poly = " << (Polynomial<8>{}) << nl
+        << "null poly = " << (polynomialFunction{8}) << nl
         << endl;
 
-    Polynomial<8> poly(coeff);
-    Polynomial<9> intPoly(poly.integral(0.0));
+    Polynomial<8> poly{coeffs};
+    Polynomial<9> intPoly{poly.integral(0)};
 
-    IStringStream is(polyDef);
-    polynomialFunction polyfunc(is);
+    polynomialFunction polyfunc;
+
+    // Could profit from a bi-directional stream
+    {
+        OTstream os;
+        os << poly;
+
+        ITstream is("input", std::move(os.tokens()));
+        is >> polyfunc;
+    }
 
     Info<< "poly = " << poly << nl
         << "intPoly = " << intPoly << nl
@@ -152,7 +166,7 @@ int main(int argc, char *argv[])
     Info<< "2.5*poly = " << polyCopy << nl << endl;
 
     Random rnd(123456);
-    for (int i=0; i<10; i++)
+    for (label i=0; i<10; ++i)
     {
         scalar x = rnd.sample01<scalar>()*100;
 
@@ -162,18 +176,18 @@ int main(int argc, char *argv[])
         scalar pxTest = poly.value(x);
         scalar ipxTest = intPoly.value(x);
 
-        Info<<"\nx = " << x << endl;
-        Info<< "    px, pxTest = " << px << ", " << pxTest << endl;
-        Info<< "    ipx, ipxTest = " << ipx << ", " << ipxTest << endl;
+        Info<<"\nx = " << x << nl
+            << "    px, pxTest = " << px << ", " << pxTest << nl
+            << "    ipx, ipxTest = " << ipx << ", " << ipxTest << nl;
 
         if (mag(px - pxTest) > SMALL)
         {
-            Info<< "    *** WARNING: px != pxTest: " << px - pxTest << endl;
+            Info<< "    *** WARNING: px != pxTest: " << px - pxTest << nl;
         }
 
         if (mag(ipx - ipxTest) > SMALL)
         {
-            Info<< "    *** WARNING: ipx != ipxTest: " << ipx - ipxTest << endl;
+            Info<< "    *** WARNING: ipx != ipxTest: " << ipx - ipxTest << nl;
         }
 
         Info<< endl;
@@ -188,7 +202,7 @@ int main(int argc, char *argv[])
 
     cpuTime timer;
 
-    for (int loop = 0; loop < n; ++loop)
+    for (label loop = 0; loop < n; ++loop)
     {
         sum = 0.0;
         for (label iter = 0; iter < nIters; ++iter)
@@ -199,7 +213,7 @@ int main(int argc, char *argv[])
     Info<< "value:        " << sum
         << " in " << timer.cpuTimeIncrement() << " s\n";
 
-    for (int loop = 0; loop < n; ++loop)
+    for (label loop = 0; loop < n; ++loop)
     {
         sum = 0.0;
         for (label iter = 0; iter < nIters; ++iter)
@@ -211,7 +225,7 @@ int main(int argc, char *argv[])
         << " in " << timer.cpuTimeIncrement() << " s\n";
 
 
-    for (int loop = 0; loop < n; ++loop)
+    for (label loop = 0; loop < n; ++loop)
     {
         sum = 0.0;
         for (label iter = 0; iter < nIters; ++iter)
@@ -223,7 +237,7 @@ int main(int argc, char *argv[])
         << " in " << timer.cpuTimeIncrement() << " s\n";
 
 
-    for (int loop = 0; loop < n; ++loop)
+    for (label loop = 0; loop < n; ++loop)
     {
         sum = 0.0;
         for (label iter = 0; iter < nIters; ++iter)
@@ -234,7 +248,7 @@ int main(int argc, char *argv[])
     Info<< "hard-coded 1: " << sum
         << " in " << timer.cpuTimeIncrement() << " s\n";
 
-    for (int loop = 0; loop < n; ++loop)
+    for (label loop = 0; loop < n; ++loop)
     {
         sum = 0.0;
         for (label iter = 0; iter < nIters; ++iter)
diff --git a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C
index 97576ba7b4..8570e7514a 100644
--- a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C
+++ b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,14 +32,33 @@ License
 
 template<int PolySize>
 Foam::Polynomial<PolySize>::Polynomial()
+:
+    VectorSpace<Polynomial<PolySize>, scalar, PolySize>(Zero),
+    logActive_(false),
+    logCoeff_(0)
+{}
+
+
+template<int PolySize>
+Foam::Polynomial<PolySize>::Polynomial(std::initializer_list<scalar> coeffs)
 :
     VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
     logActive_(false),
     logCoeff_(0)
 {
-    for (int i = 0; i < PolySize; ++i)
+    if (coeffs.size() != PolySize)
     {
-        this->v_[i] = 0;
+        FatalErrorInFunction
+            << "Size mismatch: Needed " << PolySize
+            << " but given " << label(coeffs.size())
+            << nl << exit(FatalError);
+    }
+
+    auto iter = coeffs.begin();
+    for (int i=0; i<PolySize; ++i)
+    {
+        this->v_[i] = *iter;
+        ++iter;
     }
 }
 
@@ -50,7 +70,7 @@ Foam::Polynomial<PolySize>::Polynomial(const scalar coeffs[PolySize])
     logActive_(false),
     logCoeff_(0)
 {
-    for (int i=0; i<PolySize; i++)
+    for (int i=0; i<PolySize; ++i)
     {
         this->v_[i] = coeffs[i];
     }
@@ -95,7 +115,7 @@ Foam::Polynomial<PolySize>::Polynomial(const word& name, Istream& is)
     logActive_(false),
     logCoeff_(0)
 {
-    word isName(is);
+    const word isName(is);
 
     if (isName != name)
     {
@@ -104,15 +124,11 @@ Foam::Polynomial<PolySize>::Polynomial(const word& name, Istream& is)
             << nl << exit(FatalError);
     }
 
-    VectorSpace<Polynomial<PolySize>, scalar, PolySize>::
-        operator=(VectorSpace<Polynomial<PolySize>, scalar, PolySize>(is));
-
-    if (this->size() == 0)
-    {
-        FatalErrorInFunction
-            << "Polynomial coefficients for entry " << isName
-            << " are invalid (empty)" << nl << exit(FatalError);
-    }
+    is >>
+        static_cast
+        <
+            VectorSpace<Polynomial<PolySize>, scalar, PolySize>&
+        >(*this);
 }
 
 
@@ -137,12 +153,12 @@ Foam::scalar Foam::Polynomial<PolySize>::value(const scalar x) const
 {
     scalar val = this->v_[0];
 
-    // avoid costly pow() in calculation
-    scalar powX = 1;
+    // Avoid costly pow() in calculation
+    scalar powX = x;
     for (label i=1; i<PolySize; ++i)
     {
-        powX *= x;
         val += this->v_[i]*powX;
+        powX *= x;
     }
 
     if (logActive_)
@@ -161,14 +177,14 @@ Foam::scalar Foam::Polynomial<PolySize>::derivative(const scalar x) const
 
     if (PolySize > 1)
     {
-        // avoid costly pow() in calculation
+        // Avoid costly pow() in calculation
         deriv += this->v_[1];
 
-        scalar powX = 1;
+        scalar powX = x;
         for (label i=2; i<PolySize; ++i)
         {
-            powX *= x;
             deriv += i*this->v_[i]*powX;
+            powX *= x;
         }
     }
 
@@ -188,7 +204,7 @@ Foam::scalar Foam::Polynomial<PolySize>::integral
     const scalar x2
 ) const
 {
-    // avoid costly pow() in calculation
+    // Avoid costly pow() in calculation
     scalar powX1 = x1;
     scalar powX2 = x2;
 
diff --git a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H
index 2acea6d55c..3e0bc47574 100644
--- a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H
+++ b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H
@@ -57,6 +57,7 @@ SourceFiles
 #include "scalar.H"
 #include "Ostream.H"
 #include "VectorSpace.H"
+#include <initializer_list>
 #include <type_traits>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -67,9 +68,6 @@ namespace Foam
 // Forward Declarations
 template<int PolySize> class Polynomial;
 
-template<int PolySize>
-Ostream& operator<<( Ostream&, const Polynomial<PolySize>&);
-
 /*---------------------------------------------------------------------------*\
                         Class Polynomial Declaration
 \*---------------------------------------------------------------------------*/
@@ -106,6 +104,9 @@ public:
         //- Default construct, with all coefficients = 0
         Polynomial();
 
+        //- Construct from an initializer list of coefficients
+        Polynomial(std::initializer_list<scalar> coeffs);
+
         //- Construct from C-array of coefficients
         explicit Polynomial(const scalar coeffs[PolySize]);
 
@@ -150,12 +151,7 @@ public:
             polyType integralMinus1(const scalar intConstant = 0.0) const;
 
 
-    //- Ostream Operator
-    friend Ostream& operator<< <PolySize>
-    (
-        Ostream&,
-        const Polynomial&
-    );
+    // IOstream Operators - uses VectorSpace operators
 };
 
 
@@ -167,7 +163,6 @@ public:
 
 #ifdef NoRepository
     #include "Polynomial.C"
-    #include "PolynomialIO.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/primitives/functions/Polynomial/PolynomialIO.C b/src/OpenFOAM/primitives/functions/Polynomial/PolynomialIO.C
index 2a75df9ec8..7cec1b4541 100644
--- a/src/OpenFOAM/primitives/functions/Polynomial/PolynomialIO.C
+++ b/src/OpenFOAM/primitives/functions/Polynomial/PolynomialIO.C
@@ -1,49 +1 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "Polynomial.H"
-
-// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
-
-template<int PolySize>
-Foam::Ostream& Foam::operator<<
-(
-    Ostream& os,
-    const Polynomial<PolySize>& poly
-)
-{
-    return
-        os <<
-        static_cast
-        <
-            VectorSpace<Polynomial<PolySize>, scalar, PolySize>
-        >(poly);
-}
-
-
-// ************************************************************************* //
+#warning File removed - left for old dependency check only
diff --git a/src/OpenFOAM/primitives/functions/Polynomial/polynomialFunction.C b/src/OpenFOAM/primitives/functions/Polynomial/polynomialFunction.C
index 4054f9f321..a6a979e8ab 100644
--- a/src/OpenFOAM/primitives/functions/Polynomial/polynomialFunction.C
+++ b/src/OpenFOAM/primitives/functions/Polynomial/polynomialFunction.C
@@ -39,7 +39,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
-
 Foam::polynomialFunction Foam::polynomialFunction::cloneIntegral
 (
     const polynomialFunction& poly,
@@ -82,13 +81,9 @@ Foam::polynomialFunction Foam::polynomialFunction::cloneIntegralMinus1
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::polynomialFunction::polynomialFunction(const label order)
-:
-    scalarList(order, Zero),
-    logActive_(false),
-    logCoeff_(0)
+void Foam::polynomialFunction::checkSize() const
 {
     if (this->empty())
     {
@@ -99,26 +94,46 @@ Foam::polynomialFunction::polynomialFunction(const label order)
 }
 
 
-Foam::polynomialFunction::polynomialFunction(const polynomialFunction& poly)
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::polynomialFunction::polynomialFunction()
 :
-    scalarList(poly),
-    logActive_(poly.logActive_),
-    logCoeff_(poly.logCoeff_)
+    scalarList(1, Zero),
+    logActive_(false),
+    logCoeff_(0)
 {}
 
 
+Foam::polynomialFunction::polynomialFunction(const label order)
+:
+    scalarList(order, Zero),
+    logActive_(false),
+    logCoeff_(0)
+{
+    checkSize();
+}
+
+
+Foam::polynomialFunction::polynomialFunction
+(
+    const std::initializer_list<scalar> coeffs
+)
+:
+    scalarList(coeffs),
+    logActive_(false),
+    logCoeff_(0)
+{
+    checkSize();
+}
+
+
 Foam::polynomialFunction::polynomialFunction(const UList<scalar>& coeffs)
 :
     scalarList(coeffs),
     logActive_(false),
     logCoeff_(0)
 {
-    if (this->empty())
-    {
-        FatalErrorInFunction
-            << "polynomialFunction coefficients are invalid (empty)"
-            << nl << exit(FatalError);
-    }
+    checkSize();
 }
 
 
@@ -128,12 +143,7 @@ Foam::polynomialFunction::polynomialFunction(Istream& is)
     logActive_(false),
     logCoeff_(0)
 {
-    if (this->empty())
-    {
-        FatalErrorInFunction
-            << "polynomialFunction coefficients are invalid (empty)"
-            << nl << exit(FatalError);
-    }
+    checkSize();
 }
 
 
@@ -156,7 +166,7 @@ Foam::scalar Foam::polynomialFunction::value(const scalar x) const
     const scalarList& coeffs = *this;
     scalar val = coeffs[0];
 
-    // avoid costly pow() in calculation
+    // Avoid costly pow() in calculation
     scalar powX = x;
     for (label i=1; i<coeffs.size(); ++i)
     {
@@ -188,7 +198,7 @@ Foam::scalar Foam::polynomialFunction::integrate
             << nl << abort(FatalError);
     }
 
-    // avoid costly pow() in calculation
+    // Avoid costly pow() in calculation
     scalar powX1 = x1;
     scalar powX2 = x2;
 
@@ -225,21 +235,14 @@ Foam::polynomialFunction::operator+=(const polynomialFunction& poly)
 {
     scalarList& coeffs = *this;
 
-    if (coeffs.size() > poly.size())
+    if (coeffs.size() < poly.size())
     {
-        forAll(poly, i)
-        {
-            coeffs[i] += poly[i];
-        }
+        coeffs.resize(poly.size(), Zero);
     }
-    else
-    {
-        coeffs.setSize(poly.size(), 0.0);
 
-        forAll(coeffs, i)
-        {
-            coeffs[i] += poly[i];
-        }
+    forAll(poly, i)
+    {
+        coeffs[i] += poly[i];
     }
 
     return *this;
@@ -251,21 +254,14 @@ Foam::polynomialFunction::operator-=(const polynomialFunction& poly)
 {
     scalarList& coeffs = *this;
 
-    if (coeffs.size() > poly.size())
+    if (coeffs.size() < poly.size())
     {
-        forAll(poly, i)
-        {
-            coeffs[i] -= poly[i];
-        }
+        coeffs.resize(poly.size(), Zero);
     }
-    else
-    {
-        coeffs.setSize(poly.size(), 0.0);
 
-        forAll(coeffs, i)
-        {
-            coeffs[i] -= poly[i];
-        }
+    forAll(poly, i)
+    {
+        coeffs[i] -= poly[i];
     }
 
     return *this;
@@ -300,6 +296,20 @@ Foam::polynomialFunction::operator/=(const scalar s)
 
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
+Foam::Istream& Foam::operator>>(Istream& is, polynomialFunction& poly)
+{
+    // Log handling may be unreliable
+    poly.logActive_ = false;
+    poly.logCoeff_ = 0;
+
+    is >> static_cast<scalarList&>(poly);
+
+    poly.checkSize();
+
+    return is;
+}
+
+
 Foam::Ostream& Foam::operator<<(Ostream& os, const polynomialFunction& poly)
 {
     // output like VectorSpace
diff --git a/src/OpenFOAM/primitives/functions/Polynomial/polynomialFunction.H b/src/OpenFOAM/primitives/functions/Polynomial/polynomialFunction.H
index a4fcf8691c..411260f4be 100644
--- a/src/OpenFOAM/primitives/functions/Polynomial/polynomialFunction.H
+++ b/src/OpenFOAM/primitives/functions/Polynomial/polynomialFunction.H
@@ -66,7 +66,9 @@ namespace Foam
 
 // Forward Declarations
 class polynomialFunction;
-Ostream& operator<<(Ostream&, const polynomialFunction&);
+
+Istream& operator>>(Istream&, polynomialFunction& poly);
+Ostream& operator<<(Ostream&, const polynomialFunction& poly);
 
 /*---------------------------------------------------------------------------*\
                      Class polynomialFunction Declaration
@@ -103,6 +105,9 @@ class polynomialFunction
             const scalar intConstant = 0
         );
 
+        //- Check size is non-zero or trigger FatalErrot
+        void checkSize() const;
+
         //- No copy assignment
         void operator=(const polynomialFunction&) = delete;
 
@@ -113,13 +118,19 @@ public:
     TypeName("polynomialFunction");
 
 
+    // Generated Methods: copy construct, copy assignment
+
+
     // Constructors
 
-        //- Construct a particular size, with all coefficients as 0
+        //- Default construct as size 1 with coefficient == 0
+        polynomialFunction();
+
+        //- Construct a particular size, with all coefficients == 0
         explicit polynomialFunction(const label order);
 
-        //- Copy constructor
-        polynomialFunction(const polynomialFunction& rhs);
+        //- Construct from an initializer list of coefficients
+        explicit polynomialFunction(std::initializer_list<scalar> coeffs);
 
         //- Construct from a list of coefficients
         explicit polynomialFunction(const UList<scalar>& coeffs);
@@ -134,20 +145,22 @@ public:
 
     // Member Functions
 
-            //- The number of coefficients
-            using scalarList::size;
+        //- Non-zero number of coefficients
+        using scalarList::empty;
 
-            //- Return coefficient
-            using scalarList::operator[];
+        //- The number of coefficients
+        using scalarList::size;
+
+        //- Return coefficient at given index
+        using scalarList::operator[];
 
 
         // Access
 
-
-            //- Return true if the log term is active
+            //- True if the log term is active
             bool logActive() const;
 
-            //- Return the log coefficient
+            //- The log coefficient
             scalar logCoeff() const;
 
 
@@ -164,14 +177,14 @@ public:
             //  Argument becomes zeroth element (constant of integration)
             polynomialFunction integral
             (
-                const scalar intConstant = 0.0
+                const scalar intConstant = 0
             ) const;
 
             //- Return integral coefficients when lowest order is -1.
             //  Argument becomes zeroth element (constant of integration)
             polynomialFunction integralMinus1
             (
-                const scalar intConstant = 0.0
+                const scalar intConstant = 0
             ) const;
 
 
@@ -184,13 +197,13 @@ public:
         polynomialFunction& operator/=(const scalar);
 
 
-    //- Ostream Operator
-    friend Ostream& operator<<(Ostream&, const polynomialFunction&);
+    // IOstream Operators
+
+    friend Istream& operator>>(Istream&, polynomialFunction& poly);
+    friend Ostream& operator<<(Ostream&, const polynomialFunction& poly);
 };
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
 polynomialFunction operator+