Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
function_cspline.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/function_cspline.h>
17 #include <deal.II/base/point.h>
18 
19 #ifdef DEAL_II_WITH_GSL
20 # include <algorithm>
21 # include <cmath>
22 
23 DEAL_II_NAMESPACE_OPEN
24 namespace Functions
25 {
26  template <int dim>
27  CSpline<dim>::CSpline(const std::vector<double> &x_,
28  const std::vector<double> &y_)
29  : interpolation_points(x_)
30  , interpolation_values(y_)
31  {
32  Assert(interpolation_points.size() > 0,
34 
37  interpolation_values.size()));
38 
39  // check that input vector @p interpolation_points is provided in ascending order:
40  for (unsigned int i = 0; i < interpolation_points.size() - 1; i++)
44  interpolation_points[i + 1]));
45 
46  acc = gsl_interp_accel_alloc();
47  const unsigned int n = interpolation_points.size();
48  cspline = gsl_spline_alloc(gsl_interp_cspline, n);
49  // gsl_spline_init returns something but it seems nobody knows what
50  gsl_spline_init(cspline,
51  interpolation_points.data(),
52  interpolation_values.data(),
53  n);
54  }
55 
56 
57 
58  template <int dim>
60  {
61  gsl_interp_accel_free(acc);
62  gsl_spline_free(cspline);
63  acc = nullptr;
64  cspline = nullptr;
65  }
66 
67 
68 
69  template <int dim>
70  double
71  CSpline<dim>::value(const Point<dim> &p, const unsigned int) const
72  {
73  // GSL functions may modify gsl_interp_accel *acc object (last argument).
74  // This can only work in multithreaded applications if we lock the data
75  // structures via a mutex.
76  std::lock_guard<std::mutex> lock(acc_mutex);
77 
78  const double x = p[0];
79  Assert(x >= interpolation_points.front() &&
80  x <= interpolation_points.back(),
82  interpolation_points.front(),
83  interpolation_points.back()));
84 
85  return gsl_spline_eval(cspline, x, acc);
86  }
87 
88 
89 
90  template <int dim>
92  CSpline<dim>::gradient(const Point<dim> &p, const unsigned int) const
93  {
94  // GSL functions may modify gsl_interp_accel *acc object (last argument).
95  // This can only work in multithreaded applications if we lock the data
96  // structures via a mutex.
97  std::lock_guard<std::mutex> lock(acc_mutex);
98 
99  const double x = p[0];
100  Assert(x >= interpolation_points.front() &&
101  x <= interpolation_points.back(),
102  ExcCSplineRange(x,
103  interpolation_points.front(),
104  interpolation_points.back()));
105 
106  const double deriv = gsl_spline_eval_deriv(cspline, x, acc);
107  Tensor<1, dim> res;
108  res[0] = deriv;
109  return res;
110  }
111 
112 
113 
114  template <int dim>
115  double
116  CSpline<dim>::laplacian(const Point<dim> &p, const unsigned int) const
117  {
118  // GSL functions may modify gsl_interp_accel *acc object (last argument).
119  // This can only work in multithreaded applications if we lock the data
120  // structures via a mutex.
121  std::lock_guard<std::mutex> lock(acc_mutex);
122 
123  const double x = p[0];
124  Assert(x >= interpolation_points.front() &&
125  x <= interpolation_points.back(),
126  ExcCSplineRange(x,
127  interpolation_points.front(),
128  interpolation_points.back()));
129 
130  return gsl_spline_eval_deriv2(cspline, x, acc);
131  }
132 
133 
134 
135  template <int dim>
137  CSpline<dim>::hessian(const Point<dim> &p, const unsigned int) const
138  {
140  res[0][0] = laplacian(p);
141  return res;
142  }
143 
144 
145 
146  template <int dim>
147  std::size_t
149  {
150  // only simple data elements, so
151  // use sizeof operator
152  return sizeof(*this) + 2 * sizeof(double) * interpolation_values.size();
153  }
154 
155 
156  // explicit instantiations
157  template class CSpline<1>;
158 
159 } // namespace Functions
160 
161 DEAL_II_NAMESPACE_CLOSE
162 
163 #endif
virtual SymmetricTensor< 2, dim > hessian(const Point< dim > &p, const unsigned int component=0) const override
virtual ~CSpline() override
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
virtual double value(const Point< dim > &point, const unsigned int component=0) const override
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
static ::ExceptionBase & ExcCSplineEmpty(int arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
const std::vector< double > interpolation_points
static ::ExceptionBase & ExcCSplineOrder(int arg1, double arg2, double arg3)
gsl_interp_accel * acc
const std::vector< double > interpolation_values
static ::ExceptionBase & ExcCSplineRange(double arg1, double arg2, double arg3)
static ::ExceptionBase & ExcCSplineSizeMismatch(int arg1, int arg2)
CSpline(const std::vector< double > &interpolation_points, const std::vector< double > &interpolation_values)