zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Function1s::Polynomial: Simplification

Browse Source 
The coefficients in the polynomial are now specified in order of
increasing exponent, starting from the constant coefficient (i.e., zero
exponent). Exponents are no longer specified. This better fits the
definition of a polynomial, and it prevents the strange scenario when if
exponents are given as a vector or tensor or similar then the units of
the coefficients are not the same across the different components.

For example, polynomial y = -1 + x^2 + 2x^3 can be specified as:

    <name>  polynomial (-1 0 1 2);

Or, alternatively:

    <name>
    {
        type    polynomial;
        coeffs  (-1 0 1 2);
    }
This commit is contained in:
Will Bainbridge
2024-04-23 15:17:48 +01:00
parent f28110ae94
commit c3f131e816
5 changed files with 57 additions and 102 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
applications/test/Function1/polynomial Normal file
View File

@ -0,0 +1,5 @@
function polynomial (0 4 -4 1);
x0 -0.5;
x1 2.5;
nX 51;

110
src/OpenFOAM/primitives/functions/Function1/Polynomial1/Polynomial1.C
View File

@ -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<class Type>
Foam::Function1s::Polynomial<Type>::Polynomial
(
const word& name,
const List<Tuple2<Type, Type>>& coeffs
)
:
FieldFunction1<Type, Polynomial<Type>>(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<Type>::one) < rootVSmall)
{
integrable_ = false;
break;
}
}
if (debug)
{
if (!integrable_)
{
WarningInFunction
<< "Polynomial " << this->name_ << " cannot be integrald"
<< endl;
}
}
}
template<class Type>
Foam::Function1s::Polynomial<Type>::Polynomial
(
@ -73,8 +34,17 @@ Foam::Function1s::Polynomial<Type>::Polynomial
const dictionary& dict
)
:
Polynomial<Type>(name, dict.lookup("coeffs"))
{}
FieldFunction1<Type, Polynomial<Type>>(name),
coeffs_(dict.lookup("coeffs"))
{
if (!coeffs_.size())
{
FatalIOErrorInFunction(dict)
<< typeName.capitalise() << ' ' << name
<< " must have at least one coefficient"
<< exit(FatalError);
}
}
template<class Type>
@ -84,16 +54,24 @@ Foam::Function1s::Polynomial<Type>::Polynomial
Istream& is
)
:
Polynomial<Type>(name, List<Tuple2<Type, Type>>(is))
{}
FieldFunction1<Type, Polynomial<Type>>(name),
coeffs_(is)
{
if (!coeffs_.size())
{
FatalIOErrorInFunction(is)
<< typeName.capitalise() << ' ' << name
<< " must have at least one coefficient"
<< exit(FatalError);
}
}
template<class Type>
Foam::Function1s::Polynomial<Type>::Polynomial(const Polynomial& poly)
:
FieldFunction1<Type, Polynomial<Type>>(poly),
coeffs_(poly.coeffs_),
integrable_(poly.integrable_)
coeffs_(poly.coeffs_)
{}
@ -109,14 +87,11 @@ Foam::Function1s::Polynomial<Type>::~Polynomial()
template<class Type>
Type Foam::Function1s::Polynomial<Type>::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<Type>::one*x, coeffs_[i].second())
);
y = y*x + coeffs_[i];
}
return y;
@ -130,34 +105,15 @@ Type Foam::Function1s::Polynomial<Type>::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<Type>::one
),
cmptPow
(
pTraits<Type>::one*x2,
coeffs_[i].second() + pTraits<Type>::one
)
- cmptPow
(
pTraits<Type>::one*x1,
coeffs_[i].second() + pTraits<Type>::one
)
);
}
Sy1 = Sy1*x1 + coeffs_[i]/(i + 1);
Sy2 = Sy2*x2 + coeffs_[i]/(i + 1);
}
return intx;
return x2*Sy2 - x1*Sy1;
}

40
src/OpenFOAM/primitives/functions/Function1/Polynomial1/Polynomial1.H
View File

@ -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
<name> polynomial
(
(1 2)
(2 3)
);
<name> polynomial (-1 0 1 2);
\endverbatim
Or, alternatively:
\verbatim
<name>
{
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<Tuple2<Type, Type>> coeffs_;
//- Flag to indicate whether the function can be integrated
bool integrable_;
// Private member functions
//- Polynomial coefficients
const List<Type> coeffs_;
public:
@ -85,13 +86,6 @@ public:
// Constructors
//- Construct from components
Polynomial
(
const word& name,
const List<Tuple2<Type, Type>>&
);
//- Construct from name and dictionary
Polynomial(const word& name, const dictionary& dict);

2
tutorials/incompressibleFluid/TJunctionFan/system/createBafflesDict
View File

@ -99,7 +99,7 @@ baffles
patchType cyclic;
jump uniform 0;
value uniform 0;
jumpTable polynomial 1((100 0));
jumpTable polynomial (100);
}
}
}

2
tutorials/incompressibleFluid/TJunctionFan/system/createNonConformalCouplesDict
View File

@ -32,7 +32,7 @@ nonConformalCouples
patchType nonConformalCyclic;
jump uniform 0;
value uniform 0;
jumpTable polynomial 1((100 0));
jumpTable polynomial (100);
}
}
}