diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index b7c86534..e85e4fbe 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -139,6 +139,7 @@ list (APPEND MAIN_SOURCE_FILES # originally generated with the command: # find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort list (APPEND TEST_SOURCE_FILES + tests/test_spline.cpp tests/test_dgbasis.cpp tests/test_cartgrid.cpp tests/test_cubic.cpp @@ -331,12 +332,16 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/utility/Factory.hpp opm/core/utility/MonotCubicInterpolator.hpp opm/core/utility/NonuniformTableLinear.hpp + opm/core/utility/PolynomialUtils.hpp opm/core/utility/RootFinders.hpp opm/core/utility/SparseTable.hpp opm/core/utility/SparseVector.hpp + opm/core/utility/Spline.hpp opm/core/utility/StopWatch.hpp + opm/core/utility/TridiagonalMatrix.hpp opm/core/utility/UniformTableLinear.hpp opm/core/utility/Units.hpp + opm/core/utility/Unused.hpp opm/core/utility/VelocityInterpolation.hpp opm/core/utility/WachspressCoord.hpp opm/core/utility/buildUniformMonotoneTable.hpp diff --git a/opm/core/utility/PolynomialUtils.hpp b/opm/core/utility/PolynomialUtils.hpp new file mode 100644 index 00000000..e0338427 --- /dev/null +++ b/opm/core/utility/PolynomialUtils.hpp @@ -0,0 +1,310 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2013 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \brief Define some often used mathematical functions + */ +#ifndef OPM_MATH_HH +#define OPM_MATH_HH + +#include +#include + +namespace Opm +{ +/*! + * \ingroup Math + * \brief Invert a linear polynomial analytically + * + * The polynomial is defined as + * \f[ p(x) = a\; x + b \f] + * + * This method Returns the number of solutions which are in the real + * numbers, i.e. 1 except if the slope of the line is 0. + * + * \param sol Container into which the solutions are written + * \param a The coefficient for the linear term + * \param b The coefficient for the constant term + */ +template +int invertLinearPolynomial(SolContainer &sol, + Scalar a, + Scalar b) +{ + if (a == 0.0) + return 0; + + sol[0] = -b/a; + return 1; +} + +/*! + * \ingroup Math + * \brief Invert a quadratic polynomial analytically + * + * The polynomial is defined as + * \f[ p(x) = a\; x^2 + + b\;x + c \f] + * + * This method teturns the number of solutions which are in the real + * numbers. The "sol" argument contains the real roots of the parabola + * in order with the smallest root first. + * + * \param sol Container into which the solutions are written + * \param a The coefficient for the quadratic term + * \param b The coefficient for the linear term + * \param c The coefficient for the constant term + */ +template +int invertQuadraticPolynomial(SolContainer &sol, + Scalar a, + Scalar b, + Scalar c) +{ + // check for a line + if (a == 0.0) + return invertLinearPolynomial(sol, b, c); + + // discriminant + Scalar Delta = b*b - 4*a*c; + if (Delta < 0) + return 0; // no real roots + + Delta = std::sqrt(Delta); + sol[0] = (- b + Delta)/(2*a); + sol[1] = (- b - Delta)/(2*a); + + // sort the result + if (sol[0] > sol[1]) + std::swap(sol[0], sol[1]); + return 2; // two real roots +} + +//! \cond SKIP_THIS +template +void invertCubicPolynomialPostProcess_(SolContainer &sol, + int numSol, + Scalar a, + Scalar b, + Scalar c, + Scalar d) +{ + // do one Newton iteration on the analytic solution if the + // precision is increased + for (int i = 0; i < numSol; ++i) { + Scalar x = sol[i]; + Scalar fOld = d + x*(c + x*(b + x*a)); + + Scalar fPrime = c + x*(2*b + x*3*a); + if (fPrime == 0.0) + continue; + x -= fOld/fPrime; + + Scalar fNew = d + x*(c + x*(b + x*a)); + if (std::abs(fNew) < std::abs(fOld)) + sol[i] = x; + } +} +//! \endcond + +/*! + * \ingroup Math + * \brief Invert a cubic polynomial analytically + * + * The polynomial is defined as + * \f[ p(x) = a\; x^3 + + b\;x^3 + c\;x + d \f] + * + * This method teturns the number of solutions which are in the real + * numbers. The "sol" argument contains the real roots of the cubic + * polynomial in order with the smallest root first. + * + * \param sol Container into which the solutions are written + * \param a The coefficient for the cubic term + * \param b The coefficient for the quadratic term + * \param c The coefficient for the linear term + * \param d The coefficient for the constant term + */ +template +int invertCubicPolynomial(SolContainer *sol, + Scalar a, + Scalar b, + Scalar c, + Scalar d) +{ + // reduces to a quadratic polynomial + if (a == 0) + return invertQuadraticPolynomial(sol, b, c, d); + + // normalize the polynomial + b /= a; + c /= a; + d /= a; + a = 1; + + // get rid of the quadratic term by subsituting x = t - b/3 + Scalar p = c - b*b/3; + Scalar q = d + (2*b*b*b - 9*b*c)/27; + + if (p != 0.0 && q != 0.0) { + // At this point + // + // t^3 + p*t + q = 0 + // + // with p != 0 and q != 0 holds. Introducing the variables u and v + // with the properties + // + // u + v = t and 3*u*v + p = 0 + // + // leads to + // + // u^3 + v^3 + q = 0 . + // + // multiplying both sides with u^3 and taking advantage of the + // fact that u*v = -p/3 leads to + // + // u^6 + q*u^3 - p^3/27 = 0 + // + // Now, substituting u^3 = w yields + // + // w^2 + q*w - p^3/27 = 0 + // + // This is a quadratic equation with the solutions + // + // w = -q/2 + sqrt(q^2/4 + p^3/27) and + // w = -q/2 - sqrt(q^2/4 + p^3/27) + // + // Since w is equivalent to u^3 it is sufficient to only look at + // one of the two cases. Then, there are still 2 cases: positive + // and negative discriminant. + Scalar wDisc = q*q/4 + p*p*p/27; + if (wDisc >= 0) { // the positive discriminant case: + // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27) + Scalar u = - q/2 + std::sqrt(wDisc); + if (u < 0) u = - std::pow(-u, 1.0/3); + else u = std::pow(u, 1.0/3); + + // at this point, u != 0 since p^3 = 0 is necessary in order + // for u = 0 to hold, so + sol[0] = u - p/(3*u) - b/3; + // does not produce a division by zero. the remaining two + // roots of u are rotated by +- 2/3*pi in the complex plane + // and thus not considered here + invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d); + return 1; + } + else { // the negative discriminant case: + Scalar uCubedRe = - q/2; + Scalar uCubedIm = std::sqrt(-wDisc); + // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27) + Scalar uAbs = std::pow(std::sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm), 1.0/3); + Scalar phi = std::atan2(uCubedIm, uCubedRe)/3; + + // calculate the length and the angle of the primitive root + + // with the definitions from above it follows that + // + // x = u - p/(3*u) - b/3 + // + // where x and u are complex numbers. Rewritten in polar form + // this is equivalent to + // + // x = |u|*e^(i*phi) - p*e^(-i*phi)/(3*|u|) - b/3 . + // + // Factoring out the e^ terms and subtracting the additional + // terms, yields + // + // x = (e^(i*phi) + e^(-i*phi))*(|u| - p/(3*|u|)) - y - b/3 + // + // with + // + // y = - |u|*e^(-i*phi) + p*e^(i*phi)/(3*|u|) . + // + // The crucial observation is the fact that y is the conjugate + // of - x + b/3. This means that after taking advantage of the + // relation + // + // e^(i*phi) + e^(-i*phi) = 2*cos(phi) + // + // the equation + // + // x = 2*cos(phi)*(|u| - p / (3*|u|)) - conj(x) - 2*b/3 + // + // holds. Since |u|, p, b and cos(phi) are real numbers, it + // follows that Im(x) = - Im(x) and thus Im(x) = 0. This + // implies + // + // Re(x) = x = cos(phi)*(|u| - p / (3*|u|)) - b/3 . + // + // Considering the fact that u is a cubic root, we have three + // values for phi which differ by 2/3*pi. This allows to + // calculate the three real roots of the polynomial: + for (int i = 0; i < 3; ++i) { + sol[i] = std::cos(phi)*(uAbs - p/(3*uAbs)) - b/3; + phi += 2*M_PI/3; + } + + // post process the obtained solution to increase numerical + // precision + invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d); + + // sort the result + std::sort(sol, sol + 3); + + return 3; + } + } + // Handle some (pretty unlikely) special cases required to avoid + // divisions by zero in the code above... + else if (p == 0.0 && q == 0.0) { + // t^3 = 0, i.e. triple root at t = 0 + sol[0] = sol[1] = sol[2] = 0.0 - b/3; + return 3; + } + else if (p == 0.0 && q != 0.0) { + // t^3 + q = 0, + // + // i. e. single real root at t=curt(q) + Scalar t; + if (-q > 0) t = std::pow(-q, 1./3); + else t = - std::pow(q, 1./3); + sol[0] = t - b/3; + + return 1; + } + + assert(p != 0.0 && q == 0.0); + + // t^3 + p*t = 0 = t*(t^2 + p), + // + // i. e. roots at t = 0, t^2 + p = 0 + if (p > 0) { + sol[0] = 0.0 - b/3; + return 1; // only a single real root at t=0 + } + + // two additional real roots at t = sqrt(-p) and t = -sqrt(-p) + sol[0] = -std::sqrt(-p) - b/3;; + sol[1] = 0.0 - b/3; + sol[2] = std::sqrt(-p) - b/3; + + return 3; +} +} + + +#endif diff --git a/opm/core/utility/Spline.hpp b/opm/core/utility/Spline.hpp new file mode 100644 index 00000000..b3e07882 --- /dev/null +++ b/opm/core/utility/Spline.hpp @@ -0,0 +1,1703 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2013 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \copydoc Opm::Spline + */ +#ifndef OPM_SPLINE_HH +#define OPM_SPLINE_HH + +#include +#include +#include + +#include +#include + +namespace Opm +{ +/*! + * \brief Class implementing cubic splines. + * + * This class supports full, natural, periodic and monotonic cubic + * splines. + * + * Full a splines \f$s(x)\f$ are splines which, given \f$n\f$ sampling + * points \f$x_1, \dots, x_n\f$, fulfill the following conditions + * \f{align*}{ + * s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \} \\ + * s'(x_1) & = m_1 \\ + * s'(x_n) & = m_n + * \f} + * for any given boundary slopes \f$m_1\f$ and \f$m_n\f$. + * + * Natural splines which are defined by + *\f{align*}{ + * s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \} \\ + * s''(x_1) & = 0 \\ + * s''(x_n) & = 0 + *\f} + * + * For periodic splines of splines the slopes at the boundaries are identical: + *\f{align*}{ + * s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \} \\ + * s'(x_1) & = s'(x_n) \\ + * s''(x_1) & = s''(x_n) \;. + *\f} + * + * Finally, there are monotonic splines which guarantee that the curve + * is confined by its sampling points, i.e., + * \f[ + y_i \leq s(x) \leq y_{i+1} \quad \text{for} x_i \leq x \leq x_{i+1} \;. + * \f] + * For more information on monotonic splines, see + * http://en.wikipedia.org/wiki/Monotone_cubic_interpolation + * + * Full, natural and periodic splines are continuous in their first + * and second derivatives, i.e., + * \f[ + s \in \mathcal{C}^2 + * \f] + * holds for such splines. Monotonic splines are only continuous up to + * their first derivative, i.e., for these only + * \f[ + s \in \mathcal{C}^1 + * \f] + * is true. + */ +template +class Spline +{ + typedef Opm::TridiagonalMatrix Matrix; + typedef std::vector Vector; + +public: + /*! + * \brief The type of the spline to be created. + * + * \copydetails Spline + */ + enum SplineType { + Full, + Natural, + Periodic, + Monotonic + }; + + /*! + * \brief Default constructor for a spline. + * + * To specfiy the acutal curve, use one of the set() methods. + */ + Spline() + { } + + /*! + * \brief Convenience constructor for a full spline with just two sampling points + * + * \param x0 The \f$x\f$ value of the first sampling point + * \param x1 The \f$x\f$ value of the second sampling point + * \param y0 The \f$y\f$ value of the first sampling point + * \param y1 The \f$y\f$ value of the second sampling point + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_1\f$ + */ + Spline(Scalar x0, Scalar x1, + Scalar y0, Scalar y1, + Scalar m0, Scalar m1) + { set(x0, x1, y0, y1, m0, m1); } + + /*! + * \brief Convenience constructor for a natural or a periodic spline + * + * \param nSamples The number of sampling points (must be > 2) + * \param x An array containing the \f$x\f$ values of the spline's sampling points + * \param y An array containing the \f$y\f$ values of the spline's sampling points + * \param periodic Indicates whether a natural or a periodic spline should be created + */ + template + Spline(int nSamples, + const ScalarArrayX &x, + const ScalarArrayY &y, + SplineType splineType = Natural, + bool sortInputs = false) + { this->setXYArrays(nSamples, x, y, splineType, sortInputs); } + + /*! + * \brief Convenience constructor for a natural or a periodic spline + * + * \param nSamples The number of sampling points (must be > 2) + * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points + * \param periodic Indicates whether a natural or a periodic spline should be created + */ + template + Spline(int nSamples, + const PointArray &points, + SplineType splineType = Natural, + bool sortInputs = false) + { this->setArrayOfPoints(nSamples, points, splineType, sortInputs); } + + /*! + * \brief Convenience constructor for a natural or a periodic spline + * + * \param x An array containing the \f$x\f$ values of the spline's sampling points (must have a size() method) + * \param y An array containing the \f$y\f$ values of the spline's sampling points (must have a size() method) + * \param periodic Indicates whether a natural or a periodic spline should be created + */ + template + Spline(const ScalarContainer &x, + const ScalarContainer &y, + SplineType splineType = Natural, + bool sortInputs = false) + { this->setXYContainers(x, y, splineType, sortInputs); } + + /*! + * \brief Convenience constructor for a natural or a periodic spline + * + * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points (must have a size() method) + * \param periodic Indicates whether a natural or a periodic spline should be created + */ + template + Spline(const PointContainer &points, + SplineType splineType = Natural, + bool sortInputs = false) + { this->setContainerOfPoints(points, splineType, sortInputs); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param nSamples The number of sampling points (must be >= 2) + * \param x An array containing the \f$x\f$ values of the spline's sampling points + * \param y An array containing the \f$y\f$ values of the spline's sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + * \param sortInputs Indicates whether the sample points should be sorted (this is not necessary if they are already sorted in ascending or descending order) + */ + template + Spline(int nSamples, + const ScalarArray &x, + const ScalarArray &y, + Scalar m0, + Scalar m1, + bool sortInputs = false) + { this->setXYArrays(nSamples, x, y, m0, m1, sortInputs); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param nSamples The number of sampling points (must be >= 2) + * \param points An array containing the \f$x\f$ and \f$x\f$ values of the spline's sampling points + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + * \param sortInputs Indicates whether the sample points should be sorted (this is not necessary if they are already sorted in ascending or descending order) + */ + template + Spline(int nSamples, + const PointArray &points, + Scalar m0, + Scalar m1, + bool sortInputs = false) + { this->setArrayOfPoints(nSamples, points, m0, m1, sortInputs); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param x An array containing the \f$x\f$ values of the spline's sampling points (must have a size() method) + * \param y An array containing the \f$y\f$ values of the spline's sampling points (must have a size() method) + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + * \param sortInputs Indicates whether the sample points should be sorted (this is not necessary if they are already sorted in ascending or descending order) + */ + template + Spline(const ScalarContainerX &x, + const ScalarContainerY &y, + Scalar m0, + Scalar m1, + bool sortInputs = false) + { this->setXYContainers(x, y, m0, m1, sortInputs); } + + /*! + * \brief Convenience constructor for a full spline + * + * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points (must have a size() method) + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_n\f$ + * \param sortInputs Indicates whether the sample points should be sorted (this is not necessary if they are already sorted in ascending or descending order) + */ + template + Spline(const PointContainer &points, + Scalar m0, + Scalar m1, + bool sortInputs = false) + { this->setContainerOfPoints(points, m0, m1, sortInputs); } + + /*! + * \brief Returns the number of sampling points. + */ + int numSamples() const + { return xPos_.size(); } + + /*! + * \brief Set the sampling points and the boundary slopes of the + * spline with two sampling points. + * + * \param x0 The \f$x\f$ value of the first sampling point + * \param x1 The \f$x\f$ value of the second sampling point + * \param y0 The \f$y\f$ value of the first sampling point + * \param y1 The \f$y\f$ value of the second sampling point + * \param m0 The slope of the spline at \f$x_0\f$ + * \param m1 The slope of the spline at \f$x_1\f$ + */ + void set(Scalar x0, Scalar x1, + Scalar y0, Scalar y1, + Scalar m0, Scalar m1) + { + slopeVec_.resize(2); + xPos_.resize(2); + yPos_.resize(2); + + if (x0 > x1) { + xPos_[0] = x1; + xPos_[1] = x0; + yPos_[0] = y1; + yPos_[1] = y0; + } + else { + xPos_[0] = x0; + xPos_[1] = x1; + yPos_[0] = y0; + yPos_[1] = y1; + } + + slopeVec_[0] = m0; + slopeVec_[1] = m1; + + Matrix M(numSamples()); + Vector d(numSamples()); + Vector moments(numSamples()); + this->makeFullSystem_(M, d, m0, m1); + + // solve for the moments + M.solve(moments, d); + + this->setSlopesFromMoments_(slopeVec_, moments); + } + + + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + // Full splines // + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using C-style arrays. + * + * This method uses separate array-like objects for the values of + * the X and Y coordinates. In this context 'array-like' means + * that an access to the members is provided via the [] + * operator. (e.g. C arrays, std::vector, std::array, etc.) Each + * array must be of size 'nSamples' at least. Also, the number of + * sampling points must be larger than 1. + */ + template + void setXYArrays(int nSamples, + const ScalarArrayX &x, + const ScalarArrayY &y, + Scalar m0, Scalar m1, + bool sortInputs = false) + { + assert(nSamples > 1); + + setNumSamples_(nSamples); + for (int i = 0; i < nSamples; ++i) { + xPos_[i] = x[i]; + yPos_[i] = y[i]; + } + + if (sortInputs) + sortInput_(); + else if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using STL-compatible containers. + * + * This method uses separate STL-compatible containers for the + * values of the sampling points' X and Y + * coordinates. "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setXYContainers(const ScalarContainerX &x, + const ScalarContainerY &y, + Scalar m0, Scalar m1, + bool sortInputs = false) + { + assert(x.size() == y.size()); + assert(x.size() > 1); + + setNumSamples_(x.size()); + + std::copy(x.begin(), x.end(), xPos_.begin()); + std::copy(y.begin(), y.end(), yPos_.begin()); + + if (sortInputs) + sortInput_(); + else if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using a C-style array. + * + * This method uses a single array of sampling points, which are + * seen as an array-like object which provides access to the X and + * Y coordinates. In this context 'array-like' means that an + * access to the members is provided via the [] operator. (e.g. C + * arrays, std::vector, std::array, etc.) The array containing + * the sampling points must be of size 'nSamples' at least. Also, + * the number of sampling points must be larger than 1. + */ + template + void setArrayOfPoints(int nSamples, + const PointArray &points, + Scalar m0, + Scalar m1, + bool sortInputs = false) + { + // a spline with no or just one sampling points? what an + // incredible bad idea! + assert(nSamples > 1); + + setNumSamples_(nSamples); + for (int i = 0; i < nSamples; ++i) { + xPos_[i] = points[i][0]; + yPos_[i] = points[i][1]; + } + + if (sortInputs) + sortInput_(); + else if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using a STL-compatible container of + * array-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be array-like objects storing the + * X and Y coordinates. "STL-compatible" means that the container + * provides access to iterators using the begin(), end() methods + * and also provides a size() method. Also, the number of entries + * in the X and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfPoints(const XYContainer &points, + Scalar m0, + Scalar m1, + bool sortInputs = false) + { + // a spline with no or just one sampling points? what an + // incredible bad idea! + assert(points.size() > 1); + + setNumSamples_(points.size()); + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++i, ++it) { + xPos_[i] = (*it)[0]; + yPos_[i] = (*it)[1]; + } + + if (sortInputs) + sortInput_(); + else if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + // make a full spline + makeFullSpline_(m0, m1); + } + + /*! + * \brief Set the sampling points and the boundary slopes of a + * full spline using a STL-compatible container of + * tuple-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be tuple-like objects storing the + * X and Y coordinates. "tuple-like" means that the objects + * provide access to the x values via std::get<0>(obj) and to the + * y value via std::get<1>(obj) (e.g. std::tuple or + * std::pair). "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfTuples(const XYContainer &points, + Scalar m0, + Scalar m1, + bool sortInputs = false) + { + // resize internal arrays + setNumSamples_(points.size()); + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++i, ++it) { + xPos_[i] = std::get<0>(*it); + yPos_[i] = std::get<1>(*it); + } + + if (sortInputs) + sortInput_(); + else if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + // make a full spline + makeFullSpline_(m0, m1); + } + + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + // Natural/Periodic splines // + /////////////////////////////////////// + /////////////////////////////////////// + /////////////////////////////////////// + /*! + * \brief Set the sampling points natural spline using C-style arrays. + * + * This method uses separate array-like objects for the values of + * the X and Y coordinates. In this context 'array-like' means + * that an access to the members is provided via the [] + * operator. (e.g. C arrays, std::vector, std::array, etc.) Each + * array must be of size 'nSamples' at least. Also, the number of + * sampling points must be larger than 1. + */ + template + void setXYArrays(int nSamples, + const ScalarArrayX &x, + const ScalarArrayY &y, + SplineType splineType = Natural, + bool sortInputs = false) + { + assert(nSamples > 1); + + setNumSamples_(nSamples); + for (int i = 0; i < nSamples; ++i) { + xPos_[i] = x[i]; + yPos_[i] = y[i]; + } + + if (sortInputs) + sortInput_(); + else if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + if (splineType == Periodic) + makePeriodicSpline_(); + else if (splineType == Natural) + makeNaturalSpline_(); + else if (splineType == Monotonic) + this->makeMonotonicSpline_(slopeVec_); + else + OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place"); + } + + /*! + * \brief Set the sampling points of a natural spline using + * STL-compatible containers. + * + * This method uses separate STL-compatible containers for the + * values of the sampling points' X and Y + * coordinates. "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setXYContainers(const ScalarContainerX &x, + const ScalarContainerY &y, + SplineType splineType = Natural, + bool sortInputs = false) + { + assert(x.size() == y.size()); + assert(x.size() > 1); + + setNumSamples_(x.size()); + std::copy(x.begin(), x.end(), xPos_.begin()); + std::copy(y.begin(), y.end(), yPos_.begin()); + + if (sortInputs) + sortInput_(); + else if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + if (splineType == Periodic) + makePeriodicSpline_(); + else if (splineType == Natural) + makeNaturalSpline_(); + else if (splineType == Monotonic) + this->makeMonotonicSpline_(slopeVec_); + else + OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place"); + } + + /*! + * \brief Set the sampling points of a natural spline using a + * C-style array. + * + * This method uses a single array of sampling points, which are + * seen as an array-like object which provides access to the X and + * Y coordinates. In this context 'array-like' means that an + * access to the members is provided via the [] operator. (e.g. C + * arrays, std::vector, std::array, etc.) The array containing + * the sampling points must be of size 'nSamples' at least. Also, + * the number of sampling points must be larger than 1. + */ + template + void setArrayOfPoints(int nSamples, + const PointArray &points, + SplineType splineType = Natural, + bool sortInputs = false) + { + // a spline with no or just one sampling points? what an + // incredible bad idea! + assert(nSamples > 1); + + setNumSamples_(nSamples); + for (int i = 0; i < nSamples; ++i) { + xPos_[i] = points[i][0]; + yPos_[i] = points[i][1]; + } + + if (sortInputs) + sortInput_(); + else if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + if (splineType == Periodic) + makePeriodicSpline_(); + else if (splineType == Natural) + makeNaturalSpline_(); + else if (splineType == Monotonic) + this->makeMonotonicSpline_(slopeVec_); + else + OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place"); + } + + /*! + * \brief Set the sampling points of a natural spline using a + * STL-compatible container of array-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be array-like objects storing the + * X and Y coordinates. "STL-compatible" means that the container + * provides access to iterators using the begin(), end() methods + * and also provides a size() method. Also, the number of entries + * in the X and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfPoints(const XYContainer &points, + SplineType splineType = Natural, + bool sortInputs = false) + { + // a spline with no or just one sampling points? what an + // incredible bad idea! + assert(points.size() > 1); + + setNumSamples_(points.size()); + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++ i, ++it) { + xPos_[i] = (*it)[0]; + yPos_[i] = (*it)[1]; + } + + if (sortInputs) + sortInput_(); + else if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + if (splineType == Periodic) + makePeriodicSpline_(); + else if (splineType == Natural) + makeNaturalSpline_(); + else if (splineType == Monotonic) + this->makeMonotonicSpline_(slopeVec_); + else + OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place"); + } + + /*! + * \brief Set the sampling points of a natural spline using a + * STL-compatible container of tuple-like objects. + * + * This method uses a single STL-compatible container of sampling + * points, which are assumed to be tuple-like objects storing the + * X and Y coordinates. "tuple-like" means that the objects + * provide access to the x values via std::get<0>(obj) and to the + * y value via std::get<1>(obj) (e.g. std::tuple or + * std::pair). "STL-compatible" means that the container provides + * access to iterators using the begin(), end() methods and also + * provides a size() method. Also, the number of entries in the X + * and the Y containers must be equal and larger than 1. + */ + template + void setContainerOfTuples(const XYContainer &points, + SplineType splineType = Natural, + bool sortInputs = false) + { + // resize internal arrays + setNumSamples_(points.size()); + typename XYContainer::const_iterator it = points.begin(); + typename XYContainer::const_iterator endIt = points.end(); + for (int i = 0; it != endIt; ++i, ++it) { + xPos_[i] = std::get<0>(*it); + yPos_[i] = std::get<1>(*it); + } + + if (sortInputs) + sortInput_(); + else if (xPos_[0] > xPos_[numSamples() - 1]) + reverseSamplingPoints_(); + + if (splineType == Periodic) + makePeriodicSpline_(); + else if (splineType == Natural) + makeNaturalSpline_(); + else if (splineType == Monotonic) + this->makeMonotonicSpline_(slopeVec_); + else + OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place"); + } + + /*! + * \brief Return true iff the given x is in range [x1, xn]. + */ + bool applies(Scalar x) const + { + return x_(0) <= x && x <= x_(numSamples() - 1); + } + + /*! + * \brief Return the x value of the leftmost sampling point. + */ + Scalar xMin() const + { return x_(0); } + + /*! + * \brief Return the x value of the rightmost sampling point. + */ + Scalar xMax() const + { return x_(numSamples() - 1); } + + /*! + * \brief Prints k tuples of the format (x, y, dx/dy, isMonotonic) + * to stdout. + * + * If the spline does not apply for parts of [x0, x1] it is + * extrapolated using a straight line. The result can be inspected + * using the following commands: + * + ----------- snip ----------- + ./yourProgramm > spline.csv + gnuplot + + gnuplot> plot "spline.csv" using 1:2 w l ti "Curve", \ + "spline.csv" using 1:3 w l ti "Derivative", \ + "spline.csv" using 1:4 w p ti "Monotonic" + ----------- snap ----------- + */ + void printCSV(Scalar xi0, Scalar xi1, int k, std::ostream &os = std::cout) const + { + Scalar x0 = std::min(xi0, xi1); + Scalar x1 = std::max(xi0, xi1); + const int n = numSamples() - 1; + for (int i = 0; i <= k; ++i) { + double x = i*(x1 - x0)/k + x0; + double x_p1 = x + (x1 - x0)/k; + double y; + double dy_dx; + double mono = 1; + if (!applies(x)) { + if (x < x_(0)) { + dy_dx = evalDerivative(x_(0)); + y = (x - x_(0))*dy_dx + y_(0); + mono = (dy_dx>0)?1:-1; + } + else if (x > x_(n)) { + dy_dx = evalDerivative(x_(n)); + y = (x - x_(n))*dy_dx + y_(n); + mono = (dy_dx>0)?1:-1; + } + else { + OPM_THROW(std::runtime_error, + "The sampling points given to a spline must be sorted by their x value!"); + } + } + else { + y = eval(x); + dy_dx = evalDerivative(x); + mono = monotonic(std::max(x_(0), x), std::min(x_(n), x_p1)); + } + + os << x << " " << y << " " << dy_dx << " " << mono << "\n"; + } + } + + /*! + * \brief Evaluate the spline at a given position. + * + * \param x The value on the abscissa where the spline ought to be + * evaluated + * \param extrapolate If this parameter is set to true, the spline + * will be extended beyond its range by + * straight lines, if false calling extrapolate + * for \f$ x \not [x_{min}, x_{max}]\f$ will + * cause a failed assertation. + */ + Scalar eval(Scalar x, bool extrapolate=false) const + { + assert(extrapolate || applies(x)); + + // handle extrapolation + if (extrapolate) { + if (x < xMin()) { + Scalar m = evalDerivative(xMin(), /*segmentIdx=*/0); + Scalar y0 = y_(0); + return y0 + m*(x - xMin()); + } + else if (x > xMax()) { + Scalar m = evalDerivative(xMax(), /*segmentIdx=*/numSamples()-2); + Scalar y0 = y_(numSamples() - 1); + return y0 + m*(x - xMax()); + } + } + + return eval_(x, segmentIdx_(x)); + } + + /*! + * \brief Evaluate the spline's derivative at a given position. + * + * \param x The value on the abscissa where the spline's + * derivative ought to be evaluated + * + * \param extrapolate If this parameter is set to true, the spline + * will be extended beyond its range by + * straight lines, if false calling extrapolate + * for \f$ x \not [x_{min}, x_{max}]\f$ will + * cause a failed assertation. + */ + Scalar evalDerivative(Scalar x, bool extrapolate=false) const + { + assert(extrapolate || applies(x)); + if (extrapolate) { + if (x < xMin()) + evalDerivative_(xMin(), 0); + else if (x > xMax()) + evalDerivative_(xMax(), numSamples() - 2); + } + + return evalDerivative_(x, segmentIdx_(x)); + } + + /*! + * \brief Evaluate the spline's second derivative at a given position. + * + * \param x The value on the abscissa where the spline's + * derivative ought to be evaluated + * + * \param extrapolate If this parameter is set to true, the spline + * will be extended beyond its range by + * straight lines, if false calling extrapolate + * for \f$ x \not [x_{min}, x_{max}]\f$ will + * cause a failed assertation. + */ + Scalar evalSecondDerivative(Scalar x, bool extrapolate=false) const + { + assert(extrapolate || applies(x)); + if (extrapolate) + return 0.0; + + return evalDerivative2_(x, segmentIdx_(x)); + } + + /*! + * \brief Find the intersections of the spline with a cubic + * polynomial in the whole intervall, throws + * Opm::MathError exception if there is more or less than + * one solution. + */ + Scalar intersect(Scalar a, Scalar b, Scalar c, Scalar d) const + { + return intersectIntervall(xMin(), xMax(), a, b, c, d); + } + + /*! + * \brief Find the intersections of the spline with a cubic + * polynomial in a sub-intervall of the spline, throws + * Opm::MathError exception if there is more or less than + * one solution. + */ + Scalar intersectInterval(Scalar x0, Scalar x1, + Scalar a, Scalar b, Scalar c, Scalar d) const + { + assert(applies(x0) && applies(x1)); + + Scalar tmpSol[3]; + int nSol = 0; + int iFirst = segmentIdx_(x0); + int iLast = segmentIdx_(x1); + for (int i = iFirst; i <= iLast; ++i) + { + nSol += intersectSegment_(tmpSol, i, a, b, c, d, x0, x1); + + if (nSol > 1) { + OPM_THROW(std::runtime_error, + "Spline has more than one intersection"); //< x1) + std::swap(x0, x1); + + assert(x0 < x1); + + int i = segmentIdx_(x0); + if (x_(i + 1) < x1) + // interval is fully contained within a single spline + // segment + return monotonic_(i, x0, x1); + + int iEnd = segmentIdx_(x1); + + // make sure that the segments which are completly in the + // interval [x0, x1] all exhibit the same monotonicity. + int r = monotonic_(i, x0, x_(i + 1)); + for (++i; i < iEnd - 1; ++i) { + int nextR = monotonic_(i, x_(i), x_(i + 1)); + if (nextR == 3) // spline is constant + continue; + if (r == 3) + r = nextR; + if (r != nextR) + return 0; + } + + // check for the last segment + if (x_(iEnd) < x1) { + int lastR = monotonic_(iEnd, x_(iEnd), x1); + if (lastR != 3 && r != 3 && r != lastR) + return 0; + } + + return r; + } + + /*! + * \brief Same as monotonic(x0, x1), but with the entire range of the + * spline as interval. + */ + int monotonic() const + { return monotonic(xMin(), xMax()); } + +protected: + /*! + * \brief Helper class needed to sort the input sampling points. + */ + struct ComparatorX_ + { + ComparatorX_(const std::vector &x) + : x_(x) + {}; + + bool operator ()(int idxA, int idxB) const + { return x_.at(idxA) < x_.at(idxB); } + + const std::vector &x_; + }; + + /*! + * \brief Sort the sample points in ascending order of their x value. + */ + void sortInput_() + { + int n = numSamples(); + + // create a vector containing 0...n-1 + std::vector idxVector(n); + for (int i = 0; i < n; ++i) + idxVector[i] = i; + + // sort the indices according to the x values of the sample + // points + ComparatorX_ cmp(xPos_); + std::sort(idxVector.begin(), idxVector.end(), cmp); + + // reorder the sample points + std::vector tmpX(n), tmpY(n); + for (int i = 0; i < idxVector.size(); ++ i) { + tmpX[i] = xPos_[idxVector[i]]; + tmpY[i] = yPos_[idxVector[i]]; + } + xPos_ = tmpX; + yPos_ = tmpY; + } + + /*! + * \brief Reverse order of the elements in the arrays which + * contain the sampling points. + */ + void reverseSamplingPoints_() + { + // reverse the arrays + int n = numSamples(); + for (int i = 0; i <= (n - 1)/2; ++i) { + std::swap(xPos_[i], xPos_[n - i - 1]); + std::swap(yPos_[i], yPos_[n - i - 1]); + } + } + + /*! + * \brief Resizes the internal vectors to store the sample points. + */ + void setNumSamples_(int nSamples) + { + xPos_.resize(nSamples); + yPos_.resize(nSamples); + slopeVec_.resize(nSamples); + } + + /*! + * \brief Create a natural spline from the already set sampling points. + * + * This creates a temporary matrix and right hand side vector. + */ + void makeFullSpline_(Scalar m0, Scalar m1) + { + Matrix M(numSamples()); + std::vector d(numSamples()); + std::vector moments(numSamples()); + + // create linear system of equations + this->makeFullSystem_(M, d, m0, m1); + + // solve for the moments (-> second derivatives) + M.solve(moments, d); + + // convert the moments to slopes at the sample points + this->setSlopesFromMoments_(slopeVec_, moments); + } + + /*! + * \brief Create a natural spline from the already set sampling points. + * + * This creates a temporary matrix and right hand side vector. + */ + void makeNaturalSpline_() + { + Matrix M(numSamples(), numSamples()); + Vector d(numSamples()); + Vector moments(numSamples()); + + // create linear system of equations + this->makeNaturalSystem_(M, d); + + // solve for the moments (-> second derivatives) + M.solve(moments, d); + + // convert the moments to slopes at the sample points + this->setSlopesFromMoments_(slopeVec_, moments); + } + + /*! + * \brief Create a periodic spline from the already set sampling points. + * + * This creates a temporary matrix and right hand side vector. + */ + void makePeriodicSpline_() + { + Matrix M(numSamples() - 1); + Vector d(numSamples() - 1); + Vector moments(numSamples() - 1); + + // create linear system of equations. This is a bit hacky, + // because it assumes that std::vector internally stores its + // data as a big C-style array, but it saves us from yet + // another copy operation + this->makePeriodicSystem_(M, d); + + // solve for the moments (-> second derivatives) + M.solve(moments, d); + + moments.resize(numSamples()); + for (int i = numSamples() - 2; i >= 0; --i) + moments[i+1] = moments[i]; + moments[0] = moments[numSamples() - 1]; + + // convert the moments to slopes at the sample points + this->setSlopesFromMoments_(slopeVec_, moments); + } + + /*! + * \brief Set the sampling point vectors. + * + * This takes care that the order of the x-values is ascending, + * although the input must be ordered already! + */ + template + void assignSamplingPoints_(DestVector &destX, + DestVector &destY, + const SourceVector &srcX, + const SourceVector &srcY, + int numSamples) + { + assert(numSamples >= 2); + + // copy sample points, make sure that the first x value is + // smaller than the last one + for (int i = 0; i < numSamples; ++i) { + int idx = i; + if (srcX[0] > srcX[numSamples - 1]) + idx = numSamples - i - 1; + destX[i] = srcX[idx]; + destY[i] = srcY[idx]; + } + } + + template + void assignFromArrayList_(DestVector &destX, + DestVector &destY, + const ListIterator &srcBegin, + const ListIterator &srcEnd, + int numSamples) + { + assert(numSamples >= 2); + + // find out wether the x values are in reverse order + ListIterator it = srcBegin; + ++it; + bool reverse = false; + if ((*srcBegin)[0] > (*it)[0]) + reverse = true; + --it; + + // loop over all sampling points + for (int i = 0; it != srcEnd; ++i, ++it) { + int idx = i; + if (reverse) + idx = numSamples - i - 1; + destX[i] = (*it)[0]; + destY[i] = (*it)[1]; + } + } + + /*! + * \brief Set the sampling points. + * + * Here we assume that the elements of the source vector have an + * [] operator where v[0] is the x value and v[1] is the y value + * if the sampling point. + */ + template + void assignFromTupleList_(DestVector &destX, + DestVector &destY, + ListIterator srcBegin, + ListIterator srcEnd, + int numSamples) + { + assert(numSamples >= 2); + + // copy sample points, make sure that the first x value is + // smaller than the last one + + // find out wether the x values are in reverse order + ListIterator it = srcBegin; + ++it; + bool reverse = false; + if (std::get<0>(*srcBegin) > std::get<0>(*it)) + reverse = true; + --it; + + // loop over all sampling points + for (int i = 0; it != srcEnd; ++i, ++it) { + int idx = i; + if (reverse) + idx = numSamples - i - 1; + destX[i] = std::get<0>(*it); + destY[i] = std::get<1>(*it); + } + } + + /*! + * \brief Make the linear system of equations Mx = d which results + * in the moments of the full spline. + */ + template + void makeFullSystem_(Matrix &M, Vector &d, Scalar m0, Scalar m1) + { + makeNaturalSystem_(M, d); + + int n = numSamples() - 1; + // first row + M[0][1] = 1; + d[0] = 6/h_(1) * ( (y_(1) - y_(0))/h_(1) - m0); + + // last row + M[n][n - 1] = 1; + + // right hand side + d[n] = + 6/h_(n) + * + (m1 - (y_(n) - y_(n - 1))/h_(n)); + } + + /*! + * \brief Make the linear system of equations Mx = d which results + * in the moments of the natural spline. + */ + template + void makeNaturalSystem_(Matrix &M, Vector &d) + { + M = 0.0; + + // See: J. Stoer: "Numerische Mathematik 1", 9th edition, + // Springer, 2005, p. 111 + const int n = numSamples() - 1; + + // second to next to last rows + for (int i = 1; i < n; ++i) { + Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1)); + Scalar mu_i = 1 - lambda_i; + Scalar d_i = + 6 / (h_(i) + h_(i + 1)) + * + ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i)); + + M[i][i-1] = mu_i; + M[i][i] = 2; + M[i][i + 1] = lambda_i; + d[i] = d_i; + }; + + // See Stroer, equation (2.5.2.7) + Scalar lambda_0 = 0; + Scalar d_0 = 0; + + Scalar mu_n = 0; + Scalar d_n = 0; + + // first row + M[0][0] = 2; + M[0][1] = lambda_0; + d[0] = d_0; + + // last row + M[n][n-1] = mu_n; + M[n][n] = 2; + d[n] = d_n; + } + + /*! + * \brief Make the linear system of equations Mx = d which results + * in the moments of the periodic spline. + */ + template + void makePeriodicSystem_(Matrix &M, Vector &d) + { + M = 0.0; + + // See: J. Stoer: "Numerische Mathematik 1", 9th edition, + // Springer, 2005, p. 111 + const int n = numSamples() - 1; + + assert(M.rows() == n); + + // second to next to last rows + for (int i = 2; i < n; ++i) { + Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i + 1)); + Scalar mu_i = 1 - lambda_i; + Scalar d_i = + 6 / (h_(i) + h_(i + 1)) + * + ( (y_(i + 1) - y_(i))/h_(i + 1) - (y_(i) - y_(i - 1))/h_(i)); + + M[i-1][i-2] = mu_i; + M[i-1][i-1] = 2; + M[i-1][i] = lambda_i; + d[i-1] = d_i; + }; + + Scalar lambda_n = h_(1) / (h_(n) + h_(1)); + Scalar lambda_1 = h_(2) / (h_(1) + h_(2));; + Scalar mu_1 = 1 - lambda_1; + Scalar mu_n = 1 - lambda_n; + + Scalar d_1 = + 6 / (h_(1) + h_(2)) + * + ( (y_(2) - y_(1))/h_(2) - (y_(1) - y_(0))/h_(1)); + Scalar d_n = + 6 / (h_(n) + h_(1)) + * + ( (y_(1) - y_(n))/h_(1) - (y_(n) - y_(n-1))/h_(n)); + + + // first row + M[0][0] = 2; + M[0][1] = lambda_1; + M[0][n-1] = mu_1; + d[0] = d_1; + + // last row + M[n-1][0] = lambda_n; + M[n-1][n-2] = mu_n; + M[n-1][n-1] = 2; + d[n-1] = d_n; + } + + /*! + * \brief Create a monotonic spline from the already set sampling points. + * + * This code is inspired by opm-core's "MonotCubicInterpolator" + * class and also uses the approach by Fritsch and Carlson, see + * + * http://en.wikipedia.org/wiki/Monotone_cubic_interpolation + */ + template + void makeMonotonicSpline_(Vector &slopes) + { + auto n = numSamples(); + std::vector delta(n); + for (int k = 0; k < n - 1; ++k) { + delta[k] = (y_(k + 1) - y_(k))/(x_(k + 1) - x_(k)); + } + delta[n - 1] = delta[n - 2]; + + // calculate the "raw" slopes at the sample points + for (int k = 0; k < n - 1; ++k) + slopes[k] = (delta[k] + delta[k + 1])/2; + slopes[n - 1] = delta[n - 2]; + + // post-process the "raw" slopes at the sample points + for (int k = 0; k < n - 1; ++k) { + if (std::abs(delta[k]) < 1e-20) { + // make the spline flat if the inputs are equal + slopes[k] = 0; + slopes[k + 1] = 0; + ++ k; + continue; + } + + Scalar alpha = slopes[k] / delta[k]; + Scalar beta = slopes[k + 1] / delta[k]; + + if (k > 0) { + // check if the inputs are not montonous. if yes, make + // x[k] a local extremum. + if (delta[k]*delta[k - 1] < 0) { + slopes[k] = 0; + continue; + } + } + + // limit (alpha, beta) to a circle of radius 3 + if (alpha*alpha + beta*beta > 3*3) { + Scalar tau = 3.0/std::sqrt(alpha*alpha + beta*beta); + slopes[k] = tau*alpha*delta[k]; + slopes[k + 1] = tau*beta*delta[k]; + } + } + } + + /*! + * \brief Convert the moments at the sample points to slopes. + * + * This requires to use cubic Hermite interpolation, but it is + * required because for monotonic splines the second derivative is + * not continuous. + */ + template + void setSlopesFromMoments_(SlopeVector &slopes, const MomentsVector &moments) + { + int n = numSamples(); + + // evaluate slope at the rightmost point. + // See: J. Stoer: "Numerische Mathematik 1", 9th edition, + // Springer, 2005, p. 109 + Scalar mRight; + + { + Scalar h = this->h_(n - 1); + Scalar x = h; + //Scalar x_1 = 0; + + Scalar A = + (y_(n - 1) - y_(n - 2))/h + - + h/6*(moments[n-1] - moments[n - 2]); + + mRight = + //- moments[n - 2] * x_1*x_1 / (2 * h) + //+ + moments[n - 1] * x*x / (2 * h) + + + A; + } + + // evaluate the slope for the first n-1 sample points + for (int i = 0; i < n - 1; ++ i) { + // See: J. Stoer: "Numerische Mathematik 1", 9th edition, + // Springer, 2005, p. 109 + Scalar h_i = this->h_(i + 1); + //Scalar x_i = 0; + Scalar x_i1 = h_i; + + Scalar A_i = + (y_(i+1) - y_(i))/h_i + - + h_i/6*(moments[i+1] - moments[i]); + + slopes[i] = + - moments[i] * x_i1*x_i1 / (2 * h_i) + + + //moments[i + 1] * x_i*x_i / (2 * h_i) + //+ + A_i; + + } + slopes[n - 1] = mRight; + } + + + // evaluate the spline at a given the position and given the + // segment index + Scalar eval_(Scalar x, int i) const + { + // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline + Scalar delta = h_(i + 1); + Scalar t = (x - x_(i))/delta; + + return + h00_(t) * y_(i) + + h10_(t) * slope_(i)*delta + + h01_(t) * y_(i + 1) + + h11_(t) * slope_(i + 1)*delta; + } + + // evaluate the derivative of a spline given the actual position + // and the segment index + Scalar evalDerivative_(Scalar x, int i) const + { + // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline + Scalar delta = h_(i + 1); + Scalar t = (x - x_(i))/delta; + Scalar alpha = 1 / delta; + + return + alpha * + (h00_prime_(t) * y_(i) + + h10_prime_(t) * slope_(i)*delta + + h01_prime_(t) * y_(i + 1) + + h11_prime_(t) * slope_(i + 1)*delta); + } + + // evaluate the second derivative of a spline given the actual + // position and the segment index + Scalar evalDerivative2_(Scalar x, int i) const + { + // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline + Scalar delta = h_(i + 1); + Scalar t = (x - x_(i))/delta; + Scalar alpha = 1 / delta; + + return + alpha + *(h00_prime2_(t) * y_(i) + + h10_prime2_(t) * slope_(i)*delta + + h01_prime2_(t) * y_(i + 1) + + h11_prime2_(t) * slope_(i + 1)*delta); + } + + // evaluate the third derivative of a spline given the actual + // position and the segment index + Scalar evalDerivative3_(Scalar x, int i) const + { + // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline + Scalar t = (x - x_(i))/h_(i + 1); + Scalar alpha = 1 / h_(i + 1); + + return + alpha + *(h00_prime3_(t)*y_(i) + + h10_prime3_(t)*slope_(i) + + h01_prime3_(t)*y_(i + 1) + + h11_prime3_(t)*slope_(i + 1)); + } + + // hermite basis functions + Scalar h00_(Scalar t) const + { return (2*t - 3)*t*t + 1; } + + Scalar h10_(Scalar t) const + { return ((t - 2)*t + 1)*t; } + + Scalar h01_(Scalar t) const + { return (-2*t + 3)*t*t; } + + Scalar h11_(Scalar t) const + { return (t - 1)*t*t; } + + // first derivative of the hermite basis functions + Scalar h00_prime_(Scalar t) const + { return (3*2*t - 2*3)*t; } + + Scalar h10_prime_(Scalar t) const + { return (3*t - 2*2)*t + 1; } + + Scalar h01_prime_(Scalar t) const + { return (-3*2*t + 2*3)*t; } + + Scalar h11_prime_(Scalar t) const + { return (3*t - 2)*t; } + + // second derivative of the hermite basis functions + Scalar h00_prime2_(Scalar t) const + { return 2*3*2*t - 2*3; } + + Scalar h10_prime2_(Scalar t) const + { return 2*3*t - 2*2; } + + Scalar h01_prime2_(Scalar t) const + { return -2*3*2*t + 2*3; } + + Scalar h11_prime2_(Scalar t) const + { return 2*3*t - 2; } + + // third derivative of the hermite basis functions + Scalar h00_prime3_(Scalar t) const + { return 2*3*2; } + + Scalar h10_prime3_(Scalar t) const + { return 2*3; } + + Scalar h01_prime3_(Scalar t) const + { return -2*3*2; } + + Scalar h11_prime3_(Scalar t) const + { return 2*3; } + + // returns the monotonicality of an interval of a spline segment + // + // The return value have the following meaning: + // + // 3: spline is constant within interval [x0, x1] + // 1: spline is monotonously increasing in the specified interval + // 0: spline is not monotonic (or constant) in the specified interval + // -1: spline is monotonously decreasing in the specified interval + int monotonic_(int i, Scalar x0, Scalar x1) const + { + // shift the interval so that it is consistent with the + // definitions by Stoer + x0 = x0 - x_(i); + x1 = x1 - x_(i); + + Scalar a = a_(i); + Scalar b = b_(i); + Scalar c = c_(i); + + if (std::abs(a) < 1e-20 && std::abs(b) < 1e-20 && std::abs(c) < 1e-20) + return 3; // constant in interval + + Scalar disc = 4*b*b - 12*a*c; + if (disc < 0) { + // discriminant is smaller than 0, i.e. the segment does + // not exhibit any extrema. + return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1; + } + disc = std::sqrt(disc); + Scalar xE1 = (-2*b + disc)/(6*a); + Scalar xE2 = (-2*b - disc)/(6*a); + + if (disc == 0) { + // saddle point -> no extrema + if (xE1 == x0) + // make sure that we're not picking the saddle point + // to determine whether we're monotonically increasing + // or decreasing + x0 = x1; + return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1; + }; + if ((x0 < xE1 && xE1 < x1) || + (x0 < xE2 && xE2 < x1)) + { + // there is an extremum in the range (x0, x1) + return 0; + } + // no extremum in range (x0, x1) + x0 = (x0 + x1)/2; // pick point in the middle of the interval + // to avoid extrema on the boundaries + return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1; + } + + /*! + * \brief Find all the intersections of a segment of the spline + * with a cubic polynomial within a specified interval. + */ + int intersectSegment_(Scalar *sol, + int segIdx, + Scalar a, Scalar b, Scalar c, Scalar d, + Scalar x0 = -1e100, Scalar x1 = 1e100) const + { + int n = Opm::invertCubicPolynomial(sol, + a_(segIdx) - a, + b_(segIdx) - b, + c_(segIdx) - c, + d_(segIdx) - d); + x0 = std::max(x_(segIdx), x0); + x1 = std::max(x_(segIdx+1), x1); + + // filter the intersections outside of the specified intervall + int k = 0; + for (int j = 0; j < n; ++j) { + sol[j] += x_(segIdx); // add the offset of the intervall. For details see Stoer + if (x0 <= sol[j] && sol[j] <= x1) { + sol[k] = sol[j]; + ++k; + } + } + return k; + } + + // find the segment index for a given x coordinate + int segmentIdx_(Scalar x) const + { + // bisection + int iLow = 0; + int iHigh = numSamples() - 1; + + while (iLow + 1 < iHigh) { + int i = (iLow + iHigh) / 2; + if (x_(i) > x) + iHigh = i; + else + iLow = i; + }; + return iLow; + } + + /*! + * \brief Returns x[i] - x[i - 1] + */ + Scalar h_(int i) const + { + assert(x_(i) > x_(i-1)); // the sampling points must be given + // in ascending order + return x_(i) - x_(i - 1); + } + + /*! + * \brief Returns the y coordinate of the i-th sampling point. + */ + Scalar x_(int i) const + { return xPos_[i]; } + + /*! + * \brief Returns the y coordinate of the i-th sampling point. + */ + Scalar y_(int i) const + { return yPos_[i]; } + + /*! + * \brief Returns the slope (i.e. first derivative) of the spline at + * the i-th sampling point. + */ + Scalar slope_(int i) const + { return slopeVec_[i]; } + + // returns the coefficient in front of the x^0 term. In Stoer this + // is delta. + Scalar a_(int i) const + { return evalDerivative3_(/*x=*/0, i); } + + // returns the coefficient in front of the x^2 term In Stoer this + // is gamma. + Scalar b_(int i) const + { return evalDerivative2_(/*x=*/0, i); } + + // returns the coefficient in front of the x^1 term. In Stoer this + // is beta. + Scalar c_(int i) const + { return evalDerivative_(/*x=*/0, i); } + + // returns the coefficient in front of the x^0 term. In Stoer this + // is alpha. + Scalar d_(int i) const + { return eval_(/*x=*/0, i); } + + Vector xPos_; + Vector yPos_; + Vector slopeVec_; +}; + +} + +#endif diff --git a/opm/core/utility/TridiagonalMatrix.hpp b/opm/core/utility/TridiagonalMatrix.hpp new file mode 100644 index 00000000..3b39a118 --- /dev/null +++ b/opm/core/utility/TridiagonalMatrix.hpp @@ -0,0 +1,849 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2013 by Andreas Lauser * + * * + * This program 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. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTBILITY 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \copydetails Opm::TridiagonalMatrix + */ +#ifndef OPM_TRIDIAGONAL_MATRIX_HH +#define OPM_TRIDIAGONAL_MATRIX_HH + +#include +#include +#include +#include + +#include + +namespace Opm { + +/*! + * \brief Provides a tridiagonal matrix that also supports non-zero + * entries in the upper right and lower left + * + * The entries in the lower left and upper right are supported to make + * implementing periodic systems easy. + * + * The API of this class is designed to be close to the one used by + * the DUNE matrix classes. + */ +template +class TridiagonalMatrix +{ + struct TridiagRow_ { + TridiagRow_(TridiagonalMatrix &m, size_t rowIdx) + : matrix_(m) + , rowIdx_(rowIdx) + {}; + + Scalar &operator[](size_t colIdx) + { return matrix_.at(rowIdx_, colIdx); } + + Scalar operator[](size_t colIdx) const + { return matrix_.at(rowIdx_, colIdx); } + + /*! + * \brief Prefix increment operator + */ + TridiagRow_ &operator++() + { ++ rowIdx_; return *this; } + + /*! + * \brief Prefix decrement operator + */ + TridiagRow_ &operator--() + { -- rowIdx_; return *this; } + + /*! + * \brief Comparision operator + */ + bool operator==(const TridiagRow_ &other) const + { return other.rowIdx_ == rowIdx_ && &other.matrix_ == &matrix_; } + + /*! + * \brief Comparision operator + */ + bool operator!=(const TridiagRow_ &other) const + { return !operator==(other); } + + /*! + * \brief Dereference operator + */ + TridiagRow_ &operator*() + { return *this; } + + /*! + * \brief Return the row index of the this row. + * + * 0 is the first row. + */ + size_t index() const + { return rowIdx_; } + + private: + TridiagonalMatrix &matrix_; + mutable size_t rowIdx_; + }; + +public: + typedef Scalar FieldType; + typedef TridiagRow_ RowType; + typedef size_t SizeType; + typedef TridiagRow_ iterator; + typedef TridiagRow_ const_iterator; + + explicit TridiagonalMatrix(int numRows = 0) + { + resize(numRows); + }; + + TridiagonalMatrix(int numRows, Scalar value) + { + resize(numRows); + this->operator=(value); + }; + + /*! + * \brief Copy constructor. + */ + TridiagonalMatrix(const TridiagonalMatrix &source) + { *this = source; }; + + /*! + * \brief Return the number of rows/columns of the matrix. + */ + size_t size() const + { return diag_[0].size(); } + + /*! + * \brief Return the number of rows of the matrix. + */ + size_t rows() const + { return size(); } + + /*! + * \brief Return the number of columns of the matrix. + */ + size_t cols() const + { return size(); } + + /*! + * \brief Change the number of rows of the matrix. + */ + void resize(size_t n) + { + if (n == size()) + return; + + for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) + diag_[diagIdx].resize(n); + } + + /*! + * \brief Access an entry. + */ + Scalar &at(size_t rowIdx, size_t colIdx) + { + size_t n = size(); + + // special cases + if (n > 2) { + if (rowIdx == 0 && colIdx == n - 1) + return diag_[2][0]; + if (rowIdx == n - 1 && colIdx == 0) + return diag_[0][n - 1]; + } + + int diagIdx = 1 + colIdx - rowIdx; + // make sure that the requested column is in range + assert(0 <= diagIdx && diagIdx < 3); + return diag_[diagIdx][colIdx]; + } + + /*! + * \brief Access an entry. + */ + Scalar at(size_t rowIdx, size_t colIdx) const + { + int n = size(); + + // special cases + if (rowIdx == 0 && colIdx == n - 1) + return diag_[2][0]; + if (rowIdx == n - 1 && colIdx == 0) + return diag_[0][n - 1]; + + int diagIdx = 1 + colIdx - rowIdx; + // make sure that the requested column is in range + assert(0 <= diagIdx && diagIdx < 3); + return diag_[diagIdx][colIdx]; + } + + /*! + * \brief Assignment operator from another tridiagonal matrix. + */ + TridiagonalMatrix &operator=(const TridiagonalMatrix &source) + { + for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) + diag_[diagIdx] = source.diag_[diagIdx]; + + return *this; + } + + /*! + * \brief Assignment operator from a Scalar. + */ + TridiagonalMatrix &operator=(Scalar value) + { + for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) + diag_[diagIdx].assign(size(), value); + + return *this; + } + + /*! + * \begin Iterator for the first row + */ + iterator begin() + { return TridiagRow_(*this, 0); } + + /*! + * \begin Const iterator for the first row + */ + const_iterator begin() const + { return TridiagRow_(const_cast(*this), 0); } + + /*! + * \begin Const iterator for the next-to-last row + */ + const_iterator end() const + { return TridiagRow_(const_cast(*this), size()); } + + /*! + * \brief Row access operator. + */ + TridiagRow_ operator[](size_t rowIdx) + { return TridiagRow_(*this, rowIdx); } + + /*! + * \brief Row access operator. + */ + const TridiagRow_ operator[](size_t rowIdx) const + { return TridiagRow_(*this, rowIdx); } + + /*! + * \brief Multiplication with a Scalar + */ + TridiagonalMatrix &operator*=(Scalar alpha) + { + int n = size(); + for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) { + for (int i = 0; i < n; ++i) { + diag_[diagIdx][i] *= alpha; + } + } + + return *this; + } + + /*! + * \brief Division by a Scalar + */ + TridiagonalMatrix &operator/=(Scalar alpha) + { + alpha = 1.0/alpha; + int n = size(); + for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) { + for (int i = 0; i < n; ++i) { + diag_[diagIdx][i] *= alpha; + } + } + + return *this; + } + + /*! + * \brief Subtraction operator + */ + TridiagonalMatrix &operator-=(const TridiagonalMatrix &other) + { return axpy(-1.0, other); } + + /*! + * \brief Addition operator + */ + TridiagonalMatrix &operator+=(const TridiagonalMatrix &other) + { return axpy(1.0, other); } + + + /*! + * \brief Multiply and add the matrix entries of another + * tridiagonal matrix. + * + * This means that + * \code + * A.axpy(alpha, B) + * \endcode + * is equivalent to + * \code + * A += alpha*C + * \endcode + */ + TridiagonalMatrix &axpy(Scalar alpha, const TridiagonalMatrix &other) + { + assert(size() == other.size()); + + int n = size(); + for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) + for (int i = 0; i < n; ++ i) + diag_[diagIdx][i] += alpha * other[diagIdx][i]; + + return *this; + } + + /*! + * \brief Matrix-vector product + * + * This means that + * \code + * A.mv(x, y) + * \endcode + * is equivalent to + * \code + * y = A*x + * \endcode + */ + template + void mv(const Vector &source, Vector &dest) const + { + assert(source.size() == size()); + assert(dest.size() == size()); + assert(size() > 1); + + // deal with rows 1 .. n-2 + int n = size(); + for (int i = 1; i < n - 1; ++ i) { + dest[i] = + diag_[0][i - 1]*source[i-1] + + diag_[1][i]*source[i] + + diag_[2][i + 1]*source[i + 1]; + } + + // rows 0 and n-1 + dest[0] = + diag_[1][0]*source[0] + + diag_[2][1]*source[1] + + diag_[2][0]*source[n - 1]; + + dest[n - 1] = + diag_[0][n-1]*source[0] + + diag_[0][n-2]*source[n-2] + + diag_[1][n-1]*source[n-1]; + } + + /*! + * \brief Additive matrix-vector product + * + * This means that + * \code + * A.umv(x, y) + * \endcode + * is equivalent to + * \code + * y += A*x + * \endcode + */ + template + void umv(const Vector &source, Vector &dest) const + { + assert(source.size() == size()); + assert(dest.size() == size()); + assert(size() > 1); + + // deal with rows 1 .. n-2 + int n = size(); + for (int i = 1; i < n - 1; ++ i) { + dest[i] += + diag_[0][i - 1]*source[i-1] + + diag_[1][i]*source[i] + + diag_[2][i + 1]*source[i + 1]; + } + + // rows 0 and n-1 + dest[0] += + diag_[1][0]*source[0] + + diag_[2][1]*source[1] + + diag_[2][0]*source[n - 1]; + + dest[n - 1] += + diag_[0][n-1]*source[0] + + diag_[0][n-2]*source[n-2] + + diag_[1][n-1]*source[n-1]; + } + + /*! + * \brief Subtractive matrix-vector product + * + * This means that + * \code + * A.mmv(x, y) + * \endcode + * is equivalent to + * \code + * y -= A*x + * \endcode + */ + template + void mmv(const Vector &source, Vector &dest) const + { + assert(source.size() == size()); + assert(dest.size() == size()); + assert(size() > 1); + + // deal with rows 1 .. n-2 + int n = size(); + for (int i = 1; i < n - 1; ++ i) { + dest[i] -= + diag_[0][i - 1]*source[i-1] + + diag_[1][i]*source[i] + + diag_[2][i + 1]*source[i + 1]; + } + + // rows 0 and n-1 + dest[0] -= + diag_[1][0]*source[0] + + diag_[2][1]*source[1] + + diag_[2][0]*source[n - 1]; + + dest[n - 1] -= + diag_[0][n-1]*source[0] + + diag_[0][n-2]*source[n-2] + + diag_[1][n-1]*source[n-1]; + } + + /*! + * \brief Scaled additive matrix-vector product + * + * This means that + * \code + * A.usmv(x, y) + * \endcode + * is equivalent to + * \code + * y += alpha*(A*x) + * \endcode + */ + template + void usmv(Scalar alpha, const Vector &source, Vector &dest) const + { + assert(source.size() == size()); + assert(dest.size() == size()); + assert(size() > 1); + + // deal with rows 1 .. n-2 + int n = size(); + for (int i = 1; i < n - 1; ++ i) { + dest[i] += + alpha*( + diag_[0][i - 1]*source[i-1] + + diag_[1][i]*source[i] + + diag_[2][i + 1]*source[i + 1]); + } + + // rows 0 and n-1 + dest[0] += + alpha*( + diag_[1][0]*source[0] + + diag_[2][1]*source[1] + + diag_[2][0]*source[n - 1]); + + dest[n - 1] += + alpha*( + diag_[0][n-1]*source[0] + + diag_[0][n-2]*source[n-2] + + diag_[1][n-1]*source[n-1]); + } + + /*! + * \brief Transposed matrix-vector product + * + * This means that + * \code + * A.mtv(x, y) + * \endcode + * is equivalent to + * \code + * y = A^T*x + * \endcode + */ + template + void mtv(const Vector &source, Vector &dest) const + { + assert(source.size() == size()); + assert(dest.size() == size()); + assert(size() > 1); + + // deal with rows 1 .. n-2 + int n = size(); + for (int i = 1; i < n - 1; ++ i) { + dest[i] = + diag_[2][i + 1]*source[i-1] + + diag_[1][i]*source[i] + + diag_[0][i - 1]*source[i + 1]; + } + + // rows 0 and n-1 + dest[0] = + diag_[1][0]*source[0] + + diag_[0][1]*source[1] + + diag_[0][n-1]*source[n - 1]; + + dest[n - 1] = + diag_[2][0]*source[0] + + diag_[2][n-1]*source[n-2] + + diag_[1][n-1]*source[n-1]; + } + + /*! + * \brief Transposed additive matrix-vector product + * + * This means that + * \code + * A.umtv(x, y) + * \endcode + * is equivalent to + * \code + * y += A^T*x + * \endcode + */ + template + void umtv(const Vector &source, Vector &dest) const + { + assert(source.size() == size()); + assert(dest.size() == size()); + assert(size() > 1); + + // deal with rows 1 .. n-2 + int n = size(); + for (int i = 1; i < n - 1; ++ i) { + dest[i] += + diag_[2][i + 1]*source[i-1] + + diag_[1][i]*source[i] + + diag_[0][i - 1]*source[i + 1]; + } + + // rows 0 and n-1 + dest[0] += + diag_[1][0]*source[0] + + diag_[0][1]*source[1] + + diag_[0][n-1]*source[n - 1]; + + dest[n - 1] += + diag_[2][0]*source[0] + + diag_[2][n-1]*source[n-2] + + diag_[1][n-1]*source[n-1]; + } + + /*! + * \brief Transposed subtractive matrix-vector product + * + * This means that + * \code + * A.mmtv(x, y) + * \endcode + * is equivalent to + * \code + * y -= A^T*x + * \endcode + */ + template + void mmtv (const Vector &source, Vector &dest) const + { + assert(source.size() == size()); + assert(dest.size() == size()); + assert(size() > 1); + + // deal with rows 1 .. n-2 + int n = size(); + for (int i = 1; i < n - 1; ++ i) { + dest[i] -= + diag_[2][i + 1]*source[i-1] + + diag_[1][i]*source[i] + + diag_[0][i - 1]*source[i + 1]; + } + + // rows 0 and n-1 + dest[0] -= + diag_[1][0]*source[0] + + diag_[0][1]*source[1] + + diag_[0][n-1]*source[n - 1]; + + dest[n - 1] -= + diag_[2][0]*source[0] + + diag_[2][n-1]*source[n-2] + + diag_[1][n-1]*source[n-1]; + } + + /*! + * \brief Transposed scaled additive matrix-vector product + * + * This means that + * \code + * A.umtv(alpha, x, y) + * \endcode + * is equivalent to + * \code + * y += alpha*A^T*x + * \endcode + */ + template + void usmtv(Scalar alpha, const Vector &source, Vector &dest) const + { + assert(source.size() == size()); + assert(dest.size() == size()); + assert(size() > 1); + + // deal with rows 1 .. n-2 + int n = size(); + for (int i = 1; i < n - 1; ++ i) { + dest[i] += + alpha*( + diag_[2][i + 1]*source[i-1] + + diag_[1][i]*source[i] + + diag_[0][i - 1]*source[i + 1]); + } + + // rows 0 and n-1 + dest[0] += + alpha*( + diag_[1][0]*source[0] + + diag_[0][1]*source[1] + + diag_[0][n-1]*source[n - 1]); + + dest[n - 1] += + alpha*( + diag_[2][0]*source[0] + + diag_[2][n-1]*source[n-2] + + diag_[1][n-1]*source[n-1]); + } + + /*! + * \brief Calculate the frobenius norm + * + * i.e., the square root of the sum of all squared entries. This + * corresponds to the euclidean norm for vectors. + */ + Scalar frobeniusNorm() const + { return std::sqrt(frobeniusNormSquared()); } + + /*! + * \brief Calculate the squared frobenius norm + * + * i.e., the sum of all squared entries. + */ + Scalar frobeniusNormSquared() const + { + Scalar result = 0; + + int n = size(); + for (int i = 0; i < n; ++ i) + for (int diagIdx = 0; diagIdx < 3; ++ diagIdx) + result += diag_[diagIdx][i]; + + return result; + } + + /*! + * \brief Calculate the infinity norm + * + * i.e., the maximum of the sum of the absolute values of all rows. + */ + Scalar infinityNorm() const + { + Scalar result=0; + + // deal with rows 1 .. n-2 + int n = size(); + for (int i = 1; i < n - 1; ++ i) { + result = std::max(result, + std::abs(diag_[0][i - 1]) + + std::abs(diag_[1][i]) + + std::abs(diag_[2][i + 1])); + } + + // rows 0 and n-1 + result = std::max(result, + std::abs(diag_[1][0]) + + std::abs(diag_[2][1]) + + std::abs(diag_[2][0])); + + + result = std::max(result, + std::abs(diag_[0][n-1]) + + std::abs(diag_[0][n-2]) + + std::abs(diag_[1][n-2])); + + return result; + } + + /*! + * \brief Calculate the solution for a linear system of equations + * + * i.e., calculate x, so that it solves Ax = b, where A is a + * tridiagonal matrix. + */ + template + void solve(XVector &x, const BVector &b) const + { + if (size() > 2 && diag_[2][0] != 0) + solveWithUpperRight_(x, b); + else + solveWithoutUpperRight_(x, b); + } + + /*! + * \brief Print the matrix to a given output stream. + */ + void print(std::ostream &os = std::cout) const + { + int n = size(); + + // row 0 + os << at(0, 0) << "\t" + << at(0, 1) << "\t"; + + if (n > 3) + os << "\t"; + if (n > 2) + os << at(0, n-1); + os << "\n"; + + // row 1 .. n - 2 + for (int rowIdx = 1; rowIdx < n-1; ++rowIdx) { + if (rowIdx > 1) + os << "\t"; + if (rowIdx == n - 2) + os << "\t"; + + os << at(rowIdx, rowIdx - 1) << "\t" + << at(rowIdx, rowIdx) << "\t" + << at(rowIdx, rowIdx + 1) << "\n"; + } + + // row n - 1 + if (n > 2) + os << at(n-1, 0) << "\t"; + if (n > 3) + os << "\t"; + if (n > 4) + os << "\t"; + os << at(n-1, n-2) << "\t" + << at(n-1, n-1) << "\n"; + } + +private: + template + void solveWithUpperRight_(XVector &x, const BVector &b) const + { + size_t n = size(); + + std::vector lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]), lastColumn(n); + std::vector bStar(n); + std::copy(b.begin(), b.end(), bStar.begin()); + + lastColumn[0] = upperDiag[0]; + + // forward elimination + for (size_t i = 1; i < n; ++i) { + Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1]; + + lowerDiag[i - 1] -= alpha * mainDiag[i - 1]; + mainDiag[i] -= alpha * upperDiag[i]; + + bStar[i] -= alpha * bStar[i - 1]; + }; + + // deal with the last row if the entry on the lower left is not zero + if (lowerDiag[n - 1] != 0.0 && size() > 2) { + Scalar lastRow = lowerDiag[n - 1]; + for (size_t i = 0; i < n - 1; ++i) { + Scalar alpha = lastRow/mainDiag[i]; + lastRow = - alpha*upperDiag[i + 1]; + bStar[n - 1] -= alpha * bStar[i]; + } + + mainDiag[n-1] += lastRow; + } + + // backward elimination + x[n - 1] = bStar[n - 1]/mainDiag[n-1]; + for (int i = n - 2; i >= 0; --i) { + x[i] = (bStar[i] - x[i + 1]*upperDiag[i+1] - x[n-1]*lastColumn[i])/mainDiag[i]; + } + } + + template + void solveWithoutUpperRight_(XVector &x, const BVector &b) const + { + size_t n = size(); + + std::vector lowerDiag(diag_[0]), mainDiag(diag_[1]), upperDiag(diag_[2]); + std::vector bStar(n); + std::copy(b.begin(), b.end(), bStar.begin()); + + // forward elimination + for (size_t i = 1; i < n; ++i) { + Scalar alpha = lowerDiag[i - 1]/mainDiag[i - 1]; + + lowerDiag[i - 1] -= alpha * mainDiag[i - 1]; + mainDiag[i] -= alpha * upperDiag[i]; + + bStar[i] -= alpha * bStar[i - 1]; + }; + + // deal with the last row if the entry on the lower left is not zero + if (lowerDiag[n - 1] != 0.0 && size() > 2) { + Scalar lastRow = lowerDiag[n - 1]; + for (size_t i = 0; i < n - 1; ++i) { + Scalar alpha = lastRow/mainDiag[i]; + lastRow = - alpha*upperDiag[i + 1]; + bStar[n - 1] -= alpha * bStar[i]; + } + + mainDiag[n-1] += lastRow; + } + + // backward elimination + x[n - 1] = bStar[n - 1]/mainDiag[n-1]; + for (int i = n - 2; i >= 0; --i) { + x[i] = (bStar[i] - x[i + 1]*upperDiag[i+1])/mainDiag[i]; + } + } + + mutable std::vector diag_[3]; +}; + +} // namespace Opm + +template +std::ostream &operator<<(std::ostream &os, const Opm::TridiagonalMatrix &mat) +{ + mat.print(os); + return os; +} + +#endif diff --git a/opm/core/utility/Unused.hpp b/opm/core/utility/Unused.hpp new file mode 100644 index 00000000..56ddeb7a --- /dev/null +++ b/opm/core/utility/Unused.hpp @@ -0,0 +1,28 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2013 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +#ifndef OPM_CORE_UNUSED_HH +#define OPM_CORE_UNUSED_HH + +#ifndef HAS_ATTRIBUTE_UNUSED +#define OPM_UNUSED +#else +#define OPM_UNUSED __attribute__((unused)) +#endif + +#endif diff --git a/tests/test_spline.cpp b/tests/test_spline.cpp new file mode 100644 index 00000000..f62a567a --- /dev/null +++ b/tests/test_spline.cpp @@ -0,0 +1,254 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * Copyright (C) 2010-2012 by Andreas Lauser * + * * + * This program 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. * + * * + * This program 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 this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * + * \brief This is a program to test the polynomial spline interpolation. + * + * It just prints some function to stdout. You can look at the result + * using the following commands: + * +----------- snip ----------- +./test_spline > spline.csv +gnuplot + +gnuplot> plot "spline.csv" using 1:2 w l ti "Curve", \ + "spline.csv" using 1:3 w l ti "Derivative", \ + "spline.csv" using 1:4 w p ti "Monotonical" +----------- snap ----------- + + +*/ +#include "config.h" + +#include + +#define GCC_VERSION (__GNUC__ * 10000 \ + + __GNUC_MINOR__ * 100 \ + + __GNUC_PATCHLEVEL__) + +template +void testCommon(const Spline &sp, + const double *x, + const double *y) +{ + static double eps = 1e-10; + + int n = sp.numSamples(); + for (int i = 0; i < n; ++i) { + // sure that we hit all sampling points + double y0 = (i>0)?sp.eval(x[i]-eps):y[0]; + double y1 = sp.eval(x[i]); + double y2 = (i 100*eps || std::abs(y2 - y[i]) > 100*eps) + OPM_THROW(std::runtime_error, + "Spline seems to be discontinuous at sampling point " << i << "!"); + if (std::abs(y1 - y[i]) > eps) + OPM_THROW(std::runtime_error, + "Spline does not capture sampling point " << i << "!"); + + // make sure the derivative is continuous (assuming that the + // second derivative is smaller than 1000) + double d1 = sp.evalDerivative(x[i]); + double d0 = (i>0)?sp.evalDerivative(x[i]-eps):d1; + double d2 = (i 1000*eps || std::abs(d2 - d0) > 1000*eps) + OPM_THROW(std::runtime_error, + "Spline seems to exhibit a discontinuous derivative at sampling point " << i << "!"); + } +} + +template +void testFull(const Spline &sp, + const double *x, + const double *y, + double m0, + double m1) +{ + // test the common properties of splines + testCommon(sp, x, y); + + static double eps = 1e-5; + int n = sp.numSamples(); + + // make sure the derivative at both end points is correct + double d0 = sp.evalDerivative(x[0]); + double d1 = sp.evalDerivative(x[n-1]); + if (std::abs(d0 - m0) > eps) + OPM_THROW(std::runtime_error, + "Invalid derivative at beginning of interval: is " + << d0 << " ought to be " << m0); + if (std::abs(d1 - m1) > eps) + OPM_THROW(std::runtime_error, + "Invalid derivative at end of interval: is " + << d1 << " ought to be " << m1); +} + +template +void testNatural(const Spline &sp, + const double *x, + const double *y) +{ + // test the common properties of splines + testCommon(sp, x, y); + + static double eps = 1e-5; + int n = sp.numSamples(); + + // make sure the second derivatives at both end points are 0 + double d0 = sp.evalDerivative(x[0]); + double d1 = sp.evalDerivative(x[0] + eps); + + double d2 = sp.evalDerivative(x[n-1] - eps); + double d3 = sp.evalDerivative(x[n-1]); + + if (std::abs(d1 - d0)/eps > 1000*eps) + OPM_THROW(std::runtime_error, + "Invalid second derivative at beginning of interval: is " + << (d1 - d0)/eps << " ought to be 0"); + + if (std::abs(d3 - d2)/eps > 1000*eps) + OPM_THROW(std::runtime_error, + "Invalid second derivative at end of interval: is " + << (d3 - d2)/eps << " ought to be 0"); +} + +void testAll() +{ + double x[] = { 0, 5, 7.5, 8.75, 9.375 }; + double y[] = { 10, 0, 10, 0, 10 }; + double m0 = 10; + double m1 = -10; + double points[][2] = + { + {x[0], y[0]}, + {x[1], y[1]}, + {x[2], y[2]}, + {x[3], y[3]}, + {x[4], y[4]}, + }; + + +#if GCC_VERSION >= 40500 + std::initializer_list > pointsInitList = + { + {x[0], y[0]}, + {x[1], y[1]}, + {x[2], y[2]}, + {x[3], y[3]}, + {x[4], y[4]}, + }; +#endif + + std::vector xVec; + std::vector yVec; + std::vector pointVec; + for (int i = 0; i < 5; ++i) { + xVec.push_back(x[i]); + yVec.push_back(y[i]); + pointVec.push_back(points[i]); + } + + ///////// + // test spline with two sampling points + ///////// + + // full spline + { Opm::Spline sp(x[0], x[1], y[0], y[1], m0, m1); sp.set(x[0],x[1],y[0],y[1],m0, m1); testFull(sp, x, y, m0, m1); }; + { Opm::Spline sp(2, x, y, m0, m1); sp.setXYArrays(2, x, y, m0, m1); testFull(sp, x, y, m0, m1); }; + { Opm::Spline sp(2, points, m0, m1); sp.setArrayOfPoints(2, points, m0, m1); testFull(sp, x, y, m0, m1); }; + + ///////// + // test variable length splines + ///////// + + // full spline + { Opm::Spline sp(5, x, y, m0, m1); sp.setXYArrays(5,x,y,m0, m1); testFull(sp, x, y, m0, m1); }; + { Opm::Spline sp(xVec, yVec, m0, m1); sp.setXYContainers(xVec,yVec,m0, m1); testFull(sp, x, y, m0, m1); }; + { Opm::Spline sp; sp.setArrayOfPoints(5,points,m0, m1); testFull(sp, x, y, m0, m1); }; + { Opm::Spline sp; sp.setContainerOfPoints(pointVec,m0, m1); testFull(sp, x, y, m0, m1); }; +#if GCC_VERSION >= 40500 + { Opm::Spline sp; sp.setContainerOfTuples(pointsInitList,m0, m1); testFull(sp, x, y, m0, m1); }; +#endif + + // natural spline + { Opm::Spline sp(5, x, y); sp.setXYArrays(5,x,y); testNatural(sp, x, y); }; + { Opm::Spline sp(xVec, yVec); sp.setXYContainers(xVec,yVec); testNatural(sp, x, y); }; + { Opm::Spline sp; sp.setArrayOfPoints(5,points); testNatural(sp, x, y); }; + { Opm::Spline sp; sp.setContainerOfPoints(pointVec); testNatural(sp, x, y); }; +#if GCC_VERSION >= 40500 + { Opm::Spline sp; sp.setContainerOfTuples(pointsInitList); testNatural(sp, x, y); }; +#endif +} + +void plot() +{ + const int numSamples = 5; + const int n = numSamples - 1; + typedef std::array FV; + + double x_[] = { 0, 5, 7.5, 8.75, 10 }; + double y_[] = { 10, 0, 10, 0, 10 }; + double m1 = 10; + double m2 = -10; + FV &xs = *reinterpret_cast(x_); + FV &ys = *reinterpret_cast(y_); + + Opm::Spline spFull(xs, ys, m1, m2); + Opm::Spline spNatural(xs, ys); + Opm::Spline spPeriodic(xs, ys, /*type=*/Opm::Spline::Periodic); + Opm::Spline spMonotonic(xs, ys, /*type=*/Opm::Spline::Monotonic); + + spFull.printCSV(x_[0] - 1.00001, + x_[n] + 1.00001, + 1000); + std::cout << "\n"; + spNatural.printCSV(x_[0] - 1.00001, + x_[n] + 1.00001, + 1000); + std::cout << "\n"; + + spPeriodic.printCSV(x_[0] - 1.00001, + x_[n] + 1.00001, + 1000); + std::cout << "\n"; + + spMonotonic.printCSV(x_[0] - 1.00001, + x_[n] + 1.00001, + 1000); + std::cout << "\n"; + + std::cerr << "Spline is monotonic: " << spFull.monotonic(x_[0], x_[n]) << "\n"; +} + +int main(int argc, char** argv) +{ + try { + testAll(); + + plot(); + } + catch (const std::exception &e) { + std::cout << "Caught OPM exception: " << e.what() << "\n"; + } + + return 0; +}