OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

introduce classes for monotonic, full, periodic and natural cubic splines

Browse Source
This commit is contained in:
Andreas Lauser 2013-08-22 18:44:06 +02:00
parent 2302fc6ce5
commit 8ef3b6e6c7
6 changed files with 2404 additions and 310 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
CMakeLists_files.cmake
View File

@ -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
@ -208,6 +209,13 @@ list (APPEND PROGRAM_SOURCE_FILES
# originally generated with the command:
# find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort
list (APPEND PUBLIC_HEADER_FILES
opm/core/FixedLengthSpline_.hpp
opm/core/PolynomialUtils.hpp
opm/core/Spline.hpp
opm/core/SplineCommon_.hpp
opm/core/TridiagonalMatrix.hpp
opm/core/Unused.hpp
opm/core/VariableLengthSpline_.hpp
opm/core/doxygen_main.hpp
opm/core/grid.h
opm/core/grid/CellQuadrature.hpp

310
opm/core/PolynomialUtils.hpp
View File

@ -1,310 +0,0 @@
// -*- 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Define some often used mathematical functions
*/
#ifndef OPM_MATH_HH
#define OPM_MATH_HH
#include <cmath>
#include <algorithm>
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 <class Scalar, class SolContainer>
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 <class Scalar, class SolContainer>
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 <class Scalar, class SolContainer>
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 <class Scalar, class SolContainer>
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

577
opm/core/utility/Spline.hpp Normal file
View File

@ -0,0 +1,577 @@
// -*- 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \copydoc Opm::Spline
*/
#ifndef OPM_SPLINE_HH
#define OPM_SPLINE_HH
#include "VariableLengthSpline_.hpp"
#include "SplineCommon_.hpp"
namespace Opm
{
/*!
* \ingroup Spline
* \brief A 3rd order polynomial spline.
*
* This class implements a spline \f$s(x)\f$ for which, given \f$n\f$ sampling
* points \f$x_1, \dots, x_n\f$, the following conditions hold
* \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}
* are supported. The final supported type of splines is periodic splines, i.e.,
*\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}
*/
template<class Scalar, int numSamples = 2>
class Spline : public VariableLengthSpline_<Scalar>
{
public:
/*!
* \brief Default constructor for a spline.
*
* To specfiy the acutal curve, use one of the set() methods.
*/
Spline()
{ }
/*!
* \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
* \param y An array containing the \f$y\f$ values of the spline's sampling points
* \param sortInputs Indicates whether the sample points should be sorted (this is not necessary if they are already sorted in ascending or descending order)
* \param splineType Indicate whether a periodic, natural or monotonic spline is desired
*/
template <class ScalarArray>
Spline(const ScalarArray &x,
const ScalarArray &y,
SplineType splineType = NaturalSpline,
bool sortInputs = false)
{ this->setXYArrays(numSamples, 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
* \param periodic Indicates whether a natural or a periodic spline should be created
* \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 <class PointArray>
Spline(const PointArray &points,
SplineType splineType = NaturalSpline,
bool sortInputs = false)
{ this->setArrayOfPoints(numSamples, points, splineType, sortInputs); }
/*!
* \brief Convenience constructor for a full spline
*
* \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 <class ScalarArray>
Spline(const ScalarArray &x,
const ScalarArray &y,
Scalar m0,
Scalar m1,
bool sortInputs = false)
{ this->setXYArrays(numSamples, 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
* \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 <class PointArray>
Spline(const PointArray &points,
Scalar m0,
Scalar m1,
bool sortInputs = false)
{ this->setArrayOfPoints(numSamples, points, m0, m1, sortInputs); }
};
/*!
* \brief Specialization of a spline with the number of sampling
* points only known at run time.
*
* This class implements a spline \f$s(x)\f$ for which, given \f$n\f$ sampling
* points \f$x_1, \dots, x_n\f$, the following conditions hold
* \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}
* are supported. The final supported type of splines is periodic splines, i.e.,
*\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}
*/
template<class Scalar>
class Spline<Scalar, /*numSamples=*/-1> : public VariableLengthSpline_<Scalar>
{
public:
/*!
* \brief Default constructor for a spline.
*
* To specfiy the acutal curve, use one of the set() methods.
*/
Spline()
{ }
/*!
* \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 <class ScalarArrayX, class ScalarArrayY>
Spline(int nSamples,
const ScalarArrayX &x,
const ScalarArrayY &y,
SplineType splineType = NaturalSpline,
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 <class PointArray>
Spline(int nSamples,
const PointArray &points,
SplineType splineType = NaturalSpline,
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 <class ScalarContainer>
Spline(const ScalarContainer &x,
const ScalarContainer &y,
SplineType splineType = NaturalSpline,
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 <class PointContainer>
Spline(const PointContainer &points,
SplineType splineType = NaturalSpline,
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 <class ScalarArray>
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 <class PointArray>
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 <class ScalarContainerX, class ScalarContainerY>
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 <class PointContainer>
Spline(const PointContainer &points,
Scalar m0,
Scalar m1,
bool sortInputs = false)
{ this->setContainerOfPoints(points, m0, m1, sortInputs); }
};
/*!
* \brief Do not allow splines with zero sampling points
*/
template<class Scalar>
class Spline<Scalar, /*numSamples=*/0>
// Splines with zero sampling points do not make sense!
{ private: Spline() { } };
/*!
* \brief Do not allow splines with one sampling point
*/
template<class Scalar>
class Spline<Scalar, /*numSamples=*/1>
// Splines with one sampling point do not make sense!
{ private: Spline() { } };
/*!
* \brief Spline for two sampling points.
*
* For this type of spline there is no natural and periodic spline.
*/
template<class Scalar>
class Spline<Scalar, 2> : public SplineCommon_<Scalar, Spline<Scalar, 2> >
{
friend class SplineCommon_<Scalar, Spline<Scalar, 2> >;
typedef std::array<Scalar, 2> Vector;
typedef Opm::TridiagonalMatrix<Scalar> Matrix;
public:
Spline()
{}
/*!
* \brief Convenience constructor for a full spline
*
* \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$
*/
template <class ScalarArrayX, class ScalarArrayY>
Spline(const ScalarArrayX &x,
const ScalarArrayY &y,
Scalar m0, Scalar m1)
{ setXYArrays(2, x, y, m0, m1); }
/*!
* \brief Convenience constructor for a full spline
*
* \param points An array of \f$(x,y)\f$ tuples 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$
*/
template <class PointArray>
Spline(const PointArray &points,
Scalar m0,
Scalar m1)
{ this->setArrayOfPoints(2, points, m0, m1); }
/*!
* \brief Convenience constructor for a full spline
*
* \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_n\f$
*/
Spline(Scalar x0, Scalar x1,
Scalar y0, Scalar y1,
Scalar m0, Scalar m1)
{
set(x0, x1,
y0, y1,
m0, m1);
}
/*!
* \brief Returns the number of sampling points.
*/
int numSamples() const
{ return 2; }
/*!
* \brief Set the sampling points and the boundary slopes of the
* spline.
*
* \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_[0] = m0;
slopeVec_[1] = m1;
Matrix M;
Vector d;
Vector moments;
M.resize(numSamples());
assignXY_(x0, x1, y0, y1);
this->makeFullSystem_(M, d, m0, m1);
// solve for the moments
M.solve(moments, d);
this->setSlopesFromMoments_(slopeVec_, moments);
}
/*!
* \brief Set the sampling points and the boundary slopes of the
* spline.
*
* \param nSamples The number of sampling points (must be >= 2)
* \param x An array containing the \f$x\f$ values of the sampling points
* \param y An array containing the \f$y\f$ values of the sampling points
* \param m0 The slope of the spline at \f$x_0\f$
* \param m1 The slope of the spline at \f$x_1\f$
*/
template <class ScalarContainer>
void setXYArrays(int nSamples,
const ScalarContainer &x,
const ScalarContainer &y,
Scalar m0, Scalar m1)
{
assert(nSamples == 2);
set(x[0], x[1], y[0], y[1], m0, m1);
}
/*!
* \brief Set the sampling points and the boundary slopes of the
* spline.
*
* \param x An array containing the \f$x\f$ values of the sampling points
* \param y An array containing the \f$y\f$ values of the sampling points
* \param m0 The slope of the spline at \f$x_0\f$
* \param m1 The slope of the spline at \f$x_1\f$
*/
template <class ScalarContainerX, class ScalarContainerY>
void setXYContainers(const ScalarContainerX &x,
const ScalarContainerY &y,
Scalar m0, Scalar m1)
{
assert(x.size() == y.size());
assert(x.size() == 2);
Matrix M(numSamples());
Vector d;
typename ScalarContainerX::const_iterator xIt0 = x.begin();
typename ScalarContainerX::const_iterator xIt1 = xIt0;
++xIt1;
typename ScalarContainerY::const_iterator yIt0 = y.begin();
typename ScalarContainerY::const_iterator yIt1 = yIt0;
++yIt1;
set(*xIt0, *xIt1, *yIt0, *yIt1);
}
/*!
* \brief Set the sampling points and the boundary slopes of the
* 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 m0 The slope of the spline at \f$x_0\f$
* \param m1 The slope of the spline at \f$x_1\f$
*/
template <class PointArray>
void setArrayOfPoints(int nSamples,
const PointArray &points,
Scalar m0,
Scalar m1)
{
assert(nSamples == 2);
set(points[0][0],
points[1][0],
points[0][1],
points[1][1],
m0, m1);
}
/*!
* \brief Set the sampling points and the boundary slopes from an
* STL-like container of points.
*
* \param points An array of \f$(x,y)\f$ tuples 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_1\f$
*/
template <class PointContainer>
void setContainerOfPoints(const PointContainer &points,
Scalar m0,
Scalar m1)
{
assert(points.size() == 2);
Matrix M;
Vector d;
typename PointContainer::const_iterator it0 = points.begin();
typename PointContainer::const_iterator it1 = it0;
++it1;
set((*it0)[0],
(*it0)[1],
(*it1)[0],
(*it1)[1],
m0, m1);
}
/*!
* \brief Set the sampling points and the boundary slopes from an
* STL-like container of tuples.
*
* \param tuples An array of \f$(x,y)\f$ tuples 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_1\f$
*/
template <class TupleContainer>
void setContainerOfTuples(const TupleContainer &tuples,
Scalar m0,
Scalar m1)
{
assert(tuples.size() == 2);
typename TupleContainer::const_iterator it0 = tuples.begin();
typename TupleContainer::const_iterator it1 = it0;
++it1;
set(std::get<0>(*it0),
std::get<1>(*it0),
std::get<0>(*it1),
std::get<1>(*it1),
m0, m1);
}
protected:
void assignXY_(Scalar x0, Scalar x1,
Scalar y0, Scalar y1)
{
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;
}
}
/*!
* \brief Returns the x 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 moment (i.e. second derivative) of the
* spline at the i-th sampling point.
*/
Scalar slope_(int i) const
{ return slopeVec_[i]; }
Vector xPos_;
Vector yPos_;
Vector slopeVec_;
};
}
#endif

923
opm/core/utility/SplineCommon_.hpp Normal file
View File

@ -0,0 +1,923 @@
// -*- 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \copydoc Opm::SplineCommon_
*/
#ifndef OPM_SPLINE_COMMON__HH
#define OPM_SPLINE_COMMON__HH
#include <opm/core/utility/PolynomialUtils.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/Exceptions.hpp>
#include <tuple>
#include <algorithm>
#include <iostream>
#include <cassert>
namespace Opm
{
enum SplineType {
FullSpline,
NaturalSpline,
PeriodicSpline,
MonotonicSpline
};
/*!
* \brief The common code for all 3rd order polynomial splines.
*/
template<class ScalarT, class ImplementationT>
class SplineCommon_
{
typedef ScalarT Scalar;
typedef ImplementationT Implementation;
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
const Implementation &asImp_() const
{ return *static_cast<const Implementation*>(this); }
public:
/*!
* \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<Scalar>(x_(0), x), std::min<Scalar>(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"); //<<a<<"x^3 + "<<b<"x^2 + "<<c<"x + "<<d);
}
}
if (nSol != 1)
OPM_THROW(std::runtime_error,
"Spline has no intersection"); //<<a<"x^3 + " <<b<"x^2 + "<<c<"x + "<<d<<"!");
return tmpSol[0];
}
/*!
* \brief Returns 1 if the spline is monotonically increasing, -1
* if the spline is mononously decreasing and 0 if the
* spline is not monotonous in the interval (x0, x1).
*
* In the corner case that the spline is constant within the given
* interval, this method returns 3.
*/
int monotonic(Scalar x0, Scalar x1) const
{
assert(applies(x0));
assert(applies(x1));
assert(x0 != x1);
// make sure that x0 is smaller than x1
if (x0 > 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:
// this is an internal class, so everything is protected!
SplineCommon_()
{ }
/*!
* \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 <class DestVector, class SourceVector>
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 <class DestVector, class ListIterator>
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 <class DestVector, class ListIterator>
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 <class Vector, class Matrix>
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 <class Vector, class Matrix>
void makeNaturalSystem_(Matrix &M, Vector &d)
{
M = 0.0;
// See: J. Stoer: "Numerische Mathematik 1", 9th edition,
// Springer, 2005, p. 111
const int n = asImp_().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 <class Matrix, class Vector>
void makePeriodicSystem_(Matrix &M, Vector &d)
{
M = 0.0;
// See: J. Stoer: "Numerische Mathematik 1", 9th edition,
// Springer, 2005, p. 111
const int n = asImp_().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 <class Vector>
void makeMonotonicSpline_(Vector &slopes)
{
auto n = asImp_().numSamples();
std::vector<Scalar> 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 <class MomentsVector, class SlopeVector>
void setSlopesFromMoments_(SlopeVector &slopes, const MomentsVector &moments)
{
int n = asImp_().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 asImp_().x_(i); }
/*!
* \brief Returns the y coordinate of the i-th sampling point.
*/
Scalar y_(int i) const
{ return asImp_().y_(i); }
/*!
* \brief Returns the slope (i.e. first derivative) of the spline at
* the i-th sampling point.
*/
Scalar slope_(int i) const
{ return asImp_().slope_(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); }
/*!
* \brief Returns the number of sampling points.
*/
int numSamples_() const
{ return asImp_().numSamples(); }
};
}
#endif

623
opm/core/utility/VariableLengthSpline_.hpp Normal file
View File

@ -0,0 +1,623 @@
// -*- 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \copydoc Opm::VariableLengthSpline_
*/
#ifndef OPM_VARIABLE_LENGTH_SPLINE_HH
#define OPM_VARIABLE_LENGTH_SPLINE_HH
#include "SplineCommon_.hpp"
#include "TridiagonalMatrix.hpp"
#include <array>
#include <vector>
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
namespace Opm {
/*!
* \brief The common code for all 3rd order polynomial splines with
* where the number of sampling points only known at run-time.
*/
template<class Scalar>
class VariableLengthSpline_
: public SplineCommon_<Scalar,
VariableLengthSpline_<Scalar> >
{
struct ComparatorX_
{
ComparatorX_(const std::vector<Scalar> &x)
: x_(x)
{};
bool operator ()(int idxA, int idxB) const
{ return x_.at(idxA) < x_.at(idxB); }
const std::vector<Scalar> &x_;
};
friend class SplineCommon_<Scalar, VariableLengthSpline_<Scalar> >;
typedef std::vector<Scalar> Vector;
typedef Opm::TridiagonalMatrix<Scalar> Matrix;
public:
/*!
* \brief Returns the number of sampling points.
*/
int numSamples() const
{ return xPos_.size(); }
///////////////////////////////////////
///////////////////////////////////////
///////////////////////////////////////
// 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 <class ScalarArrayX, class ScalarArrayY>
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 <class ScalarContainerX, class ScalarContainerY>
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 <class PointArray>
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 <class XYContainer>
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 <class XYContainer>
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 <class ScalarArrayX, class ScalarArrayY>
void setXYArrays(int nSamples,
const ScalarArrayX &x,
const ScalarArrayY &y,
SplineType splineType = NaturalSpline,
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 == PeriodicSpline)
makePeriodicSpline_();
else if (splineType == NaturalSpline)
makeNaturalSpline_();
else if (splineType == MonotonicSpline)
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 <class ScalarContainerX, class ScalarContainerY>
void setXYContainers(const ScalarContainerX &x,
const ScalarContainerY &y,
SplineType splineType = NaturalSpline,
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 == PeriodicSpline)
makePeriodicSpline_();
else if (splineType == NaturalSpline)
makeNaturalSpline_();
else if (splineType == MonotonicSpline)
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 <class PointArray>
void setArrayOfPoints(int nSamples,
const PointArray &points,
SplineType splineType = NaturalSpline,
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 == PeriodicSpline)
makePeriodicSpline_();
else if (splineType == NaturalSpline)
makeNaturalSpline_();
else if (splineType == MonotonicSpline)
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 <class XYContainer>
void setContainerOfPoints(const XYContainer &points,
SplineType splineType = NaturalSpline,
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 == PeriodicSpline)
makePeriodicSpline_();
else if (splineType == NaturalSpline)
makeNaturalSpline_();
else if (splineType == MonotonicSpline)
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 <class XYContainer>
void setContainerOfTuples(const XYContainer &points,
SplineType splineType = NaturalSpline,
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 == PeriodicSpline)
makePeriodicSpline_();
else if (splineType == NaturalSpline)
makeNaturalSpline_();
else if (splineType == MonotonicSpline)
this->makeMonotonicSpline_(slopeVec_);
else
OPM_THROW(std::runtime_error, "Spline type " << splineType << " not supported at this place");
}
protected:
/*!
* \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<int> 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<Scalar> 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<Scalar> d(numSamples());
std::vector<Scalar> 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_()
{
Dune::DynamicMatrix<Scalar> M(numSamples(), numSamples());
Dune::DynamicVector<Scalar> d(numSamples());
Dune::DynamicVector<Scalar> 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 Returns the x 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 slope (i.e. first derivative) of the spline at
* the i-th sampling point.
*/
Scalar slope_(int i) const
{ return slopeVec_[i]; }
Vector xPos_;
Vector yPos_;
Vector slopeVec_;
};
}
#endif

273
tests/test_spline.cpp Normal file
View File

@ -0,0 +1,273 @@
// -*- 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 <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \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 <opm/core/utility/Spline.hpp>
#define GCC_VERSION (__GNUC__ * 10000 \
+ __GNUC_MINOR__ * 100 \
+ __GNUC_PATCHLEVEL__)
template <class Spline>
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<n-1)?sp.eval(x[i]+eps):y[n-1];
if (std::abs(y0 - y[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<n-1)?sp.evalDerivative(x[i]+eps):d1;
if (std::abs(d1 - d0) > 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 <class Spline>
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 <class Spline>
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<const std::pair<double, double> > pointsInitList =
{
{x[0], y[0]},
{x[1], y[1]},
{x[2], y[2]},
{x[3], y[3]},
{x[4], y[4]},
};
#endif
std::vector<double> xVec;
std::vector<double> yVec;
std::vector<double*> pointVec;
for (int i = 0; i < 5; ++i) {
xVec.push_back(x[i]);
yVec.push_back(y[i]);
pointVec.push_back(points[i]);
}
/////////
// test fixed length spline, n = 2
/////////
// full spline
{ Opm::Spline<double, 2> 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<double, 2> sp(x, y, m0, m1); sp.setXYArrays(2, x, y, m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double, 2> sp(points, m0, m1); sp.setArrayOfPoints(2, points, m0, m1); testFull(sp, x, y, m0, m1); };
/////////
// test fixed length spline, n > 2
/////////
// full spline
{ Opm::Spline<double, 5> sp(x, y, m0, m1); sp.setXYArrays(5, x, y, m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double, 5> sp; sp.setArrayOfPoints(5, points, m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double, 5> sp; sp.setContainerOfPoints(pointVec,m0, m1); testFull(sp, x, y, m0, m1); };
#if GCC_VERSION >= 40500
{ Opm::Spline<double, 5> sp; sp.setContainerOfTuples(pointsInitList, m0, m1); testFull(sp, x, y, m0, m1); };
#endif
// natural spline
{ Opm::Spline<double, 5> sp(x, y); sp.setXYArrays(5, x, y); testNatural(sp, x, y); };
{ Opm::Spline<double, 5> sp; sp.setContainerOfPoints(pointVec); testNatural(sp, x, y); };
#if GCC_VERSION >= 40500
{ Opm::Spline<double, 5> sp; sp.setContainerOfTuples(pointsInitList); testNatural(sp, x, y); };
#endif
/////////
// test variable length splines
/////////
// full spline
{ Opm::Spline<double, -1> sp(5, x, y, m0, m1); sp.setXYArrays(5,x,y,m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double, -1> sp(xVec, yVec, m0, m1); sp.setXYContainers(xVec,yVec,m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double, -1> sp; sp.setArrayOfPoints(5,points,m0, m1); testFull(sp, x, y, m0, m1); };
{ Opm::Spline<double, -1> sp; sp.setContainerOfPoints(pointVec,m0, m1); testFull(sp, x, y, m0, m1); };
#if GCC_VERSION >= 40500
{ Opm::Spline<double, -1> sp; sp.setContainerOfTuples(pointsInitList,m0, m1); testFull(sp, x, y, m0, m1); };
#endif
// natural spline
{ Opm::Spline<double, -1> sp(5, x, y); sp.setXYArrays(5,x,y); testNatural(sp, x, y); };
{ Opm::Spline<double, -1> sp(xVec, yVec); sp.setXYContainers(xVec,yVec); testNatural(sp, x, y); };
{ Opm::Spline<double, -1> sp; sp.setArrayOfPoints(5,points); testNatural(sp, x, y); };
{ Opm::Spline<double, -1> sp; sp.setContainerOfPoints(pointVec); testNatural(sp, x, y); };
#if GCC_VERSION >= 40500
{ Opm::Spline<double, -1> sp; sp.setContainerOfTuples(pointsInitList); testNatural(sp, x, y); };
#endif
}
void plot()
{
const int numSamples = 5;
const int n = numSamples - 1;
typedef std::array<double, numSamples> 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<FV*>(x_);
FV &ys = *reinterpret_cast<FV*>(y_);
Opm::Spline<double, numSamples> spFull(xs, ys, m1, m2);
Opm::Spline<double, numSamples> spNatural(xs, ys);
Opm::Spline<double, numSamples> spPeriodic(xs, ys, /*type=*/Opm::PeriodicSpline);
Opm::Spline<double, numSamples> spMonotonic(xs, ys, /*type=*/Opm::MonotonicSpline);
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;
}