CSpline: add hessian and laplacian

diff --git a/include/deal.II/base/function_cspline.h b/include/deal.II/base/function_cspline.h
index 0d324460dc3dab75d5de7158fc25607b21a6b73d..1653ba0fd29fcbea770bec038122d7ae2b8c7f60 100644 (file)
--- 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 f06dc8b5625b721405a9ad4244f3ecb5ee974eb1..abb8d87afcebeb5c67db37e2fc4669b80c251e2a 100644 (file)
--- 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