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