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 "BilinearInterpolation.h"
11 : #include "libmesh/int_range.h"
12 :
13 : int BilinearInterpolation::_file_number = 0;
14 :
15 158 : BilinearInterpolation::BilinearInterpolation(const std::vector<Real> & xaxis,
16 : const std::vector<Real> & yaxis,
17 158 : const ColumnMajorMatrix & zsurface)
18 158 : : BidimensionalInterpolation(xaxis, yaxis), _z_surface(zsurface)
19 : {
20 158 : }
21 :
22 : void
23 6256 : BilinearInterpolation::getNeighborIndices(const std::vector<Real> & inArr,
24 : Real x,
25 : unsigned int & lowerX,
26 : unsigned int & upperX) const
27 : {
28 6256 : int N = inArr.size();
29 6256 : if (x <= inArr[0])
30 : {
31 24 : lowerX = 0;
32 24 : upperX = 0;
33 : }
34 6232 : else if (x >= inArr[N - 1])
35 : {
36 24 : lowerX = N - 1;
37 24 : upperX = N - 1;
38 : }
39 : else
40 : {
41 81885 : for (const auto i : make_range(1, N))
42 : {
43 81885 : if (x < inArr[i])
44 : {
45 3368 : lowerX = i - 1;
46 3368 : upperX = i;
47 3368 : break;
48 : }
49 78517 : else if (x == inArr[i])
50 : {
51 2840 : lowerX = i;
52 2840 : upperX = i;
53 2840 : break;
54 : }
55 : }
56 : }
57 6256 : }
58 :
59 : Real
60 3090 : BilinearInterpolation::sample(Real s1, Real s2) const
61 : {
62 3090 : return sampleInternal(s1, s2);
63 : }
64 :
65 : ADReal
66 1 : BilinearInterpolation::sample(const ADReal & s1, const ADReal & s2) const
67 : {
68 1 : return sampleInternal(s1, s2);
69 : }
70 :
71 : ChainedReal
72 1 : BilinearInterpolation::sample(const ChainedReal & s1, const ChainedReal & s2) const
73 : {
74 1 : return sampleInternal(s1, s2);
75 : }
76 :
77 : template <typename T>
78 : T
79 3092 : BilinearInterpolation::sampleInternal(const T & s1, const T & s2) const
80 : {
81 : // first find 4 neighboring points
82 3092 : unsigned int lx = 0; // index of x coordinate of adjacent grid point to left of P
83 3092 : unsigned int ux = 0; // index of x coordinate of adjacent grid point to right of P
84 3092 : getNeighborIndices(_x1, MetaPhysicL::raw_value(s1), lx, ux);
85 :
86 3092 : unsigned int ly = 0; // index of y coordinate of adjacent grid point below P
87 3092 : unsigned int uy = 0; // index of y coordinate of adjacent grid point above P
88 3092 : getNeighborIndices(_x2, MetaPhysicL::raw_value(s2), ly, uy);
89 :
90 3092 : const Real & fQ11 = _z_surface(ly, lx);
91 3092 : const Real & fQ21 = _z_surface(ly, ux);
92 3092 : const Real & fQ12 = _z_surface(uy, lx);
93 3092 : const Real & fQ22 = _z_surface(uy, ux);
94 :
95 : // if point exactly found on a node do not interpolate
96 3092 : if ((lx == ux) && (ly == uy))
97 : {
98 649 : return fQ11;
99 : }
100 :
101 2443 : const auto & x = s1;
102 2443 : const auto & y = s2;
103 2443 : const Real & x1 = _x1[lx];
104 2443 : const Real & x2 = _x1[ux];
105 2443 : const Real & y1 = _x2[ly];
106 2443 : const Real & y2 = _x2[uy];
107 :
108 : // if x1 lies exactly on an xAxis node do linear interpolation
109 2443 : if (lx == ux)
110 643 : return fQ11 + (fQ12 - fQ11) * (y - y1) / (y2 - y1);
111 :
112 : // if x2 lies exactly on an yAxis node do linear interpolation
113 1800 : if (ly == uy)
114 899 : return fQ11 + (fQ21 - fQ11) * (x - x1) / (x2 - x1);
115 :
116 901 : auto fxy = fQ11 * (x2 - x) * (y2 - y);
117 901 : fxy += fQ21 * (x - x1) * (y2 - y);
118 901 : fxy += fQ12 * (x2 - x) * (y - y1);
119 901 : fxy += fQ22 * (x - x1) * (y - y1);
120 901 : fxy /= ((x2 - x1) * (y2 - y1));
121 901 : return fxy;
122 2 : }
123 :
124 : Real
125 36 : BilinearInterpolation::sampleDerivative(Real s1, Real s2, const unsigned int deriv_var) const
126 : {
127 36 : return sampleDerivativeInternal(s1, s2, deriv_var);
128 : }
129 :
130 : ADReal
131 0 : BilinearInterpolation::sampleDerivative(const ADReal & s1,
132 : const ADReal & s2,
133 : const unsigned int deriv_var) const
134 : {
135 0 : return sampleDerivativeInternal(s1, s2, deriv_var);
136 : }
137 :
138 : ChainedReal
139 0 : BilinearInterpolation::sampleDerivative(const ChainedReal & s1,
140 : const ChainedReal & s2,
141 : const unsigned int deriv_var) const
142 : {
143 0 : return sampleDerivativeInternal(s1, s2, deriv_var);
144 : }
145 :
146 : template <typename T>
147 : T
148 36 : BilinearInterpolation::sampleDerivativeInternal(const T s1,
149 : const T s2,
150 : const unsigned int deriv_var) const
151 : {
152 : // check input
153 36 : if (deriv_var != 1 && deriv_var != 2)
154 0 : mooseError("deriv_var must equal 1 or 2");
155 :
156 : // first find 4 neighboring points
157 36 : unsigned int lx = 0; // index of x coordinate of adjacent grid point to left of P
158 36 : unsigned int ux = 0; // index of x coordinate of adjacent grid point to right of P
159 36 : getNeighborIndices(_x1, MetaPhysicL::raw_value(s1), lx, ux);
160 :
161 36 : unsigned int ly = 0; // index of y coordinate of adjacent grid point below P
162 36 : unsigned int uy = 0; // index of y coordinate of adjacent grid point above P
163 36 : getNeighborIndices(_x2, MetaPhysicL::raw_value(s2), ly, uy);
164 :
165 : // xy indexing
166 : // 11 is point on lower-left (x low, y low), 22 on upper-right (x high, y high)
167 36 : const Real & fQ11 = _z_surface(ly, lx);
168 36 : const Real & fQ21 = _z_surface(ly, ux);
169 36 : const Real & fQ12 = _z_surface(uy, lx);
170 36 : const Real & fQ22 = _z_surface(uy, ux);
171 : // when on a grid bound or node, lx can be equal to ux or ly be equal to uy
172 : // this would lead to a 0 slope, which isn't desirable
173 : // we then refer to 0 as the coordinate before, so 00 is lx-1, ly-1
174 : // and 3 as the coordinate after, so 33 is ux+1, uy+1
175 :
176 36 : const auto & x = s1;
177 36 : const auto & y = s2;
178 36 : const Real & x1 = _x1[lx];
179 36 : const Real & x2 = _x1[ux];
180 36 : const Real & y1 = _x2[ly];
181 36 : const Real & y2 = _x2[uy];
182 :
183 : // Grid sizes
184 36 : const auto nx1 = _x1.size();
185 36 : const auto nx2 = _x2.size();
186 :
187 : // On interior grid lines, the equal weighting of both sides of both neighbor
188 : // cells slopes is an implementation choice
189 : // Alternatively on interior grid nodes, we instead use the super cell around the grid node
190 : // instead of a weighting of the four neighbor cells
191 :
192 : // Check all four grid corners, use a quarter of the slope in the cell next to the corner
193 36 : if (ux == 0 && uy == 0) // bottom left node
194 : {
195 2 : const auto & fQ13 = _z_surface(ly + 1, lx); // fQ at (x1,y3)
196 2 : const auto & fQ31 = _z_surface(ly, lx + 1); // fQ at (x3,y1)
197 2 : const auto & fQ33 = _z_surface(ly + 1, lx + 1); // fQ at (x3,y3)
198 2 : const Real & x3 = _x1[lx + 1]; // ux value
199 2 : const Real & y3 = _x2[ly + 1]; // uy value
200 :
201 2 : if (deriv_var == 1)
202 : {
203 1 : auto dfdx = fQ11 * (y - y3);
204 1 : dfdx += fQ31 * (y3 - y);
205 1 : dfdx += fQ13 * (y1 - y);
206 1 : dfdx += fQ33 * (y - y1);
207 1 : dfdx /= ((x3 - x1) * (y3 - y1));
208 2 : return dfdx;
209 0 : }
210 1 : else if (deriv_var == 2)
211 : {
212 1 : auto dfdy = fQ11 * (x - x3);
213 1 : dfdy += fQ31 * (x1 - x);
214 1 : dfdy += fQ13 * (x3 - x);
215 1 : dfdy += fQ33 * (x - x1);
216 1 : dfdy /= ((x3 - x1) * (y3 - y1));
217 1 : return dfdy;
218 0 : }
219 0 : }
220 34 : else if (uy == 0 && lx == nx1 - 1) // bottom right node
221 : {
222 2 : const auto & fQ01 = _z_surface(ly, lx - 1); // fQ at (x0,y1)
223 2 : const auto & fQ03 = _z_surface(ly + 1, lx - 1); // fQ at (x0,y3)
224 2 : const auto & fQ23 = _z_surface(ly + 1, lx); // fQ at (x2,y3)
225 2 : const Real & x0 = _x1[lx - 1]; // lx value
226 2 : const Real & y3 = _x2[ly + 1]; // uy value
227 :
228 2 : if (deriv_var == 1)
229 : {
230 1 : auto dfdx = fQ01 * (y - y3);
231 1 : dfdx += fQ21 * (y3 - y);
232 1 : dfdx += fQ03 * (y1 - y);
233 1 : dfdx += fQ23 * (y - y1);
234 1 : dfdx /= ((x2 - x0) * (y3 - y1));
235 2 : return dfdx;
236 0 : }
237 1 : else if (deriv_var == 2)
238 : {
239 1 : auto dfdy = fQ01 * (x - x2);
240 1 : dfdy += fQ21 * (x0 - x);
241 1 : dfdy += fQ03 * (x2 - x);
242 1 : dfdy += fQ23 * (x - x0);
243 1 : dfdy /= ((x2 - x0) * (y3 - y1));
244 1 : return dfdy;
245 0 : }
246 0 : }
247 32 : else if (ly == nx2 - 1 && ux == 0) // top left node
248 : {
249 2 : const auto & fQ10 = _z_surface(ly - 1, lx); // fQ at (x1,y0)
250 2 : const auto & fQ30 = _z_surface(ly - 1, lx + 1); // fQ at (x3,y0)
251 2 : const auto & fQ32 = _z_surface(ly, lx + 1); // fQ at (x3,y2)
252 2 : const Real & x3 = _x1[lx + 1]; // ux value
253 2 : const Real & y0 = _x2[ly - 1]; // ly value
254 :
255 2 : if (deriv_var == 1)
256 : {
257 1 : auto dfdx = fQ10 * (y - y2);
258 1 : dfdx += fQ30 * (y2 - y);
259 1 : dfdx += fQ12 * (y0 - y);
260 1 : dfdx += fQ32 * (y - y0);
261 1 : dfdx /= ((x3 - x1) * (y2 - y0));
262 2 : return dfdx;
263 0 : }
264 1 : else if (deriv_var == 2)
265 : {
266 1 : auto dfdy = fQ10 * (x - x3);
267 1 : dfdy += fQ30 * (x1 - x);
268 1 : dfdy += fQ12 * (x3 - x);
269 1 : dfdy += fQ32 * (x - x1);
270 1 : dfdy /= ((x3 - x1) * (y2 - y0));
271 1 : return dfdy;
272 0 : }
273 0 : }
274 30 : else if (ly == nx2 - 1 && lx == nx1 - 1) // top right node
275 : {
276 2 : const auto & fQ00 = _z_surface(ly - 1, lx - 1); // fQ at (x0,y0)
277 2 : const auto & fQ20 = _z_surface(ly - 1, lx); // fQ at (x2,y0)
278 2 : const auto & fQ02 = _z_surface(ly, lx - 1); // fQ at (x0,y2)
279 2 : const Real & x0 = _x1[lx - 1]; // lx value
280 2 : const Real & y0 = _x2[ly - 1]; // ly value
281 :
282 2 : if (deriv_var == 1)
283 : {
284 1 : auto dfdx = fQ00 * (y - y2);
285 1 : dfdx += fQ20 * (y2 - y);
286 1 : dfdx += fQ02 * (y0 - y);
287 1 : dfdx += fQ22 * (y - y0);
288 1 : dfdx /= ((x2 - x0) * (y2 - y0));
289 2 : return dfdx;
290 0 : }
291 1 : else if (deriv_var == 2)
292 : {
293 1 : auto dfdy = fQ00 * (x - x2);
294 1 : dfdy += fQ20 * (x0 - x);
295 1 : dfdy += fQ02 * (x2 - x);
296 1 : dfdy += fQ22 * (x - x0);
297 1 : dfdy /= ((x2 - x0) * (y2 - y0));
298 1 : return dfdy;
299 0 : }
300 0 : }
301 :
302 : // Nodes along the four grid bounds, use the two grid cells adjacent to the node
303 28 : else if (ux == 0 && ly == uy && ux == lx) // when along left boundary, at a grid node
304 : {
305 2 : const auto & fQ10 = _z_surface(uy - 1, lx); // fQ at (x1,y0)
306 2 : const auto & fQ30 = _z_surface(uy, ux + 1); // fQ at (x3,y0)
307 2 : const auto & fQ13 = _z_surface(uy + 1, lx); // fQ at (x1,y3)
308 2 : const auto & fQ33 = _z_surface(uy + 1, ux + 1); // fQ at (x3,y3)
309 :
310 2 : const Real & x3 = _x1[lx + 1]; // ux value
311 2 : const Real & y0 = _x2[ly - 1]; // ly value
312 2 : const Real & y3 = _x2[ly + 1]; // uy value
313 :
314 2 : if (deriv_var == 1)
315 : {
316 1 : auto dfdx = fQ10 * (y - y3);
317 1 : dfdx += fQ30 * (y3 - y);
318 1 : dfdx += fQ13 * (y0 - y);
319 1 : dfdx += fQ33 * (y - y0);
320 1 : dfdx /= ((x3 - x1) * (y3 - y0));
321 2 : return dfdx;
322 0 : }
323 1 : else if (deriv_var == 2)
324 : {
325 1 : auto dfdy = fQ10 * (x - x3);
326 1 : dfdy += fQ30 * (x1 - x);
327 1 : dfdy += fQ13 * (x3 - x);
328 1 : dfdy += fQ33 * (x - x1);
329 1 : dfdy /= ((x3 - x1) * (y3 - y0));
330 1 : return dfdy;
331 0 : }
332 0 : }
333 26 : else if (lx == nx1 - 1 && ly == uy && ux == lx) // when along right boundary, at a grid node
334 : {
335 2 : const auto & fQ00 = _z_surface(uy - 1, lx - 1); // fQ at (x0,y0)
336 2 : const auto & fQ10 = _z_surface(uy - 1, lx); // fQ at (x1,y0)
337 2 : const auto & fQ03 = _z_surface(uy + 1, lx - 1); // fQ at (x0,y3)
338 2 : const auto & fQ13 = _z_surface(uy + 1, lx); // fQ at (x1,y3)
339 :
340 2 : const Real & x0 = _x1[lx - 1]; // lx value
341 2 : const Real & y0 = _x2[ly - 1]; // ly value
342 2 : const Real & y3 = _x2[ly + 1]; // uy value
343 :
344 2 : if (deriv_var == 1)
345 : {
346 1 : auto dfdx = fQ00 * (y - y3);
347 1 : dfdx += fQ10 * (y3 - y);
348 1 : dfdx += fQ03 * (y0 - y);
349 1 : dfdx += fQ13 * (y - y0);
350 1 : dfdx /= ((x1 - x0) * (y3 - y0));
351 2 : return dfdx;
352 0 : }
353 1 : else if (deriv_var == 2)
354 : {
355 1 : auto dfdy = fQ00 * (x - x1);
356 1 : dfdy += fQ10 * (x0 - x);
357 1 : dfdy += fQ03 * (x1 - x);
358 1 : dfdy += fQ13 * (x - x0);
359 1 : dfdy /= ((x1 - x0) * (y3 - y0));
360 1 : return dfdy;
361 0 : }
362 0 : }
363 24 : else if (uy == nx2 - 1 && ly == uy && ux == lx) // when along top boundary, at a grid node
364 : {
365 2 : const auto & fQ00 = _z_surface(uy - 1, lx - 1); // fQ at (x0,y0)
366 2 : const auto & fQ01 = _z_surface(uy, lx - 1); // fQ at (x0,y1)
367 2 : const auto & fQ30 = _z_surface(uy - 1, lx + 1); // fQ at (x3,y0)
368 2 : const auto & fQ31 = _z_surface(uy, lx + 1); // fQ at (x3,y1)
369 :
370 2 : const Real & x0 = _x1[lx - 1]; // lx value
371 2 : const Real & x3 = _x1[lx + 1]; // ux value
372 2 : const Real & y0 = _x2[ly - 1]; // ly value
373 :
374 2 : if (deriv_var == 1)
375 : {
376 1 : auto dfdx = fQ00 * (y - y1);
377 1 : dfdx += fQ30 * (y1 - y);
378 1 : dfdx += fQ01 * (y0 - y);
379 1 : dfdx += fQ31 * (y - y0);
380 1 : dfdx /= ((x3 - x0) * (y1 - y0));
381 2 : return dfdx;
382 0 : }
383 1 : else if (deriv_var == 2)
384 : {
385 1 : auto dfdy = fQ00 * (x - x3);
386 1 : dfdy += fQ30 * (x0 - x);
387 1 : dfdy += fQ01 * (x3 - x);
388 1 : dfdy += fQ31 * (x - x0);
389 1 : dfdy /= ((x3 - x0) * (y1 - y0));
390 1 : return dfdy;
391 0 : }
392 0 : }
393 22 : else if (uy == 0 && ly == uy && ux == lx) // when along bottom boundary, at a grid node
394 : {
395 2 : const auto & fQ01 = _z_surface(ly, lx - 1); // fQ at (x0,y1)
396 2 : const auto & fQ03 = _z_surface(ly + 1, lx - 1); // fQ at (x0,y3)
397 2 : const auto & fQ31 = _z_surface(ly, lx + 1); // fQ at (x3,y1)
398 2 : const auto & fQ33 = _z_surface(ly + 1, lx + 1); // fQ at (x3,y3)
399 :
400 2 : const Real & x0 = _x1[lx - 1]; // lx value
401 2 : const Real & x3 = _x1[lx + 1]; // ux value
402 2 : const Real & y3 = _x2[ly + 1]; // uy value
403 :
404 2 : if (deriv_var == 1)
405 : {
406 1 : auto dfdx = fQ01 * (y - y3);
407 1 : dfdx += fQ31 * (y3 - y);
408 1 : dfdx += fQ03 * (y1 - y);
409 1 : dfdx += fQ33 * (y - y1);
410 1 : dfdx /= ((x3 - x0) * (y3 - y1));
411 2 : return dfdx;
412 0 : }
413 1 : else if (deriv_var == 2)
414 : {
415 1 : auto dfdy = fQ01 * (x - x3);
416 1 : dfdy += fQ31 * (x0 - x);
417 1 : dfdy += fQ03 * (x3 - x);
418 1 : dfdy += fQ33 * (x - x0);
419 1 : dfdy /= ((x3 - x0) * (y3 - y1));
420 1 : return dfdy;
421 0 : }
422 0 : }
423 :
424 : // at a grid node inside the domain, use the super box around
425 20 : else if (ux == lx && uy == ly)
426 : {
427 2 : const auto & fQ00 = _z_surface(ly - 1, lx - 1); // fQ at (x0,y0)
428 2 : const auto & fQ03 = _z_surface(ly + 1, lx - 1); // fQ at (x0,y3)
429 2 : const auto & fQ30 = _z_surface(ly - 1, lx + 1); // fQ at (x3,y0)
430 2 : const auto & fQ33 = _z_surface(ly + 1, lx + 1); // fQ at (x3,y3)
431 :
432 2 : const Real & x0 = _x1[lx - 1]; // lx value
433 2 : const Real & x3 = _x1[lx + 1]; // ux value
434 2 : const Real & y0 = _x2[ly - 1]; // ly value
435 2 : const Real & y3 = _x2[ly + 1]; // uy value
436 :
437 2 : if (deriv_var == 1)
438 : {
439 1 : auto dfdx = fQ00 * (y - y3);
440 1 : dfdx += fQ30 * (y3 - y);
441 1 : dfdx += fQ03 * (y0 - y);
442 1 : dfdx += fQ33 * (y - y0);
443 1 : dfdx /= ((x3 - x0) * (y3 - y0));
444 2 : return dfdx;
445 0 : }
446 1 : else if (deriv_var == 2)
447 : {
448 1 : auto dfdy = fQ00 * (x - x3);
449 1 : dfdy += fQ30 * (x0 - x);
450 1 : dfdy += fQ03 * (x3 - x);
451 1 : dfdy += fQ33 * (x - x0);
452 1 : dfdy /= ((x3 - x0) * (y3 - y0));
453 1 : return dfdy;
454 0 : }
455 : }
456 :
457 : // Inside the domain, not at a grid node
458 : // calculate derivative wrt to x
459 18 : if (deriv_var == 1)
460 : {
461 : // Find derivative when on interval between two nodes
462 9 : if (y == y1)
463 3 : return (fQ21 - fQ11) / (x2 - x1);
464 : // interior grid line
465 6 : else if (lx == ux && lx > 0 && ux < nx1 - 1)
466 : {
467 : // expand grid size so x1 does not equal x2
468 1 : const auto & fQ01 = _z_surface(ly, lx - 1); // new lx at ly
469 1 : const auto & fQ31 = _z_surface(ly, lx + 1); // new ux at ly
470 1 : const auto & fQ02 = _z_surface(uy, lx - 1); // new lx at uy
471 1 : const auto & fQ32 = _z_surface(uy, lx + 1); // new ux at uy
472 :
473 1 : const Real & x0 = _x1[lx - 1]; // lx value
474 1 : const Real & x3 = _x1[lx + 1]; // ux value
475 :
476 1 : auto dfdx_a = fQ01 * (y - y2);
477 1 : dfdx_a += fQ11 * (y2 - y);
478 1 : dfdx_a += fQ02 * (y1 - y);
479 1 : dfdx_a += fQ12 * (y - y1);
480 1 : dfdx_a /= ((x1 - x0) * (y2 - y1));
481 :
482 1 : auto dfdx_b = fQ11 * (y - y2);
483 1 : dfdx_b += fQ31 * (y2 - y);
484 1 : dfdx_b += fQ12 * (y1 - y);
485 1 : dfdx_b += fQ32 * (y - y1);
486 1 : dfdx_b /= ((x3 - x1) * (y2 - y1));
487 1 : return 0.5 * (dfdx_a + dfdx_b);
488 0 : }
489 : // left boundary
490 5 : else if (lx == ux && lx == 0)
491 : {
492 1 : const auto & fQ31 = _z_surface(ly, lx + 1); // new ux at ly
493 1 : const auto & fQ32 = _z_surface(uy, lx + 1); // new ux at uy
494 :
495 1 : const Real & x3 = _x1[lx + 1]; // ux value
496 :
497 1 : auto dfdx = fQ11 * (y - y2);
498 1 : dfdx += fQ31 * (y2 - y);
499 1 : dfdx += fQ12 * (y1 - y);
500 1 : dfdx += fQ32 * (y - y1);
501 1 : dfdx /= ((x3 - x1) * (y2 - y1));
502 1 : return dfdx;
503 0 : }
504 : // right boundary
505 4 : else if (lx == ux && ux == nx1 - 1)
506 : {
507 1 : const auto & fQ01 = _z_surface(ly, ux - 1); // new lx at ly
508 1 : const auto & fQ02 = _z_surface(uy, ux - 1); // new lx at uy
509 1 : const Real & x0 = _x1[ux - 1]; // lx value
510 :
511 1 : auto dfdx = fQ01 * (y - y2);
512 1 : dfdx += fQ11 * (y2 - y);
513 1 : dfdx += fQ02 * (y1 - y);
514 1 : dfdx += fQ12 * (y - y1);
515 1 : dfdx /= ((x1 - x0) * (y2 - y1));
516 1 : return dfdx;
517 0 : }
518 : // Derivative (w/ respect to x) for some point inside box
519 : else
520 : {
521 3 : auto dfdx_xy = fQ11 * (y - y2);
522 3 : dfdx_xy += fQ21 * (y2 - y);
523 3 : dfdx_xy += fQ12 * (y1 - y);
524 3 : dfdx_xy += fQ22 * (y - y1);
525 3 : dfdx_xy /= ((x2 - x1) * (y2 - y1));
526 3 : return dfdx_xy;
527 0 : }
528 : }
529 :
530 9 : else if (deriv_var == 2)
531 : {
532 9 : if (x == x1) // if x equal to x1 node
533 3 : return (fQ12 - fQ11) / (y2 - y1);
534 : // interior grid line
535 6 : else if (ly == uy && ly > 0 && uy < nx2 - 1)
536 : {
537 : // expand grid size so x1 does not equal x2
538 1 : const auto & fQ10 = _z_surface(ly - 1, lx); // new ly at lx
539 1 : const auto & fQ20 = _z_surface(ly - 1, ux); // new uy at lx
540 1 : const auto & fQ13 = _z_surface(ly + 1, lx); // new ly at ux
541 1 : const auto & fQ23 = _z_surface(ly + 1, ux); // new uy at ux
542 1 : const Real & y0 = _x2[ly - 1];
543 1 : const Real & y3 = _x2[ly + 1];
544 :
545 1 : auto dfdy_a = fQ10 * (x - x2);
546 1 : dfdy_a += fQ20 * (x1 - x);
547 1 : dfdy_a += fQ11 * (x2 - x);
548 1 : dfdy_a += fQ21 * (x - x1);
549 1 : dfdy_a /= ((x2 - x1) * (y1 - y0));
550 :
551 1 : auto dfdy_b = fQ11 * (x - x2);
552 1 : dfdy_b += fQ21 * (x1 - x);
553 1 : dfdy_b += fQ13 * (x2 - x);
554 1 : dfdy_b += fQ23 * (x - x1);
555 1 : dfdy_b /= ((x2 - x1) * (y3 - y1));
556 1 : return 0.5 * (dfdy_a + dfdy_b);
557 0 : }
558 : // bottom boundary
559 5 : else if (ly == uy && ly == 0)
560 : {
561 1 : const auto & fQ13 = _z_surface(ly + 1, lx); // new uy at lx
562 1 : const auto & fQ23 = _z_surface(ly + 1, ux); // new uy at ux
563 :
564 1 : const Real & y3 = _x2[ly + 1]; // lx value
565 :
566 1 : auto dfdy = fQ11 * (x - x2);
567 1 : dfdy += fQ21 * (x1 - x);
568 1 : dfdy += fQ13 * (x2 - x);
569 1 : dfdy += fQ23 * (x - x1);
570 1 : dfdy /= ((x2 - x1) * (y3 - y1));
571 1 : return dfdy;
572 0 : }
573 : // top boundary
574 4 : else if (ly == uy && uy == nx2 - 1)
575 : {
576 1 : const auto & fQ10 = _z_surface(ly - 1, lx); // new ly at lx
577 1 : const auto & fQ20 = _z_surface(ly - 1, ux); // new ly at ux
578 1 : const Real & y0 = _x2[ly - 1]; // lx value
579 :
580 1 : auto dfdy = fQ10 * (x - x2);
581 1 : dfdy += fQ20 * (x1 - x);
582 1 : dfdy += fQ11 * (x2 - x);
583 1 : dfdy += fQ21 * (x - x1);
584 1 : dfdy /= ((x2 - x1) * (y1 - y0));
585 1 : return dfdy;
586 0 : }
587 : else
588 : {
589 : // Derivative (w/ respect to y) for any point inside box
590 3 : auto dfdy_xy = fQ11 * (x - x2);
591 3 : dfdy_xy += fQ21 * (x1 - x);
592 3 : dfdy_xy += fQ12 * (x2 - x);
593 3 : dfdy_xy += fQ22 * (x - x1);
594 3 : dfdy_xy /= ((x2 - x1) * (y2 - y1));
595 3 : return dfdy_xy;
596 0 : }
597 : }
598 0 : mooseError("Bilinear interpolation failed to select a case for x1= ", s1, " x2= ", s2);
599 : }
600 :
601 : void
602 1 : BilinearInterpolation::sampleValueAndDerivatives(
603 : Real s1, Real s2, Real & y, Real & dy_ds1, Real & dy_ds2) const
604 : {
605 1 : y = sample(s1, s2);
606 1 : dy_ds1 = sampleDerivative(s1, s2, 1);
607 1 : dy_ds2 = sampleDerivative(s1, s2, 2);
608 1 : }
|
