Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
function_cspline.cc
Go to the documentation of this file.
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
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 {
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,
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(),
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(),
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
162
163#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
const std::vector< double > interpolation_points
virtual double laplacian(const Point< dim > &p, const unsigned int component=0) const override
const std::vector< double > interpolation_values
virtual Tensor< 1, dim > gradient(const Point< dim > &p, const unsigned int component=0) const override
gsl_interp_accel * acc
virtual ~CSpline() override
Definition: point.h:111
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
static ::ExceptionBase & ExcCSplineSizeMismatch(int arg1, int arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1465
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)
Definition: exceptions.h:1575