GteIntpAkimaUniform3.h
Go to the documentation of this file.
1 // David Eberly, Geometric Tools, Redmond WA 98052
2 // Copyright (c) 1998-2017
3 // Distributed under the Boost Software License, Version 1.0.
4 // http://www.boost.org/LICENSE_1_0.txt
5 // http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
6 // File Version: 3.0.0 (2016/06/19)
7 
8 #pragma once
9 
10 #include <LowLevel/GteArray3.h>
11 #include <LowLevel/GteLogger.h>
12 #include <algorithm>
13 #include <cmath>
14 #include <cstring>
15 
16 // The interpolator is for uniformly spaced(x,y z)-values. The input samples
17 // must be stored in lexicographical order to represent f(x,y,z); that is,
18 // F[c + xBound*(r + yBound*s)] corresponds to f(x,y,z), where c is the index
19 // corresponding to x, r is the index corresponding to y, and s is the index
20 // corresponding to z.
21 
22 namespace gte
23 {
24 
25 template <typename Real>
26 class IntpAkimaUniform3
27 {
28 public:
29  // Construction and destruction.
30  ~IntpAkimaUniform3();
31  IntpAkimaUniform3(int xBound, int yBound, int zBound, Real xMin,
32  Real xSpacing, Real yMin, Real ySpacing, Real zMin, Real zSpacing,
33  Real const* F);
34 
35  // Member access.
36  inline int GetXBound() const;
37  inline int GetYBound() const;
38  inline int GetZBound() const;
39  inline int GetQuantity() const;
40  inline Real const* GetF() const;
41  inline Real GetXMin() const;
42  inline Real GetXMax() const;
43  inline Real GetXSpacing() const;
44  inline Real GetYMin() const;
45  inline Real GetYMax() const;
46  inline Real GetYSpacing() const;
47  inline Real GetZMin() const;
48  inline Real GetZMax() const;
49  inline Real GetZSpacing() const;
50 
51  // Evaluate the function and its derivatives. The functions clamp the
52  // inputs to xmin <= x <= xmax, ymin <= y <= ymax, and zmin <= z <= zmax.
53  // The first operator is for function evaluation. The second operator is
54  // for function or derivative evaluations. The xOrder argument is the
55  // order of the x-derivative, the yOrder argument is the order of the
56  // y-derivative, and the zOrder argument is the order of the z-derivative.
57  // All orders are zero to get the function value itself.
58  Real operator()(Real x, Real y, Real z) const;
59  Real operator()(int xOrder, int yOrder, int zOrder, Real x, Real y, Real z) const;
60 
61 private:
62  class Polynomial
63  {
64  public:
65  Polynomial();
66 
67  // P(x,y,z) = sum_{i=0}^3 sum_{j=0}^3 sum_{k=0}^3 a_{ijk} x^i y^j z^k.
68  // The tensor term A[ix][iy][iz] corresponds to the polynomial term
69  // x^{ix} y^{iy} z^{iz}.
70  Real& A(int ix, int iy, int iz);
71  Real operator()(Real x, Real y, Real z) const;
72  Real operator()(int xOrder, int yOrder, int zOrder, Real x, Real y, Real z) const;
73 
74  private:
75  Real mCoeff[4][4][4];
76  };
77 
78  // Support for construction.
79  void GetFX(Array3<Real> const& F, Array3<Real>& FX);
80  void GetFY(Array3<Real> const& F, Array3<Real>& FY);
81  void GetFZ(Array3<Real> const& F, Array3<Real>& FZ);
82  void GetFXY(Array3<Real> const& F, Array3<Real>& FXY);
83  void GetFXZ(Array3<Real> const& F, Array3<Real>& FXZ);
84  void GetFYZ(Array3<Real> const& F, Array3<Real>& FYZ);
85  void GetFXYZ(Array3<Real> const& F, Array3<Real>& FXYZ);
86  void GetPolynomials(Array3<Real> const& F, Array3<Real> const& FX,
87  Array3<Real> const& FY, Array3<Real> const& FZ, Array3<Real> const& FXY,
88  Array3<Real> const& FXZ, Array3<Real> const& FYZ, Array3<Real> const& FXYZ);
89 
90  Real ComputeDerivative(Real const* slope) const;
91  void Construct(Polynomial& poly, Real const F[2][2][2],
92  Real const FX[2][2][2], Real const FY[2][2][2],
93  Real const FZ[2][2][2], Real const FXY[2][2][2],
94  Real const FXZ[2][2][2], Real const FYZ[2][2][2],
95  Real const FXYZ[2][2][2]);
96 
97  void XLookup(Real x, int& xIndex, Real& dx) const;
98  void YLookup(Real y, int& yIndex, Real& dy) const;
99  void ZLookup(Real z, int& zIndex, Real& dz) const;
100 
101  int mXBound, mYBound, mZBound, mQuantity;
102  Real mXMin, mXMax, mXSpacing;
103  Real mYMin, mYMax, mYSpacing;
104  Real mZMin, mZMax, mZSpacing;
105  Real const* mF;
106  Array3<Polynomial> mPoly;
107 };
108 
109 
110 template <typename Real>
111 IntpAkimaUniform3<Real>::~IntpAkimaUniform3()
112 {
113 }
114 
115 template <typename Real>
116 IntpAkimaUniform3<Real>::IntpAkimaUniform3(int xBound, int yBound, int zBound,
117  Real xMin, Real xSpacing, Real yMin, Real ySpacing, Real zMin,
118  Real zSpacing, Real const* F)
119  :
120  mXBound(xBound),
121  mYBound(yBound),
122  mZBound(zBound),
123  mQuantity(xBound * yBound * zBound),
124  mXMin(xMin),
125  mXSpacing(xSpacing),
126  mYMin(yMin),
127  mYSpacing(ySpacing),
128  mZMin(zMin),
129  mZSpacing(zSpacing),
130  mF(F),
131  mPoly(xBound - 1, yBound - 1, zBound - 1)
132 {
133  // At least a 3x3x3 block of data points is needed to construct the
134  // estimates of the boundary derivatives.
135  LogAssert(mXBound >= 3 && mYBound >= 3 && mZBound >= 3 && mF,
136  "Invalid input.");
137  LogAssert(mXSpacing > (Real)0 && mYSpacing > (Real)0 &&
138  mZSpacing > (Real)0, "Invalid input.");
139 
140  mXMax = mXMin + mXSpacing * static_cast<Real>(mXBound - 1);
141  mYMax = mYMin + mYSpacing * static_cast<Real>(mYBound - 1);
142  mZMax = mZMin + mZSpacing * static_cast<Real>(mZBound - 1);
143 
144  // Create a 3D wrapper for the 1D samples.
145  Array3<Real> Fmap(mXBound, mYBound, mZBound, const_cast<Real*>(mF));
146 
147  // Construct first-order derivatives.
148  Array3<Real> FX(mXBound, mYBound, mZBound);
149  Array3<Real> FY(mXBound, mYBound, mZBound);
150  Array3<Real> FZ(mXBound, mYBound, mZBound);
151  GetFX(Fmap, FX);
152  GetFX(Fmap, FY);
153  GetFX(Fmap, FZ);
154 
155  // Construct second-order derivatives.
156  Array3<Real> FXY(mXBound, mYBound, mZBound);
157  Array3<Real> FXZ(mXBound, mYBound, mZBound);
158  Array3<Real> FYZ(mXBound, mYBound, mZBound);
159  GetFX(Fmap, FXY);
160  GetFX(Fmap, FXZ);
161  GetFX(Fmap, FYZ);
162 
163  // Construct third-order derivatives.
164  Array3<Real> FXYZ(mXBound, mYBound, mZBound);
165  GetFXYZ(Fmap, FXYZ);
166 
167  // Construct polynomials.
168  GetPolynomials(Fmap, FX, FY, FZ, FXY, FXZ, FYZ, FXYZ);
169 }
170 
171 template <typename Real> inline
172 int IntpAkimaUniform3<Real>::GetXBound() const
173 {
174  return mXBound;
175 }
176 
177 template <typename Real> inline
178 int IntpAkimaUniform3<Real>::GetYBound() const
179 {
180  return mYBound;
181 }
182 
183 template <typename Real> inline
184 int IntpAkimaUniform3<Real>::GetZBound() const
185 {
186  return mZBound;
187 }
188 
189 template <typename Real> inline
190 int IntpAkimaUniform3<Real>::GetQuantity() const
191 {
192  return mQuantity;
193 }
194 
195 template <typename Real> inline
196 Real const* IntpAkimaUniform3<Real>::GetF() const
197 {
198  return mF;
199 }
200 
201 template <typename Real> inline
202 Real IntpAkimaUniform3<Real>::GetXMin() const
203 {
204  return mXMin;
205 }
206 
207 template <typename Real> inline
208 Real IntpAkimaUniform3<Real>::GetXMax() const
209 {
210  return mXMax;
211 }
212 
213 template <typename Real> inline
214 Real IntpAkimaUniform3<Real>::GetXSpacing() const
215 {
216  return mXSpacing;
217 }
218 
219 template <typename Real> inline
220 Real IntpAkimaUniform3<Real>::GetYMin() const
221 {
222  return mYMin;
223 }
224 
225 template <typename Real> inline
226 Real IntpAkimaUniform3<Real>::GetYMax() const
227 {
228  return mYMax;
229 }
230 
231 template <typename Real> inline
232 Real IntpAkimaUniform3<Real>::GetYSpacing() const
233 {
234  return mYSpacing;
235 }
236 
237 template <typename Real> inline
238 Real IntpAkimaUniform3<Real>::GetZMin() const
239 {
240  return mZMin;
241 }
242 
243 template <typename Real> inline
244 Real IntpAkimaUniform3<Real>::GetZMax() const
245 {
246  return mZMax;
247 }
248 
249 template <typename Real> inline
250 Real IntpAkimaUniform3<Real>::GetZSpacing() const
251 {
252  return mZSpacing;
253 }
254 
255 template <typename Real>
256 Real IntpAkimaUniform3<Real>::operator()(Real x, Real y, Real z) const
257 {
258  x = std::min(std::max(x, mXMin), mXMax);
259  y = std::min(std::max(y, mYMin), mYMax);
260  z = std::min(std::max(z, mZMin), mZMax);
261  int ix, iy, iz;
262  Real dx, dy, dz;
263  XLookup(x, ix, dx);
264  YLookup(y, iy, dy);
265  ZLookup(z, iz, dz);
266  return mPoly[iz][iy][ix](dx, dy, dz);
267 }
268 
269 template <typename Real>
270 Real IntpAkimaUniform3<Real>::operator()(int xOrder, int yOrder, int zOrder,
271  Real x, Real y, Real z) const
272 {
273  x = std::min(std::max(x, mXMin), mXMax);
274  y = std::min(std::max(y, mYMin), mYMax);
275  z = std::min(std::max(z, mZMin), mZMax);
276  int ix, iy, iz;
277  Real dx, dy, dz;
278  XLookup(x, ix, dx);
279  YLookup(y, iy, dy);
280  ZLookup(z, iz, dz);
281  return mPoly[iz][iy][ix](xOrder, yOrder, zOrder, dx, dy, dz);
282 }
283 
284 template <typename Real>
285 void IntpAkimaUniform3<Real>::GetFX(Array3<Real> const& F, Array3<Real>& FX)
286 {
287  Array3<Real> slope(mXBound + 3, mYBound, mZBound);
288  Real invDX = ((Real)1) / mXSpacing;
289  int ix, iy, iz;
290  for (iz = 0; iz < mZBound; ++iz)
291  {
292  for (iy = 0; iy < mYBound; ++iy)
293  {
294  for (ix = 0; ix < mXBound - 1; ++ix)
295  {
296  slope[iz][iy][ix + 2] = (F[iz][iy][ix + 1] - F[iz][iy][ix]) * invDX;
297  }
298 
299  slope[iz][iy][1] = ((Real)2) * slope[iz][iy][2] - slope[iz][iy][3];
300  slope[iz][iy][0] = ((Real)2) * slope[iz][iy][1] - slope[iz][iy][2];
301  slope[iz][iy][mXBound + 1] = ((Real)2) * slope[iz][iy][mXBound] - slope[iz][iy][mXBound - 1];
302  slope[iz][iy][mXBound + 2] = ((Real)2) * slope[iz][iy][mXBound + 1] - slope[iz][iy][mXBound];
303  }
304  }
305 
306  for (iz = 0; iz < mZBound; ++iz)
307  {
308  for (iy = 0; iy < mYBound; ++iy)
309  {
310  for (ix = 0; ix < mXBound; ++ix)
311  {
312  FX[iz][iy][ix] = ComputeDerivative(slope[iz][iy] + ix);
313  }
314  }
315  }
316 }
317 
318 template <typename Real>
319 void IntpAkimaUniform3<Real>::GetFY(Array3<Real> const& F, Array3<Real>& FY)
320 {
321  Array3<Real> slope(mYBound + 3, mXBound, mZBound);
322  Real invDY = ((Real)1) / mYSpacing;
323  int ix, iy, iz;
324  for (iz = 0; iz < mZBound; ++iz)
325  {
326  for (ix = 0; ix < mXBound; ++ix)
327  {
328  for (iy = 0; iy < mYBound - 1; ++iy)
329  {
330  slope[iz][ix][iy + 2] = (F[iz][iy + 1][ix] - F[iz][iy][ix]) * invDY;
331  }
332 
333  slope[iz][ix][1] = ((Real)2) * slope[iz][ix][2] - slope[iz][ix][3];
334  slope[iz][ix][0] = ((Real)2) * slope[iz][ix][1] - slope[iz][ix][2];
335  slope[iz][ix][mYBound + 1] = ((Real)2) * slope[iz][ix][mYBound] - slope[iz][ix][mYBound - 1];
336  slope[iz][ix][mYBound + 2] = ((Real)2) * slope[iz][ix][mYBound + 1] - slope[iz][ix][mYBound];
337  }
338  }
339 
340  for (iz = 0; iz < mZBound; ++iz)
341  {
342  for (ix = 0; ix < mXBound; ++ix)
343  {
344  for (iy = 0; iy < mYBound; ++iy)
345  {
346  FY[iz][iy][ix] = ComputeDerivative(slope[iz][ix] + iy);
347  }
348  }
349  }
350 }
351 
352 template <typename Real>
353 void IntpAkimaUniform3<Real>::GetFZ(Array3<Real> const& F, Array3<Real>& FZ)
354 {
355  Array3<Real> slope(mZBound + 3, mXBound, mYBound);
356  Real invDZ = ((Real)1) / mZSpacing;
357  int ix, iy, iz;
358  for (iy = 0; iy < mYBound; ++iy)
359  {
360  for (ix = 0; ix < mXBound; ++ix)
361  {
362  for (iz = 0; iz < mZBound - 1; ++iz)
363  {
364  slope[iy][ix][iz + 2] = (F[iz + 1][iy][ix] - F[iz][iy][ix]) * invDZ;
365  }
366 
367  slope[iy][ix][1] = ((Real)2) * slope[iy][ix][2] - slope[iy][ix][3];
368  slope[iy][ix][0] = ((Real)2) * slope[iy][ix][1] - slope[iy][ix][2];
369  slope[iy][ix][mZBound + 1] = ((Real)2) * slope[iy][ix][mZBound] - slope[iy][ix][mZBound - 1];
370  slope[iy][ix][mZBound + 2] = ((Real)2) * slope[iy][ix][mZBound + 1] - slope[iy][ix][mZBound];
371  }
372  }
373 
374  for (iy = 0; iy < mYBound; ++iy)
375  {
376  for (ix = 0; ix < mXBound; ++ix)
377  {
378  for (iz = 0; iz < mZBound; ++iz)
379  {
380  FZ[iz][iy][ix] = ComputeDerivative(slope[iy][ix] + iz);
381  }
382  }
383  }
384 }
385 
386 template <typename Real>
387 void IntpAkimaUniform3<Real>::GetFXY(Array3<Real> const& F, Array3<Real>& FXY)
388 {
389  int xBoundM1 = mXBound - 1;
390  int yBoundM1 = mYBound - 1;
391  int ix0 = xBoundM1, ix1 = ix0 - 1, ix2 = ix1 - 1;
392  int iy0 = yBoundM1, iy1 = iy0 - 1, iy2 = iy1 - 1;
393  int ix, iy, iz;
394 
395  Real invDXDY = ((Real)1) / (mXSpacing * mYSpacing);
396  for (iz = 0; iz < mZBound; ++iz)
397  {
398  // corners of z-slice
399  FXY[iz][0][0] = ((Real)0.25)*invDXDY*(
400  ((Real)9)*F[iz][0][0]
401  - ((Real)12)*F[iz][0][1]
402  + ((Real)3)*F[iz][0][2]
403  - ((Real)12)*F[iz][1][0]
404  + ((Real)16)*F[iz][1][1]
405  - ((Real)4)*F[iz][1][2]
406  + ((Real)3)*F[iz][2][0]
407  - ((Real)4)*F[iz][2][1]
408  + F[iz][2][2]);
409 
410  FXY[iz][0][xBoundM1] = ((Real)0.25)*invDXDY*(
411  ((Real)9)*F[iz][0][ix0]
412  - ((Real)12)*F[iz][0][ix1]
413  + ((Real)3)*F[iz][0][ix2]
414  - ((Real)12)*F[iz][1][ix0]
415  + ((Real)16)*F[iz][1][ix1]
416  - ((Real)4)*F[iz][1][ix2]
417  + ((Real)3)*F[iz][2][ix0]
418  - ((Real)4)*F[iz][2][ix1]
419  + F[iz][2][ix2]);
420 
421  FXY[iz][yBoundM1][0] = ((Real)0.25)*invDXDY*(
422  ((Real)9)*F[iz][iy0][0]
423  - ((Real)12)*F[iz][iy0][1]
424  + ((Real)3)*F[iz][iy0][2]
425  - ((Real)12)*F[iz][iy1][0]
426  + ((Real)16)*F[iz][iy1][1]
427  - ((Real)4)*F[iz][iy1][2]
428  + ((Real)3)*F[iz][iy2][0]
429  - ((Real)4)*F[iz][iy2][1]
430  + F[iz][iy2][2]);
431 
432  FXY[iz][yBoundM1][xBoundM1] = ((Real)0.25)*invDXDY*(
433  ((Real)9)*F[iz][iy0][ix0]
434  - ((Real)12)*F[iz][iy0][ix1]
435  + ((Real)3)*F[iz][iy0][ix2]
436  - ((Real)12)*F[iz][iy1][ix0]
437  + ((Real)16)*F[iz][iy1][ix1]
438  - ((Real)4)*F[iz][iy1][ix2]
439  + ((Real)3)*F[iz][iy2][ix0]
440  - ((Real)4)*F[iz][iy2][ix1]
441  + F[iz][iy2][ix2]);
442 
443  // x-edges of z-slice
444  for (ix = 1; ix < xBoundM1; ++ix)
445  {
446  FXY[iz][0][ix] = ((Real)0.25)*invDXDY*(
447  ((Real)3)*(F[iz][0][ix - 1] - F[iz][0][ix + 1]) -
448  ((Real)4)*(F[iz][1][ix - 1] - F[iz][1][ix + 1]) +
449  (F[iz][2][ix - 1] - F[iz][2][ix + 1]));
450 
451  FXY[iz][yBoundM1][ix] = ((Real)0.25)*invDXDY*(
452  ((Real)3)*(F[iz][iy0][ix - 1] - F[iz][iy0][ix + 1])
453  - ((Real)4)*(F[iz][iy1][ix - 1] - F[iz][iy1][ix + 1]) +
454  (F[iz][iy2][ix - 1] - F[iz][iy2][ix + 1]));
455  }
456 
457  // y-edges of z-slice
458  for (iy = 1; iy < yBoundM1; ++iy)
459  {
460  FXY[iz][iy][0] = ((Real)0.25)*invDXDY*(
461  ((Real)3)*(F[iz][iy - 1][0] - F[iz][iy + 1][0]) -
462  ((Real)4)*(F[iz][iy - 1][1] - F[iz][iy + 1][1]) +
463  (F[iz][iy - 1][2] - F[iz][iy + 1][2]));
464 
465  FXY[iz][iy][xBoundM1] = ((Real)0.25)*invDXDY*(
466  ((Real)3)*(F[iz][iy - 1][ix0] - F[iz][iy + 1][ix0])
467  - ((Real)4)*(F[iz][iy - 1][ix1] - F[iz][iy + 1][ix1]) +
468  (F[iz][iy - 1][ix2] - F[iz][iy + 1][ix2]));
469  }
470 
471  // interior of z-slice
472  for (iy = 1; iy < yBoundM1; ++iy)
473  {
474  for (ix = 1; ix < xBoundM1; ++ix)
475  {
476  FXY[iz][iy][ix] = ((Real)0.25)*invDXDY*(
477  F[iz][iy - 1][ix - 1] - F[iz][iy - 1][ix + 1] -
478  F[iz][iy + 1][ix - 1] + F[iz][iy + 1][ix + 1]);
479  }
480  }
481  }
482 }
483 
484 template <typename Real>
485 void IntpAkimaUniform3<Real>::GetFXZ(Array3<Real> const& F, Array3<Real>& FXZ)
486 {
487  int xBoundM1 = mXBound - 1;
488  int zBoundM1 = mZBound - 1;
489  int ix0 = xBoundM1, ix1 = ix0 - 1, ix2 = ix1 - 1;
490  int iz0 = zBoundM1, iz1 = iz0 - 1, iz2 = iz1 - 1;
491  int ix, iy, iz;
492 
493  Real invDXDZ = ((Real)1) / (mXSpacing * mZSpacing);
494  for (iy = 0; iy < mYBound; ++iy)
495  {
496  // corners of z-slice
497  FXZ[0][iy][0] = ((Real)0.25)*invDXDZ*(
498  ((Real)9)*F[0][iy][0]
499  - ((Real)12)*F[0][iy][1]
500  + ((Real)3)*F[0][iy][2]
501  - ((Real)12)*F[1][iy][0]
502  + ((Real)16)*F[1][iy][1]
503  - ((Real)4)*F[1][iy][2]
504  + ((Real)3)*F[2][iy][0]
505  - ((Real)4)*F[2][iy][1]
506  + F[2][iy][2]);
507 
508  FXZ[0][iy][xBoundM1] = ((Real)0.25)*invDXDZ*(
509  ((Real)9)*F[0][iy][ix0]
510  - ((Real)12)*F[0][iy][ix1]
511  + ((Real)3)*F[0][iy][ix2]
512  - ((Real)12)*F[1][iy][ix0]
513  + ((Real)16)*F[1][iy][ix1]
514  - ((Real)4)*F[1][iy][ix2]
515  + ((Real)3)*F[2][iy][ix0]
516  - ((Real)4)*F[2][iy][ix1]
517  + F[2][iy][ix2]);
518 
519  FXZ[zBoundM1][iy][0] = ((Real)0.25)*invDXDZ*(
520  ((Real)9)*F[iz0][iy][0]
521  - ((Real)12)*F[iz0][iy][1]
522  + ((Real)3)*F[iz0][iy][2]
523  - ((Real)12)*F[iz1][iy][0]
524  + ((Real)16)*F[iz1][iy][1]
525  - ((Real)4)*F[iz1][iy][2]
526  + ((Real)3)*F[iz2][iy][0]
527  - ((Real)4)*F[iz2][iy][1]
528  + F[iz2][iy][2]);
529 
530  FXZ[zBoundM1][iy][xBoundM1] = ((Real)0.25)*invDXDZ*(
531  ((Real)9)*F[iz0][iy][ix0]
532  - ((Real)12)*F[iz0][iy][ix1]
533  + ((Real)3)*F[iz0][iy][ix2]
534  - ((Real)12)*F[iz1][iy][ix0]
535  + ((Real)16)*F[iz1][iy][ix1]
536  - ((Real)4)*F[iz1][iy][ix2]
537  + ((Real)3)*F[iz2][iy][ix0]
538  - ((Real)4)*F[iz2][iy][ix1]
539  + F[iz2][iy][ix2]);
540 
541  // x-edges of y-slice
542  for (ix = 1; ix < xBoundM1; ++ix)
543  {
544  FXZ[0][iy][ix] = ((Real)0.25)*invDXDZ*(
545  ((Real)3)*(F[0][iy][ix - 1] - F[0][iy][ix + 1]) -
546  ((Real)4)*(F[1][iy][ix - 1] - F[1][iy][ix + 1]) +
547  (F[2][iy][ix - 1] - F[2][iy][ix + 1]));
548 
549  FXZ[zBoundM1][iy][ix] = ((Real)0.25)*invDXDZ*(
550  ((Real)3)*(F[iz0][iy][ix - 1] - F[iz0][iy][ix + 1])
551  - ((Real)4)*(F[iz1][iy][ix - 1] - F[iz1][iy][ix + 1]) +
552  (F[iz2][iy][ix - 1] - F[iz2][iy][ix + 1]));
553  }
554 
555  // z-edges of y-slice
556  for (iz = 1; iz < zBoundM1; ++iz)
557  {
558  FXZ[iz][iy][0] = ((Real)0.25)*invDXDZ*(
559  ((Real)3)*(F[iz - 1][iy][0] - F[iz + 1][iy][0]) -
560  ((Real)4)*(F[iz - 1][iy][1] - F[iz + 1][iy][1]) +
561  (F[iz - 1][iy][2] - F[iz + 1][iy][2]));
562 
563  FXZ[iz][iy][xBoundM1] = ((Real)0.25)*invDXDZ*(
564  ((Real)3)*(F[iz - 1][iy][ix0] - F[iz + 1][iy][ix0])
565  - ((Real)4)*(F[iz - 1][iy][ix1] - F[iz + 1][iy][ix1]) +
566  (F[iz - 1][iy][ix2] - F[iz + 1][iy][ix2]));
567  }
568 
569  // interior of y-slice
570  for (iz = 1; iz < zBoundM1; ++iz)
571  {
572  for (ix = 1; ix < xBoundM1; ++ix)
573  {
574  FXZ[iz][iy][ix] = ((Real)0.25)*invDXDZ*(
575  F[iz - 1][iy][ix - 1] - F[iz - 1][iy][ix + 1] -
576  F[iz + 1][iy][ix - 1] + F[iz + 1][iy][ix + 1]);
577  }
578  }
579  }
580 }
581 
582 template <typename Real>
583 void IntpAkimaUniform3<Real>::GetFYZ(Array3<Real> const& F, Array3<Real>& FYZ)
584 {
585  int yBoundM1 = mYBound - 1;
586  int zBoundM1 = mZBound - 1;
587  int iy0 = yBoundM1, iy1 = iy0 - 1, iy2 = iy1 - 1;
588  int iz0 = zBoundM1, iz1 = iz0 - 1, iz2 = iz1 - 1;
589  int ix, iy, iz;
590 
591  Real invDYDZ = ((Real)1) / (mYSpacing * mZSpacing);
592  for (ix = 0; ix < mXBound; ++ix)
593  {
594  // corners of x-slice
595  FYZ[0][0][ix] = ((Real)0.25)*invDYDZ*(
596  ((Real)9)*F[0][0][ix]
597  - ((Real)12)*F[0][1][ix]
598  + ((Real)3)*F[0][2][ix]
599  - ((Real)12)*F[1][0][ix]
600  + ((Real)16)*F[1][1][ix]
601  - ((Real)4)*F[1][2][ix]
602  + ((Real)3)*F[2][0][ix]
603  - ((Real)4)*F[2][1][ix]
604  + F[2][2][ix]);
605 
606  FYZ[0][yBoundM1][ix] = ((Real)0.25)*invDYDZ*(
607  ((Real)9)*F[0][iy0][ix]
608  - ((Real)12)*F[0][iy1][ix]
609  + ((Real)3)*F[0][iy2][ix]
610  - ((Real)12)*F[1][iy0][ix]
611  + ((Real)16)*F[1][iy1][ix]
612  - ((Real)4)*F[1][iy2][ix]
613  + ((Real)3)*F[2][iy0][ix]
614  - ((Real)4)*F[2][iy1][ix]
615  + F[2][iy2][ix]);
616 
617  FYZ[zBoundM1][0][ix] = ((Real)0.25)*invDYDZ*(
618  ((Real)9)*F[iz0][0][ix]
619  - ((Real)12)*F[iz0][1][ix]
620  + ((Real)3)*F[iz0][2][ix]
621  - ((Real)12)*F[iz1][0][ix]
622  + ((Real)16)*F[iz1][1][ix]
623  - ((Real)4)*F[iz1][2][ix]
624  + ((Real)3)*F[iz2][0][ix]
625  - ((Real)4)*F[iz2][1][ix]
626  + F[iz2][2][ix]);
627 
628  FYZ[zBoundM1][yBoundM1][ix] = ((Real)0.25)*invDYDZ*(
629  ((Real)9)*F[iz0][iy0][ix]
630  - ((Real)12)*F[iz0][iy1][ix]
631  + ((Real)3)*F[iz0][iy2][ix]
632  - ((Real)12)*F[iz1][iy0][ix]
633  + ((Real)16)*F[iz1][iy1][ix]
634  - ((Real)4)*F[iz1][iy2][ix]
635  + ((Real)3)*F[iz2][iy0][ix]
636  - ((Real)4)*F[iz2][iy1][ix]
637  + F[iz2][iy2][ix]);
638 
639  // y-edges of x-slice
640  for (iy = 1; iy < yBoundM1; ++iy)
641  {
642  FYZ[0][iy][ix] = ((Real)0.25)*invDYDZ*(
643  ((Real)3)*(F[0][iy - 1][ix] - F[0][iy + 1][ix]) -
644  ((Real)4)*(F[1][iy - 1][ix] - F[1][iy + 1][ix]) +
645  (F[2][iy - 1][ix] - F[2][iy + 1][ix]));
646 
647  FYZ[zBoundM1][iy][ix] = ((Real)0.25)*invDYDZ*(
648  ((Real)3)*(F[iz0][iy - 1][ix] - F[iz0][iy + 1][ix])
649  - ((Real)4)*(F[iz1][iy - 1][ix] - F[iz1][iy + 1][ix]) +
650  (F[iz2][iy - 1][ix] - F[iz2][iy + 1][ix]));
651  }
652 
653  // z-edges of x-slice
654  for (iz = 1; iz < zBoundM1; ++iz)
655  {
656  FYZ[iz][0][ix] = ((Real)0.25)*invDYDZ*(
657  ((Real)3)*(F[iz - 1][0][ix] - F[iz + 1][0][ix]) -
658  ((Real)4)*(F[iz - 1][1][ix] - F[iz + 1][1][ix]) +
659  (F[iz - 1][2][ix] - F[iz + 1][2][ix]));
660 
661  FYZ[iz][yBoundM1][ix] = ((Real)0.25)*invDYDZ*(
662  ((Real)3)*(F[iz - 1][iy0][ix] - F[iz + 1][iy0][ix])
663  - ((Real)4)*(F[iz - 1][iy1][ix] - F[iz + 1][iy1][ix]) +
664  (F[iz - 1][iy2][ix] - F[iz + 1][iy2][ix]));
665  }
666 
667  // interior of x-slice
668  for (iz = 1; iz < zBoundM1; ++iz)
669  {
670  for (iy = 1; iy < yBoundM1; ++iy)
671  {
672  FYZ[iz][iy][ix] = ((Real)0.25)*invDYDZ*(
673  F[iz - 1][iy - 1][ix] - F[iz - 1][iy + 1][ix] -
674  F[iz + 1][iy - 1][ix] + F[iz + 1][iy + 1][ix]);
675  }
676  }
677  }
678 }
679 
680 template <typename Real>
681 void IntpAkimaUniform3<Real>::GetFXYZ(Array3<Real> const& F, Array3<Real>& FXYZ)
682 {
683  int xBoundM1 = mXBound - 1;
684  int yBoundM1 = mYBound - 1;
685  int zBoundM1 = mZBound - 1;
686  int ix, iy, iz, ix0, iy0, iz0;
687 
688  Real invDXDYDZ = ((Real)1) / (mXSpacing * mYSpacing * mZSpacing);
689 
690  // convolution masks
691  // centered difference, O(h^2)
692  Real CDer[3] = { -(Real)0.5, (Real)0, (Real)0.5 };
693  // one-sided difference, O(h^2)
694  Real ODer[3] = { -(Real)1.5, (Real)2, -(Real)0.5 };
695  Real mask;
696 
697  // corners
698  FXYZ[0][0][0] = (Real)0;
699  FXYZ[0][0][xBoundM1] = (Real)0;
700  FXYZ[0][yBoundM1][0] = (Real)0;
701  FXYZ[0][yBoundM1][xBoundM1] = (Real)0;
702  FXYZ[zBoundM1][0][0] = (Real)0;
703  FXYZ[zBoundM1][0][xBoundM1] = (Real)0;
704  FXYZ[zBoundM1][yBoundM1][0] = (Real)0;
705  FXYZ[zBoundM1][yBoundM1][xBoundM1] = (Real)0;
706  for (iz = 0; iz <= 2; ++iz)
707  {
708  for (iy = 0; iy <= 2; ++iy)
709  {
710  for (ix = 0; ix <= 2; ++ix)
711  {
712  mask = invDXDYDZ*ODer[ix] * ODer[iy] * ODer[iz];
713  FXYZ[0][0][0] += mask * F[iz][iy][ix];
714  FXYZ[0][0][xBoundM1] += mask * F[iz][iy][xBoundM1 - ix];
715  FXYZ[0][yBoundM1][0] += mask * F[iz][yBoundM1 - iy][ix];
716  FXYZ[0][yBoundM1][xBoundM1] += mask * F[iz][yBoundM1 - iy][xBoundM1 - ix];
717  FXYZ[zBoundM1][0][0] += mask * F[zBoundM1 - iz][iy][ix];
718  FXYZ[zBoundM1][0][xBoundM1] += mask * F[zBoundM1 - iz][iy][xBoundM1 - ix];
719  FXYZ[zBoundM1][yBoundM1][0] += mask * F[zBoundM1 - iz][yBoundM1 - iy][ix];
720  FXYZ[zBoundM1][yBoundM1][xBoundM1] += mask * F[zBoundM1 - iz][yBoundM1 - iy][xBoundM1 - ix];
721  }
722  }
723  }
724 
725  // x-edges
726  for (ix0 = 1; ix0 < xBoundM1; ++ix0)
727  {
728  FXYZ[0][0][ix0] = (Real)0;
729  FXYZ[0][yBoundM1][ix0] = (Real)0;
730  FXYZ[zBoundM1][0][ix0] = (Real)0;
731  FXYZ[zBoundM1][yBoundM1][ix0] = (Real)0;
732  for (iz = 0; iz <= 2; ++iz)
733  {
734  for (iy = 0; iy <= 2; ++iy)
735  {
736  for (ix = 0; ix <= 2; ++ix)
737  {
738  mask = invDXDYDZ*CDer[ix] * ODer[iy] * ODer[iz];
739  FXYZ[0][0][ix0] += mask * F[iz][iy][ix0 + ix - 1];
740  FXYZ[0][yBoundM1][ix0] += mask * F[iz][yBoundM1 - iy][ix0 + ix - 1];
741  FXYZ[zBoundM1][0][ix0] += mask * F[zBoundM1 - iz][iy][ix0 + ix - 1];
742  FXYZ[zBoundM1][yBoundM1][ix0] += mask * F[zBoundM1 - iz][yBoundM1 - iy][ix0 + ix - 1];
743  }
744  }
745  }
746  }
747 
748  // y-edges
749  for (iy0 = 1; iy0 < yBoundM1; ++iy0)
750  {
751  FXYZ[0][iy0][0] = (Real)0;
752  FXYZ[0][iy0][xBoundM1] = (Real)0;
753  FXYZ[zBoundM1][iy0][0] = (Real)0;
754  FXYZ[zBoundM1][iy0][xBoundM1] = (Real)0;
755  for (iz = 0; iz <= 2; ++iz)
756  {
757  for (iy = 0; iy <= 2; ++iy)
758  {
759  for (ix = 0; ix <= 2; ++ix)
760  {
761  mask = invDXDYDZ*ODer[ix] * CDer[iy] * ODer[iz];
762  FXYZ[0][iy0][0] += mask * F[iz][iy0 + iy - 1][ix];
763  FXYZ[0][iy0][xBoundM1] += mask * F[iz][iy0 + iy - 1][xBoundM1 - ix];
764  FXYZ[zBoundM1][iy0][0] += mask * F[zBoundM1 - iz][iy0 + iy - 1][ix];
765  FXYZ[zBoundM1][iy0][xBoundM1] += mask * F[zBoundM1 - iz][iy0 + iy - 1][xBoundM1 - ix];
766  }
767  }
768  }
769  }
770 
771  // z-edges
772  for (iz0 = 1; iz0 < zBoundM1; ++iz0)
773  {
774  FXYZ[iz0][0][0] = (Real)0;
775  FXYZ[iz0][0][xBoundM1] = (Real)0;
776  FXYZ[iz0][yBoundM1][0] = (Real)0;
777  FXYZ[iz0][yBoundM1][xBoundM1] = (Real)0;
778  for (iz = 0; iz <= 2; ++iz)
779  {
780  for (iy = 0; iy <= 2; ++iy)
781  {
782  for (ix = 0; ix <= 2; ++ix)
783  {
784  mask = invDXDYDZ*ODer[ix] * ODer[iy] * CDer[iz];
785  FXYZ[iz0][0][0] += mask * F[iz0 + iz - 1][iy][ix];
786  FXYZ[iz0][0][xBoundM1] += mask * F[iz0 + iz - 1][iy][xBoundM1 - ix];
787  FXYZ[iz0][yBoundM1][0] += mask * F[iz0 + iz - 1][yBoundM1 - iy][ix];
788  FXYZ[iz0][yBoundM1][xBoundM1] += mask * F[iz0 + iz - 1][yBoundM1 - iy][xBoundM1 - ix];
789  }
790  }
791  }
792  }
793 
794  // xy-faces
795  for (iy0 = 1; iy0 < yBoundM1; ++iy0)
796  {
797  for (ix0 = 1; ix0 < xBoundM1; ++ix0)
798  {
799  FXYZ[0][iy0][ix0] = (Real)0;
800  FXYZ[zBoundM1][iy0][ix0] = (Real)0;
801  for (iz = 0; iz <= 2; ++iz)
802  {
803  for (iy = 0; iy <= 2; ++iy)
804  {
805  for (ix = 0; ix <= 2; ++ix)
806  {
807  mask = invDXDYDZ*CDer[ix] * CDer[iy] * ODer[iz];
808  FXYZ[0][iy0][ix0] += mask * F[iz][iy0 + iy - 1][ix0 + ix - 1];
809  FXYZ[zBoundM1][iy0][ix0] += mask * F[zBoundM1 - iz][iy0 + iy - 1][ix0 + ix - 1];
810  }
811  }
812  }
813  }
814  }
815 
816  // xz-faces
817  for (iz0 = 1; iz0 < zBoundM1; ++iz0)
818  {
819  for (ix0 = 1; ix0 < xBoundM1; ++ix0)
820  {
821  FXYZ[iz0][0][ix0] = (Real)0;
822  FXYZ[iz0][yBoundM1][ix0] = (Real)0;
823  for (iz = 0; iz <= 2; ++iz)
824  {
825  for (iy = 0; iy <= 2; ++iy)
826  {
827  for (ix = 0; ix <= 2; ++ix)
828  {
829  mask = invDXDYDZ*CDer[ix] * ODer[iy] * CDer[iz];
830  FXYZ[iz0][0][ix0] += mask * F[iz0 + iz - 1][iy][ix0 + ix - 1];
831  FXYZ[iz0][yBoundM1][ix0] += mask * F[iz0 + iz - 1][yBoundM1 - iy][ix0 + ix - 1];
832  }
833  }
834  }
835  }
836  }
837 
838  // yz-faces
839  for (iz0 = 1; iz0 < zBoundM1; ++iz0)
840  {
841  for (iy0 = 1; iy0 < yBoundM1; ++iy0)
842  {
843  FXYZ[iz0][iy0][0] = (Real)0;
844  FXYZ[iz0][iy0][xBoundM1] = (Real)0;
845  for (iz = 0; iz <= 2; ++iz)
846  {
847  for (iy = 0; iy <= 2; ++iy)
848  {
849  for (ix = 0; ix <= 2; ++ix)
850  {
851  mask = invDXDYDZ*ODer[ix] * CDer[iy] * CDer[iz];
852  FXYZ[iz0][iy0][0] += mask * F[iz0 + iz - 1][iy0 + iy - 1][ix];
853  FXYZ[iz0][iy0][xBoundM1] += mask * F[iz0 + iz - 1][iy0 + iy - 1][xBoundM1 - ix];
854  }
855  }
856  }
857  }
858  }
859 
860  // interiors
861  for (iz0 = 1; iz0 < zBoundM1; ++iz0)
862  {
863  for (iy0 = 1; iy0 < yBoundM1; ++iy0)
864  {
865  for (ix0 = 1; ix0 < xBoundM1; ++ix0)
866  {
867  FXYZ[iz0][iy0][ix0] = (Real)0;
868 
869  for (iz = 0; iz <= 2; ++iz)
870  {
871  for (iy = 0; iy <= 2; ++iy)
872  {
873  for (ix = 0; ix <= 2; ++ix)
874  {
875  mask = invDXDYDZ*CDer[ix] * CDer[iy] * CDer[iz];
876  FXYZ[iz0][iy0][ix0] += mask * F[iz0 + iz - 1][iy0 + iy - 1][ix0 + ix - 1];
877  }
878  }
879  }
880  }
881  }
882  }
883 }
884 
885 template <typename Real>
886 void IntpAkimaUniform3<Real>::GetPolynomials(Array3<Real> const& F, Array3<Real> const& FX,
887  Array3<Real> const& FY, Array3<Real> const& FZ, Array3<Real> const& FXY,
888  Array3<Real> const& FXZ, Array3<Real> const& FYZ, Array3<Real> const& FXYZ)
889 {
890  int xBoundM1 = mXBound - 1;
891  int yBoundM1 = mYBound - 1;
892  int zBoundM1 = mZBound - 1;
893  for (int iz = 0; iz < zBoundM1; ++iz)
894  {
895  for (int iy = 0; iy < yBoundM1; ++iy)
896  {
897  for (int ix = 0; ix < xBoundM1; ++ix)
898  {
899  // Note the 'transposing' of the 2x2x2 blocks (to match
900  // notation used in the polynomial definition).
901  Real G[2][2][2] =
902  {
903  {
904  {
905  F[iz][iy][ix],
906  F[iz + 1][iy][ix]
907  },
908  {
909  F[iz][iy + 1][ix],
910  F[iz + 1][iy + 1][ix]
911  }
912  },
913  {
914  {
915  F[iz][iy][ix + 1],
916  F[iz + 1][iy][ix + 1]
917  },
918  {
919  F[iz][iy + 1][ix + 1],
920  F[iz + 1][iy + 1][ix + 1]
921  }
922  }
923  };
924 
925  Real GX[2][2][2] =
926  {
927  {
928  {
929  FX[iz][iy][ix],
930  FX[iz + 1][iy][ix]
931  },
932  {
933  FX[iz][iy + 1][ix],
934  FX[iz + 1][iy + 1][ix]
935  }
936  },
937  {
938  {
939  FX[iz][iy][ix + 1],
940  FX[iz + 1][iy][ix + 1]
941  },
942  {
943  FX[iz][iy + 1][ix + 1],
944  FX[iz + 1][iy + 1][ix + 1]
945  }
946  }
947  };
948 
949  Real GY[2][2][2] =
950  {
951  {
952  {
953  FY[iz][iy][ix],
954  FY[iz + 1][iy][ix]
955  },
956  {
957  FY[iz][iy + 1][ix],
958  FY[iz + 1][iy + 1][ix]
959  }
960  },
961  {
962  {
963  FY[iz][iy][ix + 1],
964  FY[iz + 1][iy][ix + 1]
965  },
966  {
967  FY[iz][iy + 1][ix + 1],
968  FY[iz + 1][iy + 1][ix + 1]
969  }
970  }
971  };
972 
973  Real GZ[2][2][2] =
974  {
975  {
976  {
977  FZ[iz][iy][ix],
978  FZ[iz + 1][iy][ix]
979  },
980  {
981  FZ[iz][iy + 1][ix],
982  FZ[iz + 1][iy + 1][ix]
983  }
984  },
985  {
986  {
987  FZ[iz][iy][ix + 1],
988  FZ[iz + 1][iy][ix + 1]
989  },
990  {
991  FZ[iz][iy + 1][ix + 1],
992  FZ[iz + 1][iy + 1][ix + 1]
993  }
994  }
995  };
996 
997  Real GXY[2][2][2] =
998  {
999  {
1000  {
1001  FXY[iz][iy][ix],
1002  FXY[iz + 1][iy][ix]
1003  },
1004  {
1005  FXY[iz][iy + 1][ix],
1006  FXY[iz + 1][iy + 1][ix]
1007  }
1008  },
1009  {
1010  {
1011  FXY[iz][iy][ix + 1],
1012  FXY[iz + 1][iy][ix + 1]
1013  },
1014  {
1015  FXY[iz][iy + 1][ix + 1],
1016  FXY[iz + 1][iy + 1][ix + 1]
1017  }
1018  }
1019  };
1020 
1021  Real GXZ[2][2][2] =
1022  {
1023  {
1024  {
1025  FXZ[iz][iy][ix],
1026  FXZ[iz + 1][iy][ix]
1027  },
1028  {
1029  FXZ[iz][iy + 1][ix],
1030  FXZ[iz + 1][iy + 1][ix]
1031  }
1032  },
1033  {
1034  {
1035  FXZ[iz][iy][ix + 1],
1036  FXZ[iz + 1][iy][ix + 1]
1037  },
1038  {
1039  FXZ[iz][iy + 1][ix + 1],
1040  FXZ[iz + 1][iy + 1][ix + 1]
1041  }
1042  }
1043  };
1044 
1045  Real GYZ[2][2][2] =
1046  {
1047  {
1048  {
1049  FYZ[iz][iy][ix],
1050  FYZ[iz + 1][iy][ix]
1051  },
1052  {
1053  FYZ[iz][iy + 1][ix],
1054  FYZ[iz + 1][iy + 1][ix]
1055  }
1056  },
1057  {
1058  {
1059  FYZ[iz][iy][ix + 1],
1060  FYZ[iz + 1][iy][ix + 1]
1061  },
1062  {
1063  FYZ[iz][iy + 1][ix + 1],
1064  FYZ[iz + 1][iy + 1][ix + 1]
1065  }
1066  }
1067  };
1068 
1069  Real GXYZ[2][2][2] =
1070  {
1071  {
1072  {
1073  FXYZ[iz][iy][ix],
1074  FXYZ[iz + 1][iy][ix]
1075  },
1076  {
1077  FXYZ[iz][iy + 1][ix],
1078  FXYZ[iz + 1][iy + 1][ix]
1079  }
1080  },
1081  {
1082  {
1083  FXYZ[iz][iy][ix + 1],
1084  FXYZ[iz + 1][iy][ix + 1]
1085  },
1086  {
1087  FXYZ[iz][iy + 1][ix + 1],
1088  FXYZ[iz + 1][iy + 1][ix + 1]
1089  }
1090  }
1091  };
1092 
1093  Construct(mPoly[iz][iy][ix], G, GX, GY, GZ, GXY, GXZ, GYZ, GXYZ);
1094  }
1095  }
1096  }
1097 }
1098 
1099 template <typename Real>
1100 Real IntpAkimaUniform3<Real>::ComputeDerivative(Real const* slope) const
1101 {
1102  if (slope[1] != slope[2])
1103  {
1104  if (slope[0] != slope[1])
1105  {
1106  if (slope[2] != slope[3])
1107  {
1108  Real ad0 = std::abs(slope[3] - slope[2]);
1109  Real ad1 = std::abs(slope[0] - slope[1]);
1110  return (ad0 * slope[1] + ad1 * slope[2]) / (ad0 + ad1);
1111  }
1112  else
1113  {
1114  return slope[2];
1115  }
1116  }
1117  else
1118  {
1119  if (slope[2] != slope[3])
1120  {
1121  return slope[1];
1122  }
1123  else
1124  {
1125  return ((Real)0.5) * (slope[1] + slope[2]);
1126  }
1127  }
1128  }
1129  else
1130  {
1131  return slope[1];
1132  }
1133 }
1134 
1135 template <typename Real>
1136 void IntpAkimaUniform3<Real>::Construct(Polynomial& poly,
1137  Real const F[2][2][2], Real const FX[2][2][2], Real const FY[2][2][2],
1138  Real const FZ[2][2][2], Real const FXY[2][2][2], Real const FXZ[2][2][2],
1139  Real const FYZ[2][2][2], Real const FXYZ[2][2][2])
1140 {
1141  Real dx = mXSpacing, dy = mYSpacing, dz = mZSpacing;
1142  Real invDX = ((Real)1) / dx, invDX2 = invDX * invDX;
1143  Real invDY = ((Real)1) / dy, invDY2 = invDY * invDY;
1144  Real invDZ = ((Real)1) / dz, invDZ2 = invDZ * invDZ;
1145  Real b0, b1, b2, b3, b4, b5, b6, b7;
1146 
1147  poly.A(0, 0, 0) = F[0][0][0];
1148  poly.A(1, 0, 0) = FX[0][0][0];
1149  poly.A(0, 1, 0) = FY[0][0][0];
1150  poly.A(0, 0, 1) = FZ[0][0][0];
1151  poly.A(1, 1, 0) = FXY[0][0][0];
1152  poly.A(1, 0, 1) = FXZ[0][0][0];
1153  poly.A(0, 1, 1) = FYZ[0][0][0];
1154  poly.A(1, 1, 1) = FXYZ[0][0][0];
1155 
1156  // solve for Aij0
1157  b0 = (F[1][0][0] - poly(0, 0, 0, dx, (Real)0, (Real)0))*invDX2;
1158  b1 = (FX[1][0][0] - poly(1, 0, 0, dx, (Real)0, (Real)0))*invDX;
1159  poly.A(2, 0, 0) = ((Real)3)*b0 - b1;
1160  poly.A(3, 0, 0) = (-((Real)2)*b0 + b1)*invDX;
1161 
1162  b0 = (F[0][1][0] - poly(0, 0, 0, (Real)0, dy, (Real)0))*invDY2;
1163  b1 = (FY[0][1][0] - poly(0, 1, 0, (Real)0, dy, (Real)0))*invDY;
1164  poly.A(0, 2, 0) = ((Real)3)*b0 - b1;
1165  poly.A(0, 3, 0) = (-((Real)2)*b0 + b1)*invDY;
1166 
1167  b0 = (FY[1][0][0] - poly(0, 1, 0, dx, (Real)0, (Real)0))*invDX2;
1168  b1 = (FXY[1][0][0] - poly(1, 1, 0, dx, (Real)0, (Real)0))*invDX;
1169  poly.A(2, 1, 0) = ((Real)3)*b0 - b1;
1170  poly.A(3, 1, 0) = (-((Real)2)*b0 + b1)*invDX;
1171 
1172  b0 = (FX[0][1][0] - poly(1, 0, 0, (Real)0, dy, (Real)0))*invDY2;
1173  b1 = (FXY[0][1][0] - poly(1, 1, 0, (Real)0, dy, (Real)0))*invDY;
1174  poly.A(1, 2, 0) = ((Real)3)*b0 - b1;
1175  poly.A(1, 3, 0) = (-((Real)2)*b0 + b1)*invDY;
1176 
1177  b0 = (F[1][1][0] - poly(0, 0, 0, dx, dy, (Real)0))*invDX2*invDY2;
1178  b1 = (FX[1][1][0] - poly(1, 0, 0, dx, dy, (Real)0))*invDX*invDY2;
1179  b2 = (FY[1][1][0] - poly(0, 1, 0, dx, dy, (Real)0))*invDX2*invDY;
1180  b3 = (FXY[1][1][0] - poly(1, 1, 0, dx, dy, (Real)0))*invDX*invDY;
1181  poly.A(2, 2, 0) = ((Real)9)*b0 - ((Real)3)*b1 - ((Real)3)*b2 + b3;
1182  poly.A(3, 2, 0) = (-((Real)6)*b0 + ((Real)3)*b1 + ((Real)2)*b2 - b3)*invDX;
1183  poly.A(2, 3, 0) = (-((Real)6)*b0 + ((Real)2)*b1 + ((Real)3)*b2 - b3)*invDY;
1184  poly.A(3, 3, 0) = (((Real)4)*b0 - ((Real)2)*b1 - ((Real)2)*b2 + b3)*invDX*invDY;
1185 
1186  // solve for Ai0k
1187  b0 = (F[0][0][1] - poly(0, 0, 0, (Real)0, (Real)0, dz))*invDZ2;
1188  b1 = (FZ[0][0][1] - poly(0, 0, 1, (Real)0, (Real)0, dz))*invDZ;
1189  poly.A(0, 0, 2) = ((Real)3)*b0 - b1;
1190  poly.A(0, 0, 3) = (-((Real)2)*b0 + b1)*invDZ;
1191 
1192  b0 = (FZ[1][0][0] - poly(0, 0, 1, dx, (Real)0, (Real)0))*invDX2;
1193  b1 = (FXZ[1][0][0] - poly(1, 0, 1, dx, (Real)0, (Real)0))*invDX;
1194  poly.A(2, 0, 1) = ((Real)3)*b0 - b1;
1195  poly.A(3, 0, 1) = (-((Real)2)*b0 + b1)*invDX;
1196 
1197  b0 = (FX[0][0][1] - poly(1, 0, 0, (Real)0, (Real)0, dz))*invDZ2;
1198  b1 = (FXZ[0][0][1] - poly(1, 0, 1, (Real)0, (Real)0, dz))*invDZ;
1199  poly.A(1, 0, 2) = ((Real)3)*b0 - b1;
1200  poly.A(1, 0, 3) = (-((Real)2)*b0 + b1)*invDZ;
1201 
1202  b0 = (F[1][0][1] - poly(0, 0, 0, dx, (Real)0, dz))*invDX2*invDZ2;
1203  b1 = (FX[1][0][1] - poly(1, 0, 0, dx, (Real)0, dz))*invDX*invDZ2;
1204  b2 = (FZ[1][0][1] - poly(0, 0, 1, dx, (Real)0, dz))*invDX2*invDZ;
1205  b3 = (FXZ[1][0][1] - poly(1, 0, 1, dx, (Real)0, dz))*invDX*invDZ;
1206  poly.A(2, 0, 2) = ((Real)9)*b0 - ((Real)3)*b1 - ((Real)3)*b2 + b3;
1207  poly.A(3, 0, 2) = (-((Real)6)*b0 + ((Real)3)*b1 + ((Real)2)*b2 - b3)*invDX;
1208  poly.A(2, 0, 3) = (-((Real)6)*b0 + ((Real)2)*b1 + ((Real)3)*b2 - b3)*invDZ;
1209  poly.A(3, 0, 3) = (((Real)4)*b0 - ((Real)2)*b1 - ((Real)2)*b2 + b3)*invDX*invDZ;
1210 
1211  // solve for A0jk
1212  b0 = (FZ[0][1][0] - poly(0, 0, 1, (Real)0, dy, (Real)0))*invDY2;
1213  b1 = (FYZ[0][1][0] - poly(0, 1, 1, (Real)0, dy, (Real)0))*invDY;
1214  poly.A(0, 2, 1) = ((Real)3)*b0 - b1;
1215  poly.A(0, 3, 1) = (-((Real)2)*b0 + b1)*invDY;
1216 
1217  b0 = (FY[0][0][1] - poly(0, 1, 0, (Real)0, (Real)0, dz))*invDZ2;
1218  b1 = (FYZ[0][0][1] - poly(0, 1, 1, (Real)0, (Real)0, dz))*invDZ;
1219  poly.A(0, 1, 2) = ((Real)3)*b0 - b1;
1220  poly.A(0, 1, 3) = (-((Real)2)*b0 + b1)*invDZ;
1221 
1222  b0 = (F[0][1][1] - poly(0, 0, 0, (Real)0, dy, dz))*invDY2*invDZ2;
1223  b1 = (FY[0][1][1] - poly(0, 1, 0, (Real)0, dy, dz))*invDY*invDZ2;
1224  b2 = (FZ[0][1][1] - poly(0, 0, 1, (Real)0, dy, dz))*invDY2*invDZ;
1225  b3 = (FYZ[0][1][1] - poly(0, 1, 1, (Real)0, dy, dz))*invDY*invDZ;
1226  poly.A(0, 2, 2) = ((Real)9)*b0 - ((Real)3)*b1 - ((Real)3)*b2 + b3;
1227  poly.A(0, 3, 2) = (-((Real)6)*b0 + ((Real)3)*b1 + ((Real)2)*b2 - b3)*invDY;
1228  poly.A(0, 2, 3) = (-((Real)6)*b0 + ((Real)2)*b1 + ((Real)3)*b2 - b3)*invDZ;
1229  poly.A(0, 3, 3) = (((Real)4)*b0 - ((Real)2)*b1 - ((Real)2)*b2 + b3)*invDY*invDZ;
1230 
1231  // solve for Aij1
1232  b0 = (FYZ[1][0][0] - poly(0, 1, 1, dx, (Real)0, (Real)0))*invDX2;
1233  b1 = (FXYZ[1][0][0] - poly(1, 1, 1, dx, (Real)0, (Real)0))*invDX;
1234  poly.A(2, 1, 1) = ((Real)3)*b0 - b1;
1235  poly.A(3, 1, 1) = (-((Real)2)*b0 + b1)*invDX;
1236 
1237  b0 = (FXZ[0][1][0] - poly(1, 0, 1, (Real)0, dy, (Real)0))*invDY2;
1238  b1 = (FXYZ[0][1][0] - poly(1, 1, 1, (Real)0, dy, (Real)0))*invDY;
1239  poly.A(1, 2, 1) = ((Real)3)*b0 - b1;
1240  poly.A(1, 3, 1) = (-((Real)2)*b0 + b1)*invDY;
1241 
1242  b0 = (FZ[1][1][0] - poly(0, 0, 1, dx, dy, (Real)0))*invDX2*invDY2;
1243  b1 = (FXZ[1][1][0] - poly(1, 0, 1, dx, dy, (Real)0))*invDX*invDY2;
1244  b2 = (FYZ[1][1][0] - poly(0, 1, 1, dx, dy, (Real)0))*invDX2*invDY;
1245  b3 = (FXYZ[1][1][0] - poly(1, 1, 1, dx, dy, (Real)0))*invDX*invDY;
1246  poly.A(2, 2, 1) = ((Real)9)*b0 - ((Real)3)*b1 - ((Real)3)*b2 + b3;
1247  poly.A(3, 2, 1) = (-((Real)6)*b0 + ((Real)3)*b1 + ((Real)2)*b2 - b3)*invDX;
1248  poly.A(2, 3, 1) = (-((Real)6)*b0 + ((Real)2)*b1 + ((Real)3)*b2 - b3)*invDY;
1249  poly.A(3, 3, 1) = (((Real)4)*b0 - ((Real)2)*b1 - ((Real)2)*b2 + b3)*invDX*invDY;
1250 
1251  // solve for Ai1k
1252  b0 = (FXY[0][0][1] - poly(1, 1, 0, (Real)0, (Real)0, dz))*invDZ2;
1253  b1 = (FXYZ[0][0][1] - poly(1, 1, 1, (Real)0, (Real)0, dz))*invDZ;
1254  poly.A(1, 1, 2) = ((Real)3)*b0 - b1;
1255  poly.A(1, 1, 3) = (-((Real)2)*b0 + b1)*invDZ;
1256 
1257  b0 = (FY[1][0][1] - poly(0, 1, 0, dx, (Real)0, dz))*invDX2*invDZ2;
1258  b1 = (FXY[1][0][1] - poly(1, 1, 0, dx, (Real)0, dz))*invDX*invDZ2;
1259  b2 = (FYZ[1][0][1] - poly(0, 1, 1, dx, (Real)0, dz))*invDX2*invDZ;
1260  b3 = (FXYZ[1][0][1] - poly(1, 1, 1, dx, (Real)0, dz))*invDX*invDZ;
1261  poly.A(2, 1, 2) = ((Real)9)*b0 - ((Real)3)*b1 - ((Real)3)*b2 + b3;
1262  poly.A(3, 1, 2) = (-((Real)6)*b0 + ((Real)3)*b1 + ((Real)2)*b2 - b3)*invDX;
1263  poly.A(2, 1, 3) = (-((Real)6)*b0 + ((Real)2)*b1 + ((Real)3)*b2 - b3)*invDZ;
1264  poly.A(3, 1, 3) = (((Real)4)*b0 - ((Real)2)*b1 - ((Real)2)*b2 + b3)*invDX*invDZ;
1265 
1266  // solve for A1jk
1267  b0 = (FX[0][1][1] - poly(1, 0, 0, (Real)0, dy, dz))*invDY2*invDZ2;
1268  b1 = (FXY[0][1][1] - poly(1, 1, 0, (Real)0, dy, dz))*invDY*invDZ2;
1269  b2 = (FXZ[0][1][1] - poly(1, 0, 1, (Real)0, dy, dz))*invDY2*invDZ;
1270  b3 = (FXYZ[0][1][1] - poly(1, 1, 1, (Real)0, dy, dz))*invDY*invDZ;
1271  poly.A(1, 2, 2) = ((Real)9)*b0 - ((Real)3)*b1 - ((Real)3)*b2 + b3;
1272  poly.A(1, 3, 2) = (-((Real)6)*b0 + ((Real)3)*b1 + ((Real)2)*b2 - b3)*invDY;
1273  poly.A(1, 2, 3) = (-((Real)6)*b0 + ((Real)2)*b1 + ((Real)3)*b2 - b3)*invDZ;
1274  poly.A(1, 3, 3) = (((Real)4)*b0 - ((Real)2)*b1 - ((Real)2)*b2 + b3)*invDY*invDZ;
1275 
1276  // solve for remaining Aijk with i >= 2, j >= 2, k >= 2
1277  b0 = (F[1][1][1] - poly(0, 0, 0, dx, dy, dz))*invDX2*invDY2*invDZ2;
1278  b1 = (FX[1][1][1] - poly(1, 0, 0, dx, dy, dz))*invDX*invDY2*invDZ2;
1279  b2 = (FY[1][1][1] - poly(0, 1, 0, dx, dy, dz))*invDX2*invDY*invDZ2;
1280  b3 = (FZ[1][1][1] - poly(0, 0, 1, dx, dy, dz))*invDX2*invDY2*invDZ;
1281  b4 = (FXY[1][1][1] - poly(1, 1, 0, dx, dy, dz))*invDX*invDY*invDZ2;
1282  b5 = (FXZ[1][1][1] - poly(1, 0, 1, dx, dy, dz))*invDX*invDY2*invDZ;
1283  b6 = (FYZ[1][1][1] - poly(0, 1, 1, dx, dy, dz))*invDX2*invDY*invDZ;
1284  b7 = (FXYZ[1][1][1] - poly(1, 1, 1, dx, dy, dz))*invDX*invDY*invDZ;
1285  poly.A(2, 2, 2) = ((Real)27)*b0 - ((Real)9)*b1 - ((Real)9)*b2 -
1286  ((Real)9)*b3 + ((Real)3)*b4 + ((Real)3)*b5 + ((Real)3)*b6 - b7;
1287  poly.A(3, 2, 2) = (((Real)-18)*b0 + ((Real)9)*b1 + ((Real)6)*b2 +
1288  ((Real)6)*b3 - ((Real)3)*b4 - ((Real)3)*b5 - ((Real)2)*b6 + b7)*invDX;
1289  poly.A(2, 3, 2) = (((Real)-18)*b0 + ((Real)6)*b1 + ((Real)9)*b2 +
1290  ((Real)6)*b3 - ((Real)3)*b4 - ((Real)2)*b5 - ((Real)3)*b6 + b7)*invDY;
1291  poly.A(2, 2, 3) = (((Real)-18)*b0 + ((Real)6)*b1 + ((Real)6)*b2 +
1292  ((Real)9)*b3 - ((Real)2)*b4 - ((Real)3)*b5 - ((Real)3)*b6 + b7)*invDZ;
1293  poly.A(3, 3, 2) = (((Real)12)*b0 - ((Real)6)*b1 - ((Real)6)*b2 -
1294  ((Real)4)*b3 + ((Real)3)*b4 + ((Real)2)*b5 + ((Real)2)*b6 - b7)*
1295  invDX*invDY;
1296  poly.A(3, 2, 3) = (((Real)12)*b0 - ((Real)6)*b1 - ((Real)4)*b2 -
1297  ((Real)6)*b3 + ((Real)2)*b4 + ((Real)3)*b5 + ((Real)2)*b6 - b7)*
1298  invDX*invDZ;
1299  poly.A(2, 3, 3) = (((Real)12)*b0 - ((Real)4)*b1 - ((Real)6)*b2 -
1300  ((Real)6)*b3 + ((Real)2)*b4 + ((Real)2)*b5 + ((Real)3)*b6 - b7)*
1301  invDY*invDZ;
1302  poly.A(3, 3, 3) = (((Real)-8)*b0 + ((Real)4)*b1 + ((Real)4)*b2 +
1303  ((Real)4)*b3 - ((Real)2)*b4 - ((Real)2)*b5 - ((Real)2)*b6 + b7)*
1304  invDX*invDY*invDZ;
1305 }
1306 
1307 template <typename Real>
1308 void IntpAkimaUniform3<Real>::XLookup(Real x, int& xIndex, Real& dx) const
1309 {
1310  for (xIndex = 0; xIndex + 1 < mXBound; ++xIndex)
1311  {
1312  if (x < mXMin + mXSpacing*(xIndex + 1))
1313  {
1314  dx = x - (mXMin + mXSpacing*xIndex);
1315  return;
1316  }
1317  }
1318 
1319  --xIndex;
1320  dx = x - (mXMin + mXSpacing*xIndex);
1321 }
1322 
1323 template <typename Real>
1324 void IntpAkimaUniform3<Real>::YLookup(Real y, int& yIndex, Real& dy) const
1325 {
1326  for (yIndex = 0; yIndex + 1 < mYBound; ++yIndex)
1327  {
1328  if (y < mYMin + mYSpacing*(yIndex + 1))
1329  {
1330  dy = y - (mYMin + mYSpacing*yIndex);
1331  return;
1332  }
1333  }
1334 
1335  --yIndex;
1336  dy = y - (mYMin + mYSpacing*yIndex);
1337 }
1338 
1339 template <typename Real>
1340 void IntpAkimaUniform3<Real>::ZLookup(Real z, int& zIndex, Real& dz) const
1341 {
1342  for (zIndex = 0; zIndex + 1 < mZBound; ++zIndex)
1343  {
1344  if (z < mZMin + mZSpacing*(zIndex + 1))
1345  {
1346  dz = z - (mZMin + mZSpacing*zIndex);
1347  return;
1348  }
1349  }
1350 
1351  --zIndex;
1352  dz = z - (mZMin + mZSpacing*zIndex);
1353 }
1354 
1355 template <typename Real>
1356 IntpAkimaUniform3<Real>::Polynomial::Polynomial()
1357 {
1358  memset(&mCoeff[0][0][0], 0, 64 * sizeof(Real));
1359 }
1360 
1361 template <typename Real>
1362 Real& IntpAkimaUniform3<Real>::Polynomial::A(int ix, int iy, int iz)
1363 {
1364  return mCoeff[ix][iy][iz];
1365 }
1366 
1367 template <typename Real>
1368 Real IntpAkimaUniform3<Real>::Polynomial::operator()(Real x, Real y, Real z)
1369 const
1370 {
1371  Real xPow[4] = { (Real)1, x, x * x, x * x * x };
1372  Real yPow[4] = { (Real)1, y, y * y, y * y * y };
1373  Real zPow[4] = { (Real)1, z, z * z, z * z * z };
1374 
1375  Real p = (Real)0;
1376  for (int iz = 0; iz <= 3; ++iz)
1377  {
1378  for (int iy = 0; iy <= 3; ++iy)
1379  {
1380  for (int ix = 0; ix <= 3; ++ix)
1381  {
1382  p += mCoeff[ix][iy][iz] * xPow[ix] * yPow[iy] * zPow[iz];
1383  }
1384  }
1385  }
1386 
1387  return p;
1388 }
1389 
1390 template <typename Real>
1391 Real IntpAkimaUniform3<Real>::Polynomial::operator()(int xOrder, int yOrder,
1392  int zOrder, Real x, Real y, Real z) const
1393 {
1394  Real xPow[4];
1395  switch (xOrder)
1396  {
1397  case 0:
1398  xPow[0] = (Real)1;
1399  xPow[1] = x;
1400  xPow[2] = x * x;
1401  xPow[3] = x * x * x;
1402  break;
1403  case 1:
1404  xPow[0] = (Real)0;
1405  xPow[1] = (Real)1;
1406  xPow[2] = ((Real)2) * x;
1407  xPow[3] = ((Real)3) * x * x;
1408  break;
1409  case 2:
1410  xPow[0] = (Real)0;
1411  xPow[1] = (Real)0;
1412  xPow[2] = (Real)2;
1413  xPow[3] = ((Real)6) * x;
1414  break;
1415  case 3:
1416  xPow[0] = (Real)0;
1417  xPow[1] = (Real)0;
1418  xPow[2] = (Real)0;
1419  xPow[3] = (Real)6;
1420  break;
1421  default:
1422  return (Real)0;
1423  }
1424 
1425  Real yPow[4];
1426  switch (yOrder)
1427  {
1428  case 0:
1429  yPow[0] = (Real)1;
1430  yPow[1] = y;
1431  yPow[2] = y * y;
1432  yPow[3] = y * y * y;
1433  break;
1434  case 1:
1435  yPow[0] = (Real)0;
1436  yPow[1] = (Real)1;
1437  yPow[2] = ((Real)2) * y;
1438  yPow[3] = ((Real)3) * y * y;
1439  break;
1440  case 2:
1441  yPow[0] = (Real)0;
1442  yPow[1] = (Real)0;
1443  yPow[2] = (Real)2;
1444  yPow[3] = ((Real)6) * y;
1445  break;
1446  case 3:
1447  yPow[0] = (Real)0;
1448  yPow[1] = (Real)0;
1449  yPow[2] = (Real)0;
1450  yPow[3] = (Real)6;
1451  break;
1452  default:
1453  return (Real)0;
1454  }
1455 
1456  Real zPow[4];
1457  switch (zOrder)
1458  {
1459  case 0:
1460  zPow[0] = (Real)1;
1461  zPow[1] = z;
1462  zPow[2] = z * z;
1463  zPow[3] = z * z * z;
1464  break;
1465  case 1:
1466  zPow[0] = (Real)0;
1467  zPow[1] = (Real)1;
1468  zPow[2] = ((Real)2) * z;
1469  zPow[3] = ((Real)3) * z * z;
1470  break;
1471  case 2:
1472  zPow[0] = (Real)0;
1473  zPow[1] = (Real)0;
1474  zPow[2] = (Real)2;
1475  zPow[3] = ((Real)6) * z;
1476  break;
1477  case 3:
1478  zPow[0] = (Real)0;
1479  zPow[1] = (Real)0;
1480  zPow[2] = (Real)0;
1481  zPow[3] = (Real)6;
1482  break;
1483  default:
1484  return (Real)0;
1485  }
1486 
1487  Real p = (Real)0;
1488 
1489  for (int iz = 0; iz <= 3; ++iz)
1490  {
1491  for (int iy = 0; iy <= 3; ++iy)
1492  {
1493  for (int ix = 0; ix <= 3; ++ix)
1494  {
1495  p += mCoeff[ix][iy][iz] * xPow[ix] * yPow[iy] * zPow[iz];
1496  }
1497  }
1498  }
1499 
1500  return p;
1501 }
1502 
1503 }
gte::IntpAkimaUniform3::GetFX
void GetFX(Array3< Real > const &F, Array3< Real > &FX)
Definition: GteIntpAkimaUniform3.h:285
gte::IntpAkimaUniform3::GetFY
void GetFY(Array3< Real > const &F, Array3< Real > &FY)
Definition: GteIntpAkimaUniform3.h:319
gte::IntpAkimaUniform3::YLookup
void YLookup(Real y, int &yIndex, Real &dy) const
Definition: GteIntpAkimaUniform3.h:1324
std::abs
gte::BSNumber< UIntegerType > abs(gte::BSNumber< UIntegerType > const &number)
Definition: GteBSNumber.h:966
gte::IntpAkimaUniform3::GetFXYZ
void GetFXYZ(Array3< Real > const &F, Array3< Real > &FXYZ)
Definition: GteIntpAkimaUniform3.h:681
LogAssert
#define LogAssert(condition, message)
Definition: GteLogger.h:86
gte::IntpAkimaUniform3::ZLookup
void ZLookup(Real z, int &zIndex, Real &dz) const
Definition: GteIntpAkimaUniform3.h:1340
gte::IntpAkimaUniform3::GetFXY
void GetFXY(Array3< Real > const &F, Array3< Real > &FXY)
Definition: GteIntpAkimaUniform3.h:387
gte::IntpAkimaUniform3::GetPolynomials
void GetPolynomials(Array3< Real > const &F, Array3< Real > const &FX, Array3< Real > const &FY, Array3< Real > const &FZ, Array3< Real > const &FXY, Array3< Real > const &FXZ, Array3< Real > const &FYZ, Array3< Real > const &FXYZ)
Definition: GteIntpAkimaUniform3.h:886
gte::IntpAkimaUniform3::ComputeDerivative
Real ComputeDerivative(Real const *slope) const
Definition: GteIntpAkimaUniform3.h:1100
gte::IntpAkimaUniform3::mPoly
Array3< Polynomial > mPoly
Definition: GteIntpAkimaUniform3.h:106
gte::IntpAkimaUniform3::GetFYZ
void GetFYZ(Array3< Real > const &F, Array3< Real > &FYZ)
Definition: GteIntpAkimaUniform3.h:583
x
GLint GLenum GLint x
Definition: glcorearb.h:404
mask
GLint GLuint mask
Definition: glcorearb.h:119
gte::IntpAkimaUniform3::GetF
Real const * GetF() const
Definition: GteIntpAkimaUniform3.h:196
gte::IntpAkimaUniform3::GetFZ
void GetFZ(Array3< Real > const &F, Array3< Real > &FZ)
Definition: GteIntpAkimaUniform3.h:353
gte::IntpAkimaUniform3::operator()
Real operator()(Real x, Real y, Real z) const
Definition: GteIntpAkimaUniform3.h:256
gte::IntpAkimaUniform3::IntpAkimaUniform3
IntpAkimaUniform3(int xBound, int yBound, int zBound, Real xMin, Real xSpacing, Real yMin, Real ySpacing, Real zMin, Real zSpacing, Real const *F)
Definition: GteIntpAkimaUniform3.h:116
gte::IntpAkimaUniform3::Polynomial::operator()
Real operator()(Real x, Real y, Real z) const
Definition: GteIntpAkimaUniform3.h:1368
z
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:843
gte::IntpAkimaUniform3::XLookup
void XLookup(Real x, int &xIndex, Real &dx) const
Definition: GteIntpAkimaUniform3.h:1308
gte::IntpAkimaUniform3::GetFXZ
void GetFXZ(Array3< Real > const &F, Array3< Real > &FXZ)
Definition: GteIntpAkimaUniform3.h:485
gte::IntpAkimaUniform3::Polynomial::A
Real & A(int ix, int iy, int iz)
Definition: GteIntpAkimaUniform3.h:1362
p
GLfloat GLfloat p
Definition: glext.h:11668
y
GLint y
Definition: glcorearb.h:98
gte::IntpAkimaUniform3::Construct
void Construct(Polynomial &poly, Real const F[2][2][2], Real const FX[2][2][2], Real const FY[2][2][2], Real const FZ[2][2][2], Real const FXY[2][2][2], Real const FXZ[2][2][2], Real const FYZ[2][2][2], Real const FXYZ[2][2][2])
Definition: GteIntpAkimaUniform3.h:1136

geometric_tools_engine
Author(s): Yijiang Huang
autogenerated on Thu Jul 18 2019 04:00:00