diff --git a/applications/test/Function1/polynomial b/applications/test/Function1/polynomial new file mode 100644 index 0000000000..e83f8a68bd --- /dev/null +++ b/applications/test/Function1/polynomial @@ -0,0 +1,5 @@ +function polynomial (0 4 -4 1); + +x0 -0.5; +x1 2.5; +nX 51; diff --git a/src/OpenFOAM/primitives/functions/Function1/Polynomial1/Polynomial1.C b/src/OpenFOAM/primitives/functions/Function1/Polynomial1/Polynomial1.C index b64e5bc933..521a135dcf 100644 --- a/src/OpenFOAM/primitives/functions/Function1/Polynomial1/Polynomial1.C +++ b/src/OpenFOAM/primitives/functions/Function1/Polynomial1/Polynomial1.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2024 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -27,45 +27,6 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -template -Foam::Function1s::Polynomial::Polynomial -( - const word& name, - const List>& coeffs -) -: - FieldFunction1>(name), - coeffs_(coeffs), - integrable_(true) -{ - if (!coeffs_.size()) - { - FatalErrorInFunction - << "Polynomial coefficients for entry " << this->name_ - << " are invalid (empty)" << nl << exit(FatalError); - } - - forAll(coeffs_, i) - { - if (mag(coeffs_[i].second() + pTraits::one) < rootVSmall) - { - integrable_ = false; - break; - } - } - - if (debug) - { - if (!integrable_) - { - WarningInFunction - << "Polynomial " << this->name_ << " cannot be integrald" - << endl; - } - } -} - - template Foam::Function1s::Polynomial::Polynomial ( @@ -73,8 +34,17 @@ Foam::Function1s::Polynomial::Polynomial const dictionary& dict ) : - Polynomial(name, dict.lookup("coeffs")) -{} + FieldFunction1>(name), + coeffs_(dict.lookup("coeffs")) +{ + if (!coeffs_.size()) + { + FatalIOErrorInFunction(dict) + << typeName.capitalise() << ' ' << name + << " must have at least one coefficient" + << exit(FatalError); + } +} template @@ -84,16 +54,24 @@ Foam::Function1s::Polynomial::Polynomial Istream& is ) : - Polynomial(name, List>(is)) -{} + FieldFunction1>(name), + coeffs_(is) +{ + if (!coeffs_.size()) + { + FatalIOErrorInFunction(is) + << typeName.capitalise() << ' ' << name + << " must have at least one coefficient" + << exit(FatalError); + } +} template Foam::Function1s::Polynomial::Polynomial(const Polynomial& poly) : FieldFunction1>(poly), - coeffs_(poly.coeffs_), - integrable_(poly.integrable_) + coeffs_(poly.coeffs_) {} @@ -109,14 +87,11 @@ Foam::Function1s::Polynomial::~Polynomial() template Type Foam::Function1s::Polynomial::value(const scalar x) const { - Type y(Zero); - forAll(coeffs_, i) + Type y = coeffs_[coeffs_.size() - 1]; + + for (label i = coeffs_.size() - 2; i >= 0; i --) { - y += cmptMultiply - ( - coeffs_[i].first(), - cmptPow(pTraits::one*x, coeffs_[i].second()) - ); + y = y*x + coeffs_[i]; } return y; @@ -130,34 +105,15 @@ Type Foam::Function1s::Polynomial::integral const scalar x2 ) const { - Type intx(Zero); + Type Sy1 = coeffs_[coeffs_.size() - 1]/coeffs_.size(), Sy2 = Sy1; - if (integrable_) + for (label i = coeffs_.size() - 2; i >= 0; i --) { - forAll(coeffs_, i) - { - intx += cmptMultiply - ( - cmptDivide - ( - coeffs_[i].first(), - coeffs_[i].second() + pTraits::one - ), - cmptPow - ( - pTraits::one*x2, - coeffs_[i].second() + pTraits::one - ) - - cmptPow - ( - pTraits::one*x1, - coeffs_[i].second() + pTraits::one - ) - ); - } + Sy1 = Sy1*x1 + coeffs_[i]/(i + 1); + Sy2 = Sy2*x2 + coeffs_[i]/(i + 1); } - return intx; + return x2*Sy2 - x1*Sy1; } diff --git a/src/OpenFOAM/primitives/functions/Function1/Polynomial1/Polynomial1.H b/src/OpenFOAM/primitives/functions/Function1/Polynomial1/Polynomial1.H index bd343c1e88..1c8b35508f 100644 --- a/src/OpenFOAM/primitives/functions/Function1/Polynomial1/Polynomial1.H +++ b/src/OpenFOAM/primitives/functions/Function1/Polynomial1/Polynomial1.H @@ -27,15 +27,23 @@ Class Description Arbitrary order polynomial Function1. - The coefficients and exponents of the terms in the polynomial are specified - as a list of Tuple2's e.g. for the polynomial y = x^2 + 2x^3 + The coefficients in the polynomial are specified in order of increasing + exponent, starting from the constant coefficient (i.e., zero exponent). + + For example, the polynomial y = -1 + x^2 + 2x^3 can be specified as: \verbatim - polynomial - ( - (1 2) - (2 3) - ); + polynomial (-1 0 1 2); + \endverbatim + + Or, alternatively: + + \verbatim + + { + type polynomial; + coeffs (-1 0 1 2); + } \endverbatim SourceFiles @@ -47,7 +55,6 @@ SourceFiles #define Polynomial1_H #include "Function1.H" -#include "Tuple2.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -67,14 +74,8 @@ class Polynomial { // Private Data - //- Polynomial coefficients - list of prefactor, exponent - const List> coeffs_; - - //- Flag to indicate whether the function can be integrated - bool integrable_; - - - // Private member functions + //- Polynomial coefficients + const List coeffs_; public: @@ -85,13 +86,6 @@ public: // Constructors - //- Construct from components - Polynomial - ( - const word& name, - const List>& - ); - //- Construct from name and dictionary Polynomial(const word& name, const dictionary& dict); diff --git a/tutorials/incompressibleFluid/TJunctionFan/system/createBafflesDict b/tutorials/incompressibleFluid/TJunctionFan/system/createBafflesDict index 5f0716b00d..e3fbb3ecbe 100644 --- a/tutorials/incompressibleFluid/TJunctionFan/system/createBafflesDict +++ b/tutorials/incompressibleFluid/TJunctionFan/system/createBafflesDict @@ -99,7 +99,7 @@ baffles patchType cyclic; jump uniform 0; value uniform 0; - jumpTable polynomial 1((100 0)); + jumpTable polynomial (100); } } } diff --git a/tutorials/incompressibleFluid/TJunctionFan/system/createNonConformalCouplesDict b/tutorials/incompressibleFluid/TJunctionFan/system/createNonConformalCouplesDict index 6bf64bcbc2..1a2d7b590d 100644 --- a/tutorials/incompressibleFluid/TJunctionFan/system/createNonConformalCouplesDict +++ b/tutorials/incompressibleFluid/TJunctionFan/system/createNonConformalCouplesDict @@ -32,7 +32,7 @@ nonConformalCouples patchType nonConformalCyclic; jump uniform 0; value uniform 0; - jumpTable polynomial 1((100 0)); + jumpTable polynomial (100); } } }