Reference documentation for deal.II version GIT c14369f203 2023-10-01 07:40:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
function_cspline.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2021 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 
17 #include <deal.II/base/point.h>
18 
19 #ifdef DEAL_II_WITH_GSL
20 # include <algorithm>
21 # include <cmath>
22 
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  , acc(gsl_interp_accel_alloc(),
32  [](gsl_interp_accel *p) { gsl_interp_accel_free(p); })
33  , cspline(gsl_spline_alloc(gsl_interp_cspline, interpolation_points.size()),
34  [](gsl_spline *p) { gsl_spline_free(p); })
35  {
36  Assert(interpolation_points.size() > 0,
37  ExcCSplineEmpty(interpolation_points.size()));
38 
39  Assert(interpolation_points.size() == interpolation_values.size(),
40  ExcCSplineSizeMismatch(interpolation_points.size(),
41  interpolation_values.size()));
42 
43  // check that input vector @p interpolation_points is provided in ascending order:
44  for (unsigned int i = 0; i < interpolation_points.size() - 1; ++i)
45  AssertThrow(interpolation_points[i] < interpolation_points[i + 1],
47  interpolation_points[i],
48  interpolation_points[i + 1]));
49 
50  const unsigned int n = interpolation_points.size();
51  // gsl_spline_init returns something but it seems nobody knows what
52  gsl_spline_init(cspline.get(),
53  interpolation_points.data(),
54  interpolation_values.data(),
55  n);
56  }
57 
58 
59 
60  template <int dim>
61  double
62  CSpline<dim>::value(const Point<dim> &p, const unsigned int) const
63  {
64  // GSL functions may modify gsl_interp_accel *acc object (last argument).
65  // This can only work in multithreaded applications if we lock the data
66  // structures via a mutex.
67  std::lock_guard<std::mutex> lock(acc_mutex);
68 
69  const double x = p[0];
70  Assert(x >= interpolation_points.front() &&
71  x <= interpolation_points.back(),
73  interpolation_points.front(),
74  interpolation_points.back()));
75 
76  return gsl_spline_eval(cspline.get(), x, acc.get());
77  }
78 
79 
80 
81  template <int dim>
83  CSpline<dim>::gradient(const Point<dim> &p, const unsigned int) const
84  {
85  // GSL functions may modify gsl_interp_accel *acc object (last argument).
86  // This can only work in multithreaded applications if we lock the data
87  // structures via a mutex.
88  std::lock_guard<std::mutex> lock(acc_mutex);
89 
90  const double x = p[0];
91  Assert(x >= interpolation_points.front() &&
92  x <= interpolation_points.back(),
94  interpolation_points.front(),
95  interpolation_points.back()));
96 
97  const double deriv = gsl_spline_eval_deriv(cspline.get(), x, acc.get());
98  Tensor<1, dim> res;
99  res[0] = deriv;
100  return res;
101  }
102 
103 
104 
105  template <int dim>
106  double
107  CSpline<dim>::laplacian(const Point<dim> &p, const unsigned int) const
108  {
109  // GSL functions may modify gsl_interp_accel *acc object (last argument).
110  // This can only work in multithreaded applications if we lock the data
111  // structures via a mutex.
112  std::lock_guard<std::mutex> lock(acc_mutex);
113 
114  const double x = p[0];
115  Assert(x >= interpolation_points.front() &&
116  x <= interpolation_points.back(),
117  ExcCSplineRange(x,
118  interpolation_points.front(),
119  interpolation_points.back()));
120 
121  return gsl_spline_eval_deriv2(cspline.get(), x, acc.get());
122  }
123 
124 
125 
126  template <int dim>
128  CSpline<dim>::hessian(const Point<dim> &p, const unsigned int) const
129  {
131  res[0][0] = laplacian(p);
132  return res;
133  }
134 
135 
136 
137  template <int dim>
138  std::size_t
140  {
141  // only simple data elements, so
142  // use sizeof operator
143  return sizeof(*this) + 2 * sizeof(double) * interpolation_values.size();
144  }
145 
146 
147  // explicit instantiations
148  template class CSpline<1>;
149 
150 } // namespace Functions
151 
153 
154 #endif
virtual SymmetricTensor< 2, dim > hessian(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
CSpline(const std::vector< double > &interpolation_points, const std::vector< double > &interpolation_values)
virtual std::size_t memory_consumption() const override
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
Definition: point.h:112
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcCSplineEmpty(int arg1)
static ::ExceptionBase & ExcCSplineRange(double arg1, double arg2, double arg3)
static ::ExceptionBase & ExcCSplineOrder(int arg1, double arg2, double arg3)
static ::ExceptionBase & ExcCSplineSizeMismatch(int arg1, int arg2)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705