From 7eb2ef301bfc1cb3e3dc7ef30580f621ccc62f32 Mon Sep 17 00:00:00 2001
From: Denis Davydov <davydden@gmail.com>
Date: Tue, 31 Jan 2017 08:43:39 +0100
Subject: [PATCH] CSpline: add hessian and laplacian

---
 include/deal.II/base/function_cspline.h |  6 +++++
 source/base/function_cspline.cc         | 31 +++++++++++++++++++++++++
 2 files changed, 37 insertions(+)

diff --git a/include/deal.II/base/function_cspline.h b/include/deal.II/base/function_cspline.h
index 0d324460dc..1653ba0fd2 100644
--- a/include/deal.II/base/function_cspline.h
+++ b/include/deal.II/base/function_cspline.h
@@ -87,6 +87,12 @@ namespace Functions
     virtual Tensor<1,dim> gradient (const Point<dim>   &p,
                                     const unsigned int  component = 0) const;
 
+    virtual SymmetricTensor<2,dim> hessian (const Point<dim>   &p,
+                                            const unsigned int  component = 0) const;
+
+    virtual double laplacian(const Point< dim > &p,
+                             const unsigned int component = 0) const;
+
     std::size_t memory_consumption () const;
 
   private:
diff --git a/source/base/function_cspline.cc b/source/base/function_cspline.cc
index f06dc8b562..abb8d87afc 100644
--- a/source/base/function_cspline.cc
+++ b/source/base/function_cspline.cc
@@ -49,6 +49,8 @@ namespace Functions
     gsl_spline_init (cspline, &interpolation_points[0], &interpolation_values[0], n);
   }
 
+
+
   template <int dim>
   CSpline<dim>::~CSpline()
   {
@@ -59,6 +61,7 @@ namespace Functions
   }
 
 
+
   template <int dim>
   double
   CSpline<dim>::value (const Point<dim>   &p,
@@ -72,6 +75,7 @@ namespace Functions
   }
 
 
+
   template <int dim>
   Tensor<1,dim>
   CSpline<dim>::gradient (const Point<dim>   &p,
@@ -88,6 +92,33 @@ namespace Functions
   }
 
 
+
+  template <int dim>
+  double
+  CSpline<dim>::laplacian (const Point<dim>   &p,
+                           const unsigned int) const
+  {
+    const double &x = p[0];
+    Assert (x >= interpolation_points.front() && x <= interpolation_points.back(),
+            ExcCSplineRange(x,interpolation_points.front(),interpolation_points.back()));
+
+    return gsl_spline_eval_deriv2 (cspline, x, acc);
+  }
+
+
+
+  template <int dim>
+  SymmetricTensor<2,dim>
+  CSpline<dim>::hessian (const Point<dim>   &p,
+                         const unsigned int) const
+  {
+    Tensor<2,dim> res;
+    res[0][0] = laplacian(p);
+    return res;
+  }
+
+
+
   template <int dim>
   std::size_t
   CSpline<dim>::memory_consumption () const
-- 
2.39.5