| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup Core | ||
| 10 | * \brief A monotone cubic spline | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_COMMON_MONOTONE_CUBIC_SPLINE_HH | ||
| 13 | #define DUMUX_COMMON_MONOTONE_CUBIC_SPLINE_HH | ||
| 14 | |||
| 15 | #include <cmath> | ||
| 16 | #include <cassert> | ||
| 17 | #include <iterator> | ||
| 18 | #include <vector> | ||
| 19 | #include <algorithm> | ||
| 20 | |||
| 21 | #include <dune/common/float_cmp.hh> | ||
| 22 | |||
| 23 | // Hermite basis functions | ||
| 24 | #include <dumux/common/cubicsplinehermitebasis.hh> | ||
| 25 | |||
| 26 | // for inversion | ||
| 27 | #include <dumux/nonlinear/findscalarroot.hh> | ||
| 28 | |||
| 29 | namespace Dumux { | ||
| 30 | |||
| 31 | /*! | ||
| 32 | * \ingroup Core | ||
| 33 | * \brief A monotone cubic spline | ||
| 34 | * \note Construction after Fritsch & Butland (1984) (see https://doi.org/10.1137/0905021) | ||
| 35 | * \note The resulting interpolation is globally monotone but only C^1 | ||
| 36 | */ | ||
| 37 | template<class Scalar = double> | ||
| 38 | class MonotoneCubicSpline | ||
| 39 | { | ||
| 40 | using Basis = CubicSplineHermiteBasis<Scalar>; | ||
| 41 | public: | ||
| 42 | /*! | ||
| 43 | * \brief Default constructor | ||
| 44 | */ | ||
| 45 |
12/24✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 1 times.
|
12 | MonotoneCubicSpline() = default; |
| 46 | |||
| 47 | /*! | ||
| 48 | * \brief Construct a monotone cubic spline from the control points (x[i], y[i]) | ||
| 49 | * \note if the data set is monotone, monotonicity is preserved | ||
| 50 | * \param x a vector of x-coordinates | ||
| 51 | * \param y a vector of y-coordinates | ||
| 52 | */ | ||
| 53 | 7 | MonotoneCubicSpline(const std::vector<Scalar>& x, const std::vector<Scalar>& y) | |
| 54 | 21 | { | |
| 55 |
1/2✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
|
7 | updatePoints(x, y); |
| 56 | 7 | } | |
| 57 | |||
| 58 | /*! | ||
| 59 | * \brief Create a monotone cubic spline from the control points (x[i], y[i]) | ||
| 60 | * \param x a vector of x-coordinates | ||
| 61 | * \param y a vector of y-coordinates | ||
| 62 | */ | ||
| 63 | 10 | void updatePoints(const std::vector<Scalar>& x, const std::vector<Scalar>& y) | |
| 64 | { | ||
| 65 | // check some requirements | ||
| 66 |
3/6✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 10 times.
|
30 | assert (x.size() == y.size()); |
| 67 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
|
10 | assert (x.size() >=2); |
| 68 |
9/16✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 7 taken 8 times.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
|
42 | assert (std::is_sorted(x.begin(), x.end()) || std::is_sorted(x.rbegin(), x.rend())); |
| 69 | |||
| 70 | // save a copy of the control points | ||
| 71 | 10 | x_ = x; | |
| 72 | 10 | y_ = y; | |
| 73 | |||
| 74 | // the number of control points | ||
| 75 | 10 | numPoints_ = x.size(); | |
| 76 | |||
| 77 | // whether we are increasing | ||
| 78 | 20 | increasingX_ = x_.back() > x_.front(); | |
| 79 | 20 | increasingY_ = y_.back() > y_.front(); | |
| 80 | |||
| 81 | // the slope at every control point | ||
| 82 | 10 | m_.resize(numPoints_); | |
| 83 | |||
| 84 | // compute slopes (see Fritsch & Butland (1984), Eq. (5)) | ||
| 85 | 10 | Scalar deltaX = (x[1]-x[0]); | |
| 86 | 10 | Scalar secant = m_.front() = (y[1]-y[0])/deltaX; | |
| 87 | 10 | Scalar prevDeltaX = deltaX; | |
| 88 | 10 | Scalar prevSecant = secant; | |
| 89 |
2/2✓ Branch 0 taken 336 times.
✓ Branch 1 taken 10 times.
|
346 | for (int i = 1; i < numPoints_-1; ++i, prevSecant = secant, prevDeltaX = deltaX) |
| 90 | { | ||
| 91 |
2/4✓ Branch 0 taken 336 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 336 times.
✗ Branch 3 not taken.
|
672 | deltaX = (x[i+1]-x[i]); |
| 92 |
2/4✓ Branch 0 taken 336 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 336 times.
✗ Branch 3 not taken.
|
672 | secant = (y[i+1]-y[i])/deltaX; |
| 93 | 336 | const auto alpha = (prevDeltaX + 2*deltaX)/(3*(prevDeltaX + deltaX)); | |
| 94 |
1/2✓ Branch 0 taken 336 times.
✗ Branch 1 not taken.
|
336 | m_[i] = prevSecant*secant > 0.0 ? prevSecant*secant/(alpha*secant + (1.0-alpha)*prevSecant) : 0.0; |
| 95 | } | ||
| 96 | 20 | m_.back() = secant; | |
| 97 | 10 | } | |
| 98 | |||
| 99 | /*! | ||
| 100 | * \brief Evaluate the y value at a given x value | ||
| 101 | * \param x the x-coordinate | ||
| 102 | * \note We extrapolate linearly if out of bounds | ||
| 103 | */ | ||
| 104 | 5015990 | Scalar eval(const Scalar x) const | |
| 105 | { | ||
| 106 |
12/12✓ Branch 0 taken 8340 times.
✓ Branch 1 taken 5007650 times.
✓ Branch 2 taken 8340 times.
✓ Branch 3 taken 5007650 times.
✓ Branch 4 taken 8000 times.
✓ Branch 5 taken 340 times.
✓ Branch 6 taken 5007652 times.
✓ Branch 7 taken 7998 times.
✓ Branch 8 taken 5007652 times.
✓ Branch 9 taken 7998 times.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 5007650 times.
|
10031980 | if ((x <= x_.front() && increasingX_) || (x >= x_.front() && !increasingX_)) |
| 107 | 1368 | return y_.front() + m_.front()*(x - x_.front()); | |
| 108 |
11/12✓ Branch 0 taken 7658 times.
✓ Branch 1 taken 5007990 times.
✓ Branch 2 taken 7658 times.
✓ Branch 3 taken 5007990 times.
✓ Branch 4 taken 7658 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 5007982 times.
✓ Branch 7 taken 7666 times.
✓ Branch 8 taken 5007982 times.
✓ Branch 9 taken 7666 times.
✓ Branch 10 taken 334 times.
✓ Branch 11 taken 5007648 times.
|
10031296 | else if ((x > x_.back() && increasingX_) || (x < x_.back() && !increasingX_)) |
| 109 | 1336 | return y_.back() + m_.back()*(x - x_.back()); | |
| 110 | |||
| 111 | 5015314 | return eval_(x); | |
| 112 | } | ||
| 113 | |||
| 114 | /*! | ||
| 115 | * \brief Evaluate the first derivative dy/dx at a given x value | ||
| 116 | * \param x the x-coordinate | ||
| 117 | * \note We extrapolate linearly if out of bounds | ||
| 118 | */ | ||
| 119 | 4007992 | Scalar evalDerivative(const Scalar x) const | |
| 120 | { | ||
| 121 |
12/12✓ Branch 0 taken 4336 times.
✓ Branch 1 taken 4003656 times.
✓ Branch 2 taken 4336 times.
✓ Branch 3 taken 4003656 times.
✓ Branch 4 taken 4000 times.
✓ Branch 5 taken 336 times.
✓ Branch 6 taken 4003658 times.
✓ Branch 7 taken 3998 times.
✓ Branch 8 taken 4003658 times.
✓ Branch 9 taken 3998 times.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 4003656 times.
|
8015984 | if ((x <= x_.front() && increasingX_) || (x >= x_.front() && !increasingX_)) |
| 122 | 676 | return m_.front(); | |
| 123 |
11/12✓ Branch 0 taken 3662 times.
✓ Branch 1 taken 4003992 times.
✓ Branch 2 taken 3662 times.
✓ Branch 3 taken 4003992 times.
✓ Branch 4 taken 3662 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 4003988 times.
✓ Branch 7 taken 3666 times.
✓ Branch 8 taken 4003988 times.
✓ Branch 9 taken 3666 times.
✓ Branch 10 taken 334 times.
✓ Branch 11 taken 4003654 times.
|
8015308 | else if ((x > x_.back() && increasingX_) || (x < x_.back() && !increasingX_)) |
| 124 | 668 | return m_.back(); | |
| 125 | |||
| 126 | 4007320 | return evalDerivative_(x); | |
| 127 | } | ||
| 128 | |||
| 129 | /*! | ||
| 130 | * \brief Evaluate the inverse function | ||
| 131 | * \param y the y-coordinate | ||
| 132 | * \note We extrapolate linearly if out of bounds | ||
| 133 | * \note Throws exception if inverse could not be found (e.g. not unique) | ||
| 134 | */ | ||
| 135 | 2007996 | Scalar evalInverse(const Scalar y) const | |
| 136 | { | ||
| 137 |
12/12✓ Branch 0 taken 2003998 times.
✓ Branch 1 taken 3998 times.
✓ Branch 2 taken 2003998 times.
✓ Branch 3 taken 3998 times.
✓ Branch 4 taken 2003996 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 4000 times.
✓ Branch 7 taken 2003994 times.
✓ Branch 8 taken 4000 times.
✓ Branch 9 taken 2003994 times.
✓ Branch 10 taken 2 times.
✓ Branch 11 taken 3998 times.
|
4015992 | if ((y <= y_.front() && increasingY_) || (y >= y_.front() && !increasingY_)) |
| 138 | 16 | return x_.front() + (y - y_.front())/m_.front(); | |
| 139 |
10/12✓ Branch 0 taken 2003992 times.
✓ Branch 1 taken 4000 times.
✓ Branch 2 taken 2003992 times.
✓ Branch 3 taken 4000 times.
✓ Branch 4 taken 2003992 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3996 times.
✓ Branch 7 taken 2003996 times.
✓ Branch 8 taken 3996 times.
✓ Branch 9 taken 2003996 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 3996 times.
|
4015984 | else if ((y > y_.back() && increasingY_) || (y < y_.back() && !increasingY_)) |
| 140 | ✗ | return x_.back() + (y - y_.back())/m_.back(); | |
| 141 | |||
| 142 | 2007992 | return evalInverse_(y); | |
| 143 | } | ||
| 144 | |||
| 145 | private: | ||
| 146 | 5015314 | Scalar eval_(const Scalar x) const | |
| 147 | { | ||
| 148 | // interpolate parametrization parameter t in [0,1] | ||
| 149 | 5015314 | const auto lookUpIndex = lookUpIndex_(x_, x, increasingX_); | |
| 150 | 15045942 | const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]); | |
| 151 | 5015314 | const auto t = (x - x_[lookUpIndex-1])/h; | |
| 152 | 25076570 | return y_[lookUpIndex-1]*Basis::h00(t) + h*m_[lookUpIndex-1]*Basis::h10(t) | |
| 153 | 25076570 | + y_[lookUpIndex]*Basis::h01(t) + h*m_[lookUpIndex]*Basis::h11(t); | |
| 154 | } | ||
| 155 | |||
| 156 | 4007320 | Scalar evalDerivative_(const Scalar x) const | |
| 157 | { | ||
| 158 | // interpolate parametrization parameter t in [0,1] | ||
| 159 | 4007320 | const auto lookUpIndex = lookUpIndex_(x_, x, increasingX_); | |
| 160 | 12021960 | const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]); | |
| 161 | 4007320 | const auto t = (x - x_[lookUpIndex-1])/h; | |
| 162 | 4007320 | const auto dtdx = 1.0/h; | |
| 163 | 20036600 | return y_[lookUpIndex-1]*Basis::dh00(t)*dtdx + m_[lookUpIndex-1]*Basis::dh10(t) | |
| 164 | 20036600 | + y_[lookUpIndex]*Basis::dh01(t)*dtdx + m_[lookUpIndex]*Basis::dh11(t); | |
| 165 | } | ||
| 166 | |||
| 167 | 2007992 | Scalar evalInverse_(const Scalar y) const | |
| 168 | { | ||
| 169 | 2007992 | const auto lookUpIndex = lookUpIndex_(y_, y, increasingY_); | |
| 170 | 49995640 | auto localPolynomial = [&](const auto x) { | |
| 171 | // interpolate parametrization parameter t in [0,1] | ||
| 172 | 143962944 | const auto h = (x_[lookUpIndex] - x_[lookUpIndex-1]); | |
| 173 | 47987648 | const auto t = (x - x_[lookUpIndex-1])/h; | |
| 174 | 239938240 | return y - (y_[lookUpIndex-1]*Basis::h00(t) + h*m_[lookUpIndex-1]*Basis::h10(t) | |
| 175 | 239938240 | + y_[lookUpIndex]*Basis::h01(t) + h*m_[lookUpIndex]*Basis::h11(t)); | |
| 176 | }; | ||
| 177 | |||
| 178 | // use an epsilon for the bracket | ||
| 179 | 4015984 | const auto eps = (x_[lookUpIndex]-x_[lookUpIndex-1])*1e-5; | |
| 180 | 2007992 | return findScalarRootBrent(x_[lookUpIndex-1]-eps, x_[lookUpIndex]+eps, localPolynomial); | |
| 181 | } | ||
| 182 | |||
| 183 | 11030626 | auto lookUpIndex_(const std::vector<Scalar>& vec, const Scalar v, bool increasing) const | |
| 184 | { | ||
| 185 |
2/2✓ Branch 0 taken 9015304 times.
✓ Branch 1 taken 2015322 times.
|
11030626 | return increasing ? lookUpIndexIncreasing_(vec, v) : lookUpIndexDecreasing_(vec, v); |
| 186 | } | ||
| 187 | |||
| 188 | 9015304 | auto lookUpIndexIncreasing_(const std::vector<Scalar>& vec, const Scalar v) const | |
| 189 | { | ||
| 190 |
2/4✓ Branch 0 taken 9015304 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 9015304 times.
✗ Branch 3 not taken.
|
45076520 | const auto lookUpIndex = std::distance(vec.begin(), std::lower_bound(vec.begin(), vec.end(), v)); |
| 191 |
3/6✓ Branch 0 taken 9015304 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 9015304 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 9015304 times.
|
9015304 | assert(lookUpIndex != 0 && lookUpIndex < vec.size()); |
| 192 | 9015304 | return lookUpIndex; | |
| 193 | } | ||
| 194 | |||
| 195 | 2015322 | auto lookUpIndexDecreasing_(const std::vector<Scalar>& vec, const Scalar v) const | |
| 196 | { | ||
| 197 |
2/4✓ Branch 4 taken 2015322 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2015322 times.
✗ Branch 7 not taken.
|
8061288 | const auto lookUpIndex = vec.size() - std::distance(vec.rbegin(), std::upper_bound(vec.rbegin(), vec.rend(), v)); |
| 198 |
3/6✓ Branch 0 taken 2015322 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2015322 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 2015322 times.
|
2015322 | assert(lookUpIndex != 0 && lookUpIndex < vec.size()); |
| 199 | 2015322 | return lookUpIndex; | |
| 200 | } | ||
| 201 | |||
| 202 | std::vector<Scalar> x_; //!< the x-coordinates | ||
| 203 | std::vector<Scalar> y_; //!< the y-coordinates | ||
| 204 | std::vector<Scalar> m_; //!< the slope for each control point | ||
| 205 | std::size_t numPoints_; //!< the number of control points | ||
| 206 | bool increasingX_; //!< if we are increasing monotone or not | ||
| 207 | bool increasingY_; //!< if we are increasing monotone or not | ||
| 208 | }; | ||
| 209 | |||
| 210 | } // end namespace Dumux | ||
| 211 | |||
| 212 | #endif | ||
| 213 |