ENH: add (non-templated) polynomialFunction

STYLE: integrateLimits() -> integrate() for consistency with
  DataEntry/polynomial.  Rename old 'integrate', 'integrateMinus1'
  methods to 'integral', 'integralMinus1' (ie, substantive instead of
  verb for clarity).
This commit is contained in:
Mark Olesen
2011-02-04 18:59:45 +01:00
parent 1a03d43b57
commit ffd20770dd
7 changed files with 809 additions and 82 deletions

View File

@ -41,6 +41,18 @@ Foam::Polynomial<PolySize>::Polynomial()
}
template<int PolySize>
Foam::Polynomial<PolySize>::Polynomial
(
const Polynomial<PolySize>& poly
)
:
VectorSpace<Polynomial<PolySize>, scalar, PolySize>(poly),
logActive_(poly.logActive_),
logCoeff_(poly.logCoeff_)
{}
template<int PolySize>
Foam::Polynomial<PolySize>::Polynomial(const scalar coeffs[PolySize])
:
@ -68,7 +80,7 @@ Foam::Polynomial<PolySize>::Polynomial(const UList<scalar>& coeffs)
(
"Polynomial<PolySize>::Polynomial(const UList<scalar>&)"
) << "Size mismatch: Needed " << PolySize
<< " but got " << coeffs.size()
<< " but given " << coeffs.size()
<< nl << exit(FatalError);
}
@ -79,6 +91,39 @@ Foam::Polynomial<PolySize>::Polynomial(const UList<scalar>& coeffs)
}
// template<int PolySize>
// Foam::Polynomial<PolySize>::Polynomial(const polynomialFunction& poly)
// :
// VectorSpace<Polynomial<PolySize>, scalar, PolySize>(),
// logActive_(poly.logActive()),
// logCoeff_(poly.logCoeff())
// {
// if (poly.size() != PolySize)
// {
// FatalErrorIn
// (
// "Polynomial<PolySize>::Polynomial(const polynomialFunction&)"
// ) << "Size mismatch: Needed " << PolySize
// << " but given " << poly.size()
// << nl << exit(FatalError);
// }
//
// for (int i = 0; i < PolySize; ++i)
// {
// this->v_[i] = poly[i];
// }
// }
template<int PolySize>
Foam::Polynomial<PolySize>::Polynomial(Istream& is)
:
VectorSpace<Polynomial<PolySize>, scalar, PolySize>(is),
logActive_(false),
logCoeff_(0.0)
{}
template<int PolySize>
Foam::Polynomial<PolySize>::Polynomial(const word& name, Istream& is)
:
@ -111,38 +156,17 @@ Foam::Polynomial<PolySize>::Polynomial(const word& name, Istream& is)
}
template<int PolySize>
Foam::Polynomial<PolySize>::Polynomial(Istream& is)
:
VectorSpace<Polynomial<PolySize>, scalar, PolySize>(is),
logActive_(false),
logCoeff_(0.0)
{}
template<int PolySize>
Foam::Polynomial<PolySize>::Polynomial
(
const Polynomial<PolySize>& poly
)
:
VectorSpace<Polynomial<PolySize>, scalar, PolySize>(poly),
logActive_(poly.logActive_),
logCoeff_(poly.logCoeff_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<int PolySize>
bool& Foam::Polynomial<PolySize>::logActive()
bool Foam::Polynomial<PolySize>::logActive() const
{
return logActive_;
}
template<int PolySize>
Foam::scalar& Foam::Polynomial<PolySize>::logCoeff()
Foam::scalar Foam::Polynomial<PolySize>::logCoeff() const
{
return logCoeff_;
}
@ -151,27 +175,27 @@ Foam::scalar& Foam::Polynomial<PolySize>::logCoeff()
template<int PolySize>
Foam::scalar Foam::Polynomial<PolySize>::value(const scalar x) const
{
scalar y = this->v_[0];
scalar val = this->v_[0];
// avoid costly pow() in calculation
scalar powX = x;
for (label i=1; i<PolySize; ++i)
{
y += this->v_[i]*powX;
val += this->v_[i]*powX;
powX *= x;
}
if (logActive_)
{
y += logCoeff_*log(x);
val += logCoeff_*log(x);
}
return y;
return val;
}
template<int PolySize>
Foam::scalar Foam::Polynomial<PolySize>::integrateLimits
Foam::scalar Foam::Polynomial<PolySize>::integrate
(
const scalar x1,
const scalar x2
@ -181,7 +205,7 @@ Foam::scalar Foam::Polynomial<PolySize>::integrateLimits
{
FatalErrorIn
(
"scalar Polynomial<PolySize>::integrateLimits"
"scalar Polynomial<PolySize>::integrate"
"("
"const scalar, "
"const scalar"
@ -190,22 +214,33 @@ Foam::scalar Foam::Polynomial<PolySize>::integrateLimits
<< nl << abort(FatalError);
}
intPolyType poly = this->integrate();
return poly.value(x2) - poly.value(x1);
// avoid costly pow() in calculation
scalar powX1 = x1;
scalar powX2 = x2;
scalar val = this->v_[0]*(powX2 - powX1);
for (label i=1; i<PolySize; ++i)
{
val += this->v_[i]/(i + 1) * (powX2 - powX1);
powX1 *= x1;
powX2 *= x2;
}
return val;
}
template<int PolySize>
typename Foam::Polynomial<PolySize>::intPolyType
Foam::Polynomial<PolySize>::integrate(const scalar intConstant)
Foam::Polynomial<PolySize>::integral(const scalar intConstant) const
{
intPolyType newCoeffs;
newCoeffs[0] = intConstant;
forAll(*this, i)
{
newCoeffs[i + 1] = this->v_[i]/(i + 1);
newCoeffs[i+1] = this->v_[i]/(i + 1);
}
return newCoeffs;
@ -214,14 +249,14 @@ Foam::Polynomial<PolySize>::integrate(const scalar intConstant)
template<int PolySize>
typename Foam::Polynomial<PolySize>::polyType
Foam::Polynomial<PolySize>::integrateMinus1(const scalar intConstant)
Foam::Polynomial<PolySize>::integralMinus1(const scalar intConstant) const
{
polyType newCoeffs;
if (this->v_[0] > VSMALL)
{
newCoeffs.logActive() = true;
newCoeffs.logCoeff() = this->v_[0];
newCoeffs.logActive_ = true;
newCoeffs.logCoeff_ = this->v_[0];
}
newCoeffs[0] = intConstant;

View File

@ -29,14 +29,14 @@ Description
poly = logCoeff*log(x) + sum(coeff_[i]*x^i)
where 0 \<= i \<= n
where 0 \<= i \<= N
- integer powers, starting at zero
- value(x) to evaluate the poly for a given value
- integrate(x1, x2) between two scalar values
- integrate() to return a new, intergated coeff polynomial
- integral() to return a new, integral coeff polynomial
- increases the size (order)
- integrateMinus1() to return a new, integrated coeff polynomial where
- integralMinus1() to return a new, integral coeff polynomial where
the base poly starts at order -1
SourceFiles
@ -85,10 +85,10 @@ class Polynomial
// Private data
//- Include the log term? - only activated using integrateMinus1()
//- Include the log term? - only activated using integralMinus1()
bool logActive_;
//- Log coefficient - only activated using integrateMinus1()
//- Log coefficient - only activated using integralMinus1()
scalar logCoeff_;
@ -104,6 +104,9 @@ public:
//- Construct null, with all coefficients = 0.0
Polynomial();
//- Copy constructor
Polynomial(const Polynomial&);
//- Construct from C-array of coefficients
explicit Polynomial(const scalar coeffs[PolySize]);
@ -111,24 +114,21 @@ public:
explicit Polynomial(const UList<scalar>& coeffs);
//- Construct from Istream
Polynomial(Istream& is);
Polynomial(Istream&);
//- Construct from name and Istream
Polynomial(const word& name, Istream& is);
//- Copy constructor
Polynomial(const Polynomial& poly);
Polynomial(const word& name, Istream&);
// Member Functions
// Access
//- Return access to the log term active flag
bool& logActive();
//- Return true if the log term is active
bool logActive() const;
//- Return access to the log coefficient
scalar& logCoeff();
//- Return the log coefficient
scalar logCoeff() const;
// Evaluation
@ -136,16 +136,17 @@ public:
//- Return polynomial value
scalar value(const scalar x) const;
//- Return integrated polynomial coefficients
// argument becomes zeroth element (constant of integration)
intPolyType integrate(const scalar intConstant = 0.0);
//- Return integrated polynomial coefficients when lowest order
// is -1. Argument added to zeroth element
polyType integrateMinus1(const scalar intConstant = 0.0);
//- Integrate between two values
scalar integrateLimits(const scalar x1, const scalar x2) const;
scalar integrate(const scalar x1, const scalar x2) const;
//- Return integral coefficients.
// Argument becomes zeroth element (constant of integration)
intPolyType integral(const scalar intConstant = 0.0) const;
//- Return integral coefficients when lowest order is -1.
// Argument becomes zeroth element (constant of integration)
polyType integralMinus1(const scalar intConstant = 0.0) const;
//- Ostream Operator

View File

@ -0,0 +1,415 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "polynomialFunction.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(polynomialFunction, 0);
}
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
Foam::polynomialFunction Foam::polynomialFunction::cloneIntegral
(
const polynomialFunction& poly,
const scalar intConstant
)
{
polynomialFunction newPoly(poly.size()+1);
newPoly[0] = intConstant;
forAll(poly, i)
{
newPoly[i+1] = poly[i]/(i + 1);
}
return newPoly;
}
Foam::polynomialFunction Foam::polynomialFunction::cloneIntegralMinus1
(
const polynomialFunction& poly,
const scalar intConstant
)
{
polynomialFunction newPoly(poly.size()+1);
if (poly[0] > VSMALL)
{
newPoly.logActive_ = true;
newPoly.logCoeff_ = poly[0];
}
newPoly[0] = intConstant;
for (label i=1; i < poly.size(); ++i)
{
newPoly[i] = poly[i]/i;
}
return newPoly;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::polynomialFunction::polynomialFunction(const label order)
:
scalarList(order, 0.0),
logActive_(false),
logCoeff_(0.0)
{
if (this->empty())
{
FatalErrorIn
(
"polynomialFunction::polynomialFunction(const label order)"
) << "polynomialFunction coefficients are invalid (empty)"
<< nl << exit(FatalError);
}
}
Foam::polynomialFunction::polynomialFunction(const polynomialFunction& poly)
:
scalarList(poly),
logActive_(poly.logActive_),
logCoeff_(poly.logCoeff_)
{}
Foam::polynomialFunction::polynomialFunction(const UList<scalar>& coeffs)
:
scalarList(coeffs),
logActive_(false),
logCoeff_(0.0)
{
if (this->empty())
{
FatalErrorIn
(
"polynomialFunction::polynomialFunction(const UList<scalar>&)"
) << "polynomialFunction coefficients are invalid (empty)"
<< nl << exit(FatalError);
}
}
Foam::polynomialFunction::polynomialFunction(Istream& is)
:
scalarList(is),
logActive_(false),
logCoeff_(0.0)
{
if (this->empty())
{
FatalErrorIn
(
"polynomialFunction::polynomialFunction(Istream&)"
) << "polynomialFunction coefficients are invalid (empty)"
<< nl << exit(FatalError);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::polynomialFunction::~polynomialFunction()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::polynomialFunction::logActive() const
{
return logActive_;
}
Foam::scalar Foam::polynomialFunction::logCoeff() const
{
return logCoeff_;
}
Foam::scalar Foam::polynomialFunction::value(const scalar x) const
{
const scalarList& coeffs = *this;
scalar val = coeffs[0];
// avoid costly pow() in calculation
scalar powX = x;
for (label i=1; i<coeffs.size(); ++i)
{
val += coeffs[i]*powX;
powX *= x;
}
if (logActive_)
{
val += this->logCoeff_*log(x);
}
return val;
}
Foam::scalar Foam::polynomialFunction::integrate
(
const scalar x1,
const scalar x2
) const
{
const scalarList& coeffs = *this;
if (logActive_)
{
FatalErrorIn
(
"scalar polynomialFunction::integrate"
"("
"const scalar, "
"const scalar"
") const"
) << "Cannot integrate polynomial with logarithmic coefficients"
<< nl << abort(FatalError);
}
// avoid costly pow() in calculation
scalar powX1 = x1;
scalar powX2 = x2;
scalar val = coeffs[0]*(powX2 - powX1);
for (label i=1; i<coeffs.size(); ++i)
{
val += coeffs[i]/(i + 1)*(powX2 - powX1);
powX1 *= x1;
powX2 *= x2;
}
return val;
}
Foam::polynomialFunction
Foam::polynomialFunction::integral(const scalar intConstant) const
{
return cloneIntegral(*this, intConstant);
}
Foam::polynomialFunction
Foam::polynomialFunction::integralMinus1(const scalar intConstant) const
{
return cloneIntegralMinus1(*this, intConstant);
}
// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
Foam::polynomialFunction&
Foam::polynomialFunction::operator+=(const polynomialFunction& poly)
{
scalarList& coeffs = *this;
if (coeffs.size() > poly.size())
{
forAll(poly, i)
{
coeffs[i] += poly[i];
}
}
else
{
coeffs.setSize(poly.size(), 0.0);
forAll(coeffs, i)
{
coeffs[i] += poly[i];
}
}
return *this;
}
Foam::polynomialFunction&
Foam::polynomialFunction::operator-=(const polynomialFunction& poly)
{
scalarList& coeffs = *this;
if (coeffs.size() > poly.size())
{
forAll(poly, i)
{
coeffs[i] -= poly[i];
}
}
else
{
coeffs.setSize(poly.size(), 0.0);
forAll(coeffs, i)
{
coeffs[i] -= poly[i];
}
}
return *this;
}
Foam::polynomialFunction&
Foam::polynomialFunction::operator*=(const scalar s)
{
scalarList& coeffs = *this;
forAll(coeffs, i)
{
coeffs[i] *= s;
}
return *this;
}
Foam::polynomialFunction&
Foam::polynomialFunction::operator/=(const scalar s)
{
scalarList& coeffs = *this;
forAll(coeffs, i)
{
coeffs[i] /= s;
}
return *this;
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<(Ostream& os, const polynomialFunction& poly)
{
// output like VectorSpace
os << token::BEGIN_LIST;
if (!poly.empty())
{
for (int i=0; i<poly.size()-1; i++)
{
os << poly[i] << token::SPACE;
}
os << poly.last();
}
os << token::END_LIST;
// Check state of Ostream
os.check("operator<<(Ostream&, const polynomialFunction&)");
return os;
}
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
Foam::polynomialFunction
Foam::operator+
(
const polynomialFunction& p1,
const polynomialFunction& p2
)
{
polynomialFunction poly(p1);
return poly += p2;
}
Foam::polynomialFunction
Foam::operator-
(
const polynomialFunction& p1,
const polynomialFunction& p2
)
{
polynomialFunction poly(p1);
return poly -= p2;
}
Foam::polynomialFunction
Foam::operator*
(
const scalar s,
const polynomialFunction& p
)
{
polynomialFunction poly(p);
return poly *= s;
}
Foam::polynomialFunction
Foam::operator/
(
const scalar s,
const polynomialFunction& p
)
{
polynomialFunction poly(p);
return poly /= s;
}
Foam::polynomialFunction
Foam::operator*
(
const polynomialFunction& p,
const scalar s
)
{
polynomialFunction poly(p);
return poly *= s;
}
Foam::polynomialFunction
Foam::operator/
(
const polynomialFunction& p,
const scalar s
)
{
polynomialFunction poly(p);
return poly /= s;
}
// ************************************************************************* //

View File

@ -0,0 +1,246 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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/>.
Class
Foam::polynomialFunction
Description
Polynomial function representation
poly = logCoeff*log(x) + sum(coeff_[i]*x^i)
where 0 \<= i \<= N
- integer powers, starting at zero
- value(x) to evaluate the poly for a given value
- integrate(x1, x2) between two scalar values
- integral() to return a new, integral coeff polynomial
- increases the size (order)
- integralMinus1() to return a new, integral coeff polynomial where
the base poly starts at order -1
SeeAlso
Foam::Polynomial for a templated implementation
SourceFiles
polynomialFunction.C
\*---------------------------------------------------------------------------*/
#ifndef polynomialFunction_H
#define polynomialFunction_H
#include "scalarList.H"
#include "Ostream.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class polynomialFunction;
// Forward declaration of friend functions
Ostream& operator<<(Ostream&, const polynomialFunction&);
/*---------------------------------------------------------------------------*\
Class polynomialFunction Declaration
\*---------------------------------------------------------------------------*/
class polynomialFunction
:
private scalarList
{
// Private data
//- Include the log term? - only activated using integralMinus1()
bool logActive_;
//- Log coefficient - only activated using integralMinus1()
scalar logCoeff_;
// Private Member Functions
//- Return integral coefficients.
// Argument becomes zeroth element (constant of integration)
static polynomialFunction cloneIntegral
(
const polynomialFunction&,
const scalar intConstant = 0.0
);
//- Return integral coefficients when lowest order is -1.
// Argument becomes zeroth element (constant of integration)
static polynomialFunction cloneIntegralMinus1
(
const polynomialFunction&,
const scalar intConstant = 0.0
);
//- Disallow default bitwise assignment
void operator=(const polynomialFunction&);
public:
//- Runtime type information
TypeName("polynomialFunction");
// Constructors
//- Construct a particular size, with all coefficients = 0.0
explicit polynomialFunction(const label);
//- Copy constructor
polynomialFunction(const polynomialFunction&);
//- Construct from a list of coefficients
explicit polynomialFunction(const UList<scalar>& coeffs);
//- Construct from Istream
polynomialFunction(Istream&);
//- Destructor
virtual ~polynomialFunction();
// Member Functions
//- Return the number of coefficients
using scalarList::size;
//- Return coefficient
using scalarList::operator[];
// Access
//- Return true if the log term is active
bool logActive() const;
//- Return the log coefficient
scalar logCoeff() const;
// Evaluation
//- Return polynomial value
scalar value(const scalar x) const;
//- Integrate between two values
scalar integrate(const scalar x1, const scalar x2) const;
//- Return integral coefficients.
// Argument becomes zeroth element (constant of integration)
polynomialFunction integral
(
const scalar intConstant = 0.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;
// Member Operators
polynomialFunction& operator+=(const polynomialFunction&);
polynomialFunction& operator-=(const polynomialFunction&);
polynomialFunction& operator*=(const scalar);
polynomialFunction& operator/=(const scalar);
//- Ostream Operator
friend Ostream& operator<<(Ostream&, const polynomialFunction&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
polynomialFunction operator+
(
const polynomialFunction&,
const polynomialFunction&
);
polynomialFunction operator-
(
const polynomialFunction&,
const polynomialFunction&
);
polynomialFunction operator*
(
const scalar,
const polynomialFunction&
);
polynomialFunction operator/
(
const scalar,
const polynomialFunction&
);
polynomialFunction operator*
(
const polynomialFunction&,
const scalar
);
polynomialFunction operator/
(
const polynomialFunction&,
const scalar
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //