function_cspline.h

// ---------------------------------------------------------------------
//
// 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.
//
// ---------------------------------------------------------------------

#ifndef dealii_function_cspline_h
#define dealii_function_cspline_h

#include <deal.II/base/config.h>

#ifdef DEAL_II_WITH_GSL
#  include <deal.II/base/function.h>
#  include <deal.II/base/point.h>
#  include <deal.II/base/thread_management.h>

#  include <gsl/gsl_spline.h>

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,
                 << "The input interpolation points are not strictly ordered : "
                 << std::endl
                 << "x[" << arg1 << "] = " << arg2 << " >= x[" << (arg1 + 1)
                 << "] = " << arg3 << ".");

  DeclException3(
    ExcCSplineRange,
    double,
    double,
    double,
    << "Spline function can not be evaluated outside of the interpolation range: "
    << std::endl
    << arg1 << " is not in [" << arg2 << ";" << arg3 << "].");

  template <int dim>
  class CSpline : public Function<dim>
  {
  public:
    CSpline(const std::vector<double> &interpolation_points,
            const std::vector<double> &interpolation_values);

    virtual ~CSpline() override;

    virtual double
    value(const Point<dim> & point,
          const unsigned int component = 0) const override;

    virtual Tensor<1, dim>
    gradient(const Point<dim> & p,
             const unsigned int component = 0) const override;

    virtual SymmetricTensor<2, dim>
    hessian(const Point<dim> & p,
            const unsigned int component = 0) const override;

    virtual double
    laplacian(const Point<dim> & p,
              const unsigned int component = 0) const override;

    std::size_t
    memory_consumption() const;

  private:
    const std::vector<double> interpolation_points;

    const std::vector<double> interpolation_values;

    gsl_interp_accel *acc;

    gsl_spline *cspline;

    mutable Threads::Mutex acc_mutex;
  };
} // namespace Functions

DEAL_II_NAMESPACE_CLOSE

#endif

#endif