Commit bc9e97cf authored by Mark Olesen's avatar Mark Olesen
Browse files

ENH: additional polynomial constructors, improved I/O

- support construct from initializer_list, which can help simplify
  code with constant coefficients.

- add default constructor for polynomialFunction and Istream reading
  to support resizable lists of polynomialFunction.

  A default constructed polynomialFunction is simply equivalent to
  a constant zero.

- no special IO handling for Polynomial required,
  it is the same as VectorSpace anyhow.
parent 95f7ed03
......@@ -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)};
polynomialFunction polyfunc;
IStringStream is(polyDef);
polynomialFunction polyfunc(is);
// 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)
......
......@@ -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;
......
......@@ -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
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
......@@ -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
......
......@@ -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");