deal.II version GIT relicensing-2076-g6b43d56a25 2024-11-02 12:30:00+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\}}\)
Loading...
Searching...
No Matches
function_cspline.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 2021 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16#include <deal.II/base/point.h>
17
18#ifdef DEAL_II_WITH_GSL
19# include <algorithm>
20# include <cmath>
21
23namespace Functions
24{
25 template <int dim>
26 CSpline<dim>::CSpline(const std::vector<double> &x_,
27 const std::vector<double> &y_)
28 : interpolation_points(x_)
29 , interpolation_values(y_)
31 [](gsl_interp_accel *p) { gsl_interp_accel_free(p); })
32 , cspline(gsl_spline_alloc(gsl_interp_cspline, interpolation_points.size()),
33 [](gsl_spline *p) { gsl_spline_free(p); })
34 {
35 Assert(interpolation_points.size() > 0,
36 ExcCSplineEmpty(interpolation_points.size()));
37
38 Assert(interpolation_points.size() == interpolation_values.size(),
39 ExcCSplineSizeMismatch(interpolation_points.size(),
40 interpolation_values.size()));
41
42 // check that input vector @p interpolation_points is provided in ascending order:
43 for (unsigned int i = 0; i < interpolation_points.size() - 1; ++i)
44 AssertThrow(interpolation_points[i] < interpolation_points[i + 1],
46 interpolation_points[i],
47 interpolation_points[i + 1]));
48
49 const unsigned int n = interpolation_points.size();
50 // gsl_spline_init returns something but it seems nobody knows what
51 gsl_spline_init(cspline.get(),
52 interpolation_points.data(),
53 interpolation_values.data(),
54 n);
55 }
56
57
58
59 template <int dim>
60 double
61 CSpline<dim>::value(const Point<dim> &p, const unsigned int) const
62 {
63 // GSL functions may modify gsl_interp_accel *acc object (last argument).
64 // This can only work in multithreaded applications if we lock the data
65 // structures via a mutex.
66 std::lock_guard<std::mutex> lock(acc_mutex);
67
68 const double x = p[0];
69 Assert(x >= interpolation_points.front() &&
70 x <= interpolation_points.back(),
72 interpolation_points.front(),
73 interpolation_points.back()));
74
75 return gsl_spline_eval(cspline.get(), x, acc.get());
76 }
77
78
79
80 template <int dim>
82 CSpline<dim>::gradient(const Point<dim> &p, const unsigned int) const
83 {
84 // GSL functions may modify gsl_interp_accel *acc object (last argument).
85 // This can only work in multithreaded applications if we lock the data
86 // structures via a mutex.
87 std::lock_guard<std::mutex> lock(acc_mutex);
88
89 const double x = p[0];
90 Assert(x >= interpolation_points.front() &&
91 x <= interpolation_points.back(),
93 interpolation_points.front(),
94 interpolation_points.back()));
95
96 const double deriv = gsl_spline_eval_deriv(cspline.get(), x, acc.get());
98 res[0] = deriv;
99 return res;
100 }
101
102
103
104 template <int dim>
105 double
106 CSpline<dim>::laplacian(const Point<dim> &p, const unsigned int) const
107 {
108 // GSL functions may modify gsl_interp_accel *acc object (last argument).
109 // This can only work in multithreaded applications if we lock the data
110 // structures via a mutex.
111 std::lock_guard<std::mutex> lock(acc_mutex);
112
113 const double x = p[0];
114 Assert(x >= interpolation_points.front() &&
115 x <= interpolation_points.back(),
117 interpolation_points.front(),
118 interpolation_points.back()));
119
120 return gsl_spline_eval_deriv2(cspline.get(), x, acc.get());
121 }
122
123
124
125 template <int dim>
127 CSpline<dim>::hessian(const Point<dim> &p, const unsigned int) const
128 {
130 res[0][0] = laplacian(p);
131 return res;
132 }
133
134
135
136 template <int dim>
137 std::size_t
139 {
140 // only simple data elements, so
141 // use sizeof operator
142 return sizeof(*this) + 2 * sizeof(double) * interpolation_values.size();
143 }
144
145
146 // explicit instantiations
147 template class CSpline<1>;
148
149} // namespace Functions
150
152
153#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
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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)