From 1593d3c31f336ce3f0a192e447826b09a8a78b0d Mon Sep 17 00:00:00 2001 From: andy Date: Thu, 18 Jun 2009 14:13:08 +0100 Subject: [PATCH] adding support for polynomials starting at order -1 --- .../functions/Polynomial/Polynomial.C | 93 ++++++++++++++++--- .../functions/Polynomial/Polynomial.H | 22 ++++- 2 files changed, 98 insertions(+), 17 deletions(-) diff --git a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C index 8c814e2acf..b842493eef 100644 --- a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C +++ b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.C @@ -25,7 +25,6 @@ License \*---------------------------------------------------------------------------*/ #include "Polynomial.H" -#include "error.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -33,7 +32,9 @@ template Foam::Polynomial::Polynomial() : VectorSpace, scalar, PolySize>(), - name_("unknownPolynomialName") + name_("unknownPolynomialName"), + logActive_(false), + logCoeff_(0.0) {} @@ -41,25 +42,27 @@ template Foam::Polynomial::Polynomial(const word& name, Istream& is) : VectorSpace, scalar, PolySize>(), - name_(is) + name_(is), + logActive_(false), + logCoeff_(0.0) { if (name_ != name) { FatalErrorIn ( - "Foam::Polynomial::Polynomial(const word&, Istream&)" + "Polynomial::Polynomial(const word&, Istream&)" ) << "Expected polynomial name " << name << " but read " << name_ << nl << exit(FatalError); } VectorSpace, scalar, PolySize>:: - operator=(polyType(is)); + operator=(VectorSpace, scalar, PolySize>(is)); if (this->size() == 0) { FatalErrorIn ( - "Foam::Polynomial::Polynomial(const word&, Istream&)" + "Polynomial::Polynomial(const word&, Istream&)" ) << "Polynomial coefficients for entry " << name_ << " are invalid (empty)" << nl << exit(FatalError); } @@ -70,7 +73,9 @@ template Foam::Polynomial::Polynomial(const Polynomial& poly) : VectorSpace, scalar, PolySize>(poly), - name_(poly.name_) + name_(poly.name_), + logActive_(poly.logActive_), + logCoeff_(poly.logCoeff_) {} @@ -82,7 +87,9 @@ Foam::Polynomial::Polynomial ) : VectorSpace, scalar, PolySize>(poly), - name_(name) + name_(name), + logActive_(poly.logActive_), + logCoeff_(poly.logCoeff_) {} @@ -102,15 +109,35 @@ const Foam::word& Foam::Polynomial::name() const } +template +bool& Foam::Polynomial::logActive() +{ + return logActive_; +} + + +template +Foam::scalar& Foam::Polynomial::logCoeff() +{ + return logCoeff_; +} + + template Foam::scalar Foam::Polynomial::evaluate(const scalar x) const { - scalar y = 0.0; - forAll(*this, i) + scalar y = this->v_[0]; + + for (label i=1; iv_[i]*pow(x, i); } + if (logActive_) + { + y += logCoeff_*log(x); + } + return y; } @@ -122,14 +149,22 @@ Foam::scalar Foam::Polynomial::integrateLimits const scalar x2 ) const { - scalar intx = 0.0; - - forAll(*this, i) + if (logActive_) { - intx += this->v_[i]/(i + 1)*(pow(x2, i + 1) - pow(x1, i + 1)); + FatalErrorIn + ( + "scalar Polynomial::integrateLimits" + "(" + "const scalar, " + "const scalar" + ") const" + ) << "Cannot integrate polynomial with logarithmic coefficients" + << nl << abort(FatalError); } - return intx; + intPolyType poly = this->integrate(); + + return poly.evaluate(x2) - poly.evaluate(x1); } @@ -149,6 +184,32 @@ Foam::Polynomial::integrate(const scalar intConstant) } +template +typename Foam::Polynomial::polyType +Foam::Polynomial::integrateMinus1(const scalar intConstant) +{ + polyType newCoeffs; + + if (this->v_[0] > VSMALL) + { + newCoeffs.logActive() = true; + newCoeffs.logCoeff() = this->v_[0]; + } + + newCoeffs[0] = intConstant; + + if (PolySize > 0) + { + for (label i=1; iv_[i]/i; + } + } + + return newCoeffs; +} + + // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // template @@ -156,6 +217,8 @@ void Foam::Polynomial::operator=(const Polynomial& poly) { name_ = poly.name_; VectorSpace, scalar, PolySize>::operator=(poly); + logActive_ = poly.logActive_; + logCoeff_ = poly.logCoeff_; } diff --git a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H index 061b4af377..89ecb9c184 100644 --- a/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H +++ b/src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H @@ -28,7 +28,7 @@ Class Description Polynomial templated on size (order): - poly = sum(coeff_[i]*x^i) + poly = logCoeff*log(x) + sum(coeff_[i]*x^i) where 0 <= i <= n @@ -37,6 +37,8 @@ Description - integrate(x1, x2) between two scalar values - integrate() to return a new, intergated coeff polynomial - increases the size (order) + - integrateMinus1() to return a new, integrated coeff polynomial where + the base poly starts at order -1 SourceFiles Polynomial.C @@ -86,10 +88,16 @@ private: //- Polynomial name word name_; + //- Include the log term? - only activated using integrateMinus1() + bool logActive_; + + //- Log coefficient - only activated using integrateMinus1() + scalar logCoeff_; + public: - typedef VectorSpace, scalar, PolySize> polyType; + typedef Polynomial polyType; typedef Polynomial intPolyType; @@ -123,6 +131,12 @@ public: //- Return const access to the polynomial name const word& name() const; + //- Return access to the log term active flag + bool& logActive(); + + //- Return access to the log coefficient + scalar& logCoeff(); + // Evaluation @@ -133,6 +147,10 @@ public: // 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;