From 489dbd160b4d320513c36aa8e0d3e37a5e14ac1f Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Tue, 31 Jan 2017 08:48:55 +0100
Subject: [PATCH] CSpline: test gradients and laplacians

---
 tests/base/functions_cspline_02.cc            | 81 +++++++++++++++++++
 .../functions_cspline_02.with_gsl=on.output   |  2 +
 2 files changed, 83 insertions(+)
 create mode 100644 tests/base/functions_cspline_02.cc
 create mode 100644 tests/base/functions_cspline_02.with_gsl=on.output

diff --git a/tests/base/functions_cspline_02.cc b/tests/base/functions_cspline_02.cc
new file mode 100644
index 0000000000..9b40148321
--- /dev/null
+++ b/tests/base/functions_cspline_02.cc
@@ -0,0 +1,81 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test CSpline value, gradient and hessian
+// 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);
+    }
+
+
+  // native:
+  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);
+
+  // dealii:
+  Functions::CSpline<dim> cspline(x,y);
+
+  for (double xi = x[0]; xi <= x.back(); xi += 0.01)
+    {
+      const double f   = gsl_spline_eval        (spline, xi, acc);
+      const double df  = gsl_spline_eval_deriv  (spline, xi, acc);
+      const double ddf = gsl_spline_eval_deriv2 (spline, xi, acc);
+
+      const double y         = cspline.value(Point<dim>(xi));
+      const Tensor<1,dim> dy = cspline.gradient(Point<dim>(xi));
+      const double ddy       = cspline.laplacian(Point<dim>(xi));
+
+      AssertThrow( std::fabs(f-y) <= std::fabs(f)*1e-10,
+                   ExcInternalError());
+
+      AssertThrow( std::fabs(df-dy[0]) <= std::fabs(df)*1e-10,
+                   ExcInternalError());
+
+      AssertThrow( std::fabs(ddf-ddy) <= std::fabs(ddf)*1e-10,
+                   ExcInternalError());
+
+    }
+  gsl_spline_free (spline);
+  gsl_interp_accel_free (acc);
+
+  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_02.with_gsl=on.output b/tests/base/functions_cspline_02.with_gsl=on.output
new file mode 100644
index 0000000000..0fd8fc12f0
--- /dev/null
+++ b/tests/base/functions_cspline_02.with_gsl=on.output
@@ -0,0 +1,2 @@
+
+DEAL::OK
-- 
2.39.5