Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
24namespace 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());
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(),
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcCSplineSizeMismatch(int arg1, int arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcCSplineRange(double arg1, double arg2, double arg3)
static ::ExceptionBase & ExcCSplineOrder(int arg1, double arg2, double arg3)
static ::ExceptionBase & ExcCSplineEmpty(int arg1)
#define AssertThrow(cond, exc)