HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FiniteDifference.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2012-2017 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
29 ///////////////////////////////////////////////////////////////////////////
30 
31 /// @file FiniteDifference.h
32 
33 #ifndef OPENVDB_MATH_FINITEDIFFERENCE_HAS_BEEN_INCLUDED
34 #define OPENVDB_MATH_FINITEDIFFERENCE_HAS_BEEN_INCLUDED
35 
36 #include <openvdb/Types.h>
37 #include "Math.h"
38 #include "Coord.h"
39 #include "Vec3.h"
40 
41 #include <hboost/algorithm/string/case_conv.hpp>
42 #include <hboost/algorithm/string/trim.hpp>
43 
44 #ifdef DWA_OPENVDB
45 #include <simd/Simd.h>
46 #endif
47 
48 namespace openvdb {
50 namespace OPENVDB_VERSION_NAME {
51 namespace math {
52 
53 
54 ////////////////////////////////////////
55 
56 
57 /// @brief Different discrete schemes used in the first derivatives.
58 // Add new items to the *end* of this list, and update NUM_DS_SCHEMES.
59 enum DScheme {
60  UNKNOWN_DS = -1,
61  CD_2NDT = 0, // center difference, 2nd order, but the result must be divided by 2
62  CD_2ND, // center difference, 2nd order
63  CD_4TH, // center difference, 4th order
64  CD_6TH, // center difference, 6th order
65  FD_1ST, // forward difference, 1st order
66  FD_2ND, // forward difference, 2nd order
67  FD_3RD, // forward difference, 3rd order
68  BD_1ST, // backward difference, 1st order
69  BD_2ND, // backward difference, 2nd order
70  BD_3RD, // backward difference, 3rd order
71  FD_WENO5, // forward difference, weno5
72  BD_WENO5, // backward difference, weno5
73  FD_HJWENO5, // forward differene, HJ-weno5
74  BD_HJWENO5 // backward difference, HJ-weno5
75 };
76 
77 enum { NUM_DS_SCHEMES = BD_HJWENO5 + 1 };
78 
79 
80 inline std::string
82 {
83  std::string ret;
84  switch (dss) {
85  case UNKNOWN_DS: ret = "unknown_ds"; break;
86  case CD_2NDT: ret = "cd_2ndt"; break;
87  case CD_2ND: ret = "cd_2nd"; break;
88  case CD_4TH: ret = "cd_4th"; break;
89  case CD_6TH: ret = "cd_6th"; break;
90  case FD_1ST: ret = "fd_1st"; break;
91  case FD_2ND: ret = "fd_2nd"; break;
92  case FD_3RD: ret = "fd_3rd"; break;
93  case BD_1ST: ret = "bd_1st"; break;
94  case BD_2ND: ret = "bd_2nd"; break;
95  case BD_3RD: ret = "bd_3rd"; break;
96  case FD_WENO5: ret = "fd_weno5"; break;
97  case BD_WENO5: ret = "bd_weno5"; break;
98  case FD_HJWENO5: ret = "fd_hjweno5"; break;
99  case BD_HJWENO5: ret = "bd_hjweno5"; break;
100  }
101  return ret;
102 }
103 
104 inline DScheme
106 {
107  DScheme ret = UNKNOWN_DS;
108 
109  std::string str = s;
110  hboost::trim(str);
111  hboost::to_lower(str);
112 
113  if (str == dsSchemeToString(CD_2NDT)) {
114  ret = CD_2NDT;
115  } else if (str == dsSchemeToString(CD_2ND)) {
116  ret = CD_2ND;
117  } else if (str == dsSchemeToString(CD_4TH)) {
118  ret = CD_4TH;
119  } else if (str == dsSchemeToString(CD_6TH)) {
120  ret = CD_6TH;
121  } else if (str == dsSchemeToString(FD_1ST)) {
122  ret = FD_1ST;
123  } else if (str == dsSchemeToString(FD_2ND)) {
124  ret = FD_2ND;
125  } else if (str == dsSchemeToString(FD_3RD)) {
126  ret = FD_3RD;
127  } else if (str == dsSchemeToString(BD_1ST)) {
128  ret = BD_1ST;
129  } else if (str == dsSchemeToString(BD_2ND)) {
130  ret = BD_2ND;
131  } else if (str == dsSchemeToString(BD_3RD)) {
132  ret = BD_3RD;
133  } else if (str == dsSchemeToString(FD_WENO5)) {
134  ret = FD_WENO5;
135  } else if (str == dsSchemeToString(BD_WENO5)) {
136  ret = BD_WENO5;
137  } else if (str == dsSchemeToString(FD_HJWENO5)) {
138  ret = FD_HJWENO5;
139  } else if (str == dsSchemeToString(BD_HJWENO5)) {
140  ret = BD_HJWENO5;
141  }
142 
143  return ret;
144 }
145 
146 inline std::string
148 {
149  std::string ret;
150  switch (dss) {
151  case UNKNOWN_DS: ret = "Unknown DS scheme"; break;
152  case CD_2NDT: ret = "Twice 2nd-order center difference"; break;
153  case CD_2ND: ret = "2nd-order center difference"; break;
154  case CD_4TH: ret = "4th-order center difference"; break;
155  case CD_6TH: ret = "6th-order center difference"; break;
156  case FD_1ST: ret = "1st-order forward difference"; break;
157  case FD_2ND: ret = "2nd-order forward difference"; break;
158  case FD_3RD: ret = "3rd-order forward difference"; break;
159  case BD_1ST: ret = "1st-order backward difference"; break;
160  case BD_2ND: ret = "2nd-order backward difference"; break;
161  case BD_3RD: ret = "3rd-order backward difference"; break;
162  case FD_WENO5: ret = "5th-order WENO forward difference"; break;
163  case BD_WENO5: ret = "5th-order WENO backward difference"; break;
164  case FD_HJWENO5: ret = "5th-order HJ-WENO forward difference"; break;
165  case BD_HJWENO5: ret = "5th-order HJ-WENO backward difference"; break;
166  }
167  return ret;
168 }
169 
170 
171 
172 ////////////////////////////////////////
173 
174 
175 /// @brief Different discrete schemes used in the second derivatives.
176 // Add new items to the *end* of this list, and update NUM_DD_SCHEMES.
177 enum DDScheme {
179  CD_SECOND = 0, // center difference, 2nd order
180  CD_FOURTH, // center difference, 4th order
181  CD_SIXTH // center difference, 6th order
182 };
183 
184 enum { NUM_DD_SCHEMES = CD_SIXTH + 1 };
185 
186 
187 ////////////////////////////////////////
188 
189 
190 /// @brief Biased Gradients are limited to non-centered differences
191 // Add new items to the *end* of this list, and update NUM_BIAS_SCHEMES.
194  FIRST_BIAS = 0, // uses FD_1ST & BD_1ST
195  SECOND_BIAS, // uses FD_2ND & BD_2ND
196  THIRD_BIAS, // uses FD_3RD & BD_3RD
197  WENO5_BIAS, // uses WENO5
198  HJWENO5_BIAS // uses HJWENO5
199 };
200 
202 
203 inline std::string
205 {
206  std::string ret;
207  switch (bgs) {
208  case UNKNOWN_BIAS: ret = "unknown_bias"; break;
209  case FIRST_BIAS: ret = "first_bias"; break;
210  case SECOND_BIAS: ret = "second_bias"; break;
211  case THIRD_BIAS: ret = "third_bias"; break;
212  case WENO5_BIAS: ret = "weno5_bias"; break;
213  case HJWENO5_BIAS: ret = "hjweno5_bias"; break;
214  }
215  return ret;
216 }
217 
220 {
222 
223  std::string str = s;
224  hboost::trim(str);
225  hboost::to_lower(str);
226 
228  ret = FIRST_BIAS;
229  } else if (str == biasedGradientSchemeToString(SECOND_BIAS)) {
230  ret = SECOND_BIAS;
231  } else if (str == biasedGradientSchemeToString(THIRD_BIAS)) {
232  ret = THIRD_BIAS;
233  } else if (str == biasedGradientSchemeToString(WENO5_BIAS)) {
234  ret = WENO5_BIAS;
235  } else if (str == biasedGradientSchemeToString(HJWENO5_BIAS)) {
236  ret = HJWENO5_BIAS;
237  }
238  return ret;
239 }
240 
241 inline std::string
243 {
244  std::string ret;
245  switch (bgs) {
246  case UNKNOWN_BIAS: ret = "Unknown biased gradient"; break;
247  case FIRST_BIAS: ret = "1st-order biased gradient"; break;
248  case SECOND_BIAS: ret = "2nd-order biased gradient"; break;
249  case THIRD_BIAS: ret = "3rd-order biased gradient"; break;
250  case WENO5_BIAS: ret = "5th-order WENO biased gradient"; break;
251  case HJWENO5_BIAS: ret = "5th-order HJ-WENO biased gradient"; break;
252  }
253  return ret;
254 }
255 
256 ////////////////////////////////////////
257 
258 
259 /// @brief Temporal integration schemes
260 // Add new items to the *end* of this list, and update NUM_TEMPORAL_SCHEMES.
263  TVD_RK1,//same as explicit Euler integration
266 };
267 
269 
270 inline std::string
272 {
273  std::string ret;
274  switch (tis) {
275  case UNKNOWN_TIS: ret = "unknown_tis"; break;
276  case TVD_RK1: ret = "tvd_rk1"; break;
277  case TVD_RK2: ret = "tvd_rk2"; break;
278  case TVD_RK3: ret = "tvd_rk3"; break;
279  }
280  return ret;
281 }
282 
285 {
287 
288  std::string str = s;
289  hboost::trim(str);
290  hboost::to_lower(str);
291 
293  ret = TVD_RK1;
294  } else if (str == temporalIntegrationSchemeToString(TVD_RK2)) {
295  ret = TVD_RK2;
296  } else if (str == temporalIntegrationSchemeToString(TVD_RK3)) {
297  ret = TVD_RK3;
298  }
299 
300  return ret;
301 }
302 
303 inline std::string
305 {
306  std::string ret;
307  switch (tis) {
308  case UNKNOWN_TIS: ret = "Unknown temporal integration"; break;
309  case TVD_RK1: ret = "Forward Euler"; break;
310  case TVD_RK2: ret = "2nd-order Runge-Kutta"; break;
311  case TVD_RK3: ret = "3rd-order Runge-Kutta"; break;
312  }
313  return ret;
314 }
315 
316 
317 //@}
318 
319 
320 /// @brief Implementation of nominally fifth-order finite-difference WENO
321 /// @details This function returns the numerical flux. See "High Order Finite Difference and
322 /// Finite Volume WENO Schemes and Discontinuous Galerkin Methods for CFD" - Chi-Wang Shu
323 /// ICASE Report No 2001-11 (page 6). Also see ICASE No 97-65 for a more complete reference
324 /// (Shu, 1997).
325 /// Given v1 = f(x-2dx), v2 = f(x-dx), v3 = f(x), v4 = f(x+dx) and v5 = f(x+2dx),
326 /// return an interpolated value f(x+dx/2) with the special property that
327 /// ( f(x+dx/2) - f(x-dx/2) ) / dx = df/dx (x) + error,
328 /// where the error is fifth-order in smooth regions: O(dx) <= error <=O(dx^5)
329 template<typename ValueType>
330 inline ValueType
331 WENO5(const ValueType& v1, const ValueType& v2, const ValueType& v3,
332  const ValueType& v4, const ValueType& v5, float scale2 = 0.01f)
333 {
334  const double C = 13.0 / 12.0;
335  // WENO is formulated for non-dimensional equations, here the optional scale2
336  // is a reference value (squared) for the function being interpolated. For
337  // example if 'v' is of order 1000, then scale2 = 10^6 is ok. But in practice
338  // leave scale2 = 1.
339  const double eps = 1.0e-6 * static_cast<double>(scale2);
340  // {\tilde \omega_k} = \gamma_k / ( \beta_k + \epsilon)^2 in Shu's ICASE report)
341  const double A1=0.1/math::Pow2(C*math::Pow2(v1-2*v2+v3)+0.25*math::Pow2(v1-4*v2+3.0*v3)+eps),
342  A2=0.6/math::Pow2(C*math::Pow2(v2-2*v3+v4)+0.25*math::Pow2(v2-v4)+eps),
343  A3=0.3/math::Pow2(C*math::Pow2(v3-2*v4+v5)+0.25*math::Pow2(3.0*v3-4*v4+v5)+eps);
344 
345  return static_cast<ValueType>(static_cast<ValueType>(
346  A1*(2.0*v1 - 7.0*v2 + 11.0*v3) +
347  A2*(5.0*v3 - v2 + 2.0*v4) +
348  A3*(2.0*v3 + 5.0*v4 - v5))/(6.0*(A1+A2+A3)));
349 }
350 
351 
352 template <typename Real>
353 inline Real GodunovsNormSqrd(bool isOutside,
354  Real dP_xm, Real dP_xp,
355  Real dP_ym, Real dP_yp,
356  Real dP_zm, Real dP_zp)
357 {
358  using math::Max;
359  using math::Min;
360  using math::Pow2;
361 
362  const Real zero(0);
363  Real dPLen2;
364  if (isOutside) { // outside
365  dPLen2 = Max(Pow2(Max(dP_xm, zero)), Pow2(Min(dP_xp,zero))); // (dP/dx)2
366  dPLen2 += Max(Pow2(Max(dP_ym, zero)), Pow2(Min(dP_yp,zero))); // (dP/dy)2
367  dPLen2 += Max(Pow2(Max(dP_zm, zero)), Pow2(Min(dP_zp,zero))); // (dP/dz)2
368  } else { // inside
369  dPLen2 = Max(Pow2(Min(dP_xm, zero)), Pow2(Max(dP_xp,zero))); // (dP/dx)2
370  dPLen2 += Max(Pow2(Min(dP_ym, zero)), Pow2(Max(dP_yp,zero))); // (dP/dy)2
371  dPLen2 += Max(Pow2(Min(dP_zm, zero)), Pow2(Max(dP_zp,zero))); // (dP/dz)2
372  }
373  return dPLen2; // |\nabla\phi|^2
374 }
375 
376 
377 template <typename Real>
379  Real dP_xm, Real dP_xp,
380  Real dP_ym, Real dP_yp,
381  Real dP_zm, Real dP_zp)
382 { return GodunovsNormSqrd(isOutside, dP_xm, dP_xp, dP_ym, dP_yp, dP_zm, dP_zp); }
383 
384 template<typename Real>
385 inline Real
386 GodunovsNormSqrd(bool isOutside, const Vec3<Real>& gradient_m, const Vec3<Real>& gradient_p)
387 {
388  return GodunovsNormSqrd<Real>(isOutside,
389  gradient_m[0], gradient_p[0],
390  gradient_m[1], gradient_p[1],
391  gradient_m[2], gradient_p[2]);
392 }
393 
394 template<typename Real>
396  const Vec3<Real>& gradient_m,
397  const Vec3<Real>& gradient_p)
398 {
399  return GodunovsNormSqrd<Real>(isOutside, gradient_m, gradient_p);
400 }
401 
402 
403 #ifdef DWA_OPENVDB
404 inline simd::Float4 simdMin(const simd::Float4& a, const simd::Float4& b) {
405  return simd::Float4(_mm_min_ps(a.base(), b.base()));
406 }
407 inline simd::Float4 simdMax(const simd::Float4& a, const simd::Float4& b) {
408  return simd::Float4(_mm_max_ps(a.base(), b.base()));
409 }
410 
411 inline float simdSum(const simd::Float4& v);
412 
413 inline simd::Float4 Pow2(const simd::Float4& v) { return v * v; }
414 
415 template<>
416 inline simd::Float4
417 WENO5<simd::Float4>(const simd::Float4& v1, const simd::Float4& v2, const simd::Float4& v3,
418  const simd::Float4& v4, const simd::Float4& v5, float scale2)
419 {
420  using math::Pow2;
421  typedef simd::Float4 F4;
422  const F4
423  C(13.f / 12.f),
424  eps(1.0e-6f * scale2),
425  two(2.0), three(3.0), four(4.0), five(5.0), fourth(0.25),
426  A1 = F4(0.1f) / Pow2(C*Pow2(v1-two*v2+v3) + fourth*Pow2(v1-four*v2+three*v3) + eps),
427  A2 = F4(0.6f) / Pow2(C*Pow2(v2-two*v3+v4) + fourth*Pow2(v2-v4) + eps),
428  A3 = F4(0.3f) / Pow2(C*Pow2(v3-two*v4+v5) + fourth*Pow2(three*v3-four*v4+v5) + eps);
429  return (A1 * (two * v1 - F4(7.0) * v2 + F4(11.0) * v3) +
430  A2 * (five * v3 - v2 + two * v4) +
431  A3 * (two * v3 + five * v4 - v5)) / (F4(6.0) * (A1 + A2 + A3));
432 }
433 
434 
435 inline float
436 simdSum(const simd::Float4& v)
437 {
438  // temp = { v3+v3, v2+v2, v1+v3, v0+v2 }
439  __m128 temp = _mm_add_ps(v.base(), _mm_movehl_ps(v.base(), v.base()));
440  // temp = { v3+v3, v2+v2, v1+v3, (v0+v2)+(v1+v3) }
441  temp = _mm_add_ss(temp, _mm_shuffle_ps(temp, temp, 1));
442  return _mm_cvtss_f32(temp);
443 }
444 
445 inline float
446 GodunovsNormSqrd(bool isOutside, const simd::Float4& dP_m, const simd::Float4& dP_p)
447 {
448  const simd::Float4 zero(0.0);
449  simd::Float4 v = isOutside
450  ? simdMax(math::Pow2(simdMax(dP_m, zero)), math::Pow2(simdMin(dP_p, zero)))
451  : simdMax(math::Pow2(simdMin(dP_m, zero)), math::Pow2(simdMax(dP_p, zero)));
452  return simdSum(v);//should be v[0]+v[1]+v[2]
453 }
454 
455 OPENVDB_DEPRECATED inline float GudonovsNormSqrd(bool isOutside,
456  const simd::Float4& dP_m,
457  const simd::Float4& dP_p)
458 {
459  return GodunovsNormSqrd(isOutside, dP_m, dP_p);
460 }
461 #endif
462 
463 template<DScheme DiffScheme>
464 struct D1
465 {
466  // random access version
467  template<typename Accessor>
468  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk);
469 
470  template<typename Accessor>
471  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk);
472 
473  template<typename Accessor>
474  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk);
475 
476  // stencil access version
477  template<typename Stencil>
478  static typename Stencil::ValueType inX(const Stencil& S);
479 
480  template<typename Stencil>
481  static typename Stencil::ValueType inY(const Stencil& S);
482 
483  template<typename Stencil>
484  static typename Stencil::ValueType inZ(const Stencil& S);
485 };
486 
487 template<>
488 struct D1<CD_2NDT>
489 {
490  // the difference opperator
491  template <typename ValueType>
492  static ValueType difference(const ValueType& xp1, const ValueType& xm1) {
493  return xp1 - xm1;
494  }
495 
496  // random access version
497  template<typename Accessor>
498  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
499  {
500  return difference(
501  grid.getValue(ijk.offsetBy(1, 0, 0)),
502  grid.getValue(ijk.offsetBy(-1, 0, 0)));
503  }
504 
505  template<typename Accessor>
506  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
507  {
508  return difference(
509  grid.getValue(ijk.offsetBy(0, 1, 0)),
510  grid.getValue(ijk.offsetBy( 0, -1, 0)));
511  }
512 
513  template<typename Accessor>
514  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
515  {
516  return difference(
517  grid.getValue(ijk.offsetBy(0, 0, 1)),
518  grid.getValue(ijk.offsetBy( 0, 0, -1)));
519  }
520 
521  // stencil access version
522  template<typename Stencil>
523  static typename Stencil::ValueType inX(const Stencil& S)
524  {
525  return difference( S.template getValue< 1, 0, 0>(), S.template getValue<-1, 0, 0>());
526  }
527 
528  template<typename Stencil>
529  static typename Stencil::ValueType inY(const Stencil& S)
530  {
531  return difference( S.template getValue< 0, 1, 0>(), S.template getValue< 0,-1, 0>());
532  }
533 
534  template<typename Stencil>
535  static typename Stencil::ValueType inZ(const Stencil& S)
536  {
537  return difference( S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0,-1>());
538  }
539 };
540 
541 template<>
542 struct D1<CD_2ND>
543 {
544 
545  // the difference opperator
546  template <typename ValueType>
547  static ValueType difference(const ValueType& xp1, const ValueType& xm1) {
548  return (xp1 - xm1)*ValueType(0.5);
549  }
550 
551 
552  // random access
553  template<typename Accessor>
554  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
555  {
556  return difference(
557  grid.getValue(ijk.offsetBy(1, 0, 0)),
558  grid.getValue(ijk.offsetBy(-1, 0, 0)));
559  }
560 
561  template<typename Accessor>
562  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
563  {
564  return difference(
565  grid.getValue(ijk.offsetBy(0, 1, 0)),
566  grid.getValue(ijk.offsetBy( 0, -1, 0)));
567  }
568 
569  template<typename Accessor>
570  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
571  {
572  return difference(
573  grid.getValue(ijk.offsetBy(0, 0, 1)),
574  grid.getValue(ijk.offsetBy( 0, 0, -1)));
575  }
576 
577 
578  // stencil access version
579  template<typename Stencil>
580  static typename Stencil::ValueType inX(const Stencil& S)
581  {
582  return difference(S.template getValue< 1, 0, 0>(), S.template getValue<-1, 0, 0>());
583  }
584  template<typename Stencil>
585  static typename Stencil::ValueType inY(const Stencil& S)
586  {
587  return difference(S.template getValue< 0, 1, 0>(), S.template getValue< 0,-1, 0>());
588  }
589 
590  template<typename Stencil>
591  static typename Stencil::ValueType inZ(const Stencil& S)
592  {
593  return difference(S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0,-1>());
594  }
595 
596 };
597 
598 template<>
599 struct D1<CD_4TH>
600 {
601 
602  // the difference opperator
603  template <typename ValueType>
604  static ValueType difference( const ValueType& xp2, const ValueType& xp1,
605  const ValueType& xm1, const ValueType& xm2 ) {
606  return ValueType(2./3.)*(xp1 - xm1) + ValueType(1./12.)*(xm2 - xp2) ;
607  }
608 
609 
610  // random access version
611  template<typename Accessor>
612  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
613  {
614  return difference(
615  grid.getValue(ijk.offsetBy( 2,0,0)), grid.getValue(ijk.offsetBy( 1,0,0)),
616  grid.getValue(ijk.offsetBy(-1,0,0)), grid.getValue(ijk.offsetBy(-2,0,0)) );
617  }
618 
619  template<typename Accessor>
620  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
621  {
622 
623  return difference(
624  grid.getValue(ijk.offsetBy( 0, 2, 0)), grid.getValue(ijk.offsetBy( 0, 1, 0)),
625  grid.getValue(ijk.offsetBy( 0,-1, 0)), grid.getValue(ijk.offsetBy( 0,-2, 0)) );
626  }
627 
628  template<typename Accessor>
629  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
630  {
631 
632  return difference(
633  grid.getValue(ijk.offsetBy( 0, 0, 2)), grid.getValue(ijk.offsetBy( 0, 0, 1)),
634  grid.getValue(ijk.offsetBy( 0, 0,-1)), grid.getValue(ijk.offsetBy( 0, 0,-2)) );
635  }
636 
637 
638  // stencil access version
639  template<typename Stencil>
640  static typename Stencil::ValueType inX(const Stencil& S)
641  {
642  return difference( S.template getValue< 2, 0, 0>(),
643  S.template getValue< 1, 0, 0>(),
644  S.template getValue<-1, 0, 0>(),
645  S.template getValue<-2, 0, 0>() );
646  }
647 
648  template<typename Stencil>
649  static typename Stencil::ValueType inY(const Stencil& S)
650  {
651  return difference( S.template getValue< 0, 2, 0>(),
652  S.template getValue< 0, 1, 0>(),
653  S.template getValue< 0,-1, 0>(),
654  S.template getValue< 0,-2, 0>() );
655  }
656 
657  template<typename Stencil>
658  static typename Stencil::ValueType inZ(const Stencil& S)
659  {
660  return difference( S.template getValue< 0, 0, 2>(),
661  S.template getValue< 0, 0, 1>(),
662  S.template getValue< 0, 0,-1>(),
663  S.template getValue< 0, 0,-2>() );
664  }
665 };
666 
667 template<>
668 struct D1<CD_6TH>
669 {
670 
671  // the difference opperator
672  template <typename ValueType>
673  static ValueType difference( const ValueType& xp3, const ValueType& xp2, const ValueType& xp1,
674  const ValueType& xm1, const ValueType& xm2, const ValueType& xm3 )
675  {
676  return ValueType(3./4.)*(xp1 - xm1) - ValueType(0.15)*(xp2 - xm2)
677  + ValueType(1./60.)*(xp3-xm3);
678  }
679 
680 
681  // random access version
682  template<typename Accessor>
683  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
684  {
685  return difference(
686  grid.getValue(ijk.offsetBy( 3,0,0)), grid.getValue(ijk.offsetBy( 2,0,0)),
687  grid.getValue(ijk.offsetBy( 1,0,0)), grid.getValue(ijk.offsetBy(-1,0,0)),
688  grid.getValue(ijk.offsetBy(-2,0,0)), grid.getValue(ijk.offsetBy(-3,0,0)));
689  }
690 
691  template<typename Accessor>
692  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
693  {
694  return difference(
695  grid.getValue(ijk.offsetBy( 0, 3, 0)), grid.getValue(ijk.offsetBy( 0, 2, 0)),
696  grid.getValue(ijk.offsetBy( 0, 1, 0)), grid.getValue(ijk.offsetBy( 0,-1, 0)),
697  grid.getValue(ijk.offsetBy( 0,-2, 0)), grid.getValue(ijk.offsetBy( 0,-3, 0)));
698  }
699 
700  template<typename Accessor>
701  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
702  {
703  return difference(
704  grid.getValue(ijk.offsetBy( 0, 0, 3)), grid.getValue(ijk.offsetBy( 0, 0, 2)),
705  grid.getValue(ijk.offsetBy( 0, 0, 1)), grid.getValue(ijk.offsetBy( 0, 0,-1)),
706  grid.getValue(ijk.offsetBy( 0, 0,-2)), grid.getValue(ijk.offsetBy( 0, 0,-3)));
707  }
708 
709  // stencil access version
710  template<typename Stencil>
711  static typename Stencil::ValueType inX(const Stencil& S)
712  {
713  return difference(S.template getValue< 3, 0, 0>(),
714  S.template getValue< 2, 0, 0>(),
715  S.template getValue< 1, 0, 0>(),
716  S.template getValue<-1, 0, 0>(),
717  S.template getValue<-2, 0, 0>(),
718  S.template getValue<-3, 0, 0>());
719  }
720 
721  template<typename Stencil>
722  static typename Stencil::ValueType inY(const Stencil& S)
723  {
724 
725  return difference( S.template getValue< 0, 3, 0>(),
726  S.template getValue< 0, 2, 0>(),
727  S.template getValue< 0, 1, 0>(),
728  S.template getValue< 0,-1, 0>(),
729  S.template getValue< 0,-2, 0>(),
730  S.template getValue< 0,-3, 0>());
731  }
732 
733  template<typename Stencil>
734  static typename Stencil::ValueType inZ(const Stencil& S)
735  {
736 
737  return difference( S.template getValue< 0, 0, 3>(),
738  S.template getValue< 0, 0, 2>(),
739  S.template getValue< 0, 0, 1>(),
740  S.template getValue< 0, 0,-1>(),
741  S.template getValue< 0, 0,-2>(),
742  S.template getValue< 0, 0,-3>());
743  }
744 };
745 
746 
747 template<>
748 struct D1<FD_1ST>
749 {
750 
751  // the difference opperator
752  template <typename ValueType>
753  static ValueType difference(const ValueType& xp1, const ValueType& xp0) {
754  return xp1 - xp0;
755  }
756 
757 
758  // random access version
759  template<typename Accessor>
760  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
761  {
762  return difference(grid.getValue(ijk.offsetBy(1, 0, 0)), grid.getValue(ijk));
763  }
764 
765  template<typename Accessor>
766  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
767  {
768  return difference(grid.getValue(ijk.offsetBy(0, 1, 0)), grid.getValue(ijk));
769  }
770 
771  template<typename Accessor>
772  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
773  {
774  return difference(grid.getValue(ijk.offsetBy(0, 0, 1)), grid.getValue(ijk));
775  }
776 
777  // stencil access version
778  template<typename Stencil>
779  static typename Stencil::ValueType inX(const Stencil& S)
780  {
781  return difference(S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>());
782  }
783 
784  template<typename Stencil>
785  static typename Stencil::ValueType inY(const Stencil& S)
786  {
787  return difference(S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>());
788  }
789 
790  template<typename Stencil>
791  static typename Stencil::ValueType inZ(const Stencil& S)
792  {
793  return difference(S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>());
794  }
795 };
796 
797 
798 template<>
799 struct D1<FD_2ND>
800 {
801  // the difference opperator
802  template <typename ValueType>
803  static ValueType difference(const ValueType& xp2, const ValueType& xp1, const ValueType& xp0)
804  {
805  return ValueType(2)*xp1 -(ValueType(0.5)*xp2 + ValueType(3./2.)*xp0);
806  }
807 
808 
809  // random access version
810  template<typename Accessor>
811  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
812  {
813  return difference(
814  grid.getValue(ijk.offsetBy(2,0,0)),
815  grid.getValue(ijk.offsetBy(1,0,0)),
816  grid.getValue(ijk));
817  }
818 
819  template<typename Accessor>
820  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
821  {
822  return difference(
823  grid.getValue(ijk.offsetBy(0,2,0)),
824  grid.getValue(ijk.offsetBy(0,1,0)),
825  grid.getValue(ijk));
826  }
827 
828  template<typename Accessor>
829  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
830  {
831  return difference(
832  grid.getValue(ijk.offsetBy(0,0,2)),
833  grid.getValue(ijk.offsetBy(0,0,1)),
834  grid.getValue(ijk));
835  }
836 
837 
838  // stencil access version
839  template<typename Stencil>
840  static typename Stencil::ValueType inX(const Stencil& S)
841  {
842  return difference( S.template getValue< 2, 0, 0>(),
843  S.template getValue< 1, 0, 0>(),
844  S.template getValue< 0, 0, 0>() );
845  }
846 
847  template<typename Stencil>
848  static typename Stencil::ValueType inY(const Stencil& S)
849  {
850  return difference( S.template getValue< 0, 2, 0>(),
851  S.template getValue< 0, 1, 0>(),
852  S.template getValue< 0, 0, 0>() );
853  }
854 
855  template<typename Stencil>
856  static typename Stencil::ValueType inZ(const Stencil& S)
857  {
858  return difference( S.template getValue< 0, 0, 2>(),
859  S.template getValue< 0, 0, 1>(),
860  S.template getValue< 0, 0, 0>() );
861  }
862 
863 };
864 
865 
866 template<>
867 struct D1<FD_3RD>
868 {
869 
870  // the difference opperator
871  template<typename ValueType>
872  static ValueType difference(const ValueType& xp3, const ValueType& xp2,
873  const ValueType& xp1, const ValueType& xp0)
874  {
875  return static_cast<ValueType>(xp3/3.0 - 1.5*xp2 + 3.0*xp1 - 11.0*xp0/6.0);
876  }
877 
878 
879  // random access version
880  template<typename Accessor>
881  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
882  {
883  return difference( grid.getValue(ijk.offsetBy(3,0,0)),
884  grid.getValue(ijk.offsetBy(2,0,0)),
885  grid.getValue(ijk.offsetBy(1,0,0)),
886  grid.getValue(ijk) );
887  }
888 
889  template<typename Accessor>
890  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
891  {
892  return difference( grid.getValue(ijk.offsetBy(0,3,0)),
893  grid.getValue(ijk.offsetBy(0,2,0)),
894  grid.getValue(ijk.offsetBy(0,1,0)),
895  grid.getValue(ijk) );
896  }
897 
898  template<typename Accessor>
899  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
900  {
901  return difference( grid.getValue(ijk.offsetBy(0,0,3)),
902  grid.getValue(ijk.offsetBy(0,0,2)),
903  grid.getValue(ijk.offsetBy(0,0,1)),
904  grid.getValue(ijk) );
905  }
906 
907 
908  // stencil access version
909  template<typename Stencil>
910  static typename Stencil::ValueType inX(const Stencil& S)
911  {
912  return difference(S.template getValue< 3, 0, 0>(),
913  S.template getValue< 2, 0, 0>(),
914  S.template getValue< 1, 0, 0>(),
915  S.template getValue< 0, 0, 0>() );
916  }
917 
918  template<typename Stencil>
919  static typename Stencil::ValueType inY(const Stencil& S)
920  {
921  return difference(S.template getValue< 0, 3, 0>(),
922  S.template getValue< 0, 2, 0>(),
923  S.template getValue< 0, 1, 0>(),
924  S.template getValue< 0, 0, 0>() );
925  }
926 
927  template<typename Stencil>
928  static typename Stencil::ValueType inZ(const Stencil& S)
929  {
930  return difference( S.template getValue< 0, 0, 3>(),
931  S.template getValue< 0, 0, 2>(),
932  S.template getValue< 0, 0, 1>(),
933  S.template getValue< 0, 0, 0>() );
934  }
935 };
936 
937 
938 template<>
939 struct D1<BD_1ST>
940 {
941 
942  // the difference opperator
943  template <typename ValueType>
944  static ValueType difference(const ValueType& xm1, const ValueType& xm0) {
945  return -D1<FD_1ST>::difference(xm1, xm0);
946  }
947 
948 
949  // random access version
950  template<typename Accessor>
951  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
952  {
953  return difference(grid.getValue(ijk.offsetBy(-1,0,0)), grid.getValue(ijk));
954  }
955 
956  template<typename Accessor>
957  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
958  {
959  return difference(grid.getValue(ijk.offsetBy(0,-1,0)), grid.getValue(ijk));
960  }
961 
962  template<typename Accessor>
963  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
964  {
965  return difference(grid.getValue(ijk.offsetBy(0, 0,-1)), grid.getValue(ijk));
966  }
967 
968 
969  // stencil access version
970  template<typename Stencil>
971  static typename Stencil::ValueType inX(const Stencil& S)
972  {
973  return difference(S.template getValue<-1, 0, 0>(), S.template getValue< 0, 0, 0>());
974  }
975 
976  template<typename Stencil>
977  static typename Stencil::ValueType inY(const Stencil& S)
978  {
979  return difference(S.template getValue< 0,-1, 0>(), S.template getValue< 0, 0, 0>());
980  }
981 
982  template<typename Stencil>
983  static typename Stencil::ValueType inZ(const Stencil& S)
984  {
985  return difference(S.template getValue< 0, 0,-1>(), S.template getValue< 0, 0, 0>());
986  }
987 };
988 
989 
990 template<>
991 struct D1<BD_2ND>
992 {
993 
994  // the difference opperator
995  template <typename ValueType>
996  static ValueType difference(const ValueType& xm2, const ValueType& xm1, const ValueType& xm0)
997  {
998  return -D1<FD_2ND>::difference(xm2, xm1, xm0);
999  }
1000 
1001 
1002  // random access version
1003  template<typename Accessor>
1004  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1005  {
1006  return difference( grid.getValue(ijk.offsetBy(-2,0,0)),
1007  grid.getValue(ijk.offsetBy(-1,0,0)),
1008  grid.getValue(ijk) );
1009  }
1010 
1011  template<typename Accessor>
1012  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1013  {
1014  return difference( grid.getValue(ijk.offsetBy(0,-2,0)),
1015  grid.getValue(ijk.offsetBy(0,-1,0)),
1016  grid.getValue(ijk) );
1017  }
1018 
1019  template<typename Accessor>
1020  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1021  {
1022  return difference( grid.getValue(ijk.offsetBy(0,0,-2)),
1023  grid.getValue(ijk.offsetBy(0,0,-1)),
1024  grid.getValue(ijk) );
1025  }
1026 
1027  // stencil access version
1028  template<typename Stencil>
1029  static typename Stencil::ValueType inX(const Stencil& S)
1030  {
1031  return difference( S.template getValue<-2, 0, 0>(),
1032  S.template getValue<-1, 0, 0>(),
1033  S.template getValue< 0, 0, 0>() );
1034  }
1035 
1036  template<typename Stencil>
1037  static typename Stencil::ValueType inY(const Stencil& S)
1038  {
1039  return difference( S.template getValue< 0,-2, 0>(),
1040  S.template getValue< 0,-1, 0>(),
1041  S.template getValue< 0, 0, 0>() );
1042  }
1043 
1044  template<typename Stencil>
1045  static typename Stencil::ValueType inZ(const Stencil& S)
1046  {
1047  return difference( S.template getValue< 0, 0,-2>(),
1048  S.template getValue< 0, 0,-1>(),
1049  S.template getValue< 0, 0, 0>() );
1050  }
1051 };
1052 
1053 
1054 template<>
1055 struct D1<BD_3RD>
1056 {
1057 
1058  // the difference opperator
1059  template <typename ValueType>
1060  static ValueType difference(const ValueType& xm3, const ValueType& xm2,
1061  const ValueType& xm1, const ValueType& xm0)
1062  {
1063  return -D1<FD_3RD>::difference(xm3, xm2, xm1, xm0);
1064  }
1065 
1066  // random access version
1067  template<typename Accessor>
1068  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1069  {
1070  return difference( grid.getValue(ijk.offsetBy(-3,0,0)),
1071  grid.getValue(ijk.offsetBy(-2,0,0)),
1072  grid.getValue(ijk.offsetBy(-1,0,0)),
1073  grid.getValue(ijk) );
1074  }
1075 
1076  template<typename Accessor>
1077  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1078  {
1079  return difference( grid.getValue(ijk.offsetBy( 0,-3,0)),
1080  grid.getValue(ijk.offsetBy( 0,-2,0)),
1081  grid.getValue(ijk.offsetBy( 0,-1,0)),
1082  grid.getValue(ijk) );
1083  }
1084 
1085  template<typename Accessor>
1086  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1087  {
1088  return difference( grid.getValue(ijk.offsetBy( 0, 0,-3)),
1089  grid.getValue(ijk.offsetBy( 0, 0,-2)),
1090  grid.getValue(ijk.offsetBy( 0, 0,-1)),
1091  grid.getValue(ijk) );
1092  }
1093 
1094  // stencil access version
1095  template<typename Stencil>
1096  static typename Stencil::ValueType inX(const Stencil& S)
1097  {
1098  return difference( S.template getValue<-3, 0, 0>(),
1099  S.template getValue<-2, 0, 0>(),
1100  S.template getValue<-1, 0, 0>(),
1101  S.template getValue< 0, 0, 0>() );
1102  }
1103 
1104  template<typename Stencil>
1105  static typename Stencil::ValueType inY(const Stencil& S)
1106  {
1107  return difference( S.template getValue< 0,-3, 0>(),
1108  S.template getValue< 0,-2, 0>(),
1109  S.template getValue< 0,-1, 0>(),
1110  S.template getValue< 0, 0, 0>() );
1111  }
1112 
1113  template<typename Stencil>
1114  static typename Stencil::ValueType inZ(const Stencil& S)
1115  {
1116  return difference( S.template getValue< 0, 0,-3>(),
1117  S.template getValue< 0, 0,-2>(),
1118  S.template getValue< 0, 0,-1>(),
1119  S.template getValue< 0, 0, 0>() );
1120  }
1121 
1122 };
1123 
1124 template<>
1125 struct D1<FD_WENO5>
1126 {
1127  // the difference operator
1128  template <typename ValueType>
1129  static ValueType difference(const ValueType& xp3, const ValueType& xp2,
1130  const ValueType& xp1, const ValueType& xp0,
1131  const ValueType& xm1, const ValueType& xm2) {
1132  return WENO5<ValueType>(xp3, xp2, xp1, xp0, xm1)
1133  - WENO5<ValueType>(xp2, xp1, xp0, xm1, xm2);
1134  }
1135 
1136 
1137  // random access version
1138  template<typename Accessor>
1139  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1140  {
1141  typedef typename Accessor::ValueType ValueType;
1142  ValueType V[6];
1143  V[0] = grid.getValue(ijk.offsetBy(3,0,0));
1144  V[1] = grid.getValue(ijk.offsetBy(2,0,0));
1145  V[2] = grid.getValue(ijk.offsetBy(1,0,0));
1146  V[3] = grid.getValue(ijk);
1147  V[4] = grid.getValue(ijk.offsetBy(-1,0,0));
1148  V[5] = grid.getValue(ijk.offsetBy(-2,0,0));
1149 
1150  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1151  }
1152 
1153  template<typename Accessor>
1154  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1155  {
1156  typedef typename Accessor::ValueType ValueType;
1157  ValueType V[6];
1158  V[0] = grid.getValue(ijk.offsetBy(0,3,0));
1159  V[1] = grid.getValue(ijk.offsetBy(0,2,0));
1160  V[2] = grid.getValue(ijk.offsetBy(0,1,0));
1161  V[3] = grid.getValue(ijk);
1162  V[4] = grid.getValue(ijk.offsetBy(0,-1,0));
1163  V[5] = grid.getValue(ijk.offsetBy(0,-2,0));
1164 
1165  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1166  }
1167 
1168  template<typename Accessor>
1169  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1170  {
1171  typedef typename Accessor::ValueType ValueType;
1172  ValueType V[6];
1173  V[0] = grid.getValue(ijk.offsetBy(0,0,3));
1174  V[1] = grid.getValue(ijk.offsetBy(0,0,2));
1175  V[2] = grid.getValue(ijk.offsetBy(0,0,1));
1176  V[3] = grid.getValue(ijk);
1177  V[4] = grid.getValue(ijk.offsetBy(0,0,-1));
1178  V[5] = grid.getValue(ijk.offsetBy(0,0,-2));
1179 
1180  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1181  }
1182 
1183  // stencil access version
1184  template<typename Stencil>
1185  static typename Stencil::ValueType inX(const Stencil& S)
1186  {
1187 
1188  return static_cast<typename Stencil::ValueType>(difference(
1189  S.template getValue< 3, 0, 0>(),
1190  S.template getValue< 2, 0, 0>(),
1191  S.template getValue< 1, 0, 0>(),
1192  S.template getValue< 0, 0, 0>(),
1193  S.template getValue<-1, 0, 0>(),
1194  S.template getValue<-2, 0, 0>() ));
1195 
1196  }
1197 
1198  template<typename Stencil>
1199  static typename Stencil::ValueType inY(const Stencil& S)
1200  {
1201  return static_cast<typename Stencil::ValueType>(difference(
1202  S.template getValue< 0, 3, 0>(),
1203  S.template getValue< 0, 2, 0>(),
1204  S.template getValue< 0, 1, 0>(),
1205  S.template getValue< 0, 0, 0>(),
1206  S.template getValue< 0,-1, 0>(),
1207  S.template getValue< 0,-2, 0>() ));
1208  }
1209 
1210  template<typename Stencil>
1211  static typename Stencil::ValueType inZ(const Stencil& S)
1212  {
1213  return static_cast<typename Stencil::ValueType>(difference(
1214  S.template getValue< 0, 0, 3>(),
1215  S.template getValue< 0, 0, 2>(),
1216  S.template getValue< 0, 0, 1>(),
1217  S.template getValue< 0, 0, 0>(),
1218  S.template getValue< 0, 0,-1>(),
1219  S.template getValue< 0, 0,-2>() ));
1220  }
1221 };
1222 
1223 template<>
1225 {
1226 
1227  // the difference opperator
1228  template <typename ValueType>
1229  static ValueType difference(const ValueType& xp3, const ValueType& xp2,
1230  const ValueType& xp1, const ValueType& xp0,
1231  const ValueType& xm1, const ValueType& xm2) {
1232  return WENO5<ValueType>(xp3 - xp2, xp2 - xp1, xp1 - xp0, xp0-xm1, xm1-xm2);
1233  }
1234 
1235  // random access version
1236  template<typename Accessor>
1237  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1238  {
1239  typedef typename Accessor::ValueType ValueType;
1240  ValueType V[6];
1241  V[0] = grid.getValue(ijk.offsetBy(3,0,0));
1242  V[1] = grid.getValue(ijk.offsetBy(2,0,0));
1243  V[2] = grid.getValue(ijk.offsetBy(1,0,0));
1244  V[3] = grid.getValue(ijk);
1245  V[4] = grid.getValue(ijk.offsetBy(-1,0,0));
1246  V[5] = grid.getValue(ijk.offsetBy(-2,0,0));
1247 
1248  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1249 
1250  }
1251 
1252  template<typename Accessor>
1253  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1254  {
1255  typedef typename Accessor::ValueType ValueType;
1256  ValueType V[6];
1257  V[0] = grid.getValue(ijk.offsetBy(0,3,0));
1258  V[1] = grid.getValue(ijk.offsetBy(0,2,0));
1259  V[2] = grid.getValue(ijk.offsetBy(0,1,0));
1260  V[3] = grid.getValue(ijk);
1261  V[4] = grid.getValue(ijk.offsetBy(0,-1,0));
1262  V[5] = grid.getValue(ijk.offsetBy(0,-2,0));
1263 
1264  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1265  }
1266 
1267  template<typename Accessor>
1268  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1269  {
1270  typedef typename Accessor::ValueType ValueType;
1271  ValueType V[6];
1272  V[0] = grid.getValue(ijk.offsetBy(0,0,3));
1273  V[1] = grid.getValue(ijk.offsetBy(0,0,2));
1274  V[2] = grid.getValue(ijk.offsetBy(0,0,1));
1275  V[3] = grid.getValue(ijk);
1276  V[4] = grid.getValue(ijk.offsetBy(0,0,-1));
1277  V[5] = grid.getValue(ijk.offsetBy(0,0,-2));
1278 
1279  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1280  }
1281 
1282  // stencil access version
1283  template<typename Stencil>
1284  static typename Stencil::ValueType inX(const Stencil& S)
1285  {
1286 
1287  return difference( S.template getValue< 3, 0, 0>(),
1288  S.template getValue< 2, 0, 0>(),
1289  S.template getValue< 1, 0, 0>(),
1290  S.template getValue< 0, 0, 0>(),
1291  S.template getValue<-1, 0, 0>(),
1292  S.template getValue<-2, 0, 0>() );
1293 
1294  }
1295 
1296  template<typename Stencil>
1297  static typename Stencil::ValueType inY(const Stencil& S)
1298  {
1299  return difference( S.template getValue< 0, 3, 0>(),
1300  S.template getValue< 0, 2, 0>(),
1301  S.template getValue< 0, 1, 0>(),
1302  S.template getValue< 0, 0, 0>(),
1303  S.template getValue< 0,-1, 0>(),
1304  S.template getValue< 0,-2, 0>() );
1305  }
1306 
1307  template<typename Stencil>
1308  static typename Stencil::ValueType inZ(const Stencil& S)
1309  {
1310 
1311  return difference( S.template getValue< 0, 0, 3>(),
1312  S.template getValue< 0, 0, 2>(),
1313  S.template getValue< 0, 0, 1>(),
1314  S.template getValue< 0, 0, 0>(),
1315  S.template getValue< 0, 0,-1>(),
1316  S.template getValue< 0, 0,-2>() );
1317  }
1318 
1319 };
1320 
1321 template<>
1322 struct D1<BD_WENO5>
1323 {
1324 
1325  template<typename ValueType>
1326  static ValueType difference(const ValueType& xm3, const ValueType& xm2, const ValueType& xm1,
1327  const ValueType& xm0, const ValueType& xp1, const ValueType& xp2)
1328  {
1329  return -D1<FD_WENO5>::difference(xm3, xm2, xm1, xm0, xp1, xp2);
1330  }
1331 
1332 
1333  // random access version
1334  template<typename Accessor>
1335  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1336  {
1337  typedef typename Accessor::ValueType ValueType;
1338  ValueType V[6];
1339  V[0] = grid.getValue(ijk.offsetBy(-3,0,0));
1340  V[1] = grid.getValue(ijk.offsetBy(-2,0,0));
1341  V[2] = grid.getValue(ijk.offsetBy(-1,0,0));
1342  V[3] = grid.getValue(ijk);
1343  V[4] = grid.getValue(ijk.offsetBy(1,0,0));
1344  V[5] = grid.getValue(ijk.offsetBy(2,0,0));
1345 
1346  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1347  }
1348 
1349  template<typename Accessor>
1350  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1351  {
1352  typedef typename Accessor::ValueType ValueType;
1353  ValueType V[6];
1354  V[0] = grid.getValue(ijk.offsetBy(0,-3,0));
1355  V[1] = grid.getValue(ijk.offsetBy(0,-2,0));
1356  V[2] = grid.getValue(ijk.offsetBy(0,-1,0));
1357  V[3] = grid.getValue(ijk);
1358  V[4] = grid.getValue(ijk.offsetBy(0,1,0));
1359  V[5] = grid.getValue(ijk.offsetBy(0,2,0));
1360 
1361  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1362  }
1363 
1364  template<typename Accessor>
1365  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1366  {
1367  typedef typename Accessor::ValueType ValueType;
1368  ValueType V[6];
1369  V[0] = grid.getValue(ijk.offsetBy(0,0,-3));
1370  V[1] = grid.getValue(ijk.offsetBy(0,0,-2));
1371  V[2] = grid.getValue(ijk.offsetBy(0,0,-1));
1372  V[3] = grid.getValue(ijk);
1373  V[4] = grid.getValue(ijk.offsetBy(0,0,1));
1374  V[5] = grid.getValue(ijk.offsetBy(0,0,2));
1375 
1376  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1377  }
1378 
1379  // stencil access version
1380  template<typename Stencil>
1381  static typename Stencil::ValueType inX(const Stencil& S)
1382  {
1383  typedef typename Stencil::ValueType ValueType;
1384  ValueType V[6];
1385  V[0] = S.template getValue<-3, 0, 0>();
1386  V[1] = S.template getValue<-2, 0, 0>();
1387  V[2] = S.template getValue<-1, 0, 0>();
1388  V[3] = S.template getValue< 0, 0, 0>();
1389  V[4] = S.template getValue< 1, 0, 0>();
1390  V[5] = S.template getValue< 2, 0, 0>();
1391 
1392  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1393  }
1394 
1395  template<typename Stencil>
1396  static typename Stencil::ValueType inY(const Stencil& S)
1397  {
1398  typedef typename Stencil::ValueType ValueType;
1399  ValueType V[6];
1400  V[0] = S.template getValue< 0,-3, 0>();
1401  V[1] = S.template getValue< 0,-2, 0>();
1402  V[2] = S.template getValue< 0,-1, 0>();
1403  V[3] = S.template getValue< 0, 0, 0>();
1404  V[4] = S.template getValue< 0, 1, 0>();
1405  V[5] = S.template getValue< 0, 2, 0>();
1406 
1407  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1408  }
1409 
1410  template<typename Stencil>
1411  static typename Stencil::ValueType inZ(const Stencil& S)
1412  {
1413  typedef typename Stencil::ValueType ValueType;
1414  ValueType V[6];
1415  V[0] = S.template getValue< 0, 0,-3>();
1416  V[1] = S.template getValue< 0, 0,-2>();
1417  V[2] = S.template getValue< 0, 0,-1>();
1418  V[3] = S.template getValue< 0, 0, 0>();
1419  V[4] = S.template getValue< 0, 0, 1>();
1420  V[5] = S.template getValue< 0, 0, 2>();
1421 
1422  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1423  }
1424 };
1425 
1426 
1427 template<>
1429 {
1430  template<typename ValueType>
1431  static ValueType difference(const ValueType& xm3, const ValueType& xm2, const ValueType& xm1,
1432  const ValueType& xm0, const ValueType& xp1, const ValueType& xp2)
1433  {
1434  return -D1<FD_HJWENO5>::difference(xm3, xm2, xm1, xm0, xp1, xp2);
1435  }
1436 
1437  // random access version
1438  template<typename Accessor>
1439  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1440  {
1441  typedef typename Accessor::ValueType ValueType;
1442  ValueType V[6];
1443  V[0] = grid.getValue(ijk.offsetBy(-3,0,0));
1444  V[1] = grid.getValue(ijk.offsetBy(-2,0,0));
1445  V[2] = grid.getValue(ijk.offsetBy(-1,0,0));
1446  V[3] = grid.getValue(ijk);
1447  V[4] = grid.getValue(ijk.offsetBy(1,0,0));
1448  V[5] = grid.getValue(ijk.offsetBy(2,0,0));
1449 
1450  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1451  }
1452 
1453  template<typename Accessor>
1454  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1455  {
1456  typedef typename Accessor::ValueType ValueType;
1457  ValueType V[6];
1458  V[0] = grid.getValue(ijk.offsetBy(0,-3,0));
1459  V[1] = grid.getValue(ijk.offsetBy(0,-2,0));
1460  V[2] = grid.getValue(ijk.offsetBy(0,-1,0));
1461  V[3] = grid.getValue(ijk);
1462  V[4] = grid.getValue(ijk.offsetBy(0,1,0));
1463  V[5] = grid.getValue(ijk.offsetBy(0,2,0));
1464 
1465  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1466  }
1467 
1468  template<typename Accessor>
1469  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1470  {
1471  typedef typename Accessor::ValueType ValueType;
1472  ValueType V[6];
1473  V[0] = grid.getValue(ijk.offsetBy(0,0,-3));
1474  V[1] = grid.getValue(ijk.offsetBy(0,0,-2));
1475  V[2] = grid.getValue(ijk.offsetBy(0,0,-1));
1476  V[3] = grid.getValue(ijk);
1477  V[4] = grid.getValue(ijk.offsetBy(0,0,1));
1478  V[5] = grid.getValue(ijk.offsetBy(0,0,2));
1479 
1480  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1481  }
1482 
1483  // stencil access version
1484  template<typename Stencil>
1485  static typename Stencil::ValueType inX(const Stencil& S)
1486  {
1487  typedef typename Stencil::ValueType ValueType;
1488  ValueType V[6];
1489  V[0] = S.template getValue<-3, 0, 0>();
1490  V[1] = S.template getValue<-2, 0, 0>();
1491  V[2] = S.template getValue<-1, 0, 0>();
1492  V[3] = S.template getValue< 0, 0, 0>();
1493  V[4] = S.template getValue< 1, 0, 0>();
1494  V[5] = S.template getValue< 2, 0, 0>();
1495 
1496  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1497  }
1498 
1499  template<typename Stencil>
1500  static typename Stencil::ValueType inY(const Stencil& S)
1501  {
1502  typedef typename Stencil::ValueType ValueType;
1503  ValueType V[6];
1504  V[0] = S.template getValue< 0,-3, 0>();
1505  V[1] = S.template getValue< 0,-2, 0>();
1506  V[2] = S.template getValue< 0,-1, 0>();
1507  V[3] = S.template getValue< 0, 0, 0>();
1508  V[4] = S.template getValue< 0, 1, 0>();
1509  V[5] = S.template getValue< 0, 2, 0>();
1510 
1511  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1512  }
1513 
1514  template<typename Stencil>
1515  static typename Stencil::ValueType inZ(const Stencil& S)
1516  {
1517  typedef typename Stencil::ValueType ValueType;
1518  ValueType V[6];
1519  V[0] = S.template getValue< 0, 0,-3>();
1520  V[1] = S.template getValue< 0, 0,-2>();
1521  V[2] = S.template getValue< 0, 0,-1>();
1522  V[3] = S.template getValue< 0, 0, 0>();
1523  V[4] = S.template getValue< 0, 0, 1>();
1524  V[5] = S.template getValue< 0, 0, 2>();
1525 
1526  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1527  }
1528 };
1529 
1530 
1531 template<DScheme DiffScheme>
1532 struct D1Vec
1533 {
1534  // random access version
1535  template<typename Accessor>
1536  static typename Accessor::ValueType::value_type
1537  inX(const Accessor& grid, const Coord& ijk, int n)
1538  {
1539  return D1<DiffScheme>::inX(grid, ijk)[n];
1540  }
1541 
1542  template<typename Accessor>
1543  static typename Accessor::ValueType::value_type
1544  inY(const Accessor& grid, const Coord& ijk, int n)
1545  {
1546  return D1<DiffScheme>::inY(grid, ijk)[n];
1547  }
1548  template<typename Accessor>
1549  static typename Accessor::ValueType::value_type
1550  inZ(const Accessor& grid, const Coord& ijk, int n)
1551  {
1552  return D1<DiffScheme>::inZ(grid, ijk)[n];
1553  }
1554 
1555 
1556  // stencil access version
1557  template<typename Stencil>
1558  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1559  {
1560  return D1<DiffScheme>::inX(S)[n];
1561  }
1562 
1563  template<typename Stencil>
1564  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1565  {
1566  return D1<DiffScheme>::inY(S)[n];
1567  }
1568 
1569  template<typename Stencil>
1570  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1571  {
1572  return D1<DiffScheme>::inZ(S)[n];
1573  }
1574 };
1575 
1576 
1577 template<>
1579 {
1580 
1581  // random access version
1582  template<typename Accessor>
1583  static typename Accessor::ValueType::value_type
1584  inX(const Accessor& grid, const Coord& ijk, int n)
1585  {
1586  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy( 1, 0, 0))[n],
1587  grid.getValue(ijk.offsetBy(-1, 0, 0))[n] );
1588  }
1589 
1590  template<typename Accessor>
1591  static typename Accessor::ValueType::value_type
1592  inY(const Accessor& grid, const Coord& ijk, int n)
1593  {
1594  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy(0, 1, 0))[n],
1595  grid.getValue(ijk.offsetBy(0,-1, 0))[n] );
1596  }
1597 
1598  template<typename Accessor>
1599  static typename Accessor::ValueType::value_type
1600  inZ(const Accessor& grid, const Coord& ijk, int n)
1601  {
1602  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy(0, 0, 1))[n],
1603  grid.getValue(ijk.offsetBy(0, 0,-1))[n] );
1604  }
1605 
1606  // stencil access version
1607  template<typename Stencil>
1608  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1609  {
1610  return D1<CD_2NDT>::difference( S.template getValue< 1, 0, 0>()[n],
1611  S.template getValue<-1, 0, 0>()[n] );
1612  }
1613 
1614  template<typename Stencil>
1615  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1616  {
1617  return D1<CD_2NDT>::difference( S.template getValue< 0, 1, 0>()[n],
1618  S.template getValue< 0,-1, 0>()[n] );
1619  }
1620 
1621  template<typename Stencil>
1622  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1623  {
1624  return D1<CD_2NDT>::difference( S.template getValue< 0, 0, 1>()[n],
1625  S.template getValue< 0, 0,-1>()[n] );
1626  }
1627 };
1628 
1629 template<>
1630 struct D1Vec<CD_2ND>
1631 {
1632 
1633  // random access version
1634  template<typename Accessor>
1635  static typename Accessor::ValueType::value_type
1636  inX(const Accessor& grid, const Coord& ijk, int n)
1637  {
1638  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy( 1, 0, 0))[n] ,
1639  grid.getValue(ijk.offsetBy(-1, 0, 0))[n] );
1640  }
1641 
1642  template<typename Accessor>
1643  static typename Accessor::ValueType::value_type
1644  inY(const Accessor& grid, const Coord& ijk, int n)
1645  {
1646  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy(0, 1, 0))[n] ,
1647  grid.getValue(ijk.offsetBy(0,-1, 0))[n] );
1648  }
1649 
1650  template<typename Accessor>
1651  static typename Accessor::ValueType::value_type
1652  inZ(const Accessor& grid, const Coord& ijk, int n)
1653  {
1654  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy(0, 0, 1))[n] ,
1655  grid.getValue(ijk.offsetBy(0, 0,-1))[n] );
1656  }
1657 
1658 
1659  // stencil access version
1660  template<typename Stencil>
1661  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1662  {
1663  return D1<CD_2ND>::difference( S.template getValue< 1, 0, 0>()[n],
1664  S.template getValue<-1, 0, 0>()[n] );
1665  }
1666 
1667  template<typename Stencil>
1668  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1669  {
1670  return D1<CD_2ND>::difference( S.template getValue< 0, 1, 0>()[n],
1671  S.template getValue< 0,-1, 0>()[n] );
1672  }
1673 
1674  template<typename Stencil>
1675  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1676  {
1677  return D1<CD_2ND>::difference( S.template getValue< 0, 0, 1>()[n],
1678  S.template getValue< 0, 0,-1>()[n] );
1679  }
1680 };
1681 
1682 
1683 template<>
1684 struct D1Vec<CD_4TH> {
1685  // typedef typename Accessor::ValueType::value_type value_type;
1686 
1687 
1688  // random access version
1689  template<typename Accessor>
1690  static typename Accessor::ValueType::value_type
1691  inX(const Accessor& grid, const Coord& ijk, int n)
1692  {
1693  return D1<CD_4TH>::difference(
1694  grid.getValue(ijk.offsetBy(2, 0, 0))[n], grid.getValue(ijk.offsetBy( 1, 0, 0))[n],
1695  grid.getValue(ijk.offsetBy(-1,0, 0))[n], grid.getValue(ijk.offsetBy(-2, 0, 0))[n]);
1696  }
1697 
1698  template<typename Accessor>
1699  static typename Accessor::ValueType::value_type
1700  inY(const Accessor& grid, const Coord& ijk, int n)
1701  {
1702  return D1<CD_4TH>::difference(
1703  grid.getValue(ijk.offsetBy( 0, 2, 0))[n], grid.getValue(ijk.offsetBy( 0, 1, 0))[n],
1704  grid.getValue(ijk.offsetBy( 0,-1, 0))[n], grid.getValue(ijk.offsetBy( 0,-2, 0))[n]);
1705  }
1706 
1707  template<typename Accessor>
1708  static typename Accessor::ValueType::value_type
1709  inZ(const Accessor& grid, const Coord& ijk, int n)
1710  {
1711  return D1<CD_4TH>::difference(
1712  grid.getValue(ijk.offsetBy(0,0, 2))[n], grid.getValue(ijk.offsetBy( 0, 0, 1))[n],
1713  grid.getValue(ijk.offsetBy(0,0,-1))[n], grid.getValue(ijk.offsetBy( 0, 0,-2))[n]);
1714  }
1715 
1716  // stencil access version
1717  template<typename Stencil>
1718  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1719  {
1720  return D1<CD_4TH>::difference(
1721  S.template getValue< 2, 0, 0>()[n], S.template getValue< 1, 0, 0>()[n],
1722  S.template getValue<-1, 0, 0>()[n], S.template getValue<-2, 0, 0>()[n] );
1723  }
1724 
1725  template<typename Stencil>
1726  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1727  {
1728  return D1<CD_4TH>::difference(
1729  S.template getValue< 0, 2, 0>()[n], S.template getValue< 0, 1, 0>()[n],
1730  S.template getValue< 0,-1, 0>()[n], S.template getValue< 0,-2, 0>()[n]);
1731  }
1732 
1733  template<typename Stencil>
1734  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1735  {
1736  return D1<CD_4TH>::difference(
1737  S.template getValue< 0, 0, 2>()[n], S.template getValue< 0, 0, 1>()[n],
1738  S.template getValue< 0, 0,-1>()[n], S.template getValue< 0, 0,-2>()[n]);
1739  }
1740 };
1741 
1742 
1743 template<>
1744 struct D1Vec<CD_6TH>
1745 {
1746  //typedef typename Accessor::ValueType::value_type::value_type ValueType;
1747 
1748  // random access version
1749  template<typename Accessor>
1750  static typename Accessor::ValueType::value_type
1751  inX(const Accessor& grid, const Coord& ijk, int n)
1752  {
1753  return D1<CD_6TH>::difference(
1754  grid.getValue(ijk.offsetBy( 3, 0, 0))[n], grid.getValue(ijk.offsetBy( 2, 0, 0))[n],
1755  grid.getValue(ijk.offsetBy( 1, 0, 0))[n], grid.getValue(ijk.offsetBy(-1, 0, 0))[n],
1756  grid.getValue(ijk.offsetBy(-2, 0, 0))[n], grid.getValue(ijk.offsetBy(-3, 0, 0))[n] );
1757  }
1758 
1759  template<typename Accessor>
1760  static typename Accessor::ValueType::value_type
1761  inY(const Accessor& grid, const Coord& ijk, int n)
1762  {
1763  return D1<CD_6TH>::difference(
1764  grid.getValue(ijk.offsetBy( 0, 3, 0))[n], grid.getValue(ijk.offsetBy( 0, 2, 0))[n],
1765  grid.getValue(ijk.offsetBy( 0, 1, 0))[n], grid.getValue(ijk.offsetBy( 0,-1, 0))[n],
1766  grid.getValue(ijk.offsetBy( 0,-2, 0))[n], grid.getValue(ijk.offsetBy( 0,-3, 0))[n] );
1767  }
1768 
1769  template<typename Accessor>
1770  static typename Accessor::ValueType::value_type
1771  inZ(const Accessor& grid, const Coord& ijk, int n)
1772  {
1773  return D1<CD_6TH>::difference(
1774  grid.getValue(ijk.offsetBy( 0, 0, 3))[n], grid.getValue(ijk.offsetBy( 0, 0, 2))[n],
1775  grid.getValue(ijk.offsetBy( 0, 0, 1))[n], grid.getValue(ijk.offsetBy( 0, 0,-1))[n],
1776  grid.getValue(ijk.offsetBy( 0, 0,-2))[n], grid.getValue(ijk.offsetBy( 0, 0,-3))[n] );
1777  }
1778 
1779 
1780  // stencil access version
1781  template<typename Stencil>
1782  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1783  {
1784  return D1<CD_6TH>::difference(
1785  S.template getValue< 3, 0, 0>()[n], S.template getValue< 2, 0, 0>()[n],
1786  S.template getValue< 1, 0, 0>()[n], S.template getValue<-1, 0, 0>()[n],
1787  S.template getValue<-2, 0, 0>()[n], S.template getValue<-3, 0, 0>()[n] );
1788  }
1789 
1790  template<typename Stencil>
1791  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1792  {
1793  return D1<CD_6TH>::difference(
1794  S.template getValue< 0, 3, 0>()[n], S.template getValue< 0, 2, 0>()[n],
1795  S.template getValue< 0, 1, 0>()[n], S.template getValue< 0,-1, 0>()[n],
1796  S.template getValue< 0,-2, 0>()[n], S.template getValue< 0,-3, 0>()[n] );
1797  }
1798 
1799  template<typename Stencil>
1800  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1801  {
1802  return D1<CD_6TH>::difference(
1803  S.template getValue< 0, 0, 3>()[n], S.template getValue< 0, 0, 2>()[n],
1804  S.template getValue< 0, 0, 1>()[n], S.template getValue< 0, 0,-1>()[n],
1805  S.template getValue< 0, 0,-2>()[n], S.template getValue< 0, 0,-3>()[n] );
1806  }
1807 };
1808 
1809 template<DDScheme DiffScheme>
1810 struct D2
1811 {
1812 
1813  template<typename Accessor>
1814  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk);
1815  template<typename Accessor>
1816  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk);
1817  template<typename Accessor>
1818  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk);
1819 
1820  // cross derivatives
1821  template<typename Accessor>
1822  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk);
1823 
1824  template<typename Accessor>
1825  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk);
1826 
1827  template<typename Accessor>
1828  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk);
1829 
1830 
1831  // stencil access version
1832  template<typename Stencil>
1833  static typename Stencil::ValueType inX(const Stencil& S);
1834  template<typename Stencil>
1835  static typename Stencil::ValueType inY(const Stencil& S);
1836  template<typename Stencil>
1837  static typename Stencil::ValueType inZ(const Stencil& S);
1838 
1839  // cross derivatives
1840  template<typename Stencil>
1841  static typename Stencil::ValueType inXandY(const Stencil& S);
1842 
1843  template<typename Stencil>
1844  static typename Stencil::ValueType inXandZ(const Stencil& S);
1845 
1846  template<typename Stencil>
1847  static typename Stencil::ValueType inYandZ(const Stencil& S);
1848 };
1849 
1850 template<>
1851 struct D2<CD_SECOND>
1852 {
1853 
1854  // the difference opperator
1855  template <typename ValueType>
1856  static ValueType difference(const ValueType& xp1, const ValueType& xp0, const ValueType& xm1)
1857  {
1858  return xp1 + xm1 - ValueType(2)*xp0;
1859  }
1860 
1861  template <typename ValueType>
1862  static ValueType crossdifference(const ValueType& xpyp, const ValueType& xpym,
1863  const ValueType& xmyp, const ValueType& xmym)
1864  {
1865  return ValueType(0.25)*(xpyp + xmym - xpym - xmyp);
1866  }
1867 
1868  // random access version
1869  template<typename Accessor>
1870  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1871  {
1872  return difference( grid.getValue(ijk.offsetBy( 1,0,0)), grid.getValue(ijk),
1873  grid.getValue(ijk.offsetBy(-1,0,0)) );
1874  }
1875 
1876  template<typename Accessor>
1877  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1878  {
1879 
1880  return difference( grid.getValue(ijk.offsetBy(0, 1,0)), grid.getValue(ijk),
1881  grid.getValue(ijk.offsetBy(0,-1,0)) );
1882  }
1883 
1884  template<typename Accessor>
1885  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1886  {
1887  return difference( grid.getValue(ijk.offsetBy( 0,0, 1)), grid.getValue(ijk),
1888  grid.getValue(ijk.offsetBy( 0,0,-1)) );
1889  }
1890 
1891  // cross derivatives
1892  template<typename Accessor>
1893  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
1894  {
1895  return crossdifference(
1896  grid.getValue(ijk.offsetBy(1, 1,0)), grid.getValue(ijk.offsetBy( 1,-1,0)),
1897  grid.getValue(ijk.offsetBy(-1,1,0)), grid.getValue(ijk.offsetBy(-1,-1,0)));
1898 
1899  }
1900 
1901  template<typename Accessor>
1902  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
1903  {
1904  return crossdifference(
1905  grid.getValue(ijk.offsetBy(1,0, 1)), grid.getValue(ijk.offsetBy(1, 0,-1)),
1906  grid.getValue(ijk.offsetBy(-1,0,1)), grid.getValue(ijk.offsetBy(-1,0,-1)) );
1907  }
1908 
1909  template<typename Accessor>
1910  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
1911  {
1912  return crossdifference(
1913  grid.getValue(ijk.offsetBy(0, 1,1)), grid.getValue(ijk.offsetBy(0, 1,-1)),
1914  grid.getValue(ijk.offsetBy(0,-1,1)), grid.getValue(ijk.offsetBy(0,-1,-1)) );
1915  }
1916 
1917 
1918  // stencil access version
1919  template<typename Stencil>
1920  static typename Stencil::ValueType inX(const Stencil& S)
1921  {
1922  return difference( S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>(),
1923  S.template getValue<-1, 0, 0>() );
1924  }
1925 
1926  template<typename Stencil>
1927  static typename Stencil::ValueType inY(const Stencil& S)
1928  {
1929  return difference( S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>(),
1930  S.template getValue< 0,-1, 0>() );
1931  }
1932 
1933  template<typename Stencil>
1934  static typename Stencil::ValueType inZ(const Stencil& S)
1935  {
1936  return difference( S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>(),
1937  S.template getValue< 0, 0,-1>() );
1938  }
1939 
1940  // cross derivatives
1941  template<typename Stencil>
1942  static typename Stencil::ValueType inXandY(const Stencil& S)
1943  {
1944  return crossdifference(S.template getValue< 1, 1, 0>(), S.template getValue< 1,-1, 0>(),
1945  S.template getValue<-1, 1, 0>(), S.template getValue<-1,-1, 0>() );
1946  }
1947 
1948  template<typename Stencil>
1949  static typename Stencil::ValueType inXandZ(const Stencil& S)
1950  {
1951  return crossdifference(S.template getValue< 1, 0, 1>(), S.template getValue< 1, 0,-1>(),
1952  S.template getValue<-1, 0, 1>(), S.template getValue<-1, 0,-1>() );
1953  }
1954 
1955  template<typename Stencil>
1956  static typename Stencil::ValueType inYandZ(const Stencil& S)
1957  {
1958  return crossdifference(S.template getValue< 0, 1, 1>(), S.template getValue< 0, 1,-1>(),
1959  S.template getValue< 0,-1, 1>(), S.template getValue< 0,-1,-1>() );
1960  }
1961 };
1962 
1963 
1964 template<>
1965 struct D2<CD_FOURTH>
1966 {
1967 
1968  // the difference opperator
1969  template <typename ValueType>
1970  static ValueType difference(const ValueType& xp2, const ValueType& xp1, const ValueType& xp0,
1971  const ValueType& xm1, const ValueType& xm2) {
1972  return ValueType(-1./12.)*(xp2 + xm2) + ValueType(4./3.)*(xp1 + xm1) -ValueType(2.5)*xp0;
1973  }
1974 
1975  template <typename ValueType>
1976  static ValueType crossdifference(const ValueType& xp2yp2, const ValueType& xp2yp1,
1977  const ValueType& xp2ym1, const ValueType& xp2ym2,
1978  const ValueType& xp1yp2, const ValueType& xp1yp1,
1979  const ValueType& xp1ym1, const ValueType& xp1ym2,
1980  const ValueType& xm2yp2, const ValueType& xm2yp1,
1981  const ValueType& xm2ym1, const ValueType& xm2ym2,
1982  const ValueType& xm1yp2, const ValueType& xm1yp1,
1983  const ValueType& xm1ym1, const ValueType& xm1ym2 ) {
1984  ValueType tmp1 =
1985  ValueType(2./3.0)*(xp1yp1 - xm1yp1 - xp1ym1 + xm1ym1)-
1986  ValueType(1./12.)*(xp2yp1 - xm2yp1 - xp2ym1 + xm2ym1);
1987  ValueType tmp2 =
1988  ValueType(2./3.0)*(xp1yp2 - xm1yp2 - xp1ym2 + xm1ym2)-
1989  ValueType(1./12.)*(xp2yp2 - xm2yp2 - xp2ym2 + xm2ym2);
1990 
1991  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
1992  }
1993 
1994 
1995 
1996  // random access version
1997  template<typename Accessor>
1998  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1999  {
2000  return difference(
2001  grid.getValue(ijk.offsetBy(2,0,0)), grid.getValue(ijk.offsetBy( 1,0,0)),
2002  grid.getValue(ijk),
2003  grid.getValue(ijk.offsetBy(-1,0,0)), grid.getValue(ijk.offsetBy(-2, 0, 0)));
2004  }
2005 
2006  template<typename Accessor>
2007  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
2008  {
2009  return difference(
2010  grid.getValue(ijk.offsetBy(0, 2,0)), grid.getValue(ijk.offsetBy(0, 1,0)),
2011  grid.getValue(ijk),
2012  grid.getValue(ijk.offsetBy(0,-1,0)), grid.getValue(ijk.offsetBy(0,-2, 0)));
2013  }
2014 
2015  template<typename Accessor>
2016  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
2017  {
2018  return difference(
2019  grid.getValue(ijk.offsetBy(0,0, 2)), grid.getValue(ijk.offsetBy(0, 0,1)),
2020  grid.getValue(ijk),
2021  grid.getValue(ijk.offsetBy(0,0,-1)), grid.getValue(ijk.offsetBy(0,0,-2)));
2022  }
2023 
2024  // cross derivatives
2025  template<typename Accessor>
2026  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
2027  {
2028  typedef typename Accessor::ValueType ValueType;
2029  typename Accessor::ValueType tmp1 =
2030  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 1, 0)) -
2031  D1<CD_4TH>::inX(grid, ijk.offsetBy(0,-1, 0));
2032  typename Accessor::ValueType tmp2 =
2033  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 2, 0)) -
2034  D1<CD_4TH>::inX(grid, ijk.offsetBy(0,-2, 0));
2035  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
2036  }
2037 
2038  template<typename Accessor>
2039  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
2040  {
2041  typedef typename Accessor::ValueType ValueType;
2042  typename Accessor::ValueType tmp1 =
2043  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0, 1)) -
2044  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0,-1));
2045  typename Accessor::ValueType tmp2 =
2046  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0, 2)) -
2047  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0,-2));
2048  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
2049  }
2050 
2051  template<typename Accessor>
2052  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
2053  {
2054  typedef typename Accessor::ValueType ValueType;
2055  typename Accessor::ValueType tmp1 =
2056  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0, 1)) -
2057  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0,-1));
2058  typename Accessor::ValueType tmp2 =
2059  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0, 2)) -
2060  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0,-2));
2061  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
2062  }
2063 
2064 
2065  // stencil access version
2066  template<typename Stencil>
2067  static typename Stencil::ValueType inX(const Stencil& S)
2068  {
2069  return difference(S.template getValue< 2, 0, 0>(), S.template getValue< 1, 0, 0>(),
2070  S.template getValue< 0, 0, 0>(),
2071  S.template getValue<-1, 0, 0>(), S.template getValue<-2, 0, 0>() );
2072  }
2073 
2074  template<typename Stencil>
2075  static typename Stencil::ValueType inY(const Stencil& S)
2076  {
2077  return difference(S.template getValue< 0, 2, 0>(), S.template getValue< 0, 1, 0>(),
2078  S.template getValue< 0, 0, 0>(),
2079  S.template getValue< 0,-1, 0>(), S.template getValue< 0,-2, 0>() );
2080  }
2081 
2082  template<typename Stencil>
2083  static typename Stencil::ValueType inZ(const Stencil& S)
2084  {
2085  return difference(S.template getValue< 0, 0, 2>(), S.template getValue< 0, 0, 1>(),
2086  S.template getValue< 0, 0, 0>(),
2087  S.template getValue< 0, 0,-1>(), S.template getValue< 0, 0,-2>() );
2088  }
2089 
2090  // cross derivatives
2091  template<typename Stencil>
2092  static typename Stencil::ValueType inXandY(const Stencil& S)
2093  {
2094  return crossdifference(
2095  S.template getValue< 2, 2, 0>(), S.template getValue< 2, 1, 0>(),
2096  S.template getValue< 2,-1, 0>(), S.template getValue< 2,-2, 0>(),
2097  S.template getValue< 1, 2, 0>(), S.template getValue< 1, 1, 0>(),
2098  S.template getValue< 1,-1, 0>(), S.template getValue< 1,-2, 0>(),
2099  S.template getValue<-2, 2, 0>(), S.template getValue<-2, 1, 0>(),
2100  S.template getValue<-2,-1, 0>(), S.template getValue<-2,-2, 0>(),
2101  S.template getValue<-1, 2, 0>(), S.template getValue<-1, 1, 0>(),
2102  S.template getValue<-1,-1, 0>(), S.template getValue<-1,-2, 0>() );
2103  }
2104 
2105  template<typename Stencil>
2106  static typename Stencil::ValueType inXandZ(const Stencil& S)
2107  {
2108  return crossdifference(
2109  S.template getValue< 2, 0, 2>(), S.template getValue< 2, 0, 1>(),
2110  S.template getValue< 2, 0,-1>(), S.template getValue< 2, 0,-2>(),
2111  S.template getValue< 1, 0, 2>(), S.template getValue< 1, 0, 1>(),
2112  S.template getValue< 1, 0,-1>(), S.template getValue< 1, 0,-2>(),
2113  S.template getValue<-2, 0, 2>(), S.template getValue<-2, 0, 1>(),
2114  S.template getValue<-2, 0,-1>(), S.template getValue<-2, 0,-2>(),
2115  S.template getValue<-1, 0, 2>(), S.template getValue<-1, 0, 1>(),
2116  S.template getValue<-1, 0,-1>(), S.template getValue<-1, 0,-2>() );
2117  }
2118 
2119  template<typename Stencil>
2120  static typename Stencil::ValueType inYandZ(const Stencil& S)
2121  {
2122  return crossdifference(
2123  S.template getValue< 0, 2, 2>(), S.template getValue< 0, 2, 1>(),
2124  S.template getValue< 0, 2,-1>(), S.template getValue< 0, 2,-2>(),
2125  S.template getValue< 0, 1, 2>(), S.template getValue< 0, 1, 1>(),
2126  S.template getValue< 0, 1,-1>(), S.template getValue< 0, 1,-2>(),
2127  S.template getValue< 0,-2, 2>(), S.template getValue< 0,-2, 1>(),
2128  S.template getValue< 0,-2,-1>(), S.template getValue< 0,-2,-2>(),
2129  S.template getValue< 0,-1, 2>(), S.template getValue< 0,-1, 1>(),
2130  S.template getValue< 0,-1,-1>(), S.template getValue< 0,-1,-2>() );
2131  }
2132 };
2133 
2134 
2135 template<>
2136 struct D2<CD_SIXTH>
2137 {
2138  // the difference opperator
2139  template <typename ValueType>
2140  static ValueType difference(const ValueType& xp3, const ValueType& xp2, const ValueType& xp1,
2141  const ValueType& xp0,
2142  const ValueType& xm1, const ValueType& xm2, const ValueType& xm3)
2143  {
2144  return ValueType(1./90.)*(xp3 + xm3) - ValueType(3./20.)*(xp2 + xm2)
2145  + ValueType(1.5)*(xp1 + xm1) - ValueType(49./18.)*xp0;
2146  }
2147 
2148  template <typename ValueType>
2149  static ValueType crossdifference( const ValueType& xp1yp1,const ValueType& xm1yp1,
2150  const ValueType& xp1ym1,const ValueType& xm1ym1,
2151  const ValueType& xp2yp1,const ValueType& xm2yp1,
2152  const ValueType& xp2ym1,const ValueType& xm2ym1,
2153  const ValueType& xp3yp1,const ValueType& xm3yp1,
2154  const ValueType& xp3ym1,const ValueType& xm3ym1,
2155  const ValueType& xp1yp2,const ValueType& xm1yp2,
2156  const ValueType& xp1ym2,const ValueType& xm1ym2,
2157  const ValueType& xp2yp2,const ValueType& xm2yp2,
2158  const ValueType& xp2ym2,const ValueType& xm2ym2,
2159  const ValueType& xp3yp2,const ValueType& xm3yp2,
2160  const ValueType& xp3ym2,const ValueType& xm3ym2,
2161  const ValueType& xp1yp3,const ValueType& xm1yp3,
2162  const ValueType& xp1ym3,const ValueType& xm1ym3,
2163  const ValueType& xp2yp3,const ValueType& xm2yp3,
2164  const ValueType& xp2ym3,const ValueType& xm2ym3,
2165  const ValueType& xp3yp3,const ValueType& xm3yp3,
2166  const ValueType& xp3ym3,const ValueType& xm3ym3 )
2167  {
2168  ValueType tmp1 =
2169  ValueType(0.7500)*(xp1yp1 - xm1yp1 - xp1ym1 + xm1ym1) -
2170  ValueType(0.1500)*(xp2yp1 - xm2yp1 - xp2ym1 + xm2ym1) +
2171  ValueType(1./60.)*(xp3yp1 - xm3yp1 - xp3ym1 + xm3ym1);
2172 
2173  ValueType tmp2 =
2174  ValueType(0.7500)*(xp1yp2 - xm1yp2 - xp1ym2 + xm1ym2) -
2175  ValueType(0.1500)*(xp2yp2 - xm2yp2 - xp2ym2 + xm2ym2) +
2176  ValueType(1./60.)*(xp3yp2 - xm3yp2 - xp3ym2 + xm3ym2);
2177 
2178  ValueType tmp3 =
2179  ValueType(0.7500)*(xp1yp3 - xm1yp3 - xp1ym3 + xm1ym3) -
2180  ValueType(0.1500)*(xp2yp3 - xm2yp3 - xp2ym3 + xm2ym3) +
2181  ValueType(1./60.)*(xp3yp3 - xm3yp3 - xp3ym3 + xm3ym3);
2182 
2183  return ValueType(0.75)*tmp1 - ValueType(0.15)*tmp2 + ValueType(1./60)*tmp3;
2184  }
2185 
2186  // random access version
2187 
2188  template<typename Accessor>
2189  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
2190  {
2191  return difference(
2192  grid.getValue(ijk.offsetBy( 3, 0, 0)), grid.getValue(ijk.offsetBy( 2, 0, 0)),
2193  grid.getValue(ijk.offsetBy( 1, 0, 0)), grid.getValue(ijk),
2194  grid.getValue(ijk.offsetBy(-1, 0, 0)), grid.getValue(ijk.offsetBy(-2, 0, 0)),
2195  grid.getValue(ijk.offsetBy(-3, 0, 0)) );
2196  }
2197 
2198  template<typename Accessor>
2199  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
2200  {
2201  return difference(
2202  grid.getValue(ijk.offsetBy( 0, 3, 0)), grid.getValue(ijk.offsetBy( 0, 2, 0)),
2203  grid.getValue(ijk.offsetBy( 0, 1, 0)), grid.getValue(ijk),
2204  grid.getValue(ijk.offsetBy( 0,-1, 0)), grid.getValue(ijk.offsetBy( 0,-2, 0)),
2205  grid.getValue(ijk.offsetBy( 0,-3, 0)) );
2206  }
2207 
2208  template<typename Accessor>
2209  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
2210  {
2211 
2212  return difference(
2213  grid.getValue(ijk.offsetBy( 0, 0, 3)), grid.getValue(ijk.offsetBy( 0, 0, 2)),
2214  grid.getValue(ijk.offsetBy( 0, 0, 1)), grid.getValue(ijk),
2215  grid.getValue(ijk.offsetBy( 0, 0,-1)), grid.getValue(ijk.offsetBy( 0, 0,-2)),
2216  grid.getValue(ijk.offsetBy( 0, 0,-3)) );
2217  }
2218 
2219  template<typename Accessor>
2220  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
2221  {
2222  typedef typename Accessor::ValueType ValueT;
2223  ValueT tmp1 =
2224  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 1, 0)) -
2225  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-1, 0));
2226  ValueT tmp2 =
2227  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 2, 0)) -
2228  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-2, 0));
2229  ValueT tmp3 =
2230  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 3, 0)) -
2231  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-3, 0));
2232  return ValueT(0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3);
2233  }
2234 
2235  template<typename Accessor>
2236  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
2237  {
2238  typedef typename Accessor::ValueType ValueT;
2239  ValueT tmp1 =
2240  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 1)) -
2241  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-1));
2242  ValueT tmp2 =
2243  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 2)) -
2244  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-2));
2245  ValueT tmp3 =
2246  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 3)) -
2247  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-3));
2248  return ValueT(0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3);
2249  }
2250 
2251  template<typename Accessor>
2252  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
2253  {
2254  typedef typename Accessor::ValueType ValueT;
2255  ValueT tmp1 =
2256  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 1)) -
2257  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-1));
2258  ValueT tmp2 =
2259  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 2)) -
2260  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-2));
2261  ValueT tmp3 =
2262  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 3)) -
2263  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-3));
2264  return ValueT(0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3);
2265  }
2266 
2267 
2268  // stencil access version
2269  template<typename Stencil>
2270  static typename Stencil::ValueType inX(const Stencil& S)
2271  {
2272  return difference( S.template getValue< 3, 0, 0>(), S.template getValue< 2, 0, 0>(),
2273  S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>(),
2274  S.template getValue<-1, 0, 0>(), S.template getValue<-2, 0, 0>(),
2275  S.template getValue<-3, 0, 0>() );
2276  }
2277 
2278  template<typename Stencil>
2279  static typename Stencil::ValueType inY(const Stencil& S)
2280  {
2281  return difference( S.template getValue< 0, 3, 0>(), S.template getValue< 0, 2, 0>(),
2282  S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>(),
2283  S.template getValue< 0,-1, 0>(), S.template getValue< 0,-2, 0>(),
2284  S.template getValue< 0,-3, 0>() );
2285 
2286  }
2287 
2288  template<typename Stencil>
2289  static typename Stencil::ValueType inZ(const Stencil& S)
2290  {
2291  return difference( S.template getValue< 0, 0, 3>(), S.template getValue< 0, 0, 2>(),
2292  S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>(),
2293  S.template getValue< 0, 0,-1>(), S.template getValue< 0, 0,-2>(),
2294  S.template getValue< 0, 0,-3>() );
2295  }
2296 
2297  template<typename Stencil>
2298  static typename Stencil::ValueType inXandY(const Stencil& S)
2299  {
2300  return crossdifference( S.template getValue< 1, 1, 0>(), S.template getValue<-1, 1, 0>(),
2301  S.template getValue< 1,-1, 0>(), S.template getValue<-1,-1, 0>(),
2302  S.template getValue< 2, 1, 0>(), S.template getValue<-2, 1, 0>(),
2303  S.template getValue< 2,-1, 0>(), S.template getValue<-2,-1, 0>(),
2304  S.template getValue< 3, 1, 0>(), S.template getValue<-3, 1, 0>(),
2305  S.template getValue< 3,-1, 0>(), S.template getValue<-3,-1, 0>(),
2306  S.template getValue< 1, 2, 0>(), S.template getValue<-1, 2, 0>(),
2307  S.template getValue< 1,-2, 0>(), S.template getValue<-1,-2, 0>(),
2308  S.template getValue< 2, 2, 0>(), S.template getValue<-2, 2, 0>(),
2309  S.template getValue< 2,-2, 0>(), S.template getValue<-2,-2, 0>(),
2310  S.template getValue< 3, 2, 0>(), S.template getValue<-3, 2, 0>(),
2311  S.template getValue< 3,-2, 0>(), S.template getValue<-3,-2, 0>(),
2312  S.template getValue< 1, 3, 0>(), S.template getValue<-1, 3, 0>(),
2313  S.template getValue< 1,-3, 0>(), S.template getValue<-1,-3, 0>(),
2314  S.template getValue< 2, 3, 0>(), S.template getValue<-2, 3, 0>(),
2315  S.template getValue< 2,-3, 0>(), S.template getValue<-2,-3, 0>(),
2316  S.template getValue< 3, 3, 0>(), S.template getValue<-3, 3, 0>(),
2317  S.template getValue< 3,-3, 0>(), S.template getValue<-3,-3, 0>() );
2318  }
2319 
2320  template<typename Stencil>
2321  static typename Stencil::ValueType inXandZ(const Stencil& S)
2322  {
2323  return crossdifference( S.template getValue< 1, 0, 1>(), S.template getValue<-1, 0, 1>(),
2324  S.template getValue< 1, 0,-1>(), S.template getValue<-1, 0,-1>(),
2325  S.template getValue< 2, 0, 1>(), S.template getValue<-2, 0, 1>(),
2326  S.template getValue< 2, 0,-1>(), S.template getValue<-2, 0,-1>(),
2327  S.template getValue< 3, 0, 1>(), S.template getValue<-3, 0, 1>(),
2328  S.template getValue< 3, 0,-1>(), S.template getValue<-3, 0,-1>(),
2329  S.template getValue< 1, 0, 2>(), S.template getValue<-1, 0, 2>(),
2330  S.template getValue< 1, 0,-2>(), S.template getValue<-1, 0,-2>(),
2331  S.template getValue< 2, 0, 2>(), S.template getValue<-2, 0, 2>(),
2332  S.template getValue< 2, 0,-2>(), S.template getValue<-2, 0,-2>(),
2333  S.template getValue< 3, 0, 2>(), S.template getValue<-3, 0, 2>(),
2334  S.template getValue< 3, 0,-2>(), S.template getValue<-3, 0,-2>(),
2335  S.template getValue< 1, 0, 3>(), S.template getValue<-1, 0, 3>(),
2336  S.template getValue< 1, 0,-3>(), S.template getValue<-1, 0,-3>(),
2337  S.template getValue< 2, 0, 3>(), S.template getValue<-2, 0, 3>(),
2338  S.template getValue< 2, 0,-3>(), S.template getValue<-2, 0,-3>(),
2339  S.template getValue< 3, 0, 3>(), S.template getValue<-3, 0, 3>(),
2340  S.template getValue< 3, 0,-3>(), S.template getValue<-3, 0,-3>() );
2341  }
2342 
2343  template<typename Stencil>
2344  static typename Stencil::ValueType inYandZ(const Stencil& S)
2345  {
2346  return crossdifference( S.template getValue< 0, 1, 1>(), S.template getValue< 0,-1, 1>(),
2347  S.template getValue< 0, 1,-1>(), S.template getValue< 0,-1,-1>(),
2348  S.template getValue< 0, 2, 1>(), S.template getValue< 0,-2, 1>(),
2349  S.template getValue< 0, 2,-1>(), S.template getValue< 0,-2,-1>(),
2350  S.template getValue< 0, 3, 1>(), S.template getValue< 0,-3, 1>(),
2351  S.template getValue< 0, 3,-1>(), S.template getValue< 0,-3,-1>(),
2352  S.template getValue< 0, 1, 2>(), S.template getValue< 0,-1, 2>(),
2353  S.template getValue< 0, 1,-2>(), S.template getValue< 0,-1,-2>(),
2354  S.template getValue< 0, 2, 2>(), S.template getValue< 0,-2, 2>(),
2355  S.template getValue< 0, 2,-2>(), S.template getValue< 0,-2,-2>(),
2356  S.template getValue< 0, 3, 2>(), S.template getValue< 0,-3, 2>(),
2357  S.template getValue< 0, 3,-2>(), S.template getValue< 0,-3,-2>(),
2358  S.template getValue< 0, 1, 3>(), S.template getValue< 0,-1, 3>(),
2359  S.template getValue< 0, 1,-3>(), S.template getValue< 0,-1,-3>(),
2360  S.template getValue< 0, 2, 3>(), S.template getValue< 0,-2, 3>(),
2361  S.template getValue< 0, 2,-3>(), S.template getValue< 0,-2,-3>(),
2362  S.template getValue< 0, 3, 3>(), S.template getValue< 0,-3, 3>(),
2363  S.template getValue< 0, 3,-3>(), S.template getValue< 0,-3,-3>() );
2364  }
2365 
2366 };
2367 
2368 } // end math namespace
2369 } // namespace OPENVDB_VERSION_NAME
2370 } // end openvdb namespace
2371 
2372 #endif // OPENVDB_MATH_FINITEDIFFERENCE_HAS_BEEN_INCLUDED
2373 
2374 // Copyright (c) 2012-2017 DreamWorks Animation LLC
2375 // All rights reserved. This software is distributed under the
2376 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inX(const Stencil &S)
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
static Accessor::ValueType inXandZ(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inZ(const Stencil &S)
static ValueType difference(const ValueType &xp1, const ValueType &xm1)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inX(const Stencil &S)
static Stencil::ValueType inXandY(const Stencil &S)
static Stencil::ValueType inY(const Stencil &S)
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inX(const Stencil &S)
static Stencil::ValueType inZ(const Stencil &S)
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
static Stencil::ValueType inY(const Stencil &S)
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
static Stencil::ValueType inZ(const Stencil &S)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
TemporalIntegrationScheme
Temporal integration schemes.
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static ValueType difference(const ValueType &xp2, const ValueType &xp1, const ValueType &xp0, const ValueType &xm1, const ValueType &xm2)
static Stencil::ValueType inY(const Stencil &S)
TemporalIntegrationScheme stringToTemporalIntegrationScheme(const std::string &s)
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
static Stencil::ValueType inY(const Stencil &S)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static ValueType difference(const ValueType &xp1, const ValueType &xp0, const ValueType &xm1)
static Accessor::ValueType inXandZ(const Accessor &grid, const Coord &ijk)
Type Pow2(Type x)
Return .
Definition: Math.h:514
static Stencil::ValueType inXandZ(const Stencil &S)
static Stencil::ValueType inY(const Stencil &S)
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xp0, const ValueType &xm1, const ValueType &xm2, const ValueType &xm3)
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
OPENVDB_DEPRECATED Real GudonovsNormSqrd(bool isOutside, Real dP_xm, Real dP_xp, Real dP_ym, Real dP_yp, Real dP_zm, Real dP_zp)
const GLdouble * v
Definition: glcorearb.h:836
static Accessor::ValueType inXandY(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inYandZ(const Accessor &grid, const Coord &ijk)
GLsizei const GLchar *const * string
Definition: glcorearb.h:813
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inX(const Stencil &S)
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
static ValueType difference(const ValueType &xm3, const ValueType &xm2, const ValueType &xm1, const ValueType &xm0)
static Stencil::ValueType inZ(const Stencil &S)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
static Stencil::ValueType inZ(const Stencil &S)
static Stencil::ValueType inXandZ(const Stencil &S)
static Stencil::ValueType inXandY(const Stencil &S)
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
static Stencil::ValueType inX(const Stencil &S)
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inZ(const Stencil &S)
static Stencil::ValueType inZ(const Stencil &S)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
GLfloat GLfloat GLfloat v2
Definition: glcorearb.h:817
static Stencil::ValueType inX(const Stencil &S)
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:48
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
GLfloat GLfloat GLfloat GLfloat v3
Definition: glcorearb.h:818
static Stencil::ValueType inX(const Stencil &S)
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
uint64 value_type
Definition: GA_PrimCompat.h:29
static ValueType difference(const ValueType &xm2, const ValueType &xm1, const ValueType &xm0)
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inZ(const Stencil &S)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inX(const Stencil &S)
static Stencil::ValueType inY(const Stencil &S)
DScheme
Different discrete schemes used in the first derivatives.
static Stencil::ValueType inZ(const Stencil &S)
static ValueType crossdifference(const ValueType &xpyp, const ValueType &xpym, const ValueType &xmyp, const ValueType &xmym)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static ValueType difference(const ValueType &xp2, const ValueType &xp1, const ValueType &xm1, const ValueType &xm2)
static ValueType difference(const ValueType &xm1, const ValueType &xm0)
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
static ValueType difference(const ValueType &xm3, const ValueType &xm2, const ValueType &xm1, const ValueType &xm0, const ValueType &xp1, const ValueType &xp2)
static Stencil::ValueType inY(const Stencil &S)
static Stencil::ValueType inY(const Stencil &S)
GLdouble n
Definition: glcorearb.h:2007
GLfloat f
Definition: glcorearb.h:1925
static Stencil::ValueType inZ(const Stencil &S)
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inX(const Stencil &S)
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
static ValueType crossdifference(const ValueType &xp1yp1, const ValueType &xm1yp1, const ValueType &xp1ym1, const ValueType &xm1ym1, const ValueType &xp2yp1, const ValueType &xm2yp1, const ValueType &xp2ym1, const ValueType &xm2ym1, const ValueType &xp3yp1, const ValueType &xm3yp1, const ValueType &xp3ym1, const ValueType &xm3ym1, const ValueType &xp1yp2, const ValueType &xm1yp2, const ValueType &xp1ym2, const ValueType &xm1ym2, const ValueType &xp2yp2, const ValueType &xm2yp2, const ValueType &xp2ym2, const ValueType &xm2ym2, const ValueType &xp3yp2, const ValueType &xm3yp2, const ValueType &xp3ym2, const ValueType &xm3ym2, const ValueType &xp1yp3, const ValueType &xm1yp3, const ValueType &xp1ym3, const ValueType &xm1ym3, const ValueType &xp2yp3, const ValueType &xm2yp3, const ValueType &xp2ym3, const ValueType &xm2ym3, const ValueType &xp3yp3, const ValueType &xm3yp3, const ValueType &xp3ym3, const ValueType &xm3ym3)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static ValueType difference(const ValueType &xm3, const ValueType &xm2, const ValueType &xm1, const ValueType &xm0, const ValueType &xp1, const ValueType &xp2)
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
#define OPENVDB_DEPRECATED
Definition: Platform.h:49
#define OPENVDB_VERSION_NAME
Definition: version.h:43
static Stencil::ValueType inY(const Stencil &S)
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
DScheme stringToDScheme(const std::string &s)
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xm1, const ValueType &xm2, const ValueType &xm3)
static Stencil::ValueType inY(const Stencil &S)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
std::string dsSchemeToString(DScheme dss)
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xp0, const ValueType &xm1, const ValueType &xm2)
static Accessor::ValueType inXandY(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inZ(const Stencil &S)
static Stencil::ValueType inX(const Stencil &S)
static Stencil::ValueType inX(const Stencil &S)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inX(const Stencil &S)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inX(const Stencil &S)
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xp0, const ValueType &xm1, const ValueType &xm2)
static Stencil::ValueType inZ(const Stencil &S)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inXandY(const Stencil &S)
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
static Stencil::ValueType inZ(const Stencil &S)
static ValueType difference(const ValueType &xp1, const ValueType &xm1)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xp0)
static Stencil::ValueType inZ(const Stencil &S)
static Stencil::ValueType inZ(const Stencil &S)
static Stencil::ValueType inY(const Stencil &S)
std::string temporalIntegrationSchemeToMenuName(TemporalIntegrationScheme tis)
static Stencil::ValueType inY(const Stencil &S)
static Stencil::ValueType inYandZ(const Stencil &S)
std::string temporalIntegrationSchemeToString(TemporalIntegrationScheme tis)
static Accessor::ValueType inXandY(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inY(const Stencil &S)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inX(const Stencil &S)
static Stencil::ValueType inY(const Stencil &S)
static Accessor::ValueType inXandZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
static Stencil::ValueType inZ(const Stencil &S)
std::string biasedGradientSchemeToMenuName(BiasedGradientScheme bgs)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inYandZ(const Stencil &S)
static Stencil::ValueType inX(const Stencil &S)
static Stencil::ValueType inY(const Stencil &S)
GLfloat GLfloat v1
Definition: glcorearb.h:816
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:622
static Accessor::ValueType inXandY(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inXandZ(const Stencil &S)
static Stencil::ValueType inYandZ(const Stencil &S)
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
static Stencil::ValueType inY(const Stencil &S)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
static Stencil::ValueType inY(const Stencil &S)
DDScheme
Different discrete schemes used in the second derivatives.
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inX(const Stencil &S)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
std::string dsSchemeToMenuName(DScheme dss)
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inXandZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static ValueType difference(const ValueType &xp2, const ValueType &xp1, const ValueType &xp0)
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
static Accessor::ValueType inYandZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static ValueType crossdifference(const ValueType &xp2yp2, const ValueType &xp2yp1, const ValueType &xp2ym1, const ValueType &xp2ym2, const ValueType &xp1yp2, const ValueType &xp1yp1, const ValueType &xp1ym1, const ValueType &xp1ym2, const ValueType &xm2yp2, const ValueType &xm2yp1, const ValueType &xm2ym1, const ValueType &xm2ym2, const ValueType &xm1yp2, const ValueType &xm1yp1, const ValueType &xm1ym1, const ValueType &xm1ym2)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
BiasedGradientScheme stringToBiasedGradientScheme(const std::string &s)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
std::string biasedGradientSchemeToString(BiasedGradientScheme bgs)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:561
static Accessor::ValueType inYandZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
static ValueType difference(const ValueType &xp1, const ValueType &xp0)
static Stencil::ValueType inZ(const Stencil &S)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Real GodunovsNormSqrd(bool isOutside, Real dP_xm, Real dP_xp, Real dP_ym, Real dP_yp, Real dP_zm, Real dP_zp)
static Stencil::ValueType inX(const Stencil &S)
ValueType WENO5(const ValueType &v1, const ValueType &v2, const ValueType &v3, const ValueType &v4, const ValueType &v5, float scale2=0.01f)
Implementation of nominally fifth-order finite-difference WENO.
static Accessor::ValueType inYandZ(const Accessor &grid, const Coord &ijk)