From a5d9e9ecb0e40d0afe96199d8586a4cf2e66fee9 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Tue, 26 Apr 2016 19:39:55 +0200 Subject: [PATCH] add Functions::CSpline --- include/deal.II/base/function_cspline.h | 113 ++++++++++++++++++++++++ source/base/CMakeLists.txt | 1 + source/base/function_cspline.cc | 99 +++++++++++++++++++++ tests/base/functions_cspline.cc | 79 +++++++++++++++++ tests/base/functions_cspline.output | 2 + 5 files changed, 294 insertions(+) create mode 100644 include/deal.II/base/function_cspline.h create mode 100644 source/base/function_cspline.cc create mode 100644 tests/base/functions_cspline.cc create mode 100644 tests/base/functions_cspline.output diff --git a/include/deal.II/base/function_cspline.h b/include/deal.II/base/function_cspline.h new file mode 100644 index 0000000000..fac6269867 --- /dev/null +++ b/include/deal.II/base/function_cspline.h @@ -0,0 +1,113 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii__function_cspline_h +#define dealii__function_cspline_h + +#include + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace Functions +{ + DeclException3 (ExcCSplineOrder, + int, + double, + double, + << "CSpline: x[" << arg1 << "] = "<< arg2 <<" >= x["<<(arg1+1)<<" = "< + class CSpline : public Function + { + 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 . + */ + CSpline(const std::vector &x, + const std::vector &f); + + /** + * Virtual destructor. + */ + virtual ~CSpline(); + + virtual double value (const Point &point, + const unsigned int component = 0) const; + + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + + std::size_t memory_consumption () const; + + private: + /** + * Points at which interpolation is provided + */ + const std::vector x; + + /** + * Values of the function at interpolation points + */ + const std::vector y; + + /** + * Maximum allowed value of x + */ + const double x_max; + + /** + * Minimum allowed value of x + */ + const double x_min; + + /** + * GSL accelerator for spline interpolation + */ + gsl_interp_accel *acc; + + /** + * GSL cubic spline object + */ + gsl_spline *cspline; + }; +} + +DEAL_II_NAMESPACE_CLOSE + +#endif + diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index d780f23651..f13a47d856 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -25,6 +25,7 @@ SET(_src exceptions.cc flow_function.cc function.cc + function_cspline.cc function_derivative.cc function_lib.cc function_lib_cutoff.cc diff --git a/source/base/function_cspline.cc b/source/base/function_cspline.cc new file mode 100644 index 0000000000..a74e311ee9 --- /dev/null +++ b/source/base/function_cspline.cc @@ -0,0 +1,99 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#include +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN +namespace Functions +{ + + template + CSpline::CSpline(const std::vector &x_, + const std::vector &y_) + : + x(x_), + y(y_), + x_max(x.size()>0? x.back() : 0.), + x_min(x.size()>0? x.front() : 0.) + { + // check that input vector @p x is provided in ascending order: + for (unsigned int i = 0; i < x.size()-1; i++) + AssertThrow(x[i] + CSpline::~CSpline() + { + gsl_interp_accel_free (acc); + gsl_spline_free (cspline); + } + + + template + double + CSpline::value (const Point &p, + const unsigned int) const + { + const double &x = p[0]; + Assert (x >= x_min && x <= x_max, + ExcCSplineRange(x,x_min,x_max)); + + return gsl_spline_eval (cspline, x, acc); + } + + + template + Tensor<1,dim> + CSpline::gradient (const Point &p, + const unsigned int) const + { + const double &x = p[0]; + Assert (x >= x_min && x <= x_max, + ExcCSplineRange(x,x_min,x_max)); + + const double deriv = gsl_spline_eval_deriv (cspline, x, acc); + Tensor<1,dim> res; + res[0] = deriv; + return res; + } + + + template + std::size_t + CSpline::memory_consumption () const + { + // only simple data elements, so + // use sizeof operator + return sizeof (*this); + } + + +// explicit instantiations + template class CSpline<1>; + +} + +DEAL_II_NAMESPACE_CLOSE diff --git a/tests/base/functions_cspline.cc b/tests/base/functions_cspline.cc new file mode 100644 index 0000000000..56338974ed --- /dev/null +++ b/tests/base/functions_cspline.cc @@ -0,0 +1,79 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +// CSpline +// take example from https://www.gnu.org/software/gsl/manual/html_node/Interpolation-Example-programs.html + +#include "../tests.h" +#include +#include + +template +void check() +{ + const unsigned int n_points = 10; + std::vector 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); + } + + + std::vector y_native; + // native version + { + 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); + + for (double xi = x[0]; xi <= x.back(); xi += 0.01) + { + const double yi = gsl_spline_eval (spline, xi, acc); + //deallog << xi << " " << yi << std::endl; + y_native.push_back(yi); + } + gsl_spline_free (spline); + gsl_interp_accel_free (acc); + } + + std::vector y_dealii; + { + Functions::CSpline cspline(x,y); + for (double xi = x[0]; xi <= x.back(); xi += 0.01) + { + const double yi = cspline.value(Point(xi)); + //deallog << xi << " " << yi << std::endl; + y_dealii.push_back(yi); + } + } + + AssertThrow(std::equal(y_native.begin(), y_native.end(), y_dealii.begin()), + ExcMessage("deal.II implementation of CSpline does not match native GSL.")); + + 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.output b/tests/base/functions_cspline.output new file mode 100644 index 0000000000..6ea1f897a5 --- /dev/null +++ b/tests/base/functions_cspline.output @@ -0,0 +1,2 @@ + +DEAL::Ok -- 2.39.5