diff --git a/dumux/common/CMakeLists.txt b/dumux/common/CMakeLists.txt
index 4ec3919ec0e3da696994947275691dae51a2f360..c743fc2994c169febabdbca4d545626b763fde93 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 0000000000000000000000000000000000000000..bf06bae262131faed7d5f72b7407fe3956b1a39b
--- /dev/null
+++ b/dumux/common/cubicspline.hh
@@ -0,0 +1,179 @@
+// -*- 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));
+ assert(lookUpIndex != 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));
+ assert(lookUpIndex != 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
diff --git a/test/common/spline/CMakeLists.txt b/test/common/spline/CMakeLists.txt
index 004a95c4e21d0e3ca120bffd70f6d79705387f33..1aafa3e3132a5010a0aea0d9f640c8405922b483 100644
--- a/test/common/spline/CMakeLists.txt
+++ b/test/common/spline/CMakeLists.txt
@@ -1,2 +1,4 @@
dune_add_test(SOURCES test_spline.cc
LABELS unit)
+dune_add_test(SOURCES test_cubicspline.cc
+ LABELS unit)
diff --git a/test/common/spline/test_cubicspline.cc b/test/common/spline/test_cubicspline.cc
new file mode 100644
index 0000000000000000000000000000000000000000..cad29e2f3fe6bb9bd23ba1e84257d39b0228df72
--- /dev/null
+++ b/test/common/spline/test_cubicspline.cc
@@ -0,0 +1,102 @@
+// -*- 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
+ * \brief Test the simple cubic spline implementation
+ */
+#include <config.h>
+#include <cmath>
+#include <vector>
+#include <algorithm>
+#include <functional>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dumux/common/cubicspline.hh>
+#include <dumux/io/gnuplotinterface.hh>
+
+std::vector<double> linspace(const double begin, const double end, const double samples)
+{
+ const double delta = (end-begin)/static_cast<double>(samples-1);
+ std::vector<double> vec(samples);
+ for (int i = 0; i < samples; ++i)
+ vec[i] = begin + i*delta;
+ return vec;
+}
+
+template<class Function>
+std::vector<double> eval(const Function& f, const std::vector<double>& x)
+{
+ auto y = x;
+ std::transform(x.begin(), x.end(), y.begin(), [&](const double x) { return f(x); });
+ return y;
+}
+
+int main(int argc, char** argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ // we test the spline interpolation against a sample function
+ const auto f = [](double x){ return 1.0 / ( 1.0 + x*x ); };
+ const auto df = [](double x){ return -2.0*x/(( 1.0 + x*x )*( 1.0 + x*x )); };
+
+ // create some test samples
+ const auto testPoints = linspace(-4.0, 4.0, 1000);
+ const auto ref = eval(f, testPoints);
+ const auto refDeriv = eval(df, testPoints);
+
+ // create the spline sample points
+ const auto samplePoints = linspace(-4.0, 4.0, 15);
+ const auto y = eval(f, samplePoints);
+
+ // create the spline
+ Dumux::CubicSpline spline(samplePoints, y);
+
+ // evaluate spline and derivative
+ const auto result = eval([&](const double x) { return spline.eval(x); }, testPoints);
+ const auto resultDeriv = eval([&](const double x) { return spline.evalDerivative(x); }, testPoints);
+
+ // compute largest difference
+ auto diff = result; auto diffDeriv = result;
+ std::transform(result.begin(), result.end(), ref.begin(), diff.begin(), [](auto a, auto b){ return std::abs(a-b); });
+ std::transform(resultDeriv.begin(), resultDeriv.end(), refDeriv.begin(), diffDeriv.begin(), [](auto a, auto b){ return std::abs(a-b); });
+ const auto maxNorm = std::accumulate(diff.begin(), diff.end(), diff[0], [](auto a, auto b){ return std::max(a, b); });
+ const auto maxNormDeriv = std::accumulate(diffDeriv.begin(), diffDeriv.end(), diffDeriv[0], [](auto a, auto b){ return std::max(a, b); });
+ std::cout << "Maximum error: " << maxNorm << "\n";
+ std::cout << "Maximum error in derivative: " << maxNormDeriv << "\n";
+
+ if (maxNorm > 0.0038 || maxNormDeriv > 0.021)
+ DUNE_THROW(Dune::Exception, "Maximum error in spline interpolation too large!");
+
+ // plot with Gnuplot (plot a bit more so we can see the linear extension)
+ const auto plotPoints = linspace(-8.0, 8.0, 1000);
+ const auto refPlot = eval(f, plotPoints);
+ const auto refDerivPlot = eval(df, plotPoints);
+ const auto resultPlot = eval([&](const double x) { return spline.eval(x); }, plotPoints);
+ const auto resultDerivPlot = eval([&](const double x) { return spline.evalDerivative(x); }, plotPoints);
+
+ Dumux::GnuplotInterface<double> gnuplot(false);
+ gnuplot.addDataSetToPlot(plotPoints, refPlot, "reference");
+ gnuplot.addDataSetToPlot(plotPoints, refDerivPlot, "reference_derivative");
+ gnuplot.addDataSetToPlot(plotPoints, resultPlot, "spline");
+ gnuplot.addDataSetToPlot(plotPoints, resultDerivPlot, "spline_derivative");
+ gnuplot.plot("spline");
+
+ return 0;
+}