From 7a4c6546b43855ddd3ed4c3407cc9d9be80a8e47 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 19 Dec 2018 14:07:16 +0100 Subject: [PATCH] fix some coding style isses in new class XYTabulated2DFunction and its test --- opm/material/common/XYTabulated2DFunction.hpp | 221 ++++--- tests/test_2dtables.cpp | 605 +++++++++--------- 2 files changed, 423 insertions(+), 403 deletions(-) diff --git a/opm/material/common/XYTabulated2DFunction.hpp b/opm/material/common/XYTabulated2DFunction.hpp index c7c93afe4..1546c2bcf 100644 --- a/opm/material/common/XYTabulated2DFunction.hpp +++ b/opm/material/common/XYTabulated2DFunction.hpp @@ -42,91 +42,69 @@ namespace Opm { /*! - * \brief Implements a scalar function that depends on two variables, which are sampled - * in X and Y directions with sampling points, respectively. - * The table might be able to be extrapolated in certian directions. + * \brief Implements a function that depends on two variables. + * + * The function is sampled in regular intervals in both directions, i.e., the + * interpolation cells are rectangles. The table can be specified to be extrapolated in + * either direction. */ - template class XYTabulated2DFunction { - public: XYTabulated2DFunction() { } - XYTabulated2DFunction(const std::vector& x_pos, - const std::vector& y_pos, - const std::vector >& data, - const bool x_extrapolate = false, - const bool y_extrapolate = false) - : xPos_(x_pos) - , yPos_(y_pos) + template + XYTabulated2DFunction(const std::vector& xPos, + const std::vector& yPos, + const DataContainer& data, + const bool xExtrapolate = false, + const bool yExtrapolate = false) + : xPos_(xPos) + , yPos_(yPos) , samples_(data) - , xExtrapolate_(x_extrapolate) - , yExtrapolate_(y_extrapolate) + , xExtrapolate_(xExtrapolate) + , yExtrapolate_(yExtrapolate) { - // make sure the size is correct - if (numX() != samples_.size()) { - throw std::runtime_error("numX() is not equal to the number of rows of the sample points"); +#ifndef NDEBUG + // in debug mode, ensure that the x and y positions arrays are strictly + // mononically increasing. + for (unsigned i = 0; i < xPos.size() - 1; ++ i) { + if (xPos[i + 1] <= xPos[i]) + throw std::runtime_error("The array for the x-positions is not strictly increasing!"); } - for (unsigned ix = 0; ix < numX(); ++ix) { - if (samples_[ix].size() != numY()) { + for (unsigned i = 0; i < yPos.size() - 1; ++ i) { + if (yPos[i + 1] <= yPos[i]) + throw std::runtime_error("The array for the y-positions is not strictly increasing!"); + } +#endif + + // make sure the size is correct + if (numX() != samples_.size()) + throw std::runtime_error("numX() is not equal to the number of rows of the sampling points"); + + for (unsigned xIdx = 0; xIdx < numX(); ++xIdx) { + if (samples_[xIdx].size() != numY()) { std::ostringstream oss; - oss << "the " << ix << "th row of the sample points has different number of data from numY() "; + oss << "The " << xIdx << "-th row of the sampling points has different size than numY() "; throw std::runtime_error(oss.str()); } } } /*! - * \brief Evaluate the function at a given (x,y) position. - * - * If this method is called for a value outside of the tabulated - * range, and extrapolation is not allowed in the corresponding direction, - * a \c Opm::NumericalIssue exception is thrown. + * \brief Returns the number of sampling points in X direction. */ - template - void eval(const DataType& x, const DataType& y, DataType& result) const - { - if ( (!xExtrapolate_ && !appliesX(x)) || - (!yExtrapolate_ && !appliesY(y)) ) { - std::ostringstream oss; - oss << "Attempt to get undefined table value (" << x << ", " << y << ")"; - throw NumericalIssue(oss.str()); - }; + size_t numX() const + { return xPos_.size(); } - // bi-linear interpolation: first, calculate the x and y indices in the lookup - // table ... - const unsigned i = xSegmentIndex(x); - const unsigned j = ySegmentIndex(y); - - // bi-linear interpolation / extrapolation - const DataType alpha = xToAlpha(x, i); - const DataType beta = yToBeta(y, j); - - const DataType s1 = valueAt(i, j) * (1.0 - beta) + valueAt(i, j + 1) * beta; - const DataType s2 = valueAt(i + 1, j) * (1.0 - beta) + valueAt(i + 1, j + 1) * beta; - - Valgrind::CheckDefined(s1); - Valgrind::CheckDefined(s2); - - // ... and combine them using the x position - result = s1 * (1.0 - alpha) + s2 * alpha; - Valgrind::CheckDefined(result); - } - -private: - // the sampling points in the x-drection - std::vector xPos_; - // the sampling points in the y-drection - std::vector yPos_; - // data at the sampling points - std::vector > samples_; - - bool xExtrapolate_ = false; - bool yExtrapolate_ = false; + /*! + * \brief Returns the number of sampling points in Y direction. + */ + size_t numY() const + { return yPos_.size(); } /*! * \brief Returns the minimum of the X coordinate of the sampling points. @@ -158,83 +136,126 @@ private: Scalar valueAt(size_t i, size_t j) const { return samples_[i][j]; } + /*! + * \brief Returns true if a coordinate lies in the tabulated range + */ + template + bool applies(const Evaluation& x, const Evaluation& y) const + { return appliesX(x) && appliesY(y); } + /*! * \brief Returns true if a coordinate lies in the tabulated range on the x direction */ template bool appliesX(const Evaluation& x) const - { - if (x < xMin() || xMax() < x) { - return false; - } else { - return true; - } - } + { return xMin() <= x && x <= xMax(); } /*! * \brief Returns true if a coordinate lies in the tabulated range on the y direction */ template bool appliesY(const Evaluation& y) const + { return yMin() <= y && y <= yMax(); } + + + /*! + * \brief Evaluate the function at a given (x,y) position. + * + * If this method is called for a value outside of the tabulated + * range, and extrapolation is not allowed in the corresponding direction, + * a \c Opm::NumericalIssue exception is thrown. + */ + template + void eval(const Evaluation& x, const Evaluation& y, Evaluation& result) const { - if (y < yMin() || yMax() < y) { - return false; - } else { - return true; - } + if ((!xExtrapolate_ && !appliesX(x)) || (!yExtrapolate_ && !appliesY(y))) { + std::ostringstream oss; + oss << "Attempt to get undefined table value (" << x << ", " << y << ")"; + throw NumericalIssue(oss.str()); + }; + + // bi-linear interpolation: first, calculate the x and y indices in the lookup + // table ... + const unsigned i = xSegmentIndex_(x); + const unsigned j = ySegmentIndex_(y); + + // bi-linear interpolation / extrapolation + const Evaluation alpha = xToAlpha(x, i); + const Evaluation beta = yToBeta(y, j); + + const Evaluation s1 = valueAt(i, j) * (1.0 - beta) + valueAt(i, j + 1) * beta; + const Evaluation s2 = valueAt(i + 1, j) * (1.0 - beta) + valueAt(i + 1, j + 1) * beta; + + Valgrind::CheckDefined(s1); + Valgrind::CheckDefined(s2); + + // ... and combine them using the x position + result = s1 * (1.0 - alpha) + s2 * alpha; + Valgrind::CheckDefined(result); } - /*! - * \brief Returns the number of sampling points in X direction. - */ - size_t numX() const - { return xPos_.size(); } +private: + // the sampling points in the x-drection + std::vector xPos_; + // the sampling points in the y-drection + std::vector yPos_; + // data at the sampling points + std::vector > samples_; - /*! - * \brief Returns the number of sampling points in Y direction. - */ - size_t numY() const - { return yPos_.size(); } + bool xExtrapolate_ = false; + bool yExtrapolate_ = false; /*! * \brief Return the interval index of a given position on the x-axis. */ template - unsigned xSegmentIndex(const Evaluation& x) const + unsigned xSegmentIndex_(const Evaluation& x) const { assert(xExtrapolate_ || appliesX(x) ); - return segmentIndex(x, xPos_); + return segmentIndex_(x, xPos_); } /*! * \brief Return the interval index of a given position on the y-axis. */ template - unsigned ySegmentIndex(const Evaluation& y) const + unsigned ySegmentIndex_(const Evaluation& y) const { assert(yExtrapolate_ || appliesY(y) ); - return segmentIndex(y, yPos_); + return segmentIndex_(y, yPos_); } template - static unsigned segmentIndex(const Evaluation& y, const std::vector& pos) + static unsigned segmentIndex_(const Evaluation& v, const std::vector& vPos) { + const unsigned n = vPos.size(); + assert(n >= 2); - const unsigned num_pos = pos.size(); - assert(num_pos >= 2); - - if (y <= pos.front() || num_pos == 2) { + if (v <= vPos.front() || n == 2) return 0; - } else if (y >= pos.back()) { - return num_pos - 2; - } else { - assert(num_pos >= 3); + else if (v >= vPos.back()) + return n - 2; - return --std::lower_bound(pos.begin(), pos.end(), y) - pos.begin(); + assert(n > 2 && v > vPos.front() && v < vPos.back()); + + // bisection. this assumes that the vPos array is strictly mononically + // increasing. + size_t lowerIdx = 0; + size_t upperIdx = vPos.size() - 1; + while (lowerIdx + 1 < upperIdx) { + size_t pivotIdx = (lowerIdx + upperIdx) / 2; + if (v < vPos[pivotIdx]) + upperIdx = pivotIdx; + else + lowerIdx = pivotIdx; } + + assert(vPos[lowerIdx] <= v); + assert(v <= vPos[lowerIdx + 1]); + return lowerIdx; } /*! diff --git a/tests/test_2dtables.cpp b/tests/test_2dtables.cpp index 8f03e461d..d3ab7fe7a 100644 --- a/tests/test_2dtables.cpp +++ b/tests/test_2dtables.cpp @@ -39,363 +39,362 @@ #include #include - template struct Test { -typedef ScalarT Scalar; + typedef ScalarT Scalar; -static Scalar testFn1(Scalar x, Scalar /* y */) -{ return x; } + static Scalar testFn1(Scalar x, Scalar /* y */) + { return x; } -static Scalar testFn2(Scalar /* x */, Scalar y) -{ return y; } + static Scalar testFn2(Scalar /* x */, Scalar y) + { return y; } -static Scalar testFn3(Scalar x, Scalar y) -{ return x*y; } + static Scalar testFn3(Scalar x, Scalar y) + { return x*y; } -template -std::shared_ptr > -createUniformTabulatedFunction(Fn& f) -{ - Scalar xMin = -2.0; - Scalar xMax = 3.0; - unsigned m = 50; + template + std::shared_ptr > + createUniformTabulatedFunction(Fn& f) + { + Scalar xMin = -2.0; + Scalar xMax = 3.0; + unsigned m = 50; - Scalar yMin = -1/2.0; - Scalar yMax = 1/3.0; - unsigned n = 40; + Scalar yMin = -1/2.0; + Scalar yMax = 1/3.0; + unsigned n = 40; - auto tab = std::make_shared>( - xMin, xMax, m, - yMin, yMax, n); - for (unsigned i = 0; i < m; ++i) { - Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin); - for (unsigned j = 0; j < n; ++j) { - Scalar y = yMin + Scalar(j)/(n - 1) * (yMax - yMin); - tab->setSamplePoint(i, j, f(x, y)); + auto tab = std::make_shared>( + xMin, xMax, m, + yMin, yMax, n); + for (unsigned i = 0; i < m; ++i) { + Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin); + for (unsigned j = 0; j < n; ++j) { + Scalar y = yMin + Scalar(j)/(n - 1) * (yMax - yMin); + tab->setSamplePoint(i, j, f(x, y)); + } } + + return tab; } - return tab; -} + template + std::shared_ptr > + createUniformXTabulatedFunction(Fn& f) + { + Scalar xMin = -2.0; + Scalar xMax = 3.0; + unsigned m = 50; -template -std::shared_ptr > -createUniformXTabulatedFunction(Fn& f) -{ - Scalar xMin = -2.0; - Scalar xMax = 3.0; - unsigned m = 50; + Scalar yMin = -1/2.0; + Scalar yMax = 1/3.0; + unsigned n = 40; - Scalar yMin = -1/2.0; - Scalar yMax = 1/3.0; - unsigned n = 40; - - auto tab = std::make_shared>(); - for (unsigned i = 0; i < m; ++i) { - Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin); - tab->appendXPos(x); - for (unsigned j = 0; j < n; ++j) { - Scalar y = yMin + Scalar(j)/(n -1) * (yMax - yMin); - tab->appendSamplePoint(i, y, f(x, y)); + auto tab = std::make_shared>(); + for (unsigned i = 0; i < m; ++i) { + Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin); + tab->appendXPos(x); + for (unsigned j = 0; j < n; ++j) { + Scalar y = yMin + Scalar(j)/(n -1) * (yMax - yMin); + tab->appendSamplePoint(i, y, f(x, y)); + } } + + return tab; } - return tab; -} + template + std::shared_ptr > + createUniformXTabulatedFunction2(Fn& f) + { + Scalar xMin = -2.0; + Scalar xMax = 3.0; + Scalar m = 50; -template -std::shared_ptr > -createUniformXTabulatedFunction2(Fn& f) -{ - Scalar xMin = -2.0; - Scalar xMax = 3.0; - Scalar m = 50; + Scalar yMin = - 4.0; + Scalar yMax = 5.0; - Scalar yMin = - 4.0; - Scalar yMax = 5.0; + auto tab = std::make_shared>(); + for (unsigned i = 0; i < m; ++i) { + Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin); + tab->appendXPos(x); - auto tab = std::make_shared>(); - for (unsigned i = 0; i < m; ++i) { - Scalar x = xMin + Scalar(i)/(m - 1) * (xMax - xMin); - tab->appendXPos(x); + Scalar n = i + 10; - Scalar n = i + 10; - - for (unsigned j = 0; j < n; ++j) { - Scalar y = yMin + Scalar(j)/(n -1) * (yMax - yMin); - tab->appendSamplePoint(i, y, f(x, y)); + for (unsigned j = 0; j < n; ++j) { + Scalar y = yMin + Scalar(j)/(n -1) * (yMax - yMin); + tab->appendSamplePoint(i, y, f(x, y)); + } } + + return tab; } - return tab; -} + template + std::shared_ptr > + createXYTabulated2DFunction(Fn& f) + { + const Scalar xMin = -2.0; + const Scalar xMax = 3.0; + const unsigned m = 50; -template -std::shared_ptr > -createXYTabulated2DFunction(Fn& f) -{ - const Scalar xMin = -2.0; - const Scalar xMax = 3.0; - const unsigned m = 50; + const Scalar yMin = -1/2.0; + const Scalar yMax = 1/3.0; + const unsigned n = 40; - const Scalar yMin = -1/2.0; - const Scalar yMax = 1/3.0; - const unsigned n = 40; + std::vector x_samples(m); + std::vector y_samples(n); + std::vector> data(m, std::vector(n)); - std::vector x_samples(m); - std::vector y_samples(n); - std::vector> data(m, std::vector(n)); + for (unsigned i = 0; i < m; ++i) { + x_samples[i] = xMin + Scalar(i)/(m - 1) * (xMax - xMin); - for (unsigned i = 0; i < m; ++i) { - x_samples[i] = xMin + Scalar(i)/(m - 1) * (xMax - xMin); - - for (unsigned j = 0; j < n; ++j) { - y_samples[j] = yMin + Scalar(j)/(n -1) * (yMax - yMin); - data[i][j] = f(x_samples[i], y_samples[j]); + for (unsigned j = 0; j < n; ++j) { + y_samples[j] = yMin + Scalar(j)/(n -1) * (yMax - yMin); + data[i][j] = f(x_samples[i], y_samples[j]); + } } + + auto tab = std::make_shared>(x_samples, y_samples, data, true, true); + + return tab; } - auto tab = std::make_shared>(x_samples, y_samples, data, true, true); + template + bool compareTableWithAnalyticFn(const Table& table, + Scalar xMin, + Scalar xMax, + unsigned numX, - return tab; -} + Scalar yMin, + Scalar yMax, + unsigned numY, -template -bool compareTableWithAnalyticFn(const Table& table, - Scalar xMin, - Scalar xMax, - unsigned numX, + Fn& f, + Scalar tolerance = 1e-8) + { + // make sure that the tabulated function evaluates to the same thing as the analytic + // one (modulo tolerance) + for (unsigned i = 1; i <= numX; ++i) { + Scalar x = xMin + Scalar(i)/numX*(xMax - xMin); - Scalar yMin, - Scalar yMax, - unsigned numY, + for (unsigned j = 0; j < numY; ++j) { + Scalar y = yMin + Scalar(j)/numY*(yMax - yMin); + if (std::abs(table->eval(x, y) - f(x, y)) > tolerance) { + std::cerr << __FILE__ << ":" << __LINE__ << ": table->eval("<eval(x,y) << " != " << f(x,y) << "\n"; + return false; + } + } + } - Fn& f, - Scalar tolerance = 1e-8) -{ - // make sure that the tabulated function evaluates to the same thing as the analytic - // one (modulo tolerance) - for (unsigned i = 1; i <= numX; ++i) { - Scalar x = xMin + Scalar(i)/numX*(xMax - xMin); + return true; + } - for (unsigned j = 0; j < numY; ++j) { - Scalar y = yMin + Scalar(j)/numY*(yMax - yMin); - if (std::abs(table->eval(x, y) - f(x, y)) > tolerance) { - std::cerr << __FILE__ << ":" << __LINE__ << ": table->eval("<eval(x,y) << " != " << f(x,y) << "\n"; + template + bool compareTableWithAnalyticFn2(const Table& table, + const Scalar xMin, + const Scalar xMax, + unsigned numX, + const Scalar yMin, + const Scalar yMax, + unsigned numY, + + Fn& f, + Scalar tolerance = 1e-8) + { + // make sure that the tabulated function evaluates to the same thing as the analytic + // one (modulo tolerance) + for (unsigned i = 1; i <= numX; ++i) { + Scalar x = xMin + Scalar(i)/numX*(xMax - xMin); + + for (unsigned j = 0; j < numY; ++j) { + Scalar y = yMin + Scalar(j)/numY*(yMax - yMin); + Scalar result; + table->eval(x, y, result); + if (std::abs(result - f(x, y)) > tolerance) { + std::cerr << __FILE__ << ":" << __LINE__ << ": table->eval("< + bool compareTables(const UniformTablePtr uTable, + const UniformXTablePtr uXTable, + Fn& f, + Scalar tolerance = 1e-8) + { + // make sure the uniform and the non-uniform tables exhibit the same dimensions + if (std::abs(uTable->xMin() - uXTable->xMin()) > tolerance) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->xMin() != uXTable->xMin(): " << uTable->xMin() << " != " << uXTable->xMin() << "\n"; + return false; + } + if (std::abs(uTable->xMax() - uXTable->xMax()) > tolerance) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->xMax() != uXTable->xMax(): " << uTable->xMax() << " != " << uXTable->xMax() << "\n"; + return false; + } + if (uTable->numX() != uXTable->numX()) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->numX() != uXTable->numX(): " << uTable->numX() << " != " << uXTable->numX() << "\n"; + return false; + } + + for (unsigned i = 0; i < uTable->numX(); ++i) { + if (std::abs(uTable->yMin() - uXTable->yMin(i)) > tolerance) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->yMin() != uXTable->yMin("<yMin() << " != " << uXTable->yMin(i) << "\n"; + return false; + } + + if (std::abs(uTable->yMax() - uXTable->yMax(i)) > tolerance) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->yMax() != uXTable->yMax("<yMax() << " != " << uXTable->yMax(i) << "\n"; + return false; + } + + if (uTable->numY() != uXTable->numY(i)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->numY() != uXTable->numY("<numY() << " != " << uXTable->numY(i) << "\n"; return false; } } - } - return true; -} - -template -bool compareTableWithAnalyticFn2(const Table& table, - const Scalar xMin, - const Scalar xMax, - unsigned numX, - const Scalar yMin, - const Scalar yMax, - unsigned numY, - - Fn& f, - Scalar tolerance = 1e-8) -{ - // make sure that the tabulated function evaluates to the same thing as the analytic - // one (modulo tolerance) - for (unsigned i = 1; i <= numX; ++i) { - Scalar x = xMin + Scalar(i)/numX*(xMax - xMin); - - for (unsigned j = 0; j < numY; ++j) { - Scalar y = yMin + Scalar(j)/numY*(yMax - yMin); - Scalar result; - table->eval(x, y, result); - if (std::abs(result - f(x, y)) > tolerance) { - std::cerr << __FILE__ << ":" << __LINE__ << ": table->eval("<numX(); ++i) { + if (std::abs(uTable->iToX(i) - uXTable->iToX(i)) > tolerance) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->iToX("<iToX("<iToX(i) << " != " << uXTable->iToX(i) << "\n"; return false; } - } - } - return true; -} - -template -bool compareTables(const UniformTablePtr uTable, - const UniformXTablePtr uXTable, - Fn& f, - Scalar tolerance = 1e-8) -{ - // make sure the uniform and the non-uniform tables exhibit the same dimensions - if (std::abs(uTable->xMin() - uXTable->xMin()) > tolerance) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->xMin() != uXTable->xMin(): " << uTable->xMin() << " != " << uXTable->xMin() << "\n"; - return false; - } - if (std::abs(uTable->xMax() - uXTable->xMax()) > tolerance) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->xMax() != uXTable->xMax(): " << uTable->xMax() << " != " << uXTable->xMax() << "\n"; - return false; - } - if (uTable->numX() != uXTable->numX()) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->numX() != uXTable->numX(): " << uTable->numX() << " != " << uXTable->numX() << "\n"; - return false; - } - - for (unsigned i = 0; i < uTable->numX(); ++i) { - if (std::abs(uTable->yMin() - uXTable->yMin(i)) > tolerance) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->yMin() != uXTable->yMin("<yMin() << " != " << uXTable->yMin(i) << "\n"; - return false; - } - - if (std::abs(uTable->yMax() - uXTable->yMax(i)) > tolerance) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->yMax() != uXTable->yMax("<yMax() << " != " << uXTable->yMax(i) << "\n"; - return false; - } - - if (uTable->numY() != uXTable->numY(i)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->numY() != uXTable->numY("<numY() << " != " << uXTable->numY(i) << "\n"; - return false; - } - } - - // make sure that the x and y values are identical - for (unsigned i = 0; i < uTable->numX(); ++i) { - if (std::abs(uTable->iToX(i) - uXTable->iToX(i)) > tolerance) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->iToX("<iToX("<iToX(i) << " != " << uXTable->iToX(i) << "\n"; - return false; - } - - for (unsigned j = 0; j < uTable->numY(); ++j) { - if (std::abs(uTable->jToY(j) - uXTable->jToY(i, j)) > tolerance) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->jToY("<jToY("<jToY(i) << " != " << uXTable->jToY(i, j) << "\n"; - return false; + for (unsigned j = 0; j < uTable->numY(); ++j) { + if (std::abs(uTable->jToY(j) - uXTable->jToY(i, j)) > tolerance) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->jToY("<jToY("<jToY(i) << " != " << uXTable->jToY(i, j) << "\n"; + return false; + } } } - } - // check that the appicable range is correct. Note that due to rounding errors it is - // undefined whether the table applies to the boundary of the tabulated domain or not - Scalar xMin = uTable->xMin(); - Scalar yMin = uTable->yMin(); - Scalar xMax = uTable->xMax(); - Scalar yMax = uTable->yMax(); + // check that the appicable range is correct. Note that due to rounding errors it is + // undefined whether the table applies to the boundary of the tabulated domain or not + Scalar xMin = uTable->xMin(); + Scalar yMin = uTable->yMin(); + Scalar xMax = uTable->xMax(); + Scalar yMax = uTable->yMax(); - Scalar x = xMin - tolerance; - Scalar y = yMin - tolerance; - if (uTable->applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": uXTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<applies(x, y)) { - std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": !uTable->applies("<applies(x, y)) { + std::cerr << __FILE__ << ":" << __LINE__ << ": !uXTable->applies("<numX()*5; - unsigned n2 = uTable->numY()*5; - if (!compareTableWithAnalyticFn(uTable, - xMin, xMax, m2, - yMin, yMax, n2, - f, - tolerance)) - return false; - if (!compareTableWithAnalyticFn(uXTable, - xMin, xMax, m2, - yMin, yMax, n2, - f, - tolerance)) - return false; + // make sure that the function values at the sampling points are identical and that + // they correspond to the analytic function + unsigned m2 = uTable->numX()*5; + unsigned n2 = uTable->numY()*5; + if (!compareTableWithAnalyticFn(uTable, + xMin, xMax, m2, + yMin, yMax, n2, + f, + tolerance)) + return false; + if (!compareTableWithAnalyticFn(uXTable, + xMin, xMax, m2, + yMin, yMax, n2, + f, + tolerance)) + return false; - return true; -} + return true; + } }; template -inline int testAll( const typename TestType::Scalar tolerance = 1e-6 ) +inline int testAll(const typename TestType::Scalar tolerance = 1e-6) { TestType test; auto uniformTab = test.createUniformTabulatedFunction(TestType::testFn1); @@ -415,10 +414,10 @@ inline int testAll( const typename TestType::Scalar tolerance = 1e-6 ) uniformXTab = test.createUniformXTabulatedFunction2(TestType::testFn3); if (!test.compareTableWithAnalyticFn(uniformXTab, - -2.0, 3.0, 100, - -4.0, 5.0, 100, - TestType::testFn3, - /*tolerance=*/1e-2)) + -2.0, 3.0, 100, + -4.0, 5.0, 100, + TestType::testFn3, + /*tolerance=*/1e-2)) return 1; { @@ -434,19 +433,19 @@ inline int testAll( const typename TestType::Scalar tolerance = 1e-6 ) const unsigned n = 170; // extrapolation and interpolation involved, the tolerance needs to be bigger - const ScalarType temp_tolerance = 1000. * tolerance; + const ScalarType tmpTolerance = 1000. * tolerance; - if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn1, temp_tolerance)) + if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn1, tmpTolerance)) return 1; xytab = test.createXYTabulated2DFunction(TestType::testFn2); - if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn2, temp_tolerance)) + if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn2, tmpTolerance)) return 1; xytab = test.createXYTabulated2DFunction(TestType::testFn3); - if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn3, temp_tolerance)) + if (!test.compareTableWithAnalyticFn2(xytab, xMin, xMax, m, yMin, yMax, n, TestType::testFn3, tmpTolerance)) return 1; }