Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "SplineInterpolationBase.h"
11 : #include "MooseError.h"
12 : #include <limits>
13 :
14 : const Real SplineInterpolationBase::_deriv_bound = std::numeric_limits<Real>::max();
15 :
16 351 : SplineInterpolationBase::SplineInterpolationBase() {}
17 :
18 : void
19 27069 : SplineInterpolationBase::spline(const std::vector<Real> & x,
20 : const std::vector<Real> & y,
21 : std::vector<Real> & y2,
22 : Real yp1 /* = _deriv_bound*/,
23 : Real ypn /* = _deriv_bound*/)
24 : {
25 27069 : auto n = x.size();
26 27069 : if (n < 2)
27 0 : mooseError("You must have at least two knots to create a spline.");
28 :
29 27069 : std::vector<Real> u(n, 0.);
30 27069 : y2.assign(n, 0.);
31 :
32 27069 : if (yp1 >= 1e30)
33 309 : y2[0] = u[0] = 0.;
34 : else
35 : {
36 26760 : y2[0] = -0.5;
37 26760 : u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1);
38 : }
39 : // decomposition of tri-diagonal algorithm (y2 and u are used for temporary storage)
40 65175 : for (decltype(n) i = 1; i < n - 1; i++)
41 : {
42 38106 : Real sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
43 38106 : Real p = sig * y2[i - 1] + 2.0;
44 38106 : y2[i] = (sig - 1.0) / p;
45 38106 : u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
46 38106 : u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
47 : }
48 :
49 : Real qn, un;
50 27069 : if (ypn >= 1e30)
51 309 : qn = un = 0.;
52 : else
53 : {
54 26760 : qn = 0.5;
55 26760 : un = (3.0 / (x[n - 1] - x[n - 2])) * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
56 : }
57 :
58 27069 : y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.);
59 : // back substitution
60 92244 : for (auto k = n - 1; k >= 1; k--)
61 65175 : y2[k - 1] = y2[k - 1] * y2[k] + u[k - 1];
62 27069 : }
63 :
64 : void
65 59509 : SplineInterpolationBase::findInterval(const std::vector<Real> & x,
66 : Real x_int,
67 : unsigned int & klo,
68 : unsigned int & khi) const
69 : {
70 59509 : klo = 0;
71 : mooseAssert(x.size() >= 2, "You must have at least two knots to create a spline.");
72 59509 : khi = x.size() - 1;
73 140981 : while (khi - klo > 1)
74 : {
75 81472 : unsigned int k = (khi + klo) >> 1;
76 81472 : if (x[k] > x_int)
77 36814 : khi = k;
78 : else
79 44658 : klo = k;
80 : }
81 59509 : }
82 :
83 : template <typename T>
84 : void
85 121937 : SplineInterpolationBase::computeCoeffs(const std::vector<Real> & x,
86 : unsigned int klo,
87 : unsigned int khi,
88 : const T & x_int,
89 : Real & h,
90 : T & a,
91 : T & b) const
92 : {
93 121937 : h = x[khi] - x[klo];
94 121937 : if (h == 0)
95 0 : mooseError("The values of x must be distinct");
96 121937 : a = (x[khi] - x_int) / h;
97 121937 : b = (x_int - x[klo]) / h;
98 121937 : }
99 :
100 : Real
101 10974 : SplineInterpolationBase::sample(const std::vector<Real> & x,
102 : const std::vector<Real> & y,
103 : const std::vector<Real> & y2,
104 : Real x_int) const
105 : {
106 : unsigned int klo, khi;
107 10974 : findInterval(x, x_int, klo, khi);
108 :
109 21948 : return sample(x, y, y2, x_int, klo, khi);
110 : }
111 :
112 : ADReal
113 0 : SplineInterpolationBase::sample(const std::vector<Real> & x,
114 : const std::vector<Real> & y,
115 : const std::vector<Real> & y2,
116 : const ADReal & x_int) const
117 : {
118 : unsigned int klo, khi;
119 0 : findInterval(x, MetaPhysicL::raw_value(x_int), klo, khi);
120 :
121 0 : return sample(x, y, y2, x_int, klo, khi);
122 : }
123 :
124 : Real
125 20924 : SplineInterpolationBase::sampleDerivative(const std::vector<Real> & x,
126 : const std::vector<Real> & y,
127 : const std::vector<Real> & y2,
128 : Real x_int) const
129 : {
130 : unsigned int klo, khi;
131 20924 : findInterval(x, x_int, klo, khi);
132 :
133 : Real h, a, b;
134 20924 : computeCoeffs(x, klo, khi, x_int, h, a, b);
135 :
136 20924 : return (y[khi] - y[klo]) / h -
137 20924 : (((3.0 * a * a - 1.0) * y2[klo] + (3.0 * b * b - 1.0) * -y2[khi]) * h / 6.0);
138 : }
139 :
140 : Real
141 1154 : SplineInterpolationBase::sample2ndDerivative(const std::vector<Real> & x,
142 : const std::vector<Real> & /*y*/,
143 : const std::vector<Real> & y2,
144 : Real x_int) const
145 : {
146 : unsigned int klo, khi;
147 1154 : findInterval(x, x_int, klo, khi);
148 :
149 : Real h, a, b;
150 1154 : computeCoeffs(x, klo, khi, x_int, h, a, b);
151 :
152 1154 : return a * y2[klo] + b * y2[khi];
153 : }
154 :
155 : template <typename T>
156 : T
157 99859 : SplineInterpolationBase::sample(const std::vector<Real> & x,
158 : const std::vector<Real> & y,
159 : const std::vector<Real> & y2,
160 : const T & x_int,
161 : unsigned int klo,
162 : unsigned int khi) const
163 : {
164 : Real h;
165 0 : T a, b;
166 99859 : computeCoeffs(x, klo, khi, x_int, h, a, b);
167 :
168 99859 : return a * y[klo] + b * y[khi] +
169 99859 : ((a * a * a - a) * y2[klo] + (b * b * b - b) * y2[khi]) * (h * h) / 6.0;
170 0 : }
171 :
172 : template Real SplineInterpolationBase::sample<Real>(const std::vector<Real> & x,
173 : const std::vector<Real> & y,
174 : const std::vector<Real> & y2,
175 : const Real & x_int,
176 : unsigned int klo,
177 : unsigned int khi) const;
178 : template ADReal SplineInterpolationBase::sample<ADReal>(const std::vector<Real> & x,
179 : const std::vector<Real> & y,
180 : const std::vector<Real> & y2,
181 : const ADReal & x_int,
182 : unsigned int klo,
183 : unsigned int khi) const;
|
