function_cspline.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2016 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_function_cspline_h
#define dealii_function_cspline_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exception_macros.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>
 
#include <cstddef>
#include <ostream>
#include <vector>
 
#ifdef DEAL_II_WITH_GSL
#  include <deal.II/base/function.h>
#  include <deal.II/base/mutex.h>
#  include <deal.II/base/point.h>
 
#  include <gsl/gsl_interp.h>
#  include <gsl/gsl_spline.h>
 
#  include <memory>
 
#endif
 
DEAL_II_NAMESPACE_OPEN
 
#ifdef DEAL_II_WITH_GSL
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 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;
 
    virtual std::size_t
    memory_consumption() const override;
 
  private:
    const std::vector<double> interpolation_points;
 
    const std::vector<double> interpolation_values;
 
    std::unique_ptr<gsl_interp_accel, void (*)(gsl_interp_accel *)> acc;
 
    std::unique_ptr<gsl_spline, void (*)(gsl_spline *)> cspline;
 
    mutable Threads::Mutex acc_mutex;
  };
} // namespace Functions
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif