deal.II version GIT relicensing-2865-ga6be64acbe 2025-03-19 21:10: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 <gsl/gsl_spline.h>
20
21# include <algorithm>
22# include <cmath>
23
24
26namespace Functions
27{
28 template <int dim>
29 CSpline<dim>::CSpline(const std::vector<double> &x_,
30 const std::vector<double> &y_)
31 : interpolation_points(x_)
32 , interpolation_values(y_)
33 , acc(gsl_interp_accel_alloc(),
34 [](gsl_interp_accel *p) { gsl_interp_accel_free(p); })
35 , cspline(gsl_spline_alloc(gsl_interp_cspline, interpolation_points.size()),
36 [](gsl_spline *p) { gsl_spline_free(p); })
37 {
38 Assert(interpolation_points.size() > 0,
39 ExcCSplineEmpty(interpolation_points.size()));
40
41 Assert(interpolation_points.size() == interpolation_values.size(),
42 ExcCSplineSizeMismatch(interpolation_points.size(),
43 interpolation_values.size()));
44
45 // check that input vector @p interpolation_points is provided in ascending order:
46 for (unsigned int i = 0; i < interpolation_points.size() - 1; ++i)
47 AssertThrow(interpolation_points[i] < interpolation_points[i + 1],
49 interpolation_points[i],
50 interpolation_points[i + 1]));
51
52 const unsigned int n = interpolation_points.size();
53 // gsl_spline_init returns something but it seems nobody knows what
54 gsl_spline_init(cspline.get(),
55 interpolation_points.data(),
56 interpolation_values.data(),
57 n);
58 }
59
60
61
62 template <int dim>
63 double
64 CSpline<dim>::value(const Point<dim> &p, const unsigned int) const
65 {
66 // GSL functions may modify gsl_interp_accel *acc object (last argument).
67 // This can only work in multithreaded applications if we lock the data
68 // structures via a mutex.
69 std::lock_guard<std::mutex> lock(acc_mutex);
70
71 const double x = p[0];
72 Assert(x >= interpolation_points.front() &&
73 x <= interpolation_points.back(),
75 interpolation_points.front(),
76 interpolation_points.back()));
77
78 return gsl_spline_eval(cspline.get(), x, acc.get());
79 }
80
81
82
83 template <int dim>
85 CSpline<dim>::gradient(const Point<dim> &p, const unsigned int) const
86 {
87 // GSL functions may modify gsl_interp_accel *acc object (last argument).
88 // This can only work in multithreaded applications if we lock the data
89 // structures via a mutex.
90 std::lock_guard<std::mutex> lock(acc_mutex);
91
92 const double x = p[0];
93 Assert(x >= interpolation_points.front() &&
94 x <= interpolation_points.back(),
96 interpolation_points.front(),
97 interpolation_points.back()));
98
99 const double deriv = gsl_spline_eval_deriv(cspline.get(), x, acc.get());
100 Tensor<1, dim> res;
101 res[0] = deriv;
102 return res;
103 }
104
105
106
107 template <int dim>
108 double
109 CSpline<dim>::laplacian(const Point<dim> &p, const unsigned int) const
110 {
111 // GSL functions may modify gsl_interp_accel *acc object (last argument).
112 // This can only work in multithreaded applications if we lock the data
113 // structures via a mutex.
114 std::lock_guard<std::mutex> lock(acc_mutex);
115
116 const double x = p[0];
117 Assert(x >= interpolation_points.front() &&
118 x <= interpolation_points.back(),
120 interpolation_points.front(),
121 interpolation_points.back()));
122
123 return gsl_spline_eval_deriv2(cspline.get(), x, acc.get());
124 }
125
126
127
128 template <int dim>
130 CSpline<dim>::hessian(const Point<dim> &p, const unsigned int) const
131 {
133 res[0][0] = laplacian(p);
134 return res;
135 }
136
137
138
139 template <int dim>
140 std::size_t
142 {
143 // only simple data elements, so
144 // use sizeof operator
145 return sizeof(*this) + 2 * sizeof(double) * interpolation_values.size();
146 }
147
148
149 // explicit instantiations
150 template class CSpline<1>;
151
152} // namespace Functions
153
155
156#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:113
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
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)