function_cspline.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2017 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#include <deal.II/base/function_cspline.h>
#include <deal.II/base/point.h>

#ifdef DEAL_II_WITH_GSL
#  include <algorithm>
#  include <cmath>

DEAL_II_NAMESPACE_OPEN
namespace Functions
{
  template <int dim>
  CSpline<dim>::CSpline(const std::vector<double> &x_,
                        const std::vector<double> &y_)
    : interpolation_points(x_)
    , interpolation_values(y_)
  {
    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 = interpolation_points.size();
    cspline              = gsl_spline_alloc(gsl_interp_cspline, n);
    // gsl_spline_init returns something but it seems nobody knows what
    gsl_spline_init(cspline,
                    interpolation_points.data(),
                    interpolation_values.data(),
                    n);
  }



  template <int dim>
  CSpline<dim>::~CSpline()
  {
    gsl_interp_accel_free(acc);
    gsl_spline_free(cspline);
    acc     = nullptr;
    cspline = nullptr;
  }



  template <int dim>
  double
  CSpline<dim>::value(const Point<dim> &p, const unsigned int) const
  {
    // GSL functions may modify gsl_interp_accel *acc object (last argument).
    // This can only work in multithreaded applications if we lock the data
    // structures via a mutex.
    Threads::Mutex::ScopedLock lock(acc_mutex);

    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(cspline, x, acc);
  }



  template <int dim>
  Tensor<1, dim>
  CSpline<dim>::gradient(const Point<dim> &p, const unsigned int) const
  {
    // GSL functions may modify gsl_interp_accel *acc object (last argument).
    // This can only work in multithreaded applications if we lock the data
    // structures via a mutex.
    Threads::Mutex::ScopedLock lock(acc_mutex);

    const double &x = p[0];
    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;
    res[0] = deriv;
    return res;
  }



  template <int dim>
  double
  CSpline<dim>::laplacian(const Point<dim> &p, const unsigned int) const
  {
    // GSL functions may modify gsl_interp_accel *acc object (last argument).
    // This can only work in multithreaded applications if we lock the data
    // structures via a mutex.
    Threads::Mutex::ScopedLock lock(acc_mutex);

    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
  {
    SymmetricTensor<2, dim> res;
    res[0][0] = laplacian(p);
    return res;
  }



  template <int dim>
  std::size_t
  CSpline<dim>::memory_consumption() const
  {
    // only simple data elements, so
    // use sizeof operator
    return sizeof(*this) + 2 * sizeof(double) * interpolation_values.size();
  }


  // explicit instantiations
  template class CSpline<1>;

} // namespace Functions

DEAL_II_NAMESPACE_CLOSE

#endif