GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/common/cubicspline.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 56 56 100.0%
Functions: 4 4 100.0%
Branches: 61 110 55.5%
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 simple implementation of a cubic spline
11 */
12 #ifndef DUMUX_COMMON_CUBIC_SPLINE_HH
13 #define DUMUX_COMMON_CUBIC_SPLINE_HH
14
15 #include <iterator>
16 #include <vector>
17 #include <algorithm>
18
19 #include <dune/common/fvector.hh>
20 #include <dune/common/fmatrix.hh>
21 #include <dune/istl/btdmatrix.hh>
22 #include <dune/istl/bvector.hh>
23
24 namespace Dumux {
25
26 /*!
27 * \ingroup Core
28 * \brief A simple implementation of a natural cubic spline
29 * \note We follow the notation at http://mathworld.wolfram.com/CubicSpline.html
30 */
31 template<class Scalar = double>
32 class CubicSpline
33 {
34 public:
35 /*!
36 * \brief Default constructor
37 */
38 CubicSpline() = default;
39
40 /*!
41 * \brief Construct a natural cubic spline from the control points (x[i], y[i])
42 * \param x a vector of x-coordinates
43 * \param y a vector of y-coordinates
44 */
45 1 CubicSpline(const std::vector<Scalar>& x, const std::vector<Scalar> y)
46 2 {
47 // check some requirements
48
3/6
✗ 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.
 3 assert (x.size() == y.size());
49
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 assert (x.size() >=2);
50
4/8
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
 4 assert (std::is_sorted(x.begin(), x.end()));
51
52
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 updatePoints(x, y);
53 1 }
54
55 /*!
56 * \brief Create a natural cubic spline from the control points (x[i], y[i])
57 * \note we enforce continuous second derivatives in the inside and zero second derivatives at the boundary
58 * \param x a vector of x-coordinates
59 * \param y a vector of y-coordinates
60 */
61 1 void updatePoints(const std::vector<Scalar>& x, const std::vector<Scalar>& y)
62 {
63 // save a copy of the control points
64 1 x_ = x;
65
66 // the number of control points
67 1 numPoints_ = x.size();
68
69 1 const auto numSegments = numPoints_-1;
70 1 coeff_.resize(numSegments*4+2); // 4 coefficients for each segment + last y-value and derivative
71
72 // construct a block-tridiagonal matrix system solving for the derivatives
73 1 Dune::BTDMatrix<Dune::FieldMatrix<double, 1, 1>> matrix(numPoints_);
74
1/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2 Dune::BlockVector<Dune::FieldVector<double, 1>> rhs(numPoints_);
75
2/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 2 Dune::BlockVector<Dune::FieldVector<double, 1>> d(numPoints_);
76
77 // assemble matrix and rhs row-wise
78
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 matrix[0][0] = 2.0;
79
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 matrix[0][1] = 1.0;
80 1 rhs[0] = 3.0*(y[1] - y[0]);
81
2/2
✓ Branch 0 taken 13 times.
✓ Branch 1 taken 1 times.
 14 for (int i = 1; i < numPoints_-1; ++i)
82 {
83
3/6
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 13 times.
✗ Branch 8 not taken.
 26 matrix[i][i-1] = 1.0;
84
3/6
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 13 times.
✗ Branch 8 not taken.
 26 matrix[i][i] = 4.0;
85
2/4
✓ Branch 1 taken 13 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 13 times.
✗ Branch 5 not taken.
 26 matrix[i][i+1] = 1.0;
86 65 rhs[i] = 3.0*(y[i+1] - y[i-1]);
87 }
88
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 2 matrix[numPoints_-1][numPoints_-1] = 2.0;
89
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 2 matrix[numPoints_-1][numPoints_-2] = 1.0;
90
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
 4 rhs[numPoints_-1] = 3.0*(y[numPoints_-1] - y[numPoints_-2]);
91
92 // solve for derivatives
93
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 matrix.solve(d, rhs);
94
95 // compute coefficients
96 std::size_t offset = 0;
97
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 1 times.
 15 for (int i = 0; i < numSegments; ++i, offset += 4)
98 {
99 28 coeff_[offset+0] = y[i];
100 28 coeff_[offset+1] = d[i];
101 84 coeff_[offset+2] = 3.0*(y[i+1]-y[i]) - 2.0*d[i] - d[i+1];
102 56 coeff_[offset+3] = 2.0*(y[i]-y[i+1]) + d[i] + d[i+1];
103 }
104
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 coeff_[offset+0] = y[numPoints_-1];
105
3/6
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 3 coeff_[offset+1] = d[numPoints_-1];
106 1 }
107
108 /*!
109 * \brief Evaluate the y value at a given x value
110 * \param x the x-coordinate
111 * \note We extrapolate linearly if out of bounds
112 */
113 2000 Scalar eval(const Scalar x) const
114 {
115
2/2
✓ Branch 0 taken 251 times.
✓ Branch 1 taken 1749 times.
 2000 if (x <= x_[0])
116 251 return coeff_[0] + evalDerivative(x_[0])*(x - x_[0]);
117
4/4
✓ Branch 0 taken 250 times.
✓ Branch 1 taken 1499 times.
✓ Branch 2 taken 250 times.
✓ Branch 3 taken 1499 times.
 3498 else if (x > x_[numPoints_-1])
118 500 return coeff_[(numPoints_-1)*4] + evalDerivative(x_[numPoints_-1])*(x - x_[numPoints_-1]);
119 else
120 {
121
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1499 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1499 times.
 7495 const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
122
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1499 times.
 1499 assert(lookUpIndex != 0);
123
124 // get coefficients
125 1499 const auto* coeff = coeff_.data() + (lookUpIndex-1)*4;
126
127 // interpolate parametrization parameter t in [0,1]
128 1499 const auto t = (x - x_[lookUpIndex-1])/(x_[lookUpIndex] - x_[lookUpIndex-1]);
129 1499 return coeff[0] + t*(coeff[1] + t*coeff[2] + t*t*coeff[3]);
130 }
131 }
132
133 /*!
134 * \brief Evaluate the first derivative dy/dx at a given x value
135 * \param x the x-coordinate
136 * \note We extrapolate linearly if out of bounds
137 */
138 2501 Scalar evalDerivative(const Scalar x) const
139 {
140
2/2
✓ Branch 0 taken 502 times.
✓ Branch 1 taken 1999 times.
 2501 if (x <= x_[0])
141 502 return coeff_[1]/(x_[1] - x_[0]);
142
4/4
✓ Branch 0 taken 250 times.
✓ Branch 1 taken 1749 times.
✓ Branch 2 taken 250 times.
✓ Branch 3 taken 1749 times.
 3998 else if (x > x_[numPoints_-1])
143 750 return coeff_[(numPoints_-1)*4 + 1]/(x_[numPoints_-1] - x_[numPoints_-2]);
144 else
145 {
146
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1749 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1749 times.
 8745 const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
147
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1749 times.
 1749 assert(lookUpIndex != 0);
148
149 // get coefficients
150 1749 const auto* coeff = coeff_.data() + (lookUpIndex-1)*4;
151
152 // interpolate parametrization parameter t in [0,1]
153 1749 const auto t = (x - x_[lookUpIndex-1])/(x_[lookUpIndex] - x_[lookUpIndex-1]);
154 1749 const auto dtdx = 1.0/(x_[lookUpIndex] - x_[lookUpIndex-1]);
155 1749 return dtdx*(coeff[1] + t*(2.0*coeff[2] + t*3.0*coeff[3]));
156 }
157 }
158
159 private:
160 std::vector<Scalar> x_; //!< the x-coordinates
161 std::size_t numPoints_; //!< the number of control points
162 std::vector<Scalar> coeff_; //!< the spline coefficients
163 };
164
165 } // end namespace Dumux
166
167 #endif
168