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