openfoam/applications/test/Polynomial/PolynomialTest.C
2009-12-10 12:50:20 +00:00

129 lines
3.6 KiB
C

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2009-2009 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
PolynomialTest
Description
Test application for the templated Polynomial class
\*---------------------------------------------------------------------------*/
#include "IFstream.H"
#include "Polynomial.H"
#include "Random.H"
using namespace Foam;
scalar polyValue(const scalar x)
{
// Hard-coded polynomial 8 coeff (7th order)
return
0.11
+ 0.45*x
- 0.94*sqr(x)
+ 1.58*pow3(x)
- 2.58*pow4(x)
+ 0.08*pow5(x)
+ 3.15*pow6(x)
- 4.78*x*pow6(x);
}
scalar intPolyValue(const scalar x)
{
// Hard-coded integrated form of above polynomial
return
0.11*x
+ 0.45/2.0*sqr(x)
- 0.94/3.0*pow3(x)
+ 1.58/4.0*pow4(x)
- 2.58/5.0*pow5(x)
+ 0.08/6.0*pow6(x)
+ 3.15/7.0*x*pow6(x)
- 4.78/8.0*x*x*pow6(x);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
IFstream is("polyTestInput");
Polynomial<8> poly("testPoly", is);
Polynomial<9> intPoly(poly.integrate(0.0));
Info<< "poly = " << poly << endl;
Info<< "intPoly = " << intPoly << nl << endl;
Info<< "2*poly = " << 2*poly << endl;
Info<< "poly+poly = " << poly + poly << nl << endl;
Info<< "3*poly = " << 3*poly << endl;
Info<< "poly+poly+poly = " << poly + poly + poly << nl << endl;
Info<< "3*poly - 2*poly = " << 3*poly - 2*poly << nl << endl;
Polynomial<8> polyCopy = poly;
Info<< "poly, polyCopy = " << poly << ", " << polyCopy << nl << endl;
polyCopy = 2.5*poly;
Info<< "2.5*polyCopy = " << polyCopy << nl << endl;
Random rnd(123456);
for (int i=0; i<10; i++)
{
scalar x = rnd.scalar01()*100;
scalar px = polyValue(x);
scalar ipx = intPolyValue(x);
scalar pxTest = poly.evaluate(x);
scalar ipxTest = intPoly.evaluate(x);
Info<<"\nx = " << x << endl;
Info<< " px, pxTest = " << px << ", " << pxTest << endl;
Info<< " ipx, ipxTest = " << ipx << ", " << ipxTest << endl;
if (mag(px - pxTest) > SMALL)
{
Info<< " *** WARNING: px != pxTest: " << px - pxTest << endl;
}
if (mag(ipx - ipxTest) > SMALL)
{
Info<< " *** WARNING: ipx != ipxTest: " << ipx - ipxTest << endl;
}
Info<< endl;
}
Info<< nl << "Done." << endl;
return 0;
}
// ************************************************************************* //