1763 lines
58 KiB
C++
1763 lines
58 KiB
C++
// -*- 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 <opm/core/utility/TridiagonalMatrix.hpp>
|
|
#include <opm/core/utility/PolynomialUtils.hpp>
|
|
#include <opm/core/utility/ErrorMacros.hpp>
|
|
|
|
#include <ostream>
|
|
#include <vector>
|
|
#include <tuple>
|
|
|
|
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 Scalar>
|
|
class Spline
|
|
{
|
|
typedef Opm::TridiagonalMatrix<Scalar> Matrix;
|
|
typedef std::vector<Scalar> 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 <class ScalarArrayX, class ScalarArrayY>
|
|
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 <class PointArray>
|
|
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 <class ScalarContainer>
|
|
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 <class PointContainer>
|
|
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 <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 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 <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 = 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 <class ScalarContainerX, class ScalarContainerY>
|
|
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 <class PointArray>
|
|
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 <class XYContainer>
|
|
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 <class XYContainer>
|
|
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(x, x_p1, /*extrapolate=*/true);
|
|
}
|
|
|
|
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())
|
|
return evalDerivative_(xMin(), /*segmentIdx=*/0);
|
|
else if (x > xMax())
|
|
return evalDerivative_(xMax(), /*segmentIdx=*/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 Evaluate the spline's third 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 evalThirdDerivative(Scalar x, bool extrapolate=false) const
|
|
{
|
|
assert(extrapolate || applies(x));
|
|
if (extrapolate)
|
|
return 0.0;
|
|
|
|
return evalDerivative3_(x, segmentIdx_(x));
|
|
}
|
|
|
|
/*!
|
|
* \brief Find the intersections of the spline with a cubic
|
|
* polynomial in the whole interval, 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 intersectInterval(xMin(), xMax(), a, b, c, d);
|
|
}
|
|
|
|
/*!
|
|
* \brief Find the intersections of the spline with a cubic
|
|
* polynomial in a sub-interval 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], sol = 0;
|
|
int nSol = 0;
|
|
int iFirst = segmentIdx_(x0);
|
|
int iLast = segmentIdx_(x1);
|
|
for (int i = iFirst; i <= iLast; ++i)
|
|
{
|
|
int nCur = intersectSegment_(tmpSol, i, a, b, c, d, x0, x1);
|
|
if (nCur == 1)
|
|
sol = tmpSol[0];
|
|
|
|
nSol += nCur;
|
|
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 sol;
|
|
}
|
|
|
|
/*!
|
|
* \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, bool extrapolate=false) const
|
|
{
|
|
assert(x0 != x1);
|
|
|
|
// make sure that x0 is smaller than x1
|
|
if (x0 > x1)
|
|
std::swap(x0, x1);
|
|
|
|
assert(x0 < x1);
|
|
|
|
int r = 3;
|
|
if (x0 < xMin()) {
|
|
static_cast<void>(extrapolate);
|
|
assert(extrapolate);
|
|
Scalar m = evalDerivative_(xMin(), /*segmentIdx=*/0);
|
|
if (std::abs(m) < 1e-20)
|
|
r = (m < 0)?-1:1;
|
|
x0 = xMin();
|
|
};
|
|
|
|
int i = segmentIdx_(x0);
|
|
if (x_(i + 1) >= x1) {
|
|
// interval is fully contained within a single spline
|
|
// segment
|
|
monotonic_(i, x0, x1, r);
|
|
return r;
|
|
}
|
|
|
|
// the first segment overlaps with the specified interval
|
|
// partially
|
|
monotonic_(i, x0, x_(i+1), r);
|
|
++ i;
|
|
|
|
// make sure that the segments which are completly in the
|
|
// interval [x0, x1] all exhibit the same monotonicity.
|
|
int iEnd = segmentIdx_(x1);
|
|
for (; i < iEnd - 1; ++i) {
|
|
monotonic_(i, x_(i), x_(i + 1), r);
|
|
if (!r)
|
|
return 0;
|
|
}
|
|
|
|
// if the user asked for a part of the spline which is
|
|
// extrapolated, we need to check the slope at the spline's
|
|
// endpoint
|
|
if (x1 > xMax()) {
|
|
assert(extrapolate);
|
|
|
|
Scalar m = evalDerivative_(xMax(), /*segmentIdx=*/numSamples() - 2);
|
|
if (m < 0)
|
|
return (r < 0 || r==3)?-1:0;
|
|
else if (m > 0)
|
|
return (r > 0 || r==3)?1:0;
|
|
|
|
return r;
|
|
}
|
|
|
|
// check for the last segment
|
|
monotonic_(iEnd, x_(iEnd), x1, r);
|
|
|
|
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<Scalar> &x)
|
|
: x_(x)
|
|
{};
|
|
|
|
bool operator ()(int idxA, int idxB) const
|
|
{ return x_.at(idxA) < x_.at(idxB); }
|
|
|
|
const std::vector<Scalar> &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<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 (size_t 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_()
|
|
{
|
|
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 <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 = 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 size_t n = numSamples() - 1;
|
|
|
|
assert(M.rows() == n);
|
|
|
|
// second to next to last rows
|
|
for (size_t 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 = numSamples();
|
|
|
|
// calculate the slopes of the secant lines
|
|
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));
|
|
}
|
|
|
|
// calculate the "raw" slopes at the sample points
|
|
for (int k = 1; k < n - 1; ++k)
|
|
slopes[k] = (delta[k - 1] + delta[k])/2;
|
|
slopes[0] = delta[0];
|
|
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-50) {
|
|
// make the spline flat if the inputs are equal
|
|
slopes[k] = 0;
|
|
slopes[k + 1] = 0;
|
|
++ k;
|
|
continue;
|
|
}
|
|
else {
|
|
Scalar alpha = slopes[k] / delta[k];
|
|
Scalar beta = slopes[k + 1] / delta[k];
|
|
|
|
if (alpha < 0 || (k > 0 && slopes[k] / delta[k - 1] < 0)) {
|
|
slopes[k] = 0;
|
|
}
|
|
// limit (alpha, beta) to a circle of radius 3
|
|
else 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 = 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*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 delta = h_(i + 1);
|
|
Scalar t = (x - x_(i))/delta;
|
|
Scalar alpha = 1 / delta;
|
|
|
|
return
|
|
alpha*alpha*alpha
|
|
*(h00_prime3_(t)*y_(i)
|
|
+ h10_prime3_(t)*slope_(i)*delta
|
|
+ h01_prime3_(t)*y_(i + 1)
|
|
+ h11_prime3_(t)*slope_(i + 1)*delta);
|
|
}
|
|
|
|
// 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, int &r) const
|
|
{
|
|
// coefficients of derivative in monomial basis
|
|
Scalar a = 3*a_(i);
|
|
Scalar b = 2*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, r stays unchanged!
|
|
|
|
Scalar disc = b*b - 4*a*c;
|
|
if (disc < 0) {
|
|
// discriminant of derivative is smaller than 0, i.e. the
|
|
// segment's derivative does not exhibit any extrema.
|
|
if (x0*(x0*a + b) + c > 0) {
|
|
r = (r==3 || r == 1)?1:0;
|
|
return 1;
|
|
}
|
|
else {
|
|
r = (r==3 || r == -1)?-1:0;
|
|
return -1;
|
|
}
|
|
}
|
|
disc = std::sqrt(disc);
|
|
Scalar xE1 = (-b + disc)/(2*a);
|
|
Scalar xE2 = (-b - disc)/(2*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;
|
|
if (x0*(x0*a + b) + c > 0) {
|
|
r = (r==3 || r == 1)?1:0;
|
|
return 1;
|
|
}
|
|
else {
|
|
r = (r==3 || r == -1)?-1:0;
|
|
return -1;
|
|
}
|
|
};
|
|
if ((x0 < xE1 && xE1 < x1) ||
|
|
(x0 < xE2 && xE2 < x1))
|
|
{
|
|
// there is an extremum in the range (x0, x1)
|
|
r = 0;
|
|
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
|
|
if (x0*(x0*a + b) + c > 0) {
|
|
r = (r==3 || r == 1)?1:0;
|
|
return 1;
|
|
}
|
|
else {
|
|
r = (r==3 || r == -1)?-1:0;
|
|
return -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::min(x_(segIdx+1), x1);
|
|
|
|
// filter the intersections outside of the specified interval
|
|
int k = 0;
|
|
for (int j = 0; j < n; ++j) {
|
|
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^3 term. In Stoer this
|
|
// is delta.
|
|
Scalar a_(int i) const
|
|
{ return evalDerivative3_(/*x=*/0, i)/6.0; }
|
|
|
|
// 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)/2.0; }
|
|
|
|
// 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
|