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 "BicubicSplineInterpolation.h"
11 : #include "MooseError.h"
12 :
13 : int BicubicSplineInterpolation::_file_number = 0;
14 :
15 42 : BicubicSplineInterpolation::BicubicSplineInterpolation() {}
16 :
17 2 : BicubicSplineInterpolation::BicubicSplineInterpolation(const std::vector<Real> & x1,
18 : const std::vector<Real> & x2,
19 : const std::vector<std::vector<Real>> & y,
20 : const std::vector<Real> & yx11,
21 : const std::vector<Real> & yx1n,
22 : const std::vector<Real> & yx21,
23 2 : const std::vector<Real> & yx2n)
24 : : SplineInterpolationBase(),
25 2 : _x1(x1),
26 2 : _x2(x2),
27 2 : _y(y),
28 2 : _yx11(yx11),
29 2 : _yx1n(yx1n),
30 2 : _yx21(yx21),
31 6 : _yx2n(yx2n)
32 : {
33 2 : auto n = _x2.size();
34 2 : _row_spline_second_derivs.resize(n);
35 2 : _column_spline_eval.resize(n);
36 :
37 2 : auto m = _x1.size();
38 2 : _column_spline_second_derivs.resize(m);
39 2 : _row_spline_eval.resize(m);
40 :
41 2 : errorCheck();
42 2 : solve();
43 2 : }
44 :
45 : void
46 42 : BicubicSplineInterpolation::setData(const std::vector<Real> & x1,
47 : const std::vector<Real> & x2,
48 : const std::vector<std::vector<Real>> & y,
49 : const std::vector<Real> & yx11,
50 : const std::vector<Real> & yx1n,
51 : const std::vector<Real> & yx21,
52 : const std::vector<Real> & yx2n)
53 : {
54 42 : _x1 = x1;
55 42 : _x2 = x2;
56 42 : _y = y;
57 42 : _yx11 = yx11;
58 42 : _yx1n = yx1n;
59 42 : _yx21 = yx21;
60 42 : _yx2n = yx2n;
61 :
62 42 : auto n = _x2.size();
63 42 : _row_spline_second_derivs.resize(n);
64 42 : _column_spline_eval.resize(n);
65 :
66 42 : auto m = _x1.size();
67 42 : _column_spline_second_derivs.resize(m);
68 42 : _row_spline_eval.resize(m);
69 :
70 42 : errorCheck();
71 42 : solve();
72 42 : }
73 :
74 : void
75 44 : BicubicSplineInterpolation::errorCheck()
76 : {
77 44 : auto m = _x1.size(), n = _x2.size();
78 :
79 44 : if (_y.size() != m)
80 0 : mooseError("y row dimension does not match the size of x1.");
81 : else
82 180 : for (decltype(m) i = 0; i < _y.size(); ++i)
83 136 : if (_y[i].size() != n)
84 0 : mooseError("y column dimension does not match the size of x2.");
85 :
86 44 : if (_yx11.empty())
87 0 : _yx11.resize(n, _deriv_bound);
88 44 : else if (_yx11.size() != n)
89 0 : mooseError("The length of the vectors holding the first derivatives of y with respect to x1 "
90 : "must match the length of x2.");
91 :
92 44 : if (_yx1n.empty())
93 0 : _yx1n.resize(n, _deriv_bound);
94 44 : else if (_yx1n.size() != n)
95 0 : mooseError("The length of the vectors holding the first derivatives of y with respect to x1 "
96 : "must match the length of x2.");
97 :
98 44 : if (_yx21.empty())
99 0 : _yx21.resize(m, _deriv_bound);
100 44 : else if (_yx21.size() != m)
101 0 : mooseError("The length of the vectors holding the first derivatives of y with respect to x2 "
102 : "must match the length of x1.");
103 :
104 44 : if (_yx2n.empty())
105 0 : _yx2n.resize(m, _deriv_bound);
106 44 : else if (_yx2n.size() != m)
107 0 : mooseError("The length of the vectors holding the first derivatives of y with respect to x2 "
108 : "must match the length of x1.");
109 44 : }
110 :
111 : void
112 44 : BicubicSplineInterpolation::constructRowSplineSecondDerivativeTable()
113 : {
114 44 : auto m = _x1.size();
115 44 : _y2_rows.resize(m);
116 :
117 44 : if (_yx21.empty())
118 0 : for (decltype(m) j = 0; j < m; ++j)
119 0 : spline(_x2, _y[j], _y2_rows[j]);
120 :
121 : else
122 180 : for (decltype(m) j = 0; j < m; ++j)
123 136 : spline(_x2, _y[j], _y2_rows[j], _yx21[j], _yx2n[j]);
124 44 : }
125 :
126 : void
127 44 : BicubicSplineInterpolation::constructColumnSplineSecondDerivativeTable()
128 : {
129 44 : auto m = _x1.size();
130 44 : auto n = _x2.size();
131 44 : _y2_columns.resize(n);
132 44 : _y_trans.resize(n);
133 :
134 222 : for (decltype(n) j = 0; j < n; ++j)
135 178 : _y_trans[j].resize(m);
136 :
137 : // transpose the _y values so the columns can be easily iterated over
138 180 : for (decltype(n) i = 0; i < _y.size(); ++i)
139 690 : for (decltype(n) j = 0; j < _y[0].size(); ++j)
140 554 : _y_trans[j][i] = _y[i][j];
141 :
142 44 : if (_yx11.empty())
143 0 : for (decltype(n) j = 0; j < n; ++j)
144 0 : spline(_x1, _y_trans[j], _y2_columns[j]);
145 :
146 : else
147 222 : for (decltype(n) j = 0; j < n; ++j)
148 178 : spline(_x1, _y_trans[j], _y2_columns[j], _yx11[j], _yx1n[j]);
149 44 : }
150 :
151 : void
152 44 : BicubicSplineInterpolation::solve()
153 : {
154 44 : constructRowSplineSecondDerivativeTable();
155 44 : constructColumnSplineSecondDerivativeTable();
156 44 : }
157 :
158 : Real
159 7452 : BicubicSplineInterpolation::sample(Real x1,
160 : Real x2,
161 : Real yx11 /* = _deriv_bound*/,
162 : Real yx1n /* = _deriv_bound*/)
163 : {
164 7452 : constructColumnSpline(x2, _row_spline_eval, _column_spline_second_derivs, yx11, yx1n);
165 :
166 : // Evaluate newly constructed column spline
167 7452 : return SplineInterpolationBase::sample(_x1, _row_spline_eval, _column_spline_second_derivs, x1);
168 : }
169 :
170 : Real
171 19004 : BicubicSplineInterpolation::sampleDerivative(Real x1,
172 : Real x2,
173 : unsigned int deriv_var,
174 : Real yp1 /* = _deriv_bound*/,
175 : Real ypn /* = _deriv_bound*/)
176 : {
177 : // Take derivative along x1 axis
178 19004 : if (deriv_var == 1)
179 : {
180 9502 : constructColumnSpline(x2, _row_spline_eval, _column_spline_second_derivs, yp1, ypn);
181 :
182 : // Evaluate derivative wrt x1 of newly constructed column spline
183 19004 : return SplineInterpolationBase::sampleDerivative(
184 9502 : _x1, _row_spline_eval, _column_spline_second_derivs, x1);
185 : }
186 :
187 : // Take derivative along x2 axis
188 9502 : else if (deriv_var == 2)
189 : {
190 9502 : constructRowSpline(x1, _column_spline_eval, _row_spline_second_derivs, yp1, ypn);
191 :
192 : // Evaluate derivative wrt x2 of newly constructed row spline
193 19004 : return SplineInterpolationBase::sampleDerivative(
194 9502 : _x2, _column_spline_eval, _row_spline_second_derivs, x2);
195 : }
196 :
197 : else
198 0 : mooseError("deriv_var must be either 1 or 2 in BicubicSplineInterpolation");
199 : }
200 :
201 : Real
202 4 : BicubicSplineInterpolation::sample2ndDerivative(Real x1,
203 : Real x2,
204 : unsigned int deriv_var,
205 : Real yp1 /* = _deriv_bound*/,
206 : Real ypn /* = _deriv_bound*/)
207 : {
208 : // Take second derivative along x1 axis
209 4 : if (deriv_var == 1)
210 : {
211 2 : constructColumnSpline(x2, _row_spline_eval, _column_spline_second_derivs, yp1, ypn);
212 :
213 : // Evaluate second derivative wrt x1 of newly constructed column spline
214 4 : return SplineInterpolationBase::sample2ndDerivative(
215 2 : _x1, _row_spline_eval, _column_spline_second_derivs, x1);
216 : }
217 :
218 : // Take second derivative along x2 axis
219 2 : else if (deriv_var == 2)
220 : {
221 2 : constructRowSpline(x1, _column_spline_eval, _row_spline_second_derivs, yp1, ypn);
222 :
223 : // Evaluate second derivative wrt x2 of newly constructed row spline
224 4 : return SplineInterpolationBase::sample2ndDerivative(
225 2 : _x2, _column_spline_eval, _row_spline_second_derivs, x2);
226 : }
227 :
228 : else
229 0 : mooseError("deriv_var must be either 1 or 2 in BicubicSplineInterpolation");
230 : }
231 :
232 : void
233 2 : BicubicSplineInterpolation::sampleValueAndDerivatives(Real x1,
234 : Real x2,
235 : Real & y,
236 : Real & dy1,
237 : Real & dy2,
238 : Real yx11 /* = _deriv_bound*/,
239 : Real yx1n /* = _deriv_bound*/,
240 : Real yx21 /* = _deriv_bound*/,
241 : Real yx2n /* = _deriv_bound*/)
242 : {
243 2 : constructColumnSpline(x2, _row_spline_eval, _column_spline_second_derivs, yx11, yx1n);
244 2 : y = SplineInterpolationBase::sample(_x1, _row_spline_eval, _column_spline_second_derivs, x1);
245 4 : dy1 = SplineInterpolationBase::sampleDerivative(
246 2 : _x1, _row_spline_eval, _column_spline_second_derivs, x1);
247 :
248 2 : constructRowSpline(x1, _column_spline_eval, _row_spline_second_derivs, yx21, yx2n);
249 4 : dy2 = SplineInterpolationBase::sampleDerivative(
250 2 : _x2, _column_spline_eval, _row_spline_second_derivs, x2);
251 2 : }
252 :
253 : void
254 9506 : BicubicSplineInterpolation::constructRowSpline(Real x1,
255 : std::vector<Real> & column_spline_eval,
256 : std::vector<Real> & row_spline_second_derivs,
257 : Real yx11 /*= _deriv_bound*/,
258 : Real yx1n /*= _deriv_bound*/)
259 : {
260 9506 : auto n = _x2.size();
261 :
262 : // Find the indices that bound the point x1
263 : unsigned int klo, khi;
264 9506 : findInterval(_x1, x1, klo, khi);
265 :
266 : // Evaluate n column-splines to get y-values for row spline construction using
267 : // the indices above to avoid computing them for each j
268 47536 : for (decltype(n) j = 0; j < n; ++j)
269 38030 : _column_spline_eval[j] =
270 38030 : SplineInterpolationBase::sample(_x1, _y_trans[j], _y2_columns[j], x1, klo, khi);
271 :
272 : // Construct single row spline; get back the second derivatives wrt x2 coord
273 : // on the x2 grid points
274 9506 : spline(_x2, column_spline_eval, row_spline_second_derivs, yx11, yx1n);
275 9506 : }
276 :
277 : void
278 16958 : BicubicSplineInterpolation::constructColumnSpline(Real x2,
279 : std::vector<Real> & row_spline_eval,
280 : std::vector<Real> & column_spline_second_derivs,
281 : Real yx21 /*= _deriv_bound*/,
282 : Real yx2n /*= _deriv_bound*/)
283 : {
284 16958 : auto m = _x1.size();
285 :
286 : // Find the indices that bound the point x2
287 : unsigned int klo, khi;
288 16958 : findInterval(_x2, x2, klo, khi);
289 :
290 : // Evaluate m row-splines to get y-values for column spline construction using
291 : // the indices above to avoid computing them for each j
292 67848 : for (decltype(m) j = 0; j < m; ++j)
293 50890 : _row_spline_eval[j] = SplineInterpolationBase::sample(_x2, _y[j], _y2_rows[j], x2, klo, khi);
294 :
295 : // Construct single column spline; get back the second derivatives wrt x1 coord
296 : // on the x1 grid points
297 16958 : spline(_x1, row_spline_eval, column_spline_second_derivs, yx21, yx2n);
298 16958 : }
|
