From 5a82bd3923571fa50e4cee2292455b2d5877426c Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Sat, 18 May 2019 02:45:58 +0200
Subject: [PATCH] [common] Implement simple natural cubic spline
---
dumux/common/CMakeLists.txt | 1 +
dumux/common/cubicspline.hh | 181 ++++++++++++++++++++++++++++++++++++
2 files changed, 182 insertions(+)
create mode 100644 dumux/common/cubicspline.hh
diff --git a/dumux/common/CMakeLists.txt b/dumux/common/CMakeLists.txt
index 4ec3919ec0..c743fc2994 100644
--- a/dumux/common/CMakeLists.txt
+++ b/dumux/common/CMakeLists.txt
@@ -8,6 +8,7 @@ boundaryconditions.hh
boundaryflag.hh
boundarytypes.hh
boundingboxtree.hh
+cubicspline.hh
defaultmappertraits.hh
defaultusagemessage.hh
dimensionlessnumbers.hh
diff --git a/dumux/common/cubicspline.hh b/dumux/common/cubicspline.hh
new file mode 100644
index 0000000000..ed6739194d
--- /dev/null
+++ b/dumux/common/cubicspline.hh
@@ -0,0 +1,181 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * 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 3 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
+ * \ingroup Common
+ * \brief A simple implementation of a cubic spline
+ */
+#ifndef DUMUX_COMMON_CUBIC_SPLINE_HH
+#define DUMUX_COMMON_CUBIC_SPLINE_HH
+
+#include <iterator>
+#include <vector>
+#include <algorithm>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/btdmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Common
+ * \brief A simple implementation of a natural cubic spline
+ * \note We follow the notation at http://mathworld.wolfram.com/CubicSpline.html
+ */
+template<class Scalar = double>
+class CubicSpline
+{
+public:
+ /*!
+ * \brief Default constructor
+ */
+ CubicSpline() = default;
+
+ /*!
+ * \brief Contruct a natural cubic spline from the control points (x[i], y[i])
+ * \param x a vector of x-coordinates
+ * \param y a vector of y-coordinates
+ */
+ CubicSpline(const std::vector<Scalar>& x, const std::vector<Scalar> y)
+ {
+ // check some requirements
+ assert (x.size() == y.size());
+ assert (x.size() >=2);
+ assert (std::is_sorted(x.begin(), x.end()));
+
+ updatePoints(x, y);
+ }
+
+ /*!
+ * \brief Create a natural cubic spline from the control points (x[i], y[i])
+ * \note we enforce continuous second derivatives in the inside and zero second derivatives at the boundary
+ * \param x a vector of x-coordinates
+ * \param y a vector of y-coordinates
+ */
+ void updatePoints(const std::vector<Scalar>& x, const std::vector<Scalar>& y)
+ {
+ // save a copy of the control points
+ x_ = x;
+
+ // the number of control points
+ numPoints_ = x.size();
+
+ const auto numSegments = numPoints_-1;
+ coeff_.resize(numSegments*4+2); // 4 coefficients for each segment + last y-value and derivative
+
+ // construct a block-tridiagonal matrix system solving for the derivatives
+ Dune::BTDMatrix<Dune::FieldMatrix<double, 1, 1>> matrix(numPoints_);
+ Dune::BlockVector<Dune::FieldVector<double, 1>> rhs(numPoints_);
+ Dune::BlockVector<Dune::FieldVector<double, 1>> d(numPoints_);
+
+ // assemble matrix and rhs row-wise
+ matrix[0][0] = 2.0;
+ matrix[0][1] = 1.0;
+ rhs[0] = 3.0*(y[1] - y[0]);
+ for (int i = 1; i < numPoints_-1; ++i)
+ {
+ matrix[i][i-1] = 1.0;
+ matrix[i][i] = 4.0;
+ matrix[i][i+1] = 1.0;
+ rhs[i] = 3.0*(y[i+1] - y[i-1]);
+ }
+ matrix[numPoints_-1][numPoints_-1] = 2.0;
+ matrix[numPoints_-1][numPoints_-2] = 1.0;
+ rhs[numPoints_-1] = 3.0*(y[numPoints_-1] - y[numPoints_-2]);
+
+ // solve for derivatives
+ matrix.solve(d, rhs);
+
+ // compute coefficients
+ std::size_t offset = 0;
+ for (int i = 0; i < numSegments; ++i, offset += 4)
+ {
+ coeff_[offset+0] = y[i];
+ coeff_[offset+1] = d[i];
+ coeff_[offset+2] = 3.0*(y[i+1]-y[i]) - 2.0*d[i] - d[i+1];
+ coeff_[offset+3] = 2.0*(y[i]-y[i+1]) + d[i] + d[i+1];
+ }
+ coeff_[offset+0] = y[numPoints_-1];
+ coeff_[offset+1] = d[numPoints_-1];
+ }
+
+ /*!
+ * \brief Evaluate the y value at a given x value
+ * \param x the x-coordinate
+ * \note We extrapolate linearly if out of bounds
+ */
+ Scalar eval(const Scalar x) const
+ {
+ if (x < x_[0])
+ return coeff_[0] + evalDerivative(x_[0])*(x - x_[0]);
+ else if (x > x_[numPoints_-1])
+ return coeff_[(numPoints_-1)*4] + evalDerivative(x_[numPoints_-1])*(x - x_[numPoints_-1]);
+ else
+ {
+ const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
+ if (lookUpIndex == 0)
+ return coeff_[0];
+
+ // get coefficients
+ const auto* coeff = coeff_.data() + (lookUpIndex-1)*4;
+
+ // interpolate parametrization parameter t in [0,1]
+ const auto t = (x - x_[lookUpIndex-1])/(x_[lookUpIndex] - x_[lookUpIndex-1]);
+ return coeff[0] + t*(coeff[1] + t*coeff[2] + t*t*coeff[3]);
+ }
+ }
+
+ /*!
+ * \brief Evaluate the first derivative dy/dx at a given x value
+ * \param x the x-coordinate
+ * \note We extrapolate linearly if out of bounds
+ */
+ Scalar evalDerivative(const Scalar x) const
+ {
+ if (x < x_[0])
+ return coeff_[1]/(x_[1] - x_[0]);
+ else if (x > x_[numPoints_-1])
+ return coeff_[(numPoints_-1)*4 + 1]/(x_[numPoints_-1] - x_[numPoints_-2]);
+ else
+ {
+ const auto lookUpIndex = std::distance(x_.begin(), std::lower_bound(x_.begin(), x_.end(), x));
+ if (lookUpIndex == 0)
+ return coeff_[1]/(x_[1] - x_[0]);
+
+ // get coefficients
+ const auto* coeff = coeff_.data() + (lookUpIndex-1)*4;
+
+ // interpolate parametrization parameter t in [0,1]
+ const auto t = (x - x_[lookUpIndex-1])/(x_[lookUpIndex] - x_[lookUpIndex-1]);
+ const auto dtdx = 1.0/(x_[lookUpIndex] - x_[lookUpIndex-1]);
+ return dtdx*(coeff[1] + t*(2.0*coeff[2] + t*3.0*coeff[3]));
+ }
+ }
+
+private:
+ std::vector<Scalar> x_; //!< the x-coordinates
+ std::size_t numPoints_; //!< the number of control points
+ std::vector<Scalar> coeff_; //!< the spline coefficients
+};
+
+} // end namespace Dumux
+
+#endif
--
GitLab