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 326 : SplineInterpolationBase::SplineInterpolationBase() {}
17 :
18 : void
19 23976 : 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 23976 : auto n = x.size();
26 23976 : if (n < 2)
27 0 : mooseError("You must have at least two knots to create a spline.");
28 :
29 23976 : std::vector<Real> u(n, 0.);
30 23976 : y2.assign(n, 0.);
31 :
32 23976 : if (yp1 >= 1e30)
33 287 : y2[0] = u[0] = 0.;
34 : else
35 : {
36 23689 : y2[0] = -0.5;
37 23689 : 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 57784 : for (decltype(n) i = 1; i < n - 1; i++)
41 : {
42 33808 : Real sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
43 33808 : Real p = sig * y2[i - 1] + 2.0;
44 33808 : y2[i] = (sig - 1.0) / p;
45 33808 : u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
46 33808 : u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
47 : }
48 :
49 : Real qn, un;
50 23976 : if (ypn >= 1e30)
51 287 : qn = un = 0.;
52 : else
53 : {
54 23689 : qn = 0.5;
55 23689 : 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 23976 : y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.);
59 : // back substitution
60 81760 : for (auto k = n - 1; k >= 1; k--)
61 57784 : y2[k - 1] = y2[k - 1] * y2[k] + u[k - 1];
62 23976 : }
63 :
64 : void
65 52679 : SplineInterpolationBase::findInterval(const std::vector<Real> & x,
66 : Real x_int,
67 : unsigned int & klo,
68 : unsigned int & khi) const
69 : {
70 52679 : klo = 0;
71 : mooseAssert(x.size() >= 2, "You must have at least two knots to create a spline.");
72 52679 : khi = x.size() - 1;
73 124855 : while (khi - klo > 1)
74 : {
75 72176 : unsigned int k = (khi + klo) >> 1;
76 72176 : if (x[k] > x_int)
77 32612 : khi = k;
78 : else
79 39564 : klo = k;
80 : }
81 52679 : }
82 :
83 : template <typename T>
84 : void
85 107907 : 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 107907 : h = x[khi] - x[klo];
94 107907 : if (h == 0)
95 0 : mooseError("The values of x must be distinct");
96 107907 : a = (x[khi] - x_int) / h;
97 107907 : b = (x_int - x[klo]) / h;
98 107907 : }
99 :
100 : Real
101 9762 : 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 9762 : findInterval(x, x_int, klo, khi);
108 :
109 19524 : 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 18484 : 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 18484 : findInterval(x, x_int, klo, khi);
132 :
133 : Real h, a, b;
134 18484 : computeCoeffs(x, klo, khi, x_int, h, a, b);
135 :
136 18484 : return (y[khi] - y[klo]) / h -
137 18484 : (((3.0 * a * a - 1.0) * y2[klo] + (3.0 * b * b - 1.0) * -y2[khi]) * h / 6.0);
138 : }
139 :
140 : Real
141 1026 : 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 1026 : findInterval(x, x_int, klo, khi);
148 :
149 : Real h, a, b;
150 1026 : computeCoeffs(x, klo, khi, x_int, h, a, b);
151 :
152 1026 : return a * y2[klo] + b * y2[khi];
153 : }
154 :
155 : template <typename T>
156 : T
157 88397 : 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 88397 : computeCoeffs(x, klo, khi, x_int, h, a, b);
167 :
168 88397 : return a * y[klo] + b * y[khi] +
169 88397 : ((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;
|
