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 "BicubicInterpolation.h"
11 : #include "MooseError.h"
12 : #include "MathUtils.h"
13 : #include "MooseUtils.h"
14 :
15 : #include "libmesh/int_range.h"
16 :
17 1 : BicubicInterpolation::BicubicInterpolation(const std::vector<Real> & x1,
18 : const std::vector<Real> & x2,
19 1 : const std::vector<std::vector<Real>> & y)
20 18 : : BidimensionalInterpolation(x1, x2), _y(y)
21 : {
22 1 : errorCheck();
23 :
24 1 : auto m = _x1.size();
25 1 : auto n = _x2.size();
26 :
27 1 : _bicubic_coeffs.resize(m - 1);
28 9 : for (const auto i : index_range(_bicubic_coeffs))
29 : {
30 8 : _bicubic_coeffs[i].resize(n - 1);
31 72 : for (const auto j : index_range(_bicubic_coeffs[i]))
32 : {
33 64 : _bicubic_coeffs[i][j].resize(4);
34 320 : for (const auto k : make_range(4))
35 256 : _bicubic_coeffs[i][j][k].resize(4);
36 : }
37 : }
38 :
39 : // Precompute the coefficients
40 1 : precomputeCoefficients();
41 3 : }
42 :
43 : Real
44 2 : BicubicInterpolation::sample(const Real x1, const Real x2) const
45 : {
46 2 : return sampleInternal(x1, x2);
47 : }
48 :
49 : ADReal
50 1 : BicubicInterpolation::sample(const ADReal & x1, const ADReal & x2) const
51 : {
52 1 : return sampleInternal(x1, x2);
53 : }
54 :
55 : ChainedReal
56 1 : BicubicInterpolation::sample(const ChainedReal & x1, const ChainedReal & x2) const
57 : {
58 1 : return sampleInternal(x1, x2);
59 : }
60 :
61 : template <typename T>
62 : T
63 4 : BicubicInterpolation::sampleInternal(const T & x1, const T & x2) const
64 : {
65 : unsigned int x1l, x1u, x2l, x2u;
66 2 : T t, u;
67 4 : findInterval(_x1, x1, x1l, x1u, t);
68 4 : findInterval(_x2, x2, x2l, x2u, u);
69 :
70 4 : T sample = 0.0;
71 20 : for (const auto i : make_range(4))
72 80 : for (const auto j : make_range(4))
73 64 : sample += _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i) * MathUtils::pow(u, j);
74 :
75 6 : return sample;
76 2 : }
77 :
78 : Real
79 4 : BicubicInterpolation::sampleDerivative(Real x1, Real x2, unsigned int deriv_var) const
80 : {
81 : unsigned int x1l, x1u, x2l, x2u;
82 : Real t, u;
83 4 : findInterval(_x1, x1, x1l, x1u, t);
84 4 : findInterval(_x2, x2, x2l, x2u, u);
85 :
86 : // Take derivative along x1 axis
87 : // Note: sum from i = 1 as the first term is zero
88 4 : if (deriv_var == 1)
89 : {
90 2 : Real sample_deriv = 0.0;
91 8 : for (const auto i : make_range(1, 4))
92 30 : for (const auto j : make_range(4))
93 24 : sample_deriv +=
94 24 : i * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i - 1) * MathUtils::pow(u, j);
95 :
96 2 : Real d = _x1[x1u] - _x1[x1l];
97 :
98 2 : if (!MooseUtils::absoluteFuzzyEqual(d, 0.0))
99 2 : sample_deriv /= d;
100 :
101 2 : return sample_deriv;
102 : }
103 :
104 : // Take derivative along x2 axis
105 : // Note: sum from j = 1 as the first term is zero
106 2 : else if (deriv_var == 2)
107 : {
108 2 : Real sample_deriv = 0.0;
109 :
110 10 : for (const auto i : make_range(4))
111 32 : for (const auto j : make_range(1, 4))
112 24 : sample_deriv +=
113 24 : j * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i) * MathUtils::pow(u, j - 1);
114 :
115 2 : Real d = _x2[x2u] - _x2[x2l];
116 :
117 2 : if (!MooseUtils::absoluteFuzzyEqual(d, Real(0.0)))
118 2 : sample_deriv /= d;
119 :
120 2 : return sample_deriv;
121 : }
122 :
123 : else
124 0 : mooseError("deriv_var must be either 1 or 2 in BicubicInterpolation");
125 : }
126 :
127 : Real
128 6 : BicubicInterpolation::sample2ndDerivative(Real x1, Real x2, unsigned int deriv_var) const
129 : {
130 : unsigned int x1l, x1u, x2l, x2u;
131 : Real t, u;
132 6 : findInterval(_x1, x1, x1l, x1u, t);
133 6 : findInterval(_x2, x2, x2l, x2u, u);
134 :
135 : // Take derivative along x1 axis
136 : // Note: sum from i = 2 as the first two terms are zero
137 6 : if (deriv_var == 1)
138 : {
139 2 : Real sample_deriv = 0.0;
140 6 : for (const auto i : make_range(2, 4))
141 20 : for (const auto j : make_range(4))
142 32 : sample_deriv += i * (i - 1) * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i - 2) *
143 16 : MathUtils::pow(u, j);
144 :
145 2 : Real d = _x1[x1u] - _x1[x1l];
146 :
147 2 : if (!MooseUtils::absoluteFuzzyEqual(d, Real(0.0)))
148 2 : sample_deriv /= (d * d);
149 :
150 2 : return sample_deriv;
151 : }
152 :
153 : // Take derivative along x2 axis
154 : // Note: sum from j = 2 as the first two terms are zero
155 4 : else if (deriv_var == 2)
156 : {
157 2 : Real sample_deriv = 0.0;
158 10 : for (const auto i : make_range(4))
159 24 : for (const auto j : make_range(2, 4))
160 32 : sample_deriv += j * (j - 1) * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i) *
161 16 : MathUtils::pow(u, j - 2);
162 :
163 2 : Real d = _x2[x2u] - _x2[x2l];
164 :
165 2 : if (!MooseUtils::absoluteFuzzyEqual(d, Real(0.0)))
166 2 : sample_deriv /= (d * d);
167 :
168 2 : return sample_deriv;
169 : }
170 :
171 : // Take mixed derivative along x1 and x2 axes
172 : // Note: sum from i = 1 and j = 1 as the first terms are zero
173 2 : else if (deriv_var == 3)
174 : {
175 2 : Real sample_deriv = 0.0;
176 8 : for (const auto i : make_range(1, 4))
177 24 : for (const auto j : make_range(1, 4))
178 36 : sample_deriv += i * j * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i - 1) *
179 18 : MathUtils::pow(u, j - 1);
180 :
181 2 : const Real d1 = _x1[x1u] - _x1[x1l];
182 2 : const Real d2 = _x2[x2u] - _x2[x2l];
183 :
184 4 : if (!MooseUtils::absoluteFuzzyEqual(d1, Real(0.0)) &&
185 4 : !MooseUtils::absoluteFuzzyEqual(d2, Real(0.0)))
186 2 : sample_deriv /= (d1 * d2);
187 :
188 2 : return sample_deriv;
189 : }
190 :
191 : else
192 0 : mooseError("deriv_var must be either 1, 2 or 3 in BicubicInterpolation");
193 : }
194 :
195 : void
196 1 : BicubicInterpolation::sampleValueAndDerivatives(
197 : Real x1, Real x2, Real & y, Real & dy1, Real & dy2) const
198 : {
199 1 : sampleValueAndDerivativesInternal(x1, x2, y, dy1, dy2);
200 1 : }
201 :
202 : void
203 1 : BicubicInterpolation::sampleValueAndDerivatives(
204 : const ADReal & x1, const ADReal & x2, ADReal & y, ADReal & dy1, ADReal & dy2) const
205 : {
206 1 : sampleValueAndDerivativesInternal(x1, x2, y, dy1, dy2);
207 1 : }
208 : void
209 1 : BicubicInterpolation::sampleValueAndDerivatives(const ChainedReal & x1,
210 : const ChainedReal & x2,
211 : ChainedReal & y,
212 : ChainedReal & dy1,
213 : ChainedReal & dy2) const
214 : {
215 1 : sampleValueAndDerivativesInternal(x1, x2, y, dy1, dy2);
216 1 : }
217 :
218 : template <typename T>
219 : void
220 3 : BicubicInterpolation::sampleValueAndDerivativesInternal(T x1, T x2, T & y, T & dy1, T & dy2) const
221 : {
222 : unsigned int x1l, x1u, x2l, x2u;
223 2 : T t, u;
224 3 : findInterval(_x1, x1, x1l, x1u, t);
225 3 : findInterval(_x2, x2, x2l, x2u, u);
226 :
227 3 : y = 0.0;
228 15 : for (const auto i : make_range(4))
229 60 : for (const auto j : make_range(4))
230 48 : y += _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i) * MathUtils::pow(u, j);
231 :
232 : // Note: sum from i = 1 as the first term is zero
233 3 : dy1 = 0.0;
234 12 : for (const auto i : make_range(1, 4))
235 45 : for (const auto j : make_range(4))
236 36 : dy1 += i * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i - 1) * MathUtils::pow(u, j);
237 :
238 : // Note: sum from j = 1 as the first term is zero
239 3 : dy2 = 0.0;
240 15 : for (const auto i : make_range(4))
241 48 : for (const auto j : make_range(1, 4))
242 36 : dy2 += j * _bicubic_coeffs[x1l][x2l][i][j] * MathUtils::pow(t, i) * MathUtils::pow(u, j - 1);
243 :
244 3 : T d1 = _x1[x1u] - _x1[x1l];
245 :
246 3 : if (!MooseUtils::absoluteFuzzyEqual(d1, 0.0))
247 3 : dy1 /= d1;
248 :
249 3 : T d2 = _x2[x2u] - _x2[x2l];
250 :
251 3 : if (!MooseUtils::absoluteFuzzyEqual(d2, 0.0))
252 3 : dy2 /= d2;
253 3 : }
254 :
255 : void
256 1 : BicubicInterpolation::precomputeCoefficients()
257 : {
258 : // Calculate the first derivatives in each direction at each point, and the
259 : // mixed second derivative
260 1 : std::vector<std::vector<Real>> dy_dx1, dy_dx2, d2y_dx1x2;
261 1 : tableDerivatives(dy_dx1, dy_dx2, d2y_dx1x2);
262 :
263 : // Now solve for the coefficients at each point in the grid
264 9 : for (const auto i : index_range(_bicubic_coeffs))
265 72 : for (const auto j : index_range(_bicubic_coeffs[i]))
266 : {
267 : // Distance between corner points in each direction
268 64 : const Real d1 = _x1[i + 1] - _x1[i];
269 64 : const Real d2 = _x2[j + 1] - _x2[j];
270 :
271 : // Values of function and derivatives of the four corner points
272 64 : std::vector<Real> y = {_y[i][j], _y[i + 1][j], _y[i + 1][j + 1], _y[i][j + 1]};
273 : std::vector<Real> y1 = {
274 64 : dy_dx1[i][j], dy_dx1[i + 1][j], dy_dx1[i + 1][j + 1], dy_dx1[i][j + 1]};
275 : std::vector<Real> y2 = {
276 64 : dy_dx2[i][j], dy_dx2[i + 1][j], dy_dx2[i + 1][j + 1], dy_dx2[i][j + 1]};
277 : std::vector<Real> y12 = {
278 64 : d2y_dx1x2[i][j], d2y_dx1x2[i + 1][j], d2y_dx1x2[i + 1][j + 1], d2y_dx1x2[i][j + 1]};
279 :
280 64 : std::vector<Real> cl(16), x(16);
281 : Real xx;
282 64 : unsigned int count = 0;
283 :
284 : // Temporary vector used in the matrix multiplication
285 320 : for (const auto k : make_range(4))
286 : {
287 256 : x[k] = y[k];
288 256 : x[k + 4] = y1[k] * d1;
289 256 : x[k + 8] = y2[k] * d2;
290 256 : x[k + 12] = y12[k] * d1 * d2;
291 : }
292 :
293 : // Multiply by the matrix of constants
294 1088 : for (const auto k : make_range(16))
295 : {
296 1024 : xx = 0.0;
297 17408 : for (const auto q : make_range(16))
298 16384 : xx += _wt[k][q] * x[q];
299 :
300 1024 : cl[k] = xx;
301 : }
302 :
303 : // Unpack results into coefficient table
304 320 : for (const auto k : make_range(4))
305 1280 : for (const auto q : make_range(4))
306 : {
307 1024 : _bicubic_coeffs[i][j][k][q] = cl[count];
308 1024 : count++;
309 : }
310 64 : }
311 1 : }
312 :
313 : void
314 1 : BicubicInterpolation::tableDerivatives(std::vector<std::vector<Real>> & dy_dx1,
315 : std::vector<std::vector<Real>> & dy_dx2,
316 : std::vector<std::vector<Real>> & d2y_dx1x2)
317 : {
318 1 : const auto m = _x1.size();
319 1 : const auto n = _x2.size();
320 1 : dy_dx1.resize(m);
321 1 : dy_dx2.resize(m);
322 1 : d2y_dx1x2.resize(m);
323 :
324 10 : for (const auto i : make_range(m))
325 : {
326 9 : dy_dx1[i].resize(n);
327 9 : dy_dx2[i].resize(n);
328 9 : d2y_dx1x2[i].resize(n);
329 : }
330 :
331 : // Central difference calculations of derivatives
332 10 : for (const auto i : make_range(m))
333 90 : for (const auto j : make_range(n))
334 : {
335 : // Derivative wrt x1
336 81 : if (i == 0)
337 : {
338 9 : dy_dx1[i][j] = (_y[i + 1][j] - _y[i][j]) / (_x1[i + 1] - _x1[i]);
339 : }
340 72 : else if (i == m - 1)
341 9 : dy_dx1[i][j] = (_y[i][j] - _y[i - 1][j]) / (_x1[i] - _x1[i - 1]);
342 : else
343 63 : dy_dx1[i][j] = (_y[i + 1][j] - _y[i - 1][j]) / (_x1[i + 1] - _x1[i - 1]);
344 :
345 : // Derivative wrt x2
346 81 : if (j == 0)
347 9 : dy_dx2[i][j] = (_y[i][j + 1] - _y[i][j]) / (_x2[j + 1] - _x2[j]);
348 72 : else if (j == n - 1)
349 9 : dy_dx2[i][j] = (_y[i][j] - _y[i][j - 1]) / (_x2[j] - _x2[j - 1]);
350 : else
351 63 : dy_dx2[i][j] = (_y[i][j + 1] - _y[i][j - 1]) / (_x2[j + 1] - _x2[j - 1]);
352 :
353 : // Mixed derivative d2y/dx1x2
354 81 : if (i == 0 && j == 0)
355 2 : d2y_dx1x2[i][j] = (_y[i + 1][j + 1] - _y[i + 1][j] - _y[i][j + 1] + _y[i][j]) /
356 1 : (_x1[i + 1] - _x1[i]) / (_x2[j + 1] - _x2[j]);
357 80 : else if (i == 0 && j == n - 1)
358 2 : d2y_dx1x2[i][j] = (_y[i + 1][j] - _y[i + 1][j - 1] - _y[i][j] + _y[i][j - 1]) /
359 1 : (_x1[i + 1] - _x1[i]) / (_x2[j] - _x2[j - 1]);
360 79 : else if (i == m - 1 && j == 0)
361 2 : d2y_dx1x2[i][j] = (_y[i][j + 1] - _y[i][j] - _y[i - 1][j + 1] + _y[i - 1][j]) /
362 1 : (_x1[i] - _x1[i - 1]) / (_x2[j + 1] - _x2[j]);
363 78 : else if (i == m - 1 && j == n - 1)
364 2 : d2y_dx1x2[i][j] = (_y[i][j] - _y[i][j - 1] - _y[i - 1][j] + _y[i - 1][j - 1]) /
365 1 : (_x1[i] - _x1[i - 1]) / (_x2[j] - _x2[j - 1]);
366 77 : else if (i == 0)
367 14 : d2y_dx1x2[i][j] = (_y[i + 1][j + 1] - _y[i + 1][j - 1] - _y[i][j + 1] + _y[i][j - 1]) /
368 7 : (_x1[i + 1] - _x1[i]) / (_x2[j + 1] - _x2[j - 1]);
369 70 : else if (i == m - 1)
370 14 : d2y_dx1x2[i][j] = (_y[i][j + 1] - _y[i][j - 1] - _y[i - 1][j + 1] + _y[i - 1][j - 1]) /
371 7 : (_x1[i] - _x1[i - 1]) / (_x2[j + 1] - _x2[j - 1]);
372 63 : else if (j == 0)
373 14 : d2y_dx1x2[i][j] = (_y[i + 1][j + 1] - _y[i + 1][j] - _y[i - 1][j + 1] + _y[i - 1][j]) /
374 7 : (_x1[i + 1] - _x1[i - 1]) / (_x2[j + 1] - _x2[j]);
375 56 : else if (j == n - 1)
376 14 : d2y_dx1x2[i][j] = (_y[i + 1][j] - _y[i + 1][j - 1] - _y[i - 1][j] + _y[i - 1][j - 1]) /
377 7 : (_x1[i + 1] - _x1[i - 1]) / (_x2[j] - _x2[j - 1]);
378 : else
379 49 : d2y_dx1x2[i][j] =
380 49 : (_y[i + 1][j + 1] - _y[i + 1][j - 1] - _y[i - 1][j + 1] + _y[i - 1][j - 1]) /
381 49 : (_x1[i + 1] - _x1[i - 1]) / (_x2[j + 1] - _x2[j - 1]);
382 : }
383 1 : }
384 :
385 : template <typename T>
386 : void
387 34 : BicubicInterpolation::findInterval(
388 : const std::vector<Real> & x, const T & xi, unsigned int & klo, unsigned int & khi, T & xs) const
389 : {
390 : // Find the indices that bracket the point xi
391 34 : klo = 0;
392 : mooseAssert(x.size() >= 2,
393 : "There must be at least two points in the table in BicubicInterpolation");
394 :
395 34 : khi = x.size() - 1;
396 136 : while (khi - klo > 1)
397 : {
398 102 : unsigned int kmid = (khi + klo) / 2;
399 102 : if (x[kmid] > xi)
400 51 : khi = kmid;
401 : else
402 51 : klo = kmid;
403 : }
404 :
405 : // Now find the scaled position, normalized to [0,1]
406 34 : Real d = x[khi] - x[klo];
407 34 : xs = xi - x[klo];
408 :
409 34 : if (!MooseUtils::absoluteFuzzyEqual(d, Real(0.0)))
410 34 : xs /= d;
411 34 : }
412 :
413 : template void BicubicInterpolation::findInterval(const std::vector<Real> & x,
414 : const Real & xi,
415 : unsigned int & klo,
416 : unsigned int & khi,
417 : Real & xs) const;
418 : template void BicubicInterpolation::findInterval(const std::vector<Real> & x,
419 : const ADReal & xi,
420 : unsigned int & klo,
421 : unsigned int & khi,
422 : ADReal & xs) const;
423 :
424 : void
425 1 : BicubicInterpolation::errorCheck()
426 : {
427 1 : auto m = _x1.size(), n = _x2.size();
428 :
429 1 : if (_y.size() != m)
430 0 : mooseError("y row dimension does not match the size of x1 in BicubicInterpolation");
431 : else
432 10 : for (const auto i : index_range(_y))
433 9 : if (_y[i].size() != n)
434 0 : mooseError("y column dimension does not match the size of x2 in BicubicInterpolation");
435 1 : }
|
