From 2e477e4b915d811e20e37b859ccbc336ba7044e1 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Mon, 6 Aug 2018 14:54:58 +0200 Subject: [PATCH] [common] Add function for generic interpolation --- dumux/common/math.hh | 73 +++++++++++++++++++++++++++++++++++ test/common/math/test_math.cc | 27 +++++++++++++ 2 files changed, 100 insertions(+) diff --git a/dumux/common/math.hh b/dumux/common/math.hh index 330a0a9059..a2ebd0f91a 100644 --- a/dumux/common/math.hh +++ b/dumux/common/math.hh @@ -26,6 +26,7 @@ #include #include +#include #include #include @@ -489,6 +490,78 @@ bool isBetween(const Dune::FieldVector &pos, return false; } +//! forward declaration of the linear interpolation policy (default) +namespace InterpolationPolicy { struct Linear; } + +/*! + * \ingroup Common + * \brief a generic function to interpolate given a set of parameters and an interpolation point + * \param params the parameters used for interpolation (depends on the policy used) + * \param ip the interpolation point + * \param policy the interpolation policy + */ +template +Scalar interpolate(Scalar ip, Parameter&& params) +{ return Policy::interpolate(ip, std::forward(params)); } + +/*! + * \ingroup Common + * \brief Interpolation policies + */ +namespace InterpolationPolicy { + +/*! + * \ingroup Common + * \brief interpolate linearly between two given values + */ +struct Linear +{ + /*! + * \brief interpolate linearly between two given values + * \param ip the interpolation point in [0,1] + * \param array with the lower and upper bound + */ + template + static constexpr Scalar interpolate(Scalar ip, const std::array& params) + { + return params[0]*(1.0 - ip) + params[1]*ip; + } +}; + +/*! + * \ingroup Common + * \brief interpolate linearly in a piecewise linear function (tabularized function) + */ +struct LinearTable +{ + /*! + * \brief interpolate linearly in a piecewise linear function (tabularized function) + * \param ip the interpolation point + * \param table the table as a pair of sorted vectors (have to be same size) + * \note if the interpolation point is out of bounds this will return the bounds + */ + template + static constexpr Scalar interpolate(Scalar ip, const std::pair& table) + { + const auto& range = table.first; + const auto& values = table.second; + + // check bounds + if (ip > range.back()) return values.back(); + if (ip < range[0]) return values[0]; + + // if we are within bounds find the index of the lower bound + const auto lookUpIndex = std::distance(range.begin(), std::lower_bound(range.begin(), range.end(), ip)); + if (lookUpIndex == 0) + return values[0]; + + const auto ipLinear = (ip - range[lookUpIndex-1])/(range[lookUpIndex] - range[lookUpIndex-1]); + return Dumux::interpolate(ipLinear, std::array{{values[lookUpIndex-1], values[lookUpIndex]}}); + } +}; + +} // end namespace InterpolationPolicy + /*! * \ingroup Common diff --git a/test/common/math/test_math.cc b/test/common/math/test_math.cc index 7c62d3e0c3..d723a8f143 100644 --- a/test/common/math/test_math.cc +++ b/test/common/math/test_math.cc @@ -31,6 +31,7 @@ #include #include +#include #include #include @@ -40,6 +41,18 @@ #include +namespace Test { + +template +void checkTableInterpolation(Scalar ip, Scalar expected, const Table& table, Scalar eps = 1e-15) +{ + using namespace Dumux; + auto interp = interpolate(ip, table); + if (!Dune::FloatCmp::eq(interp, expected, eps)) + DUNE_THROW(Dune::Exception, "Wrong interpolation, expected " << expected << " got " << interp); +} + +} // end namespace Test int main() try { @@ -129,6 +142,20 @@ int main() try static_assert(Dumux::sign(2.0) == 1, "Wrong sign!"); static_assert(Dumux::sign(-2) == -1, "Wrong sign!"); static_assert(Dumux::sign(-3.5) == -1, "Wrong sign!"); + + ////////////////////////////////////////////////////////////////// + ///// Dumux::interpolate ///////////////////////////////////////// + ////////////////////////////////////////////////////////////////// + std::vector a{0.0, 1.0, 2.0}; + std::vector b{-1.0, 1.0, 3.0}; + const auto table = std::make_pair(a, b); + + Test::checkTableInterpolation(-1.0, -1.0, table); + Test::checkTableInterpolation(+0.0, -1.0, table); + Test::checkTableInterpolation(0.001, -0.998, table); + Test::checkTableInterpolation(1.5, 2.0, table); + Test::checkTableInterpolation(2.0, 3.0, table); + Test::checkTableInterpolation(3.0, 3.0, table); } catch (Dune::Exception& e) { std::cerr << e << std::endl; -- GitLab