Skip to content
Snippets Groups Projects
Commit 43d6b16a authored by Carsten Gräser's avatar Carsten Gräser
Browse files

Remove Polynomial

This class implements the old function interface
and is no longer needed.
parent 629cfe3f
No related branches found
No related tags found
1 merge request!138Remove some Dune::VirtualFunction related stuff
...@@ -28,7 +28,8 @@ ...@@ -28,7 +28,8 @@
`dune-common` 2.7. `dune-common` 2.7.
- The following implementations of the deprecated `Dune::VirtualFunction` interface - The following implementations of the deprecated `Dune::VirtualFunction` interface
have been removed: `ConstantFunction`, `SumFunction`, `SumGridFunction`, `ComposedFunction`, `ComposedGridFunction`. have been removed: `ConstantFunction`, `SumFunction`, `SumGridFunction`,
`ComposedFunction`, `ComposedGridFunction`, `Polynomial`.
# 2.9 Release # 2.9 Release
......
...@@ -7,7 +7,6 @@ install(FILES ...@@ -7,7 +7,6 @@ install(FILES
deformationfunction.hh deformationfunction.hh
localbasisderivativefunction.hh localbasisderivativefunction.hh
localbasisjacobianfunction.hh localbasisjacobianfunction.hh
polynomial.hh
portablegreymap.hh portablegreymap.hh
portablepixelmap.hh portablepixelmap.hh
vintagebasisgridfunction.hh vintagebasisgridfunction.hh
......
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set ts=8 sw=2 et sts=2:
#ifndef POLYNOMIAL_HH
#define POLYNOMIAL_HH
#include <map>
#include <memory>
#include <dune/fufem/functions/virtualdifferentiablefunction.hh>
/** \brief A class implementing polynomials \f$ P:\mathbb{R}\rightarrow\mathbb{R} \f$ of arbitrary degree.
*
* The polynomial expects the coefficients of the monomial representation and implements its evaluation and
* its derivatives of arbitrary order.
* \tparam DT Domain Field type
* \tparam RT Range Field type
*/
template < class DT, class RT >
class Polynomial:
public VirtualDifferentiableFunction<DT,RT>
{
private:
typedef VirtualDifferentiableFunction<DT,RT> BaseType;
public:
typedef std::vector<RT> CoeffType;
typedef typename BaseType::DomainType DomainType;
typedef typename BaseType::RangeType RangeType;
typedef typename BaseType::DerivativeType DerivativeType;
/** \brief Constructor
*
* \param coeffs vector containing the coefficients of the monomial representation
*/
Polynomial(const CoeffType coeffs):
coeffs_(coeffs),
degree_(coeffs.size()-1)
{}
/** \brief polynomial degree
*/
int degree() const
{
return degree_;
}
/** \brief evaluation of the polynomial at point x
*
* The polynomial is evaluated efficiently using Horner's scheme
* \param x the point at which to evaluate
*/
RangeType evaluate(const DomainType& x) const
{
RangeType p = 0.0;
for (int k=degree_; k>=0; --k)
p = x*p + coeffs_[k];
return p;
}
/** \brief evaluation of the polynomial at point x
*
* The polynomial is evaluated efficiently using Horner's scheme
* \param x the point at which to evaluate
* \param y the object to store the function value in
*/
virtual void evaluate(const DomainType& x, RangeType& y) const
{
y = evaluate(x);
return;
}
/** \brief evaluation of derivative
*
* \param x point at which to evaluate
* \param d object to store the derivative
*/
virtual void evaluateDerivative(const DomainType& x, DerivativeType& d) const
{
// Compute result as field, because FM<K,1,1> and FV<k,1,1>
// d not interact nicely and we're having univariate polynomials
// anyway.
using Field = typename Dune::FieldTraits<DerivativeType>::field_type;
Field y = 0.0;
for (int k=degree_; k>=1; --k)
y = x*y + dCoeff(1,k)*coeffs_[k];
d = y;
return;
}
/** \brief Returns the derivative function as a Polynomial object
*
* \param d order of the derivative
*/
Polynomial<DomainType, RangeType>& D(const std::size_t d)
{
if (d==0)
return *this;
typename DerivativeContainer::iterator ret=dP.find(d);
if (ret==dP.end())
{
CoeffType newCoeffs = coeffs_;
for (std::size_t k = degree_; k>=d; --k)
newCoeffs[k] *= dCoeff(d,k);
for (std::size_t k=0; k<d; ++k)
{
typename CoeffType::iterator firstEntry=newCoeffs.begin();
newCoeffs.erase(firstEntry);
}
ret = dP.insert(std::make_pair(d,std::make_shared<Polynomial<DT,RT> >(newCoeffs))).first;
}
return *(ret->second);
}
~Polynomial()
{}
private:
const CoeffType coeffs_;
const unsigned int degree_;
typedef std::map<std::size_t, std::shared_ptr<Polynomial<DT,RT> > > DerivativeContainer;
DerivativeContainer dP;
unsigned int dCoeff(const int& d, const int& k) const
{
unsigned int r=1;
for (int i=0; i<d; ++i)
r *= k-i;
return r;
}
};
#endif
...@@ -20,7 +20,6 @@ set(GRID_BASED_TESTS ...@@ -20,7 +20,6 @@ set(GRID_BASED_TESTS
istlbackendtest istlbackendtest
laplaceassemblertest laplaceassemblertest
localassemblertest localassemblertest
polynomialtest
secondorderassemblertest secondorderassemblertest
subgridxyfunctionalassemblertest subgridxyfunctionalassemblertest
tensortest tensortest
......
#include <config.h>
#include <cmath>
#include <cstdio>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/function.hh>
#include <dune/fufem/functions/polynomial.hh>
#include "common.hh"
/** \brief TestSuite for Polynomial
*
* + compares Polynomial::evaluate to manual evaluation of polynomial (hm, well...)
* + compares Polynomial::evaluateDerivative to manual evaluateDerivative of derivative of polynomial (hm, well)
* + checks consistency of Polynomial::evaluateDerivative() and Polynomial::D().evaluate()
*
*/
struct PolynomialTestSuite
{
bool check()
{
const int n_testpoints=10,
n_testpolynomials = 5,
maxdegree = 10,
maxcoeff = 100;
typedef Dune::FieldVector<double,1> DomainType;
typedef Dune::FieldVector<double,1> RangeType;
typedef Polynomial<DomainType,RangeType> PType;
typedef PType::CoeffType CoeffType;
typedef PType::DerivativeType DerivativeType;
for (int p=0; p<n_testpolynomials; ++p)
{
CoeffType coeffs((rand() % maxdegree) + 2);
std::cout << p << ": coeffs= ";
for (size_t c=0; c<coeffs.size(); ++c)
{
coeffs[c] = (2.0*maxcoeff*rand())/RAND_MAX - maxcoeff;
std::cout << coeffs[c] << " " ;
}
std::cout << std::endl;
PType polynomial(coeffs);
DomainType x(0.0);
RangeType y(0.0),
ycheck(0.0),
Dcheck(0.0);
DerivativeType d(0.0),
dcheck(0.0);
for (int run=0; run<n_testpoints; ++run)
{
x = (4.0*rand())/RAND_MAX - 2.0;
ycheck = 0.0;
for (size_t j=0; j<coeffs.size(); ++j)
ycheck += coeffs[j]*std::pow(x,j);
polynomial.evaluate(x,y);
if (std::abs((y-ycheck)/y)>1e-12)
{
std::cout << "Polynomial::evaluate(" << x << ") = " << y << std::endl;
std::cout << "manual evaluation: " << ycheck << std::endl;
return false;
}
dcheck = 0.0;
for (size_t j=1; j<coeffs.size(); ++j)
dcheck += coeffs[j]*j*std::pow(x,j-1);
polynomial.evaluateDerivative(x,d);
polynomial.D(1).evaluate(x,Dcheck);
if (std::abs((d-dcheck)/d)>1e-12)
{
std::cout << "Polynomial::evaluateDerivative(" << x << ") = " << d << std::endl;
std::cout << "manual derivative evaluation: " << dcheck << std::endl;
return false;
}
if (std::abs((d-Dcheck)/d)>1e-12)
{
std::cout << "Polynomial::evaluateDerivative(" << x << ") = " << d << std::endl;
std::cout << "Polynomial::D().evaluate(" << x << ") = " << Dcheck << std::endl;
return false;
}
}
}
return true;
}
};
int main(int argc, char** argv)
{
Dune::MPIHelper::instance(argc, argv);
std::cout << "This is the PolynomialTest" << std::endl;
PolynomialTestSuite tests;
bool passed = tests.check();
if (passed)
std::cout << "Polynomial test passed" << std::endl;
else
std::cout << "Polynomial test failed" << std::endl;
return passed ? 0 : 1;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment