add Functions::CSpline

diff --git a/include/deal.II/base/function_cspline.h b/include/deal.II/base/function_cspline.h
new file mode 100644 (file)
index 0000000..fac6269
--- /dev/null
+++ b/include/deal.II/base/function_cspline.h
@@ -0,0 +1,113 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__function_cspline_h
+#define dealii__function_cspline_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/point.h>
+#include <gsl/gsl_spline.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace Functions
+{
+  DeclException3 (ExcCSplineOrder,
+                  int,
+                  double,
+                  double,
+                  << "CSpline: x[" << arg1 << "] = "<< arg2 <<" >= x["<<(arg1+1)<<" = "<<arg3
+                 );
+
+  DeclException3 (ExcCSplineRange,
+                  double,
+                  double,
+                  double,
+                  << "CSpline: " << arg1 << " is not in ["<< arg2<<";"<<arg3<<"]."
+                 );
+  /**
+   * The cubic spline function using GNU Scientific Library.
+   * The resulting curve is piecewise cubic on each interval, with matching first
+   * and second derivatives at the supplied data-points. The second derivative
+   * is chosen to be zero at the first point and last point.
+   *
+   * @note This function is only implemented for dim==1 .
+   *
+   * @author Denis Davydov
+   * @date 2016
+   */
+  template <int dim>
+  class CSpline : public Function<dim>
+  {
+  public:
+    /**
+     * Constructor which should be provided with a set of points at which interpolation is to be done @p x
+     * and a set of function values @p f .
+     */
+    CSpline(const std::vector<double> &x,
+            const std::vector<double> &f);
+
+    /**
+     * Virtual destructor.
+     */
+    virtual ~CSpline();
+
+    virtual double value (const Point<dim> &point,
+                          const unsigned int component = 0) const;
+
+    virtual Tensor<1,dim> gradient (const Point<dim>   &p,
+                                    const unsigned int  component = 0) const;
+
+    std::size_t memory_consumption () const;
+
+  private:
+    /**
+     * Points at which interpolation is provided
+     */
+    const std::vector<double> x;
+
+    /**
+     * Values of the function at interpolation points
+     */
+    const std::vector<double> y;
+
+    /**
+     * Maximum allowed value of x
+     */
+    const double x_max;
+
+    /**
+     * Minimum allowed value of x
+     */
+    const double x_min;
+
+    /**
+     * GSL accelerator for spline interpolation
+     */
+    gsl_interp_accel *acc;
+
+    /**
+     * GSL cubic spline object
+     */
+    gsl_spline *cspline;
+  };
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+

diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt
index d780f23651b0c3310f6020f9a4742120a11b6f6f..f13a47d856c36ef2fb5ddefbb81d43b8af4f93ac 100644 (file)
--- a/source/base/CMakeLists.txt
+++ b/source/base/CMakeLists.txt
@@ -25,6 +25,7 @@ SET(_src
   exceptions.cc
   flow_function.cc
   function.cc
+  function_cspline.cc
   function_derivative.cc
   function_lib.cc
   function_lib_cutoff.cc

diff --git a/source/base/function_cspline.cc b/source/base/function_cspline.cc
new file mode 100644 (file)
index 0000000..a74e311
--- /dev/null
+++ b/source/base/function_cspline.cc
@@ -0,0 +1,99 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/point.h>
+#include <deal.II/base/function_cspline.h>
+
+#include <cmath>
+#include <algorithm>
+
+DEAL_II_NAMESPACE_OPEN
+namespace Functions
+{
+
+  template <int dim>
+  CSpline<dim>::CSpline(const std::vector<double> &x_,
+                        const std::vector<double> &y_)
+    :
+    x(x_),
+    y(y_),
+    x_max(x.size()>0? x.back()  : 0.),
+    x_min(x.size()>0? x.front() : 0.)
+  {
+    // check that input vector @p x is provided in ascending order:
+    for (unsigned int i = 0; i < x.size()-1; i++)
+      AssertThrow(x[i]<x[i+1],
+                  ExcCSplineOrder(i,x[i],x[i+1]));
+
+    acc = gsl_interp_accel_alloc ();
+    const unsigned int n = x.size();
+    cspline = gsl_spline_alloc (gsl_interp_cspline, n);
+    // gsl_spline_init returns something but it seems nobody knows what
+    const int ierr = gsl_spline_init (cspline, &x[0], &y[0], n);
+  }
+
+  template <int dim>
+  CSpline<dim>::~CSpline()
+  {
+    gsl_interp_accel_free (acc);
+    gsl_spline_free (cspline);
+  }
+
+
+  template <int dim>
+  double
+  CSpline<dim>::value (const Point<dim>   &p,
+                       const unsigned int) const
+  {
+    const double &x = p[0];
+    Assert (x >= x_min && x <= x_max,
+            ExcCSplineRange(x,x_min,x_max));
+
+    return gsl_spline_eval (cspline, x, acc);
+  }
+
+
+  template <int dim>
+  Tensor<1,dim>
+  CSpline<dim>::gradient (const Point<dim>   &p,
+                          const unsigned int) const
+  {
+    const double &x = p[0];
+    Assert (x >= x_min && x <= x_max,
+            ExcCSplineRange(x,x_min,x_max));
+
+    const double deriv = gsl_spline_eval_deriv (cspline, x, acc);
+    Tensor<1,dim> res;
+    res[0] = deriv;
+    return res;
+  }
+
+
+  template <int dim>
+  std::size_t
+  CSpline<dim>::memory_consumption () const
+  {
+    // only simple data elements, so
+    // use sizeof operator
+    return sizeof (*this);
+  }
+
+
+// explicit instantiations
+  template class CSpline<1>;
+
+}
+
+DEAL_II_NAMESPACE_CLOSE

diff --git a/tests/base/functions_cspline.cc b/tests/base/functions_cspline.cc
new file mode 100644 (file)
index 0000000..5633897
--- /dev/null
+++ b/tests/base/functions_cspline.cc
@@ -0,0 +1,79 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// CSpline
+// take example from https://www.gnu.org/software/gsl/manual/html_node/Interpolation-Example-programs.html
+
+#include "../tests.h"
+#include <cmath>
+#include <deal.II/base/function_cspline.h>
+
+template <int dim>
+void check()
+{
+  const unsigned int n_points = 10;
+  std::vector<double> x(n_points), y(n_points);
+  for (unsigned int i = 0; i < n_points; i++)
+    {
+      x[i] = i + 0.5 * std::sin (i);
+      y[i] = i + std::cos (i * i);
+    }
+
+
+  std::vector<double> y_native;
+  // native version
+  {
+    gsl_interp_accel *acc = gsl_interp_accel_alloc ();
+    gsl_spline *spline = gsl_spline_alloc (gsl_interp_cspline, n_points);
+
+    gsl_spline_init (spline, &x[0], &y[0], n_points);
+
+    for (double xi = x[0]; xi <= x.back(); xi += 0.01)
+      {
+        const double yi = gsl_spline_eval (spline, xi, acc);
+        //deallog << xi << " " << yi << std::endl;
+        y_native.push_back(yi);
+      }
+    gsl_spline_free (spline);
+    gsl_interp_accel_free (acc);
+  }
+
+  std::vector<double> y_dealii;
+  {
+    Functions::CSpline<dim> cspline(x,y);
+    for (double xi = x[0]; xi <= x.back(); xi += 0.01)
+      {
+        const double yi = cspline.value(Point<dim>(xi));
+        //deallog << xi << " " << yi << std::endl;
+        y_dealii.push_back(yi);
+      }
+  }
+
+  AssertThrow(std::equal(y_native.begin(), y_native.end(), y_dealii.begin()),
+              ExcMessage("deal.II implementation of CSpline does not match native GSL."));
+
+  deallog << "Ok"<< std::endl;
+}
+
+int main()
+{
+  std::string logname = "output";
+  std::ofstream logfile(logname.c_str());
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  check<1>();
+}

diff --git a/tests/base/functions_cspline.output b/tests/base/functions_cspline.output
new file mode 100644 (file)
index 0000000..6ea1f89
--- /dev/null
+++ b/tests/base/functions_cspline.output
@@ -0,0 +1,2 @@
+
+DEAL::Ok