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 39 : BicubicSplineInterpolation::BicubicSplineInterpolation() {}
16 :
17 1 : 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 1 : const std::vector<Real> & yx2n)
24 : : SplineInterpolationBase(),
25 1 : _x1(x1),
26 1 : _x2(x2),
27 1 : _y(y),
28 1 : _yx11(yx11),
29 1 : _yx1n(yx1n),
30 1 : _yx21(yx21),
31 3 : _yx2n(yx2n)
32 : {
33 1 : auto n = _x2.size();
34 1 : _row_spline_second_derivs.resize(n);
35 1 : _column_spline_eval.resize(n);
36 :
37 1 : auto m = _x1.size();
38 1 : _column_spline_second_derivs.resize(m);
39 1 : _row_spline_eval.resize(m);
40 :
41 1 : errorCheck();
42 1 : solve();
43 1 : }
44 :
45 : void
46 39 : 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 39 : _x1 = x1;
55 39 : _x2 = x2;
56 39 : _y = y;
57 39 : _yx11 = yx11;
58 39 : _yx1n = yx1n;
59 39 : _yx21 = yx21;
60 39 : _yx2n = yx2n;
61 :
62 39 : auto n = _x2.size();
63 39 : _row_spline_second_derivs.resize(n);
64 39 : _column_spline_eval.resize(n);
65 :
66 39 : auto m = _x1.size();
67 39 : _column_spline_second_derivs.resize(m);
68 39 : _row_spline_eval.resize(m);
69 :
70 39 : errorCheck();
71 39 : solve();
72 39 : }
73 :
74 : void
75 40 : BicubicSplineInterpolation::errorCheck()
76 : {
77 40 : auto m = _x1.size(), n = _x2.size();
78 :
79 40 : if (_y.size() != m)
80 0 : mooseError("y row dimension does not match the size of x1.");
81 : else
82 162 : for (decltype(m) i = 0; i < _y.size(); ++i)
83 122 : if (_y[i].size() != n)
84 0 : mooseError("y column dimension does not match the size of x2.");
85 :
86 40 : if (_yx11.empty())
87 0 : _yx11.resize(n, _deriv_bound);
88 40 : 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 40 : if (_yx1n.empty())
93 0 : _yx1n.resize(n, _deriv_bound);
94 40 : 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 40 : if (_yx21.empty())
99 0 : _yx21.resize(m, _deriv_bound);
100 40 : 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 40 : if (_yx2n.empty())
105 0 : _yx2n.resize(m, _deriv_bound);
106 40 : 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 40 : }
110 :
111 : void
112 40 : BicubicSplineInterpolation::constructRowSplineSecondDerivativeTable()
113 : {
114 40 : auto m = _x1.size();
115 40 : _y2_rows.resize(m);
116 :
117 40 : 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 162 : for (decltype(m) j = 0; j < m; ++j)
123 122 : spline(_x2, _y[j], _y2_rows[j], _yx21[j], _yx2n[j]);
124 40 : }
125 :
126 : void
127 40 : BicubicSplineInterpolation::constructColumnSplineSecondDerivativeTable()
128 : {
129 40 : auto m = _x1.size();
130 40 : auto n = _x2.size();
131 40 : _y2_columns.resize(n);
132 40 : _y_trans.resize(n);
133 :
134 201 : for (decltype(n) j = 0; j < n; ++j)
135 161 : _y_trans[j].resize(m);
136 :
137 : // transpose the _y values so the columns can be easily iterated over
138 162 : for (decltype(n) i = 0; i < _y.size(); ++i)
139 615 : for (decltype(n) j = 0; j < _y[0].size(); ++j)
140 493 : _y_trans[j][i] = _y[i][j];
141 :
142 40 : 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 201 : for (decltype(n) j = 0; j < n; ++j)
148 161 : spline(_x1, _y_trans[j], _y2_columns[j], _yx11[j], _yx1n[j]);
149 40 : }
150 :
151 : void
152 40 : BicubicSplineInterpolation::solve()
153 : {
154 40 : constructRowSplineSecondDerivativeTable();
155 40 : constructColumnSplineSecondDerivativeTable();
156 40 : }
157 :
158 : Real
159 6601 : BicubicSplineInterpolation::sample(Real x1,
160 : Real x2,
161 : Real yx11 /* = _deriv_bound*/,
162 : Real yx1n /* = _deriv_bound*/)
163 : {
164 6601 : constructColumnSpline(x2, _row_spline_eval, _column_spline_second_derivs, yx11, yx1n);
165 :
166 : // Evaluate newly constructed column spline
167 6601 : return SplineInterpolationBase::sample(_x1, _row_spline_eval, _column_spline_second_derivs, x1);
168 : }
169 :
170 : Real
171 16802 : 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 16802 : if (deriv_var == 1)
179 : {
180 8401 : constructColumnSpline(x2, _row_spline_eval, _column_spline_second_derivs, yp1, ypn);
181 :
182 : // Evaluate derivative wrt x1 of newly constructed column spline
183 16802 : return SplineInterpolationBase::sampleDerivative(
184 8401 : _x1, _row_spline_eval, _column_spline_second_derivs, x1);
185 : }
186 :
187 : // Take derivative along x2 axis
188 8401 : else if (deriv_var == 2)
189 : {
190 8401 : constructRowSpline(x1, _column_spline_eval, _row_spline_second_derivs, yp1, ypn);
191 :
192 : // Evaluate derivative wrt x2 of newly constructed row spline
193 16802 : return SplineInterpolationBase::sampleDerivative(
194 8401 : _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 2 : 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 2 : if (deriv_var == 1)
210 : {
211 1 : constructColumnSpline(x2, _row_spline_eval, _column_spline_second_derivs, yp1, ypn);
212 :
213 : // Evaluate second derivative wrt x1 of newly constructed column spline
214 2 : return SplineInterpolationBase::sample2ndDerivative(
215 1 : _x1, _row_spline_eval, _column_spline_second_derivs, x1);
216 : }
217 :
218 : // Take second derivative along x2 axis
219 1 : else if (deriv_var == 2)
220 : {
221 1 : constructRowSpline(x1, _column_spline_eval, _row_spline_second_derivs, yp1, ypn);
222 :
223 : // Evaluate second derivative wrt x2 of newly constructed row spline
224 2 : return SplineInterpolationBase::sample2ndDerivative(
225 1 : _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 1 : 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 1 : constructColumnSpline(x2, _row_spline_eval, _column_spline_second_derivs, yx11, yx1n);
244 1 : y = SplineInterpolationBase::sample(_x1, _row_spline_eval, _column_spline_second_derivs, x1);
245 2 : dy1 = SplineInterpolationBase::sampleDerivative(
246 1 : _x1, _row_spline_eval, _column_spline_second_derivs, x1);
247 :
248 1 : constructRowSpline(x1, _column_spline_eval, _row_spline_second_derivs, yx21, yx2n);
249 2 : dy2 = SplineInterpolationBase::sampleDerivative(
250 1 : _x2, _column_spline_eval, _row_spline_second_derivs, x2);
251 1 : }
252 :
253 : void
254 8403 : 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 8403 : auto n = _x2.size();
261 :
262 : // Find the indices that bound the point x1
263 : unsigned int klo, khi;
264 8403 : 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 42018 : for (decltype(n) j = 0; j < n; ++j)
269 33615 : _column_spline_eval[j] =
270 33615 : 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 8403 : spline(_x2, column_spline_eval, row_spline_second_derivs, yx11, yx1n);
275 8403 : }
276 :
277 : void
278 15004 : 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 15004 : auto m = _x1.size();
285 :
286 : // Find the indices that bound the point x2
287 : unsigned int klo, khi;
288 15004 : 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 60024 : for (decltype(m) j = 0; j < m; ++j)
293 45020 : _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 15004 : spline(_x1, row_spline_eval, column_spline_second_derivs, yx21, yx2n);
298 15004 : }
|
