3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
cubicspline.hh
Go to the documentation of this file.
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 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_COMMON_CUBIC_SPLINE_HH
25#define DUMUX_COMMON_CUBIC_SPLINE_HH
26
27#include <iterator>
28#include <vector>
29#include <algorithm>
30
31#include <dune/common/fvector.hh>
32#include <dune/common/fmatrix.hh>
33#include <dune/istl/btdmatrix.hh>
34#include <dune/istl/bvector.hh>
35
36namespace Dumux {
37
43template<class Scalar = double>
45{
46public:
50 CubicSpline() = default;
51
57 CubicSpline(const std::vector<Scalar>& x, const std::vector<Scalar> y)
58 {
59 // check some requirements
60 assert (x.size() == y.size());
61 assert (x.size() >=2);
62 assert (std::is_sorted(x.begin(), x.end()));
63
64 updatePoints(x, y);
65 }
66
73 void updatePoints(const std::vector<Scalar>& x, const std::vector<Scalar>& y)
74 {
75 // save a copy of the control points
76 x_ = x;
77
78 // the number of control points
79 numPoints_ = x.size();
80
81 const auto numSegments = numPoints_-1;
82 coeff_.resize(numSegments*4+2); // 4 coefficients for each segment + last y-value and derivative
83
84 // construct a block-tridiagonal matrix system solving for the derivatives
85 Dune::BTDMatrix<Dune::FieldMatrix<double, 1, 1>> matrix(numPoints_);
86 Dune::BlockVector<Dune::FieldVector<double, 1>> rhs(numPoints_);
87 Dune::BlockVector<Dune::FieldVector<double, 1>> d(numPoints_);
88
89 // assemble matrix and rhs row-wise
90 matrix[0][0] = 2.0;
91 matrix[0][1] = 1.0;
92 rhs[0] = 3.0*(y[1] - y[0]);
93 for (int i = 1; i < numPoints_-1; ++i)
94 {
95 matrix[i][i-1] = 1.0;
96 matrix[i][i] = 4.0;
97 matrix[i][i+1] = 1.0;
98 rhs[i] = 3.0*(y[i+1] - y[i-1]);
99 }
100 matrix[numPoints_-1][numPoints_-1] = 2.0;
101 matrix[numPoints_-1][numPoints_-2] = 1.0;
102 rhs[numPoints_-1] = 3.0*(y[numPoints_-1] - y[numPoints_-2]);
103
104 // solve for derivatives
105 matrix.solve(d, rhs);
106
107 // compute coefficients
108 std::size_t offset = 0;
109 for (int i = 0; i < numSegments; ++i, offset += 4)
110 {
111 coeff_[offset+0] = y[i];
112 coeff_[offset+1] = d[i];
113 coeff_[offset+2] = 3.0*(y[i+1]-y[i]) - 2.0*d[i] - d[i+1];
114 coeff_[offset+3] = 2.0*(y[i]-y[i+1]) + d[i] + d[i+1];
115 }
116 coeff_[offset+0] = y[numPoints_-1];
117 coeff_[offset+1] = d[numPoints_-1];
118 }
119
125 Scalar eval(const Scalar x) const
126 {
127 if (x <= x_[0])
128 return coeff_[0] + evalDerivative(x_[0])*(x - x_[0]);
129 else if (x > x_[numPoints_-1])
130 return coeff_[(numPoints_-1)*4] + evalDerivative(x_[numPoints_-1])*(x - x_[numPoints_-1]);
131 else
132 {
133 const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
134 assert(lookUpIndex != 0);
135
136 // get coefficients
137 const auto* coeff = coeff_.data() + (lookUpIndex-1)*4;
138
139 // interpolate parametrization parameter t in [0,1]
140 const auto t = (x - x_[lookUpIndex-1])/(x_[lookUpIndex] - x_[lookUpIndex-1]);
141 return coeff[0] + t*(coeff[1] + t*coeff[2] + t*t*coeff[3]);
142 }
143 }
144
150 Scalar evalDerivative(const Scalar x) const
151 {
152 if (x <= x_[0])
153 return coeff_[1]/(x_[1] - x_[0]);
154 else if (x > x_[numPoints_-1])
155 return coeff_[(numPoints_-1)*4 + 1]/(x_[numPoints_-1] - x_[numPoints_-2]);
156 else
157 {
158 const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
159 assert(lookUpIndex != 0);
160
161 // get coefficients
162 const auto* coeff = coeff_.data() + (lookUpIndex-1)*4;
163
164 // interpolate parametrization parameter t in [0,1]
165 const auto t = (x - x_[lookUpIndex-1])/(x_[lookUpIndex] - x_[lookUpIndex-1]);
166 const auto dtdx = 1.0/(x_[lookUpIndex] - x_[lookUpIndex-1]);
167 return dtdx*(coeff[1] + t*(2.0*coeff[2] + t*3.0*coeff[3]));
168 }
169 }
170
171private:
172 std::vector<Scalar> x_;
173 std::size_t numPoints_;
174 std::vector<Scalar> coeff_;
175};
176
177} // end namespace Dumux
178
179#endif
ctype distance(const Dune::FieldVector< ctype, dimWorld > &a, const Dune::FieldVector< ctype, dimWorld > &b)
Compute the shortest distance between two points.
Definition: geometry/distance.hh:138
Definition: adapt.hh:29
A simple implementation of a natural cubic spline.
Definition: cubicspline.hh:45
CubicSpline(const std::vector< Scalar > &x, const std::vector< Scalar > y)
Contruct a natural cubic spline from the control points (x[i], y[i])
Definition: cubicspline.hh:57
CubicSpline()=default
Default constructor.
Scalar eval(const Scalar x) const
Evaluate the y value at a given x value.
Definition: cubicspline.hh:125
Scalar evalDerivative(const Scalar x) const
Evaluate the first derivative dy/dx at a given x value.
Definition: cubicspline.hh:150
void updatePoints(const std::vector< Scalar > &x, const std::vector< Scalar > &y)
Create a natural cubic spline from the control points (x[i], y[i])
Definition: cubicspline.hh:73