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