From e0b30acd6c24780901f48778f5c9922092a2b44c Mon Sep 17 00:00:00 2001 From: Denis Davydov <davydden@gmail.com> Date: Thu, 28 Apr 2016 17:10:16 +0200 Subject: [PATCH] minor cleanup of Functions::CSpline class --- include/deal.II/base/function_cspline.h | 42 ++++++++++++++----------- source/base/function_cspline.cc | 36 ++++++++++++--------- 2 files changed, 44 insertions(+), 34 deletions(-) diff --git a/include/deal.II/base/function_cspline.h b/include/deal.II/base/function_cspline.h index 361c7f7464..cf5577014d 100644 --- 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 af8df06dd3..f06dc8b562 100644 --- 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(); } -- 2.39.5