minor cleanup of Functions::CSpline class

diff --git a/include/deal.II/base/function_cspline.h b/include/deal.II/base/function_cspline.h
index 361c7f7464202b729d5ea0188d0003fa7f32a002..cf5577014d1799d549d16d866c37a92dc426aec1 100644 (file)
--- a/include/deal.II/base/function_cspline.h
+++ b/include/deal.II/base/function_cspline.h
@@ -27,19 +27,32 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace Functions
 {
+  DeclException1 (ExcCSplineEmpty,
+                  int,
+                  << "Interpolation points vector size can not be <"<<arg1<<">."
+                 );
+
+  DeclException2 (ExcCSplineSizeMismatch,
+                  int,
+                  int,
+                  << "The size of interpolation points <"<<arg1<<"> is different from the size of interpolation values <" << arg2 <<">."
+                 );
+
+
   DeclException3 (ExcCSplineOrder,
                   int,
                   double,
                   double,
-                  << "CSpline: x[" << arg1 << "] = "<< arg2 <<" >= x["<<(arg1+1)<<" = "<<arg3
+                  << "The input interpolation points are not strictly ordered : " << std::endl << "x[" << arg1 << "] = "<< arg2 <<" >= x["<<(arg1+1)<<" = "<<arg3 <<"."
                  );
 
   DeclException3 (ExcCSplineRange,
                   double,
                   double,
                   double,
-                  << "CSpline: " << arg1 << " is not in ["<< arg2<<";"<<arg3<<"]."
+                  << "Spline function can not be evaluated outside of the interpolation range: "<< std::endl << 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
@@ -56,11 +69,12 @@ namespace Functions
   {
   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 .
+     * Constructor which should be provided with a set of points at which
+     * interpolation is to be done @p interpolation_points and a set of function
+     * values @p interpolation_values .
      */
-    CSpline(const std::vector<double> &x,
-            const std::vector<double> &f);
+    CSpline(const std::vector<double> &interpolation_points,
+            const std::vector<double> &interpolation_values);
 
     /**
      * Virtual destructor.
@@ -79,22 +93,12 @@ namespace Functions
     /**
      * Points at which interpolation is provided
      */
-    const std::vector<double> x;
+    const std::vector<double> interpolation_points;
 
     /**
      * 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;
+    const std::vector<double> interpolation_values;
 
     /**
      * GSL accelerator for spline interpolation
@@ -102,7 +106,7 @@ namespace Functions
     gsl_interp_accel *acc;
 
     /**
-     * GSL cubic spline object
+     * GSL cubic spline interpolator
      */
     gsl_spline *cspline;
   };

diff --git a/source/base/function_cspline.cc b/source/base/function_cspline.cc
index af8df06dd33539990d0387787c810591797ded90..f06dc8b5625b721405a9ad4244f3ecb5ee974eb1 100644 (file)
--- a/source/base/function_cspline.cc
+++ b/source/base/function_cspline.cc
@@ -28,21 +28,25 @@ namespace Functions
   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.)
+    interpolation_points(x_),
+    interpolation_values(y_)
   {
-    // 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]));
+    Assert (interpolation_points.size() > 0,
+            ExcCSplineEmpty(interpolation_points.size()));
+
+    Assert (interpolation_points.size() == interpolation_values.size(),
+            ExcCSplineSizeMismatch(interpolation_points.size(),interpolation_values.size()));
+
+    // check that input vector @p interpolation_points is provided in ascending order:
+    for (unsigned int i = 0; i < interpolation_points.size()-1; i++)
+      AssertThrow(interpolation_points[i]<interpolation_points[i+1],
+                  ExcCSplineOrder(i,interpolation_points[i],interpolation_points[i+1]));
 
     acc = gsl_interp_accel_alloc ();
-    const unsigned int n = x.size();
+    const unsigned int n = interpolation_points.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);
+    gsl_spline_init (cspline, &interpolation_points[0], &interpolation_values[0], n);
   }
 
   template <int dim>
@@ -50,6 +54,8 @@ namespace Functions
   {
     gsl_interp_accel_free (acc);
     gsl_spline_free (cspline);
+    acc = NULL;
+    cspline = NULL;
   }
 
 
@@ -59,8 +65,8 @@ namespace Functions
                        const unsigned int) const
   {
     const double &x = p[0];
-    Assert (x >= x_min && x <= x_max,
-            ExcCSplineRange(x,x_min,x_max));
+    Assert (x >= interpolation_points.front() && x <= interpolation_points.back(),
+            ExcCSplineRange(x,interpolation_points.front(),interpolation_points.back()));
 
     return gsl_spline_eval (cspline, x, acc);
   }
@@ -72,8 +78,8 @@ namespace Functions
                           const unsigned int) const
   {
     const double &x = p[0];
-    Assert (x >= x_min && x <= x_max,
-            ExcCSplineRange(x,x_min,x_max));
+    Assert (x >= interpolation_points.front() && x <= interpolation_points.back(),
+            ExcCSplineRange(x,interpolation_points.front(),interpolation_points.back()));
 
     const double deriv = gsl_spline_eval_deriv (cspline, x, acc);
     Tensor<1,dim> res;
@@ -88,7 +94,7 @@ namespace Functions
   {
     // only simple data elements, so
     // use sizeof operator
-    return sizeof (*this);
+    return sizeof (*this) + 2*sizeof(double)*interpolation_values.size();
   }