HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FiniteDifference.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2012-2018 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 math/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 #include <string>
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>
378 inline Real
379 GodunovsNormSqrd(bool isOutside, const Vec3<Real>& gradient_m, const Vec3<Real>& gradient_p)
380 {
381  return GodunovsNormSqrd<Real>(isOutside,
382  gradient_m[0], gradient_p[0],
383  gradient_m[1], gradient_p[1],
384  gradient_m[2], gradient_p[2]);
385 }
386 
387 
388 #ifdef DWA_OPENVDB
389 inline simd::Float4 simdMin(const simd::Float4& a, const simd::Float4& b) {
390  return simd::Float4(_mm_min_ps(a.base(), b.base()));
391 }
392 inline simd::Float4 simdMax(const simd::Float4& a, const simd::Float4& b) {
393  return simd::Float4(_mm_max_ps(a.base(), b.base()));
394 }
395 
396 inline float simdSum(const simd::Float4& v);
397 
398 inline simd::Float4 Pow2(const simd::Float4& v) { return v * v; }
399 
400 template<>
401 inline simd::Float4
402 WENO5<simd::Float4>(const simd::Float4& v1, const simd::Float4& v2, const simd::Float4& v3,
403  const simd::Float4& v4, const simd::Float4& v5, float scale2)
404 {
405  using math::Pow2;
406  using F4 = simd::Float4;
407  const F4
408  C(13.f / 12.f),
409  eps(1.0e-6f * scale2),
410  two(2.0), three(3.0), four(4.0), five(5.0), fourth(0.25),
411  A1 = F4(0.1f) / Pow2(C*Pow2(v1-two*v2+v3) + fourth*Pow2(v1-four*v2+three*v3) + eps),
412  A2 = F4(0.6f) / Pow2(C*Pow2(v2-two*v3+v4) + fourth*Pow2(v2-v4) + eps),
413  A3 = F4(0.3f) / Pow2(C*Pow2(v3-two*v4+v5) + fourth*Pow2(three*v3-four*v4+v5) + eps);
414  return (A1 * (two * v1 - F4(7.0) * v2 + F4(11.0) * v3) +
415  A2 * (five * v3 - v2 + two * v4) +
416  A3 * (two * v3 + five * v4 - v5)) / (F4(6.0) * (A1 + A2 + A3));
417 }
418 
419 
420 inline float
421 simdSum(const simd::Float4& v)
422 {
423  // temp = { v3+v3, v2+v2, v1+v3, v0+v2 }
424  __m128 temp = _mm_add_ps(v.base(), _mm_movehl_ps(v.base(), v.base()));
425  // temp = { v3+v3, v2+v2, v1+v3, (v0+v2)+(v1+v3) }
426  temp = _mm_add_ss(temp, _mm_shuffle_ps(temp, temp, 1));
427  return _mm_cvtss_f32(temp);
428 }
429 
430 inline float
431 GodunovsNormSqrd(bool isOutside, const simd::Float4& dP_m, const simd::Float4& dP_p)
432 {
433  const simd::Float4 zero(0.0);
434  simd::Float4 v = isOutside
435  ? simdMax(math::Pow2(simdMax(dP_m, zero)), math::Pow2(simdMin(dP_p, zero)))
436  : simdMax(math::Pow2(simdMin(dP_m, zero)), math::Pow2(simdMax(dP_p, zero)));
437  return simdSum(v);//should be v[0]+v[1]+v[2]
438 }
439 #endif
440 
441 
442 template<DScheme DiffScheme>
443 struct D1
444 {
445  // random access version
446  template<typename Accessor>
447  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk);
448 
449  template<typename Accessor>
450  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk);
451 
452  template<typename Accessor>
453  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk);
454 
455  // stencil access version
456  template<typename Stencil>
457  static typename Stencil::ValueType inX(const Stencil& S);
458 
459  template<typename Stencil>
460  static typename Stencil::ValueType inY(const Stencil& S);
461 
462  template<typename Stencil>
463  static typename Stencil::ValueType inZ(const Stencil& S);
464 };
465 
466 template<>
467 struct D1<CD_2NDT>
468 {
469  // the difference opperator
470  template <typename ValueType>
471  static ValueType difference(const ValueType& xp1, const ValueType& xm1) {
472  return xp1 - xm1;
473  }
474 
475  // random access version
476  template<typename Accessor>
477  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
478  {
479  return difference(
480  grid.getValue(ijk.offsetBy(1, 0, 0)),
481  grid.getValue(ijk.offsetBy(-1, 0, 0)));
482  }
483 
484  template<typename Accessor>
485  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
486  {
487  return difference(
488  grid.getValue(ijk.offsetBy(0, 1, 0)),
489  grid.getValue(ijk.offsetBy( 0, -1, 0)));
490  }
491 
492  template<typename Accessor>
493  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
494  {
495  return difference(
496  grid.getValue(ijk.offsetBy(0, 0, 1)),
497  grid.getValue(ijk.offsetBy( 0, 0, -1)));
498  }
499 
500  // stencil access version
501  template<typename Stencil>
502  static typename Stencil::ValueType inX(const Stencil& S)
503  {
504  return difference( S.template getValue< 1, 0, 0>(), S.template getValue<-1, 0, 0>());
505  }
506 
507  template<typename Stencil>
508  static typename Stencil::ValueType inY(const Stencil& S)
509  {
510  return difference( S.template getValue< 0, 1, 0>(), S.template getValue< 0,-1, 0>());
511  }
512 
513  template<typename Stencil>
514  static typename Stencil::ValueType inZ(const Stencil& S)
515  {
516  return difference( S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0,-1>());
517  }
518 };
519 
520 template<>
521 struct D1<CD_2ND>
522 {
523 
524  // the difference opperator
525  template <typename ValueType>
526  static ValueType difference(const ValueType& xp1, const ValueType& xm1) {
527  return (xp1 - xm1)*ValueType(0.5);
528  }
529 
530 
531  // random access
532  template<typename Accessor>
533  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
534  {
535  return difference(
536  grid.getValue(ijk.offsetBy(1, 0, 0)),
537  grid.getValue(ijk.offsetBy(-1, 0, 0)));
538  }
539 
540  template<typename Accessor>
541  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
542  {
543  return difference(
544  grid.getValue(ijk.offsetBy(0, 1, 0)),
545  grid.getValue(ijk.offsetBy( 0, -1, 0)));
546  }
547 
548  template<typename Accessor>
549  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
550  {
551  return difference(
552  grid.getValue(ijk.offsetBy(0, 0, 1)),
553  grid.getValue(ijk.offsetBy( 0, 0, -1)));
554  }
555 
556 
557  // stencil access version
558  template<typename Stencil>
559  static typename Stencil::ValueType inX(const Stencil& S)
560  {
561  return difference(S.template getValue< 1, 0, 0>(), S.template getValue<-1, 0, 0>());
562  }
563  template<typename Stencil>
564  static typename Stencil::ValueType inY(const Stencil& S)
565  {
566  return difference(S.template getValue< 0, 1, 0>(), S.template getValue< 0,-1, 0>());
567  }
568 
569  template<typename Stencil>
570  static typename Stencil::ValueType inZ(const Stencil& S)
571  {
572  return difference(S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0,-1>());
573  }
574 
575 };
576 
577 template<>
578 struct D1<CD_4TH>
579 {
580 
581  // the difference opperator
582  template <typename ValueType>
583  static ValueType difference( const ValueType& xp2, const ValueType& xp1,
584  const ValueType& xm1, const ValueType& xm2 ) {
585  return ValueType(2./3.)*(xp1 - xm1) + ValueType(1./12.)*(xm2 - xp2) ;
586  }
587 
588 
589  // random access version
590  template<typename Accessor>
591  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
592  {
593  return difference(
594  grid.getValue(ijk.offsetBy( 2,0,0)), grid.getValue(ijk.offsetBy( 1,0,0)),
595  grid.getValue(ijk.offsetBy(-1,0,0)), grid.getValue(ijk.offsetBy(-2,0,0)) );
596  }
597 
598  template<typename Accessor>
599  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
600  {
601 
602  return difference(
603  grid.getValue(ijk.offsetBy( 0, 2, 0)), grid.getValue(ijk.offsetBy( 0, 1, 0)),
604  grid.getValue(ijk.offsetBy( 0,-1, 0)), grid.getValue(ijk.offsetBy( 0,-2, 0)) );
605  }
606 
607  template<typename Accessor>
608  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
609  {
610 
611  return difference(
612  grid.getValue(ijk.offsetBy( 0, 0, 2)), grid.getValue(ijk.offsetBy( 0, 0, 1)),
613  grid.getValue(ijk.offsetBy( 0, 0,-1)), grid.getValue(ijk.offsetBy( 0, 0,-2)) );
614  }
615 
616 
617  // stencil access version
618  template<typename Stencil>
619  static typename Stencil::ValueType inX(const Stencil& S)
620  {
621  return difference( S.template getValue< 2, 0, 0>(),
622  S.template getValue< 1, 0, 0>(),
623  S.template getValue<-1, 0, 0>(),
624  S.template getValue<-2, 0, 0>() );
625  }
626 
627  template<typename Stencil>
628  static typename Stencil::ValueType inY(const Stencil& S)
629  {
630  return difference( S.template getValue< 0, 2, 0>(),
631  S.template getValue< 0, 1, 0>(),
632  S.template getValue< 0,-1, 0>(),
633  S.template getValue< 0,-2, 0>() );
634  }
635 
636  template<typename Stencil>
637  static typename Stencil::ValueType inZ(const Stencil& S)
638  {
639  return difference( S.template getValue< 0, 0, 2>(),
640  S.template getValue< 0, 0, 1>(),
641  S.template getValue< 0, 0,-1>(),
642  S.template getValue< 0, 0,-2>() );
643  }
644 };
645 
646 template<>
647 struct D1<CD_6TH>
648 {
649 
650  // the difference opperator
651  template <typename ValueType>
652  static ValueType difference( const ValueType& xp3, const ValueType& xp2, const ValueType& xp1,
653  const ValueType& xm1, const ValueType& xm2, const ValueType& xm3 )
654  {
655  return ValueType(3./4.)*(xp1 - xm1) - ValueType(0.15)*(xp2 - xm2)
656  + ValueType(1./60.)*(xp3-xm3);
657  }
658 
659 
660  // random access version
661  template<typename Accessor>
662  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
663  {
664  return difference(
665  grid.getValue(ijk.offsetBy( 3,0,0)), grid.getValue(ijk.offsetBy( 2,0,0)),
666  grid.getValue(ijk.offsetBy( 1,0,0)), grid.getValue(ijk.offsetBy(-1,0,0)),
667  grid.getValue(ijk.offsetBy(-2,0,0)), grid.getValue(ijk.offsetBy(-3,0,0)));
668  }
669 
670  template<typename Accessor>
671  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
672  {
673  return difference(
674  grid.getValue(ijk.offsetBy( 0, 3, 0)), grid.getValue(ijk.offsetBy( 0, 2, 0)),
675  grid.getValue(ijk.offsetBy( 0, 1, 0)), grid.getValue(ijk.offsetBy( 0,-1, 0)),
676  grid.getValue(ijk.offsetBy( 0,-2, 0)), grid.getValue(ijk.offsetBy( 0,-3, 0)));
677  }
678 
679  template<typename Accessor>
680  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
681  {
682  return difference(
683  grid.getValue(ijk.offsetBy( 0, 0, 3)), grid.getValue(ijk.offsetBy( 0, 0, 2)),
684  grid.getValue(ijk.offsetBy( 0, 0, 1)), grid.getValue(ijk.offsetBy( 0, 0,-1)),
685  grid.getValue(ijk.offsetBy( 0, 0,-2)), grid.getValue(ijk.offsetBy( 0, 0,-3)));
686  }
687 
688  // stencil access version
689  template<typename Stencil>
690  static typename Stencil::ValueType inX(const Stencil& S)
691  {
692  return difference(S.template getValue< 3, 0, 0>(),
693  S.template getValue< 2, 0, 0>(),
694  S.template getValue< 1, 0, 0>(),
695  S.template getValue<-1, 0, 0>(),
696  S.template getValue<-2, 0, 0>(),
697  S.template getValue<-3, 0, 0>());
698  }
699 
700  template<typename Stencil>
701  static typename Stencil::ValueType inY(const Stencil& S)
702  {
703 
704  return difference( S.template getValue< 0, 3, 0>(),
705  S.template getValue< 0, 2, 0>(),
706  S.template getValue< 0, 1, 0>(),
707  S.template getValue< 0,-1, 0>(),
708  S.template getValue< 0,-2, 0>(),
709  S.template getValue< 0,-3, 0>());
710  }
711 
712  template<typename Stencil>
713  static typename Stencil::ValueType inZ(const Stencil& S)
714  {
715 
716  return difference( S.template getValue< 0, 0, 3>(),
717  S.template getValue< 0, 0, 2>(),
718  S.template getValue< 0, 0, 1>(),
719  S.template getValue< 0, 0,-1>(),
720  S.template getValue< 0, 0,-2>(),
721  S.template getValue< 0, 0,-3>());
722  }
723 };
724 
725 
726 template<>
727 struct D1<FD_1ST>
728 {
729 
730  // the difference opperator
731  template <typename ValueType>
732  static ValueType difference(const ValueType& xp1, const ValueType& xp0) {
733  return xp1 - xp0;
734  }
735 
736 
737  // random access version
738  template<typename Accessor>
739  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
740  {
741  return difference(grid.getValue(ijk.offsetBy(1, 0, 0)), grid.getValue(ijk));
742  }
743 
744  template<typename Accessor>
745  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
746  {
747  return difference(grid.getValue(ijk.offsetBy(0, 1, 0)), grid.getValue(ijk));
748  }
749 
750  template<typename Accessor>
751  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
752  {
753  return difference(grid.getValue(ijk.offsetBy(0, 0, 1)), grid.getValue(ijk));
754  }
755 
756  // stencil access version
757  template<typename Stencil>
758  static typename Stencil::ValueType inX(const Stencil& S)
759  {
760  return difference(S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>());
761  }
762 
763  template<typename Stencil>
764  static typename Stencil::ValueType inY(const Stencil& S)
765  {
766  return difference(S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>());
767  }
768 
769  template<typename Stencil>
770  static typename Stencil::ValueType inZ(const Stencil& S)
771  {
772  return difference(S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>());
773  }
774 };
775 
776 
777 template<>
778 struct D1<FD_2ND>
779 {
780  // the difference opperator
781  template <typename ValueType>
782  static ValueType difference(const ValueType& xp2, const ValueType& xp1, const ValueType& xp0)
783  {
784  return ValueType(2)*xp1 -(ValueType(0.5)*xp2 + ValueType(3./2.)*xp0);
785  }
786 
787 
788  // random access version
789  template<typename Accessor>
790  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
791  {
792  return difference(
793  grid.getValue(ijk.offsetBy(2,0,0)),
794  grid.getValue(ijk.offsetBy(1,0,0)),
795  grid.getValue(ijk));
796  }
797 
798  template<typename Accessor>
799  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
800  {
801  return difference(
802  grid.getValue(ijk.offsetBy(0,2,0)),
803  grid.getValue(ijk.offsetBy(0,1,0)),
804  grid.getValue(ijk));
805  }
806 
807  template<typename Accessor>
808  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
809  {
810  return difference(
811  grid.getValue(ijk.offsetBy(0,0,2)),
812  grid.getValue(ijk.offsetBy(0,0,1)),
813  grid.getValue(ijk));
814  }
815 
816 
817  // stencil access version
818  template<typename Stencil>
819  static typename Stencil::ValueType inX(const Stencil& S)
820  {
821  return difference( S.template getValue< 2, 0, 0>(),
822  S.template getValue< 1, 0, 0>(),
823  S.template getValue< 0, 0, 0>() );
824  }
825 
826  template<typename Stencil>
827  static typename Stencil::ValueType inY(const Stencil& S)
828  {
829  return difference( S.template getValue< 0, 2, 0>(),
830  S.template getValue< 0, 1, 0>(),
831  S.template getValue< 0, 0, 0>() );
832  }
833 
834  template<typename Stencil>
835  static typename Stencil::ValueType inZ(const Stencil& S)
836  {
837  return difference( S.template getValue< 0, 0, 2>(),
838  S.template getValue< 0, 0, 1>(),
839  S.template getValue< 0, 0, 0>() );
840  }
841 
842 };
843 
844 
845 template<>
846 struct D1<FD_3RD>
847 {
848 
849  // the difference opperator
850  template<typename ValueType>
851  static ValueType difference(const ValueType& xp3, const ValueType& xp2,
852  const ValueType& xp1, const ValueType& xp0)
853  {
854  return static_cast<ValueType>(xp3/3.0 - 1.5*xp2 + 3.0*xp1 - 11.0*xp0/6.0);
855  }
856 
857 
858  // random access version
859  template<typename Accessor>
860  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
861  {
862  return difference( grid.getValue(ijk.offsetBy(3,0,0)),
863  grid.getValue(ijk.offsetBy(2,0,0)),
864  grid.getValue(ijk.offsetBy(1,0,0)),
865  grid.getValue(ijk) );
866  }
867 
868  template<typename Accessor>
869  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
870  {
871  return difference( grid.getValue(ijk.offsetBy(0,3,0)),
872  grid.getValue(ijk.offsetBy(0,2,0)),
873  grid.getValue(ijk.offsetBy(0,1,0)),
874  grid.getValue(ijk) );
875  }
876 
877  template<typename Accessor>
878  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
879  {
880  return difference( grid.getValue(ijk.offsetBy(0,0,3)),
881  grid.getValue(ijk.offsetBy(0,0,2)),
882  grid.getValue(ijk.offsetBy(0,0,1)),
883  grid.getValue(ijk) );
884  }
885 
886 
887  // stencil access version
888  template<typename Stencil>
889  static typename Stencil::ValueType inX(const Stencil& S)
890  {
891  return difference(S.template getValue< 3, 0, 0>(),
892  S.template getValue< 2, 0, 0>(),
893  S.template getValue< 1, 0, 0>(),
894  S.template getValue< 0, 0, 0>() );
895  }
896 
897  template<typename Stencil>
898  static typename Stencil::ValueType inY(const Stencil& S)
899  {
900  return difference(S.template getValue< 0, 3, 0>(),
901  S.template getValue< 0, 2, 0>(),
902  S.template getValue< 0, 1, 0>(),
903  S.template getValue< 0, 0, 0>() );
904  }
905 
906  template<typename Stencil>
907  static typename Stencil::ValueType inZ(const Stencil& S)
908  {
909  return difference( S.template getValue< 0, 0, 3>(),
910  S.template getValue< 0, 0, 2>(),
911  S.template getValue< 0, 0, 1>(),
912  S.template getValue< 0, 0, 0>() );
913  }
914 };
915 
916 
917 template<>
918 struct D1<BD_1ST>
919 {
920 
921  // the difference opperator
922  template <typename ValueType>
923  static ValueType difference(const ValueType& xm1, const ValueType& xm0) {
924  return -D1<FD_1ST>::difference(xm1, xm0);
925  }
926 
927 
928  // random access version
929  template<typename Accessor>
930  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
931  {
932  return difference(grid.getValue(ijk.offsetBy(-1,0,0)), grid.getValue(ijk));
933  }
934 
935  template<typename Accessor>
936  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
937  {
938  return difference(grid.getValue(ijk.offsetBy(0,-1,0)), grid.getValue(ijk));
939  }
940 
941  template<typename Accessor>
942  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
943  {
944  return difference(grid.getValue(ijk.offsetBy(0, 0,-1)), grid.getValue(ijk));
945  }
946 
947 
948  // stencil access version
949  template<typename Stencil>
950  static typename Stencil::ValueType inX(const Stencil& S)
951  {
952  return difference(S.template getValue<-1, 0, 0>(), S.template getValue< 0, 0, 0>());
953  }
954 
955  template<typename Stencil>
956  static typename Stencil::ValueType inY(const Stencil& S)
957  {
958  return difference(S.template getValue< 0,-1, 0>(), S.template getValue< 0, 0, 0>());
959  }
960 
961  template<typename Stencil>
962  static typename Stencil::ValueType inZ(const Stencil& S)
963  {
964  return difference(S.template getValue< 0, 0,-1>(), S.template getValue< 0, 0, 0>());
965  }
966 };
967 
968 
969 template<>
970 struct D1<BD_2ND>
971 {
972 
973  // the difference opperator
974  template <typename ValueType>
975  static ValueType difference(const ValueType& xm2, const ValueType& xm1, const ValueType& xm0)
976  {
977  return -D1<FD_2ND>::difference(xm2, xm1, xm0);
978  }
979 
980 
981  // random access version
982  template<typename Accessor>
983  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
984  {
985  return difference( grid.getValue(ijk.offsetBy(-2,0,0)),
986  grid.getValue(ijk.offsetBy(-1,0,0)),
987  grid.getValue(ijk) );
988  }
989 
990  template<typename Accessor>
991  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
992  {
993  return difference( grid.getValue(ijk.offsetBy(0,-2,0)),
994  grid.getValue(ijk.offsetBy(0,-1,0)),
995  grid.getValue(ijk) );
996  }
997 
998  template<typename Accessor>
999  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1000  {
1001  return difference( grid.getValue(ijk.offsetBy(0,0,-2)),
1002  grid.getValue(ijk.offsetBy(0,0,-1)),
1003  grid.getValue(ijk) );
1004  }
1005 
1006  // stencil access version
1007  template<typename Stencil>
1008  static typename Stencil::ValueType inX(const Stencil& S)
1009  {
1010  return difference( S.template getValue<-2, 0, 0>(),
1011  S.template getValue<-1, 0, 0>(),
1012  S.template getValue< 0, 0, 0>() );
1013  }
1014 
1015  template<typename Stencil>
1016  static typename Stencil::ValueType inY(const Stencil& S)
1017  {
1018  return difference( S.template getValue< 0,-2, 0>(),
1019  S.template getValue< 0,-1, 0>(),
1020  S.template getValue< 0, 0, 0>() );
1021  }
1022 
1023  template<typename Stencil>
1024  static typename Stencil::ValueType inZ(const Stencil& S)
1025  {
1026  return difference( S.template getValue< 0, 0,-2>(),
1027  S.template getValue< 0, 0,-1>(),
1028  S.template getValue< 0, 0, 0>() );
1029  }
1030 };
1031 
1032 
1033 template<>
1034 struct D1<BD_3RD>
1035 {
1036 
1037  // the difference opperator
1038  template <typename ValueType>
1039  static ValueType difference(const ValueType& xm3, const ValueType& xm2,
1040  const ValueType& xm1, const ValueType& xm0)
1041  {
1042  return -D1<FD_3RD>::difference(xm3, xm2, xm1, xm0);
1043  }
1044 
1045  // random access version
1046  template<typename Accessor>
1047  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1048  {
1049  return difference( grid.getValue(ijk.offsetBy(-3,0,0)),
1050  grid.getValue(ijk.offsetBy(-2,0,0)),
1051  grid.getValue(ijk.offsetBy(-1,0,0)),
1052  grid.getValue(ijk) );
1053  }
1054 
1055  template<typename Accessor>
1056  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1057  {
1058  return difference( grid.getValue(ijk.offsetBy( 0,-3,0)),
1059  grid.getValue(ijk.offsetBy( 0,-2,0)),
1060  grid.getValue(ijk.offsetBy( 0,-1,0)),
1061  grid.getValue(ijk) );
1062  }
1063 
1064  template<typename Accessor>
1065  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1066  {
1067  return difference( grid.getValue(ijk.offsetBy( 0, 0,-3)),
1068  grid.getValue(ijk.offsetBy( 0, 0,-2)),
1069  grid.getValue(ijk.offsetBy( 0, 0,-1)),
1070  grid.getValue(ijk) );
1071  }
1072 
1073  // stencil access version
1074  template<typename Stencil>
1075  static typename Stencil::ValueType inX(const Stencil& S)
1076  {
1077  return difference( S.template getValue<-3, 0, 0>(),
1078  S.template getValue<-2, 0, 0>(),
1079  S.template getValue<-1, 0, 0>(),
1080  S.template getValue< 0, 0, 0>() );
1081  }
1082 
1083  template<typename Stencil>
1084  static typename Stencil::ValueType inY(const Stencil& S)
1085  {
1086  return difference( S.template getValue< 0,-3, 0>(),
1087  S.template getValue< 0,-2, 0>(),
1088  S.template getValue< 0,-1, 0>(),
1089  S.template getValue< 0, 0, 0>() );
1090  }
1091 
1092  template<typename Stencil>
1093  static typename Stencil::ValueType inZ(const Stencil& S)
1094  {
1095  return difference( S.template getValue< 0, 0,-3>(),
1096  S.template getValue< 0, 0,-2>(),
1097  S.template getValue< 0, 0,-1>(),
1098  S.template getValue< 0, 0, 0>() );
1099  }
1100 
1101 };
1102 
1103 template<>
1104 struct D1<FD_WENO5>
1105 {
1106  // the difference operator
1107  template <typename ValueType>
1108  static ValueType difference(const ValueType& xp3, const ValueType& xp2,
1109  const ValueType& xp1, const ValueType& xp0,
1110  const ValueType& xm1, const ValueType& xm2) {
1111  return WENO5<ValueType>(xp3, xp2, xp1, xp0, xm1)
1112  - WENO5<ValueType>(xp2, xp1, xp0, xm1, xm2);
1113  }
1114 
1115 
1116  // random access version
1117  template<typename Accessor>
1118  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1119  {
1120  using ValueType = typename Accessor::ValueType;
1121  ValueType V[6];
1122  V[0] = grid.getValue(ijk.offsetBy(3,0,0));
1123  V[1] = grid.getValue(ijk.offsetBy(2,0,0));
1124  V[2] = grid.getValue(ijk.offsetBy(1,0,0));
1125  V[3] = grid.getValue(ijk);
1126  V[4] = grid.getValue(ijk.offsetBy(-1,0,0));
1127  V[5] = grid.getValue(ijk.offsetBy(-2,0,0));
1128 
1129  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1130  }
1131 
1132  template<typename Accessor>
1133  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1134  {
1135  using ValueType = typename Accessor::ValueType;
1136  ValueType V[6];
1137  V[0] = grid.getValue(ijk.offsetBy(0,3,0));
1138  V[1] = grid.getValue(ijk.offsetBy(0,2,0));
1139  V[2] = grid.getValue(ijk.offsetBy(0,1,0));
1140  V[3] = grid.getValue(ijk);
1141  V[4] = grid.getValue(ijk.offsetBy(0,-1,0));
1142  V[5] = grid.getValue(ijk.offsetBy(0,-2,0));
1143 
1144  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1145  }
1146 
1147  template<typename Accessor>
1148  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1149  {
1150  using ValueType = typename Accessor::ValueType;
1151  ValueType V[6];
1152  V[0] = grid.getValue(ijk.offsetBy(0,0,3));
1153  V[1] = grid.getValue(ijk.offsetBy(0,0,2));
1154  V[2] = grid.getValue(ijk.offsetBy(0,0,1));
1155  V[3] = grid.getValue(ijk);
1156  V[4] = grid.getValue(ijk.offsetBy(0,0,-1));
1157  V[5] = grid.getValue(ijk.offsetBy(0,0,-2));
1158 
1159  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1160  }
1161 
1162  // stencil access version
1163  template<typename Stencil>
1164  static typename Stencil::ValueType inX(const Stencil& S)
1165  {
1166 
1167  return static_cast<typename Stencil::ValueType>(difference(
1168  S.template getValue< 3, 0, 0>(),
1169  S.template getValue< 2, 0, 0>(),
1170  S.template getValue< 1, 0, 0>(),
1171  S.template getValue< 0, 0, 0>(),
1172  S.template getValue<-1, 0, 0>(),
1173  S.template getValue<-2, 0, 0>() ));
1174 
1175  }
1176 
1177  template<typename Stencil>
1178  static typename Stencil::ValueType inY(const Stencil& S)
1179  {
1180  return static_cast<typename Stencil::ValueType>(difference(
1181  S.template getValue< 0, 3, 0>(),
1182  S.template getValue< 0, 2, 0>(),
1183  S.template getValue< 0, 1, 0>(),
1184  S.template getValue< 0, 0, 0>(),
1185  S.template getValue< 0,-1, 0>(),
1186  S.template getValue< 0,-2, 0>() ));
1187  }
1188 
1189  template<typename Stencil>
1190  static typename Stencil::ValueType inZ(const Stencil& S)
1191  {
1192  return static_cast<typename Stencil::ValueType>(difference(
1193  S.template getValue< 0, 0, 3>(),
1194  S.template getValue< 0, 0, 2>(),
1195  S.template getValue< 0, 0, 1>(),
1196  S.template getValue< 0, 0, 0>(),
1197  S.template getValue< 0, 0,-1>(),
1198  S.template getValue< 0, 0,-2>() ));
1199  }
1200 };
1201 
1202 template<>
1204 {
1205 
1206  // the difference opperator
1207  template <typename ValueType>
1208  static ValueType difference(const ValueType& xp3, const ValueType& xp2,
1209  const ValueType& xp1, const ValueType& xp0,
1210  const ValueType& xm1, const ValueType& xm2) {
1211  return WENO5<ValueType>(xp3 - xp2, xp2 - xp1, xp1 - xp0, xp0-xm1, xm1-xm2);
1212  }
1213 
1214  // random access version
1215  template<typename Accessor>
1216  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1217  {
1218  using ValueType = typename Accessor::ValueType;
1219  ValueType V[6];
1220  V[0] = grid.getValue(ijk.offsetBy(3,0,0));
1221  V[1] = grid.getValue(ijk.offsetBy(2,0,0));
1222  V[2] = grid.getValue(ijk.offsetBy(1,0,0));
1223  V[3] = grid.getValue(ijk);
1224  V[4] = grid.getValue(ijk.offsetBy(-1,0,0));
1225  V[5] = grid.getValue(ijk.offsetBy(-2,0,0));
1226 
1227  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1228 
1229  }
1230 
1231  template<typename Accessor>
1232  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1233  {
1234  using ValueType = typename Accessor::ValueType;
1235  ValueType V[6];
1236  V[0] = grid.getValue(ijk.offsetBy(0,3,0));
1237  V[1] = grid.getValue(ijk.offsetBy(0,2,0));
1238  V[2] = grid.getValue(ijk.offsetBy(0,1,0));
1239  V[3] = grid.getValue(ijk);
1240  V[4] = grid.getValue(ijk.offsetBy(0,-1,0));
1241  V[5] = grid.getValue(ijk.offsetBy(0,-2,0));
1242 
1243  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1244  }
1245 
1246  template<typename Accessor>
1247  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1248  {
1249  using ValueType = typename Accessor::ValueType;
1250  ValueType V[6];
1251  V[0] = grid.getValue(ijk.offsetBy(0,0,3));
1252  V[1] = grid.getValue(ijk.offsetBy(0,0,2));
1253  V[2] = grid.getValue(ijk.offsetBy(0,0,1));
1254  V[3] = grid.getValue(ijk);
1255  V[4] = grid.getValue(ijk.offsetBy(0,0,-1));
1256  V[5] = grid.getValue(ijk.offsetBy(0,0,-2));
1257 
1258  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1259  }
1260 
1261  // stencil access version
1262  template<typename Stencil>
1263  static typename Stencil::ValueType inX(const Stencil& S)
1264  {
1265 
1266  return difference( S.template getValue< 3, 0, 0>(),
1267  S.template getValue< 2, 0, 0>(),
1268  S.template getValue< 1, 0, 0>(),
1269  S.template getValue< 0, 0, 0>(),
1270  S.template getValue<-1, 0, 0>(),
1271  S.template getValue<-2, 0, 0>() );
1272 
1273  }
1274 
1275  template<typename Stencil>
1276  static typename Stencil::ValueType inY(const Stencil& S)
1277  {
1278  return difference( S.template getValue< 0, 3, 0>(),
1279  S.template getValue< 0, 2, 0>(),
1280  S.template getValue< 0, 1, 0>(),
1281  S.template getValue< 0, 0, 0>(),
1282  S.template getValue< 0,-1, 0>(),
1283  S.template getValue< 0,-2, 0>() );
1284  }
1285 
1286  template<typename Stencil>
1287  static typename Stencil::ValueType inZ(const Stencil& S)
1288  {
1289 
1290  return difference( S.template getValue< 0, 0, 3>(),
1291  S.template getValue< 0, 0, 2>(),
1292  S.template getValue< 0, 0, 1>(),
1293  S.template getValue< 0, 0, 0>(),
1294  S.template getValue< 0, 0,-1>(),
1295  S.template getValue< 0, 0,-2>() );
1296  }
1297 
1298 };
1299 
1300 template<>
1301 struct D1<BD_WENO5>
1302 {
1303 
1304  template<typename ValueType>
1305  static ValueType difference(const ValueType& xm3, const ValueType& xm2, const ValueType& xm1,
1306  const ValueType& xm0, const ValueType& xp1, const ValueType& xp2)
1307  {
1308  return -D1<FD_WENO5>::difference(xm3, xm2, xm1, xm0, xp1, xp2);
1309  }
1310 
1311 
1312  // random access version
1313  template<typename Accessor>
1314  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1315  {
1316  using ValueType = typename Accessor::ValueType;
1317  ValueType V[6];
1318  V[0] = grid.getValue(ijk.offsetBy(-3,0,0));
1319  V[1] = grid.getValue(ijk.offsetBy(-2,0,0));
1320  V[2] = grid.getValue(ijk.offsetBy(-1,0,0));
1321  V[3] = grid.getValue(ijk);
1322  V[4] = grid.getValue(ijk.offsetBy(1,0,0));
1323  V[5] = grid.getValue(ijk.offsetBy(2,0,0));
1324 
1325  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1326  }
1327 
1328  template<typename Accessor>
1329  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1330  {
1331  using ValueType = typename Accessor::ValueType;
1332  ValueType V[6];
1333  V[0] = grid.getValue(ijk.offsetBy(0,-3,0));
1334  V[1] = grid.getValue(ijk.offsetBy(0,-2,0));
1335  V[2] = grid.getValue(ijk.offsetBy(0,-1,0));
1336  V[3] = grid.getValue(ijk);
1337  V[4] = grid.getValue(ijk.offsetBy(0,1,0));
1338  V[5] = grid.getValue(ijk.offsetBy(0,2,0));
1339 
1340  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1341  }
1342 
1343  template<typename Accessor>
1344  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1345  {
1346  using ValueType = typename Accessor::ValueType;
1347  ValueType V[6];
1348  V[0] = grid.getValue(ijk.offsetBy(0,0,-3));
1349  V[1] = grid.getValue(ijk.offsetBy(0,0,-2));
1350  V[2] = grid.getValue(ijk.offsetBy(0,0,-1));
1351  V[3] = grid.getValue(ijk);
1352  V[4] = grid.getValue(ijk.offsetBy(0,0,1));
1353  V[5] = grid.getValue(ijk.offsetBy(0,0,2));
1354 
1355  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1356  }
1357 
1358  // stencil access version
1359  template<typename Stencil>
1360  static typename Stencil::ValueType inX(const Stencil& S)
1361  {
1362  using ValueType = typename Stencil::ValueType;
1363  ValueType V[6];
1364  V[0] = S.template getValue<-3, 0, 0>();
1365  V[1] = S.template getValue<-2, 0, 0>();
1366  V[2] = S.template getValue<-1, 0, 0>();
1367  V[3] = S.template getValue< 0, 0, 0>();
1368  V[4] = S.template getValue< 1, 0, 0>();
1369  V[5] = S.template getValue< 2, 0, 0>();
1370 
1371  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1372  }
1373 
1374  template<typename Stencil>
1375  static typename Stencil::ValueType inY(const Stencil& S)
1376  {
1377  using ValueType = typename Stencil::ValueType;
1378  ValueType V[6];
1379  V[0] = S.template getValue< 0,-3, 0>();
1380  V[1] = S.template getValue< 0,-2, 0>();
1381  V[2] = S.template getValue< 0,-1, 0>();
1382  V[3] = S.template getValue< 0, 0, 0>();
1383  V[4] = S.template getValue< 0, 1, 0>();
1384  V[5] = S.template getValue< 0, 2, 0>();
1385 
1386  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1387  }
1388 
1389  template<typename Stencil>
1390  static typename Stencil::ValueType inZ(const Stencil& S)
1391  {
1392  using ValueType = typename Stencil::ValueType;
1393  ValueType V[6];
1394  V[0] = S.template getValue< 0, 0,-3>();
1395  V[1] = S.template getValue< 0, 0,-2>();
1396  V[2] = S.template getValue< 0, 0,-1>();
1397  V[3] = S.template getValue< 0, 0, 0>();
1398  V[4] = S.template getValue< 0, 0, 1>();
1399  V[5] = S.template getValue< 0, 0, 2>();
1400 
1401  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1402  }
1403 };
1404 
1405 
1406 template<>
1408 {
1409  template<typename ValueType>
1410  static ValueType difference(const ValueType& xm3, const ValueType& xm2, const ValueType& xm1,
1411  const ValueType& xm0, const ValueType& xp1, const ValueType& xp2)
1412  {
1413  return -D1<FD_HJWENO5>::difference(xm3, xm2, xm1, xm0, xp1, xp2);
1414  }
1415 
1416  // random access version
1417  template<typename Accessor>
1418  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1419  {
1420  using ValueType = typename Accessor::ValueType;
1421  ValueType V[6];
1422  V[0] = grid.getValue(ijk.offsetBy(-3,0,0));
1423  V[1] = grid.getValue(ijk.offsetBy(-2,0,0));
1424  V[2] = grid.getValue(ijk.offsetBy(-1,0,0));
1425  V[3] = grid.getValue(ijk);
1426  V[4] = grid.getValue(ijk.offsetBy(1,0,0));
1427  V[5] = grid.getValue(ijk.offsetBy(2,0,0));
1428 
1429  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1430  }
1431 
1432  template<typename Accessor>
1433  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1434  {
1435  using ValueType = typename Accessor::ValueType;
1436  ValueType V[6];
1437  V[0] = grid.getValue(ijk.offsetBy(0,-3,0));
1438  V[1] = grid.getValue(ijk.offsetBy(0,-2,0));
1439  V[2] = grid.getValue(ijk.offsetBy(0,-1,0));
1440  V[3] = grid.getValue(ijk);
1441  V[4] = grid.getValue(ijk.offsetBy(0,1,0));
1442  V[5] = grid.getValue(ijk.offsetBy(0,2,0));
1443 
1444  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1445  }
1446 
1447  template<typename Accessor>
1448  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1449  {
1450  using ValueType = typename Accessor::ValueType;
1451  ValueType V[6];
1452  V[0] = grid.getValue(ijk.offsetBy(0,0,-3));
1453  V[1] = grid.getValue(ijk.offsetBy(0,0,-2));
1454  V[2] = grid.getValue(ijk.offsetBy(0,0,-1));
1455  V[3] = grid.getValue(ijk);
1456  V[4] = grid.getValue(ijk.offsetBy(0,0,1));
1457  V[5] = grid.getValue(ijk.offsetBy(0,0,2));
1458 
1459  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1460  }
1461 
1462  // stencil access version
1463  template<typename Stencil>
1464  static typename Stencil::ValueType inX(const Stencil& S)
1465  {
1466  using ValueType = typename Stencil::ValueType;
1467  ValueType V[6];
1468  V[0] = S.template getValue<-3, 0, 0>();
1469  V[1] = S.template getValue<-2, 0, 0>();
1470  V[2] = S.template getValue<-1, 0, 0>();
1471  V[3] = S.template getValue< 0, 0, 0>();
1472  V[4] = S.template getValue< 1, 0, 0>();
1473  V[5] = S.template getValue< 2, 0, 0>();
1474 
1475  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1476  }
1477 
1478  template<typename Stencil>
1479  static typename Stencil::ValueType inY(const Stencil& S)
1480  {
1481  using ValueType = typename Stencil::ValueType;
1482  ValueType V[6];
1483  V[0] = S.template getValue< 0,-3, 0>();
1484  V[1] = S.template getValue< 0,-2, 0>();
1485  V[2] = S.template getValue< 0,-1, 0>();
1486  V[3] = S.template getValue< 0, 0, 0>();
1487  V[4] = S.template getValue< 0, 1, 0>();
1488  V[5] = S.template getValue< 0, 2, 0>();
1489 
1490  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1491  }
1492 
1493  template<typename Stencil>
1494  static typename Stencil::ValueType inZ(const Stencil& S)
1495  {
1496  using ValueType = typename Stencil::ValueType;
1497  ValueType V[6];
1498  V[0] = S.template getValue< 0, 0,-3>();
1499  V[1] = S.template getValue< 0, 0,-2>();
1500  V[2] = S.template getValue< 0, 0,-1>();
1501  V[3] = S.template getValue< 0, 0, 0>();
1502  V[4] = S.template getValue< 0, 0, 1>();
1503  V[5] = S.template getValue< 0, 0, 2>();
1504 
1505  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1506  }
1507 };
1508 
1509 
1510 template<DScheme DiffScheme>
1511 struct D1Vec
1512 {
1513  // random access version
1514  template<typename Accessor>
1515  static typename Accessor::ValueType::value_type
1516  inX(const Accessor& grid, const Coord& ijk, int n)
1517  {
1518  return D1<DiffScheme>::inX(grid, ijk)[n];
1519  }
1520 
1521  template<typename Accessor>
1522  static typename Accessor::ValueType::value_type
1523  inY(const Accessor& grid, const Coord& ijk, int n)
1524  {
1525  return D1<DiffScheme>::inY(grid, ijk)[n];
1526  }
1527  template<typename Accessor>
1528  static typename Accessor::ValueType::value_type
1529  inZ(const Accessor& grid, const Coord& ijk, int n)
1530  {
1531  return D1<DiffScheme>::inZ(grid, ijk)[n];
1532  }
1533 
1534 
1535  // stencil access version
1536  template<typename Stencil>
1537  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1538  {
1539  return D1<DiffScheme>::inX(S)[n];
1540  }
1541 
1542  template<typename Stencil>
1543  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1544  {
1545  return D1<DiffScheme>::inY(S)[n];
1546  }
1547 
1548  template<typename Stencil>
1549  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1550  {
1551  return D1<DiffScheme>::inZ(S)[n];
1552  }
1553 };
1554 
1555 
1556 template<>
1558 {
1559 
1560  // random access version
1561  template<typename Accessor>
1562  static typename Accessor::ValueType::value_type
1563  inX(const Accessor& grid, const Coord& ijk, int n)
1564  {
1565  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy( 1, 0, 0))[n],
1566  grid.getValue(ijk.offsetBy(-1, 0, 0))[n] );
1567  }
1568 
1569  template<typename Accessor>
1570  static typename Accessor::ValueType::value_type
1571  inY(const Accessor& grid, const Coord& ijk, int n)
1572  {
1573  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy(0, 1, 0))[n],
1574  grid.getValue(ijk.offsetBy(0,-1, 0))[n] );
1575  }
1576 
1577  template<typename Accessor>
1578  static typename Accessor::ValueType::value_type
1579  inZ(const Accessor& grid, const Coord& ijk, int n)
1580  {
1581  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy(0, 0, 1))[n],
1582  grid.getValue(ijk.offsetBy(0, 0,-1))[n] );
1583  }
1584 
1585  // stencil access version
1586  template<typename Stencil>
1587  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1588  {
1589  return D1<CD_2NDT>::difference( S.template getValue< 1, 0, 0>()[n],
1590  S.template getValue<-1, 0, 0>()[n] );
1591  }
1592 
1593  template<typename Stencil>
1594  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1595  {
1596  return D1<CD_2NDT>::difference( S.template getValue< 0, 1, 0>()[n],
1597  S.template getValue< 0,-1, 0>()[n] );
1598  }
1599 
1600  template<typename Stencil>
1601  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1602  {
1603  return D1<CD_2NDT>::difference( S.template getValue< 0, 0, 1>()[n],
1604  S.template getValue< 0, 0,-1>()[n] );
1605  }
1606 };
1607 
1608 template<>
1609 struct D1Vec<CD_2ND>
1610 {
1611 
1612  // random access version
1613  template<typename Accessor>
1614  static typename Accessor::ValueType::value_type
1615  inX(const Accessor& grid, const Coord& ijk, int n)
1616  {
1617  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy( 1, 0, 0))[n] ,
1618  grid.getValue(ijk.offsetBy(-1, 0, 0))[n] );
1619  }
1620 
1621  template<typename Accessor>
1622  static typename Accessor::ValueType::value_type
1623  inY(const Accessor& grid, const Coord& ijk, int n)
1624  {
1625  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy(0, 1, 0))[n] ,
1626  grid.getValue(ijk.offsetBy(0,-1, 0))[n] );
1627  }
1628 
1629  template<typename Accessor>
1630  static typename Accessor::ValueType::value_type
1631  inZ(const Accessor& grid, const Coord& ijk, int n)
1632  {
1633  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy(0, 0, 1))[n] ,
1634  grid.getValue(ijk.offsetBy(0, 0,-1))[n] );
1635  }
1636 
1637 
1638  // stencil access version
1639  template<typename Stencil>
1640  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1641  {
1642  return D1<CD_2ND>::difference( S.template getValue< 1, 0, 0>()[n],
1643  S.template getValue<-1, 0, 0>()[n] );
1644  }
1645 
1646  template<typename Stencil>
1647  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1648  {
1649  return D1<CD_2ND>::difference( S.template getValue< 0, 1, 0>()[n],
1650  S.template getValue< 0,-1, 0>()[n] );
1651  }
1652 
1653  template<typename Stencil>
1654  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1655  {
1656  return D1<CD_2ND>::difference( S.template getValue< 0, 0, 1>()[n],
1657  S.template getValue< 0, 0,-1>()[n] );
1658  }
1659 };
1660 
1661 
1662 template<>
1663 struct D1Vec<CD_4TH> {
1664  // using value_type = typename Accessor::ValueType::value_type;
1665 
1666 
1667  // random access version
1668  template<typename Accessor>
1669  static typename Accessor::ValueType::value_type
1670  inX(const Accessor& grid, const Coord& ijk, int n)
1671  {
1672  return D1<CD_4TH>::difference(
1673  grid.getValue(ijk.offsetBy(2, 0, 0))[n], grid.getValue(ijk.offsetBy( 1, 0, 0))[n],
1674  grid.getValue(ijk.offsetBy(-1,0, 0))[n], grid.getValue(ijk.offsetBy(-2, 0, 0))[n]);
1675  }
1676 
1677  template<typename Accessor>
1678  static typename Accessor::ValueType::value_type
1679  inY(const Accessor& grid, const Coord& ijk, int n)
1680  {
1681  return D1<CD_4TH>::difference(
1682  grid.getValue(ijk.offsetBy( 0, 2, 0))[n], grid.getValue(ijk.offsetBy( 0, 1, 0))[n],
1683  grid.getValue(ijk.offsetBy( 0,-1, 0))[n], grid.getValue(ijk.offsetBy( 0,-2, 0))[n]);
1684  }
1685 
1686  template<typename Accessor>
1687  static typename Accessor::ValueType::value_type
1688  inZ(const Accessor& grid, const Coord& ijk, int n)
1689  {
1690  return D1<CD_4TH>::difference(
1691  grid.getValue(ijk.offsetBy(0,0, 2))[n], grid.getValue(ijk.offsetBy( 0, 0, 1))[n],
1692  grid.getValue(ijk.offsetBy(0,0,-1))[n], grid.getValue(ijk.offsetBy( 0, 0,-2))[n]);
1693  }
1694 
1695  // stencil access version
1696  template<typename Stencil>
1697  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1698  {
1699  return D1<CD_4TH>::difference(
1700  S.template getValue< 2, 0, 0>()[n], S.template getValue< 1, 0, 0>()[n],
1701  S.template getValue<-1, 0, 0>()[n], S.template getValue<-2, 0, 0>()[n] );
1702  }
1703 
1704  template<typename Stencil>
1705  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1706  {
1707  return D1<CD_4TH>::difference(
1708  S.template getValue< 0, 2, 0>()[n], S.template getValue< 0, 1, 0>()[n],
1709  S.template getValue< 0,-1, 0>()[n], S.template getValue< 0,-2, 0>()[n]);
1710  }
1711 
1712  template<typename Stencil>
1713  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1714  {
1715  return D1<CD_4TH>::difference(
1716  S.template getValue< 0, 0, 2>()[n], S.template getValue< 0, 0, 1>()[n],
1717  S.template getValue< 0, 0,-1>()[n], S.template getValue< 0, 0,-2>()[n]);
1718  }
1719 };
1720 
1721 
1722 template<>
1723 struct D1Vec<CD_6TH>
1724 {
1725  //using ValueType = typename Accessor::ValueType::value_type::value_type;
1726 
1727  // random access version
1728  template<typename Accessor>
1729  static typename Accessor::ValueType::value_type
1730  inX(const Accessor& grid, const Coord& ijk, int n)
1731  {
1732  return D1<CD_6TH>::difference(
1733  grid.getValue(ijk.offsetBy( 3, 0, 0))[n], grid.getValue(ijk.offsetBy( 2, 0, 0))[n],
1734  grid.getValue(ijk.offsetBy( 1, 0, 0))[n], grid.getValue(ijk.offsetBy(-1, 0, 0))[n],
1735  grid.getValue(ijk.offsetBy(-2, 0, 0))[n], grid.getValue(ijk.offsetBy(-3, 0, 0))[n] );
1736  }
1737 
1738  template<typename Accessor>
1739  static typename Accessor::ValueType::value_type
1740  inY(const Accessor& grid, const Coord& ijk, int n)
1741  {
1742  return D1<CD_6TH>::difference(
1743  grid.getValue(ijk.offsetBy( 0, 3, 0))[n], grid.getValue(ijk.offsetBy( 0, 2, 0))[n],
1744  grid.getValue(ijk.offsetBy( 0, 1, 0))[n], grid.getValue(ijk.offsetBy( 0,-1, 0))[n],
1745  grid.getValue(ijk.offsetBy( 0,-2, 0))[n], grid.getValue(ijk.offsetBy( 0,-3, 0))[n] );
1746  }
1747 
1748  template<typename Accessor>
1749  static typename Accessor::ValueType::value_type
1750  inZ(const Accessor& grid, const Coord& ijk, int n)
1751  {
1752  return D1<CD_6TH>::difference(
1753  grid.getValue(ijk.offsetBy( 0, 0, 3))[n], grid.getValue(ijk.offsetBy( 0, 0, 2))[n],
1754  grid.getValue(ijk.offsetBy( 0, 0, 1))[n], grid.getValue(ijk.offsetBy( 0, 0,-1))[n],
1755  grid.getValue(ijk.offsetBy( 0, 0,-2))[n], grid.getValue(ijk.offsetBy( 0, 0,-3))[n] );
1756  }
1757 
1758 
1759  // stencil access version
1760  template<typename Stencil>
1761  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1762  {
1763  return D1<CD_6TH>::difference(
1764  S.template getValue< 3, 0, 0>()[n], S.template getValue< 2, 0, 0>()[n],
1765  S.template getValue< 1, 0, 0>()[n], S.template getValue<-1, 0, 0>()[n],
1766  S.template getValue<-2, 0, 0>()[n], S.template getValue<-3, 0, 0>()[n] );
1767  }
1768 
1769  template<typename Stencil>
1770  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1771  {
1772  return D1<CD_6TH>::difference(
1773  S.template getValue< 0, 3, 0>()[n], S.template getValue< 0, 2, 0>()[n],
1774  S.template getValue< 0, 1, 0>()[n], S.template getValue< 0,-1, 0>()[n],
1775  S.template getValue< 0,-2, 0>()[n], S.template getValue< 0,-3, 0>()[n] );
1776  }
1777 
1778  template<typename Stencil>
1779  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1780  {
1781  return D1<CD_6TH>::difference(
1782  S.template getValue< 0, 0, 3>()[n], S.template getValue< 0, 0, 2>()[n],
1783  S.template getValue< 0, 0, 1>()[n], S.template getValue< 0, 0,-1>()[n],
1784  S.template getValue< 0, 0,-2>()[n], S.template getValue< 0, 0,-3>()[n] );
1785  }
1786 };
1787 
1788 template<DDScheme DiffScheme>
1789 struct D2
1790 {
1791 
1792  template<typename Accessor>
1793  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk);
1794  template<typename Accessor>
1795  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk);
1796  template<typename Accessor>
1797  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk);
1798 
1799  // cross derivatives
1800  template<typename Accessor>
1801  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk);
1802 
1803  template<typename Accessor>
1804  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk);
1805 
1806  template<typename Accessor>
1807  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk);
1808 
1809 
1810  // stencil access version
1811  template<typename Stencil>
1812  static typename Stencil::ValueType inX(const Stencil& S);
1813  template<typename Stencil>
1814  static typename Stencil::ValueType inY(const Stencil& S);
1815  template<typename Stencil>
1816  static typename Stencil::ValueType inZ(const Stencil& S);
1817 
1818  // cross derivatives
1819  template<typename Stencil>
1820  static typename Stencil::ValueType inXandY(const Stencil& S);
1821 
1822  template<typename Stencil>
1823  static typename Stencil::ValueType inXandZ(const Stencil& S);
1824 
1825  template<typename Stencil>
1826  static typename Stencil::ValueType inYandZ(const Stencil& S);
1827 };
1828 
1829 template<>
1830 struct D2<CD_SECOND>
1831 {
1832 
1833  // the difference opperator
1834  template <typename ValueType>
1835  static ValueType difference(const ValueType& xp1, const ValueType& xp0, const ValueType& xm1)
1836  {
1837  return xp1 + xm1 - ValueType(2)*xp0;
1838  }
1839 
1840  template <typename ValueType>
1841  static ValueType crossdifference(const ValueType& xpyp, const ValueType& xpym,
1842  const ValueType& xmyp, const ValueType& xmym)
1843  {
1844  return ValueType(0.25)*(xpyp + xmym - xpym - xmyp);
1845  }
1846 
1847  // random access version
1848  template<typename Accessor>
1849  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1850  {
1851  return difference( grid.getValue(ijk.offsetBy( 1,0,0)), grid.getValue(ijk),
1852  grid.getValue(ijk.offsetBy(-1,0,0)) );
1853  }
1854 
1855  template<typename Accessor>
1856  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1857  {
1858 
1859  return difference( grid.getValue(ijk.offsetBy(0, 1,0)), grid.getValue(ijk),
1860  grid.getValue(ijk.offsetBy(0,-1,0)) );
1861  }
1862 
1863  template<typename Accessor>
1864  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1865  {
1866  return difference( grid.getValue(ijk.offsetBy( 0,0, 1)), grid.getValue(ijk),
1867  grid.getValue(ijk.offsetBy( 0,0,-1)) );
1868  }
1869 
1870  // cross derivatives
1871  template<typename Accessor>
1872  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
1873  {
1874  return crossdifference(
1875  grid.getValue(ijk.offsetBy(1, 1,0)), grid.getValue(ijk.offsetBy( 1,-1,0)),
1876  grid.getValue(ijk.offsetBy(-1,1,0)), grid.getValue(ijk.offsetBy(-1,-1,0)));
1877 
1878  }
1879 
1880  template<typename Accessor>
1881  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
1882  {
1883  return crossdifference(
1884  grid.getValue(ijk.offsetBy(1,0, 1)), grid.getValue(ijk.offsetBy(1, 0,-1)),
1885  grid.getValue(ijk.offsetBy(-1,0,1)), grid.getValue(ijk.offsetBy(-1,0,-1)) );
1886  }
1887 
1888  template<typename Accessor>
1889  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
1890  {
1891  return crossdifference(
1892  grid.getValue(ijk.offsetBy(0, 1,1)), grid.getValue(ijk.offsetBy(0, 1,-1)),
1893  grid.getValue(ijk.offsetBy(0,-1,1)), grid.getValue(ijk.offsetBy(0,-1,-1)) );
1894  }
1895 
1896 
1897  // stencil access version
1898  template<typename Stencil>
1899  static typename Stencil::ValueType inX(const Stencil& S)
1900  {
1901  return difference( S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>(),
1902  S.template getValue<-1, 0, 0>() );
1903  }
1904 
1905  template<typename Stencil>
1906  static typename Stencil::ValueType inY(const Stencil& S)
1907  {
1908  return difference( S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>(),
1909  S.template getValue< 0,-1, 0>() );
1910  }
1911 
1912  template<typename Stencil>
1913  static typename Stencil::ValueType inZ(const Stencil& S)
1914  {
1915  return difference( S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>(),
1916  S.template getValue< 0, 0,-1>() );
1917  }
1918 
1919  // cross derivatives
1920  template<typename Stencil>
1921  static typename Stencil::ValueType inXandY(const Stencil& S)
1922  {
1923  return crossdifference(S.template getValue< 1, 1, 0>(), S.template getValue< 1,-1, 0>(),
1924  S.template getValue<-1, 1, 0>(), S.template getValue<-1,-1, 0>() );
1925  }
1926 
1927  template<typename Stencil>
1928  static typename Stencil::ValueType inXandZ(const Stencil& S)
1929  {
1930  return crossdifference(S.template getValue< 1, 0, 1>(), S.template getValue< 1, 0,-1>(),
1931  S.template getValue<-1, 0, 1>(), S.template getValue<-1, 0,-1>() );
1932  }
1933 
1934  template<typename Stencil>
1935  static typename Stencil::ValueType inYandZ(const Stencil& S)
1936  {
1937  return crossdifference(S.template getValue< 0, 1, 1>(), S.template getValue< 0, 1,-1>(),
1938  S.template getValue< 0,-1, 1>(), S.template getValue< 0,-1,-1>() );
1939  }
1940 };
1941 
1942 
1943 template<>
1944 struct D2<CD_FOURTH>
1945 {
1946 
1947  // the difference opperator
1948  template <typename ValueType>
1949  static ValueType difference(const ValueType& xp2, const ValueType& xp1, const ValueType& xp0,
1950  const ValueType& xm1, const ValueType& xm2) {
1951  return ValueType(-1./12.)*(xp2 + xm2) + ValueType(4./3.)*(xp1 + xm1) -ValueType(2.5)*xp0;
1952  }
1953 
1954  template <typename ValueType>
1955  static ValueType crossdifference(const ValueType& xp2yp2, const ValueType& xp2yp1,
1956  const ValueType& xp2ym1, const ValueType& xp2ym2,
1957  const ValueType& xp1yp2, const ValueType& xp1yp1,
1958  const ValueType& xp1ym1, const ValueType& xp1ym2,
1959  const ValueType& xm2yp2, const ValueType& xm2yp1,
1960  const ValueType& xm2ym1, const ValueType& xm2ym2,
1961  const ValueType& xm1yp2, const ValueType& xm1yp1,
1962  const ValueType& xm1ym1, const ValueType& xm1ym2 ) {
1963  ValueType tmp1 =
1964  ValueType(2./3.0)*(xp1yp1 - xm1yp1 - xp1ym1 + xm1ym1)-
1965  ValueType(1./12.)*(xp2yp1 - xm2yp1 - xp2ym1 + xm2ym1);
1966  ValueType tmp2 =
1967  ValueType(2./3.0)*(xp1yp2 - xm1yp2 - xp1ym2 + xm1ym2)-
1968  ValueType(1./12.)*(xp2yp2 - xm2yp2 - xp2ym2 + xm2ym2);
1969 
1970  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
1971  }
1972 
1973 
1974 
1975  // random access version
1976  template<typename Accessor>
1977  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1978  {
1979  return difference(
1980  grid.getValue(ijk.offsetBy(2,0,0)), grid.getValue(ijk.offsetBy( 1,0,0)),
1981  grid.getValue(ijk),
1982  grid.getValue(ijk.offsetBy(-1,0,0)), grid.getValue(ijk.offsetBy(-2, 0, 0)));
1983  }
1984 
1985  template<typename Accessor>
1986  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1987  {
1988  return difference(
1989  grid.getValue(ijk.offsetBy(0, 2,0)), grid.getValue(ijk.offsetBy(0, 1,0)),
1990  grid.getValue(ijk),
1991  grid.getValue(ijk.offsetBy(0,-1,0)), grid.getValue(ijk.offsetBy(0,-2, 0)));
1992  }
1993 
1994  template<typename Accessor>
1995  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1996  {
1997  return difference(
1998  grid.getValue(ijk.offsetBy(0,0, 2)), grid.getValue(ijk.offsetBy(0, 0,1)),
1999  grid.getValue(ijk),
2000  grid.getValue(ijk.offsetBy(0,0,-1)), grid.getValue(ijk.offsetBy(0,0,-2)));
2001  }
2002 
2003  // cross derivatives
2004  template<typename Accessor>
2005  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
2006  {
2007  using ValueType = typename Accessor::ValueType;
2008  typename Accessor::ValueType tmp1 =
2009  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 1, 0)) -
2010  D1<CD_4TH>::inX(grid, ijk.offsetBy(0,-1, 0));
2011  typename Accessor::ValueType tmp2 =
2012  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 2, 0)) -
2013  D1<CD_4TH>::inX(grid, ijk.offsetBy(0,-2, 0));
2014  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
2015  }
2016 
2017  template<typename Accessor>
2018  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
2019  {
2020  using ValueType = typename Accessor::ValueType;
2021  typename Accessor::ValueType tmp1 =
2022  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0, 1)) -
2023  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0,-1));
2024  typename Accessor::ValueType tmp2 =
2025  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0, 2)) -
2026  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0,-2));
2027  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
2028  }
2029 
2030  template<typename Accessor>
2031  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
2032  {
2033  using ValueType = typename Accessor::ValueType;
2034  typename Accessor::ValueType tmp1 =
2035  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0, 1)) -
2036  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0,-1));
2037  typename Accessor::ValueType tmp2 =
2038  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0, 2)) -
2039  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0,-2));
2040  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
2041  }
2042 
2043 
2044  // stencil access version
2045  template<typename Stencil>
2046  static typename Stencil::ValueType inX(const Stencil& S)
2047  {
2048  return difference(S.template getValue< 2, 0, 0>(), S.template getValue< 1, 0, 0>(),
2049  S.template getValue< 0, 0, 0>(),
2050  S.template getValue<-1, 0, 0>(), S.template getValue<-2, 0, 0>() );
2051  }
2052 
2053  template<typename Stencil>
2054  static typename Stencil::ValueType inY(const Stencil& S)
2055  {
2056  return difference(S.template getValue< 0, 2, 0>(), S.template getValue< 0, 1, 0>(),
2057  S.template getValue< 0, 0, 0>(),
2058  S.template getValue< 0,-1, 0>(), S.template getValue< 0,-2, 0>() );
2059  }
2060 
2061  template<typename Stencil>
2062  static typename Stencil::ValueType inZ(const Stencil& S)
2063  {
2064  return difference(S.template getValue< 0, 0, 2>(), S.template getValue< 0, 0, 1>(),
2065  S.template getValue< 0, 0, 0>(),
2066  S.template getValue< 0, 0,-1>(), S.template getValue< 0, 0,-2>() );
2067  }
2068 
2069  // cross derivatives
2070  template<typename Stencil>
2071  static typename Stencil::ValueType inXandY(const Stencil& S)
2072  {
2073  return crossdifference(
2074  S.template getValue< 2, 2, 0>(), S.template getValue< 2, 1, 0>(),
2075  S.template getValue< 2,-1, 0>(), S.template getValue< 2,-2, 0>(),
2076  S.template getValue< 1, 2, 0>(), S.template getValue< 1, 1, 0>(),
2077  S.template getValue< 1,-1, 0>(), S.template getValue< 1,-2, 0>(),
2078  S.template getValue<-2, 2, 0>(), S.template getValue<-2, 1, 0>(),
2079  S.template getValue<-2,-1, 0>(), S.template getValue<-2,-2, 0>(),
2080  S.template getValue<-1, 2, 0>(), S.template getValue<-1, 1, 0>(),
2081  S.template getValue<-1,-1, 0>(), S.template getValue<-1,-2, 0>() );
2082  }
2083 
2084  template<typename Stencil>
2085  static typename Stencil::ValueType inXandZ(const Stencil& S)
2086  {
2087  return crossdifference(
2088  S.template getValue< 2, 0, 2>(), S.template getValue< 2, 0, 1>(),
2089  S.template getValue< 2, 0,-1>(), S.template getValue< 2, 0,-2>(),
2090  S.template getValue< 1, 0, 2>(), S.template getValue< 1, 0, 1>(),
2091  S.template getValue< 1, 0,-1>(), S.template getValue< 1, 0,-2>(),
2092  S.template getValue<-2, 0, 2>(), S.template getValue<-2, 0, 1>(),
2093  S.template getValue<-2, 0,-1>(), S.template getValue<-2, 0,-2>(),
2094  S.template getValue<-1, 0, 2>(), S.template getValue<-1, 0, 1>(),
2095  S.template getValue<-1, 0,-1>(), S.template getValue<-1, 0,-2>() );
2096  }
2097 
2098  template<typename Stencil>
2099  static typename Stencil::ValueType inYandZ(const Stencil& S)
2100  {
2101  return crossdifference(
2102  S.template getValue< 0, 2, 2>(), S.template getValue< 0, 2, 1>(),
2103  S.template getValue< 0, 2,-1>(), S.template getValue< 0, 2,-2>(),
2104  S.template getValue< 0, 1, 2>(), S.template getValue< 0, 1, 1>(),
2105  S.template getValue< 0, 1,-1>(), S.template getValue< 0, 1,-2>(),
2106  S.template getValue< 0,-2, 2>(), S.template getValue< 0,-2, 1>(),
2107  S.template getValue< 0,-2,-1>(), S.template getValue< 0,-2,-2>(),
2108  S.template getValue< 0,-1, 2>(), S.template getValue< 0,-1, 1>(),
2109  S.template getValue< 0,-1,-1>(), S.template getValue< 0,-1,-2>() );
2110  }
2111 };
2112 
2113 
2114 template<>
2115 struct D2<CD_SIXTH>
2116 {
2117  // the difference opperator
2118  template <typename ValueType>
2119  static ValueType difference(const ValueType& xp3, const ValueType& xp2, const ValueType& xp1,
2120  const ValueType& xp0,
2121  const ValueType& xm1, const ValueType& xm2, const ValueType& xm3)
2122  {
2123  return ValueType(1./90.)*(xp3 + xm3) - ValueType(3./20.)*(xp2 + xm2)
2124  + ValueType(1.5)*(xp1 + xm1) - ValueType(49./18.)*xp0;
2125  }
2126 
2127  template <typename ValueType>
2128  static ValueType crossdifference( const ValueType& xp1yp1,const ValueType& xm1yp1,
2129  const ValueType& xp1ym1,const ValueType& xm1ym1,
2130  const ValueType& xp2yp1,const ValueType& xm2yp1,
2131  const ValueType& xp2ym1,const ValueType& xm2ym1,
2132  const ValueType& xp3yp1,const ValueType& xm3yp1,
2133  const ValueType& xp3ym1,const ValueType& xm3ym1,
2134  const ValueType& xp1yp2,const ValueType& xm1yp2,
2135  const ValueType& xp1ym2,const ValueType& xm1ym2,
2136  const ValueType& xp2yp2,const ValueType& xm2yp2,
2137  const ValueType& xp2ym2,const ValueType& xm2ym2,
2138  const ValueType& xp3yp2,const ValueType& xm3yp2,
2139  const ValueType& xp3ym2,const ValueType& xm3ym2,
2140  const ValueType& xp1yp3,const ValueType& xm1yp3,
2141  const ValueType& xp1ym3,const ValueType& xm1ym3,
2142  const ValueType& xp2yp3,const ValueType& xm2yp3,
2143  const ValueType& xp2ym3,const ValueType& xm2ym3,
2144  const ValueType& xp3yp3,const ValueType& xm3yp3,
2145  const ValueType& xp3ym3,const ValueType& xm3ym3 )
2146  {
2147  ValueType tmp1 =
2148  ValueType(0.7500)*(xp1yp1 - xm1yp1 - xp1ym1 + xm1ym1) -
2149  ValueType(0.1500)*(xp2yp1 - xm2yp1 - xp2ym1 + xm2ym1) +
2150  ValueType(1./60.)*(xp3yp1 - xm3yp1 - xp3ym1 + xm3ym1);
2151 
2152  ValueType tmp2 =
2153  ValueType(0.7500)*(xp1yp2 - xm1yp2 - xp1ym2 + xm1ym2) -
2154  ValueType(0.1500)*(xp2yp2 - xm2yp2 - xp2ym2 + xm2ym2) +
2155  ValueType(1./60.)*(xp3yp2 - xm3yp2 - xp3ym2 + xm3ym2);
2156 
2157  ValueType tmp3 =
2158  ValueType(0.7500)*(xp1yp3 - xm1yp3 - xp1ym3 + xm1ym3) -
2159  ValueType(0.1500)*(xp2yp3 - xm2yp3 - xp2ym3 + xm2ym3) +
2160  ValueType(1./60.)*(xp3yp3 - xm3yp3 - xp3ym3 + xm3ym3);
2161 
2162  return ValueType(0.75)*tmp1 - ValueType(0.15)*tmp2 + ValueType(1./60)*tmp3;
2163  }
2164 
2165  // random access version
2166 
2167  template<typename Accessor>
2168  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
2169  {
2170  return difference(
2171  grid.getValue(ijk.offsetBy( 3, 0, 0)), grid.getValue(ijk.offsetBy( 2, 0, 0)),
2172  grid.getValue(ijk.offsetBy( 1, 0, 0)), grid.getValue(ijk),
2173  grid.getValue(ijk.offsetBy(-1, 0, 0)), grid.getValue(ijk.offsetBy(-2, 0, 0)),
2174  grid.getValue(ijk.offsetBy(-3, 0, 0)) );
2175  }
2176 
2177  template<typename Accessor>
2178  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
2179  {
2180  return difference(
2181  grid.getValue(ijk.offsetBy( 0, 3, 0)), grid.getValue(ijk.offsetBy( 0, 2, 0)),
2182  grid.getValue(ijk.offsetBy( 0, 1, 0)), grid.getValue(ijk),
2183  grid.getValue(ijk.offsetBy( 0,-1, 0)), grid.getValue(ijk.offsetBy( 0,-2, 0)),
2184  grid.getValue(ijk.offsetBy( 0,-3, 0)) );
2185  }
2186 
2187  template<typename Accessor>
2188  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
2189  {
2190 
2191  return difference(
2192  grid.getValue(ijk.offsetBy( 0, 0, 3)), grid.getValue(ijk.offsetBy( 0, 0, 2)),
2193  grid.getValue(ijk.offsetBy( 0, 0, 1)), grid.getValue(ijk),
2194  grid.getValue(ijk.offsetBy( 0, 0,-1)), grid.getValue(ijk.offsetBy( 0, 0,-2)),
2195  grid.getValue(ijk.offsetBy( 0, 0,-3)) );
2196  }
2197 
2198  template<typename Accessor>
2199  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
2200  {
2201  using ValueT = typename Accessor::ValueType;
2202  ValueT tmp1 =
2203  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 1, 0)) -
2204  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-1, 0));
2205  ValueT tmp2 =
2206  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 2, 0)) -
2207  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-2, 0));
2208  ValueT tmp3 =
2209  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 3, 0)) -
2210  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-3, 0));
2211  return ValueT(0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3);
2212  }
2213 
2214  template<typename Accessor>
2215  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
2216  {
2217  using ValueT = typename Accessor::ValueType;
2218  ValueT tmp1 =
2219  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 1)) -
2220  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-1));
2221  ValueT tmp2 =
2222  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 2)) -
2223  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-2));
2224  ValueT tmp3 =
2225  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 3)) -
2226  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-3));
2227  return ValueT(0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3);
2228  }
2229 
2230  template<typename Accessor>
2231  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
2232  {
2233  using ValueT = typename Accessor::ValueType;
2234  ValueT tmp1 =
2235  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 1)) -
2236  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-1));
2237  ValueT tmp2 =
2238  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 2)) -
2239  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-2));
2240  ValueT tmp3 =
2241  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 3)) -
2242  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-3));
2243  return ValueT(0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3);
2244  }
2245 
2246 
2247  // stencil access version
2248  template<typename Stencil>
2249  static typename Stencil::ValueType inX(const Stencil& S)
2250  {
2251  return difference( S.template getValue< 3, 0, 0>(), S.template getValue< 2, 0, 0>(),
2252  S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>(),
2253  S.template getValue<-1, 0, 0>(), S.template getValue<-2, 0, 0>(),
2254  S.template getValue<-3, 0, 0>() );
2255  }
2256 
2257  template<typename Stencil>
2258  static typename Stencil::ValueType inY(const Stencil& S)
2259  {
2260  return difference( S.template getValue< 0, 3, 0>(), S.template getValue< 0, 2, 0>(),
2261  S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>(),
2262  S.template getValue< 0,-1, 0>(), S.template getValue< 0,-2, 0>(),
2263  S.template getValue< 0,-3, 0>() );
2264 
2265  }
2266 
2267  template<typename Stencil>
2268  static typename Stencil::ValueType inZ(const Stencil& S)
2269  {
2270  return difference( S.template getValue< 0, 0, 3>(), S.template getValue< 0, 0, 2>(),
2271  S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>(),
2272  S.template getValue< 0, 0,-1>(), S.template getValue< 0, 0,-2>(),
2273  S.template getValue< 0, 0,-3>() );
2274  }
2275 
2276  template<typename Stencil>
2277  static typename Stencil::ValueType inXandY(const Stencil& S)
2278  {
2279  return crossdifference( S.template getValue< 1, 1, 0>(), S.template getValue<-1, 1, 0>(),
2280  S.template getValue< 1,-1, 0>(), S.template getValue<-1,-1, 0>(),
2281  S.template getValue< 2, 1, 0>(), S.template getValue<-2, 1, 0>(),
2282  S.template getValue< 2,-1, 0>(), S.template getValue<-2,-1, 0>(),
2283  S.template getValue< 3, 1, 0>(), S.template getValue<-3, 1, 0>(),
2284  S.template getValue< 3,-1, 0>(), S.template getValue<-3,-1, 0>(),
2285  S.template getValue< 1, 2, 0>(), S.template getValue<-1, 2, 0>(),
2286  S.template getValue< 1,-2, 0>(), S.template getValue<-1,-2, 0>(),
2287  S.template getValue< 2, 2, 0>(), S.template getValue<-2, 2, 0>(),
2288  S.template getValue< 2,-2, 0>(), S.template getValue<-2,-2, 0>(),
2289  S.template getValue< 3, 2, 0>(), S.template getValue<-3, 2, 0>(),
2290  S.template getValue< 3,-2, 0>(), S.template getValue<-3,-2, 0>(),
2291  S.template getValue< 1, 3, 0>(), S.template getValue<-1, 3, 0>(),
2292  S.template getValue< 1,-3, 0>(), S.template getValue<-1,-3, 0>(),
2293  S.template getValue< 2, 3, 0>(), S.template getValue<-2, 3, 0>(),
2294  S.template getValue< 2,-3, 0>(), S.template getValue<-2,-3, 0>(),
2295  S.template getValue< 3, 3, 0>(), S.template getValue<-3, 3, 0>(),
2296  S.template getValue< 3,-3, 0>(), S.template getValue<-3,-3, 0>() );
2297  }
2298 
2299  template<typename Stencil>
2300  static typename Stencil::ValueType inXandZ(const Stencil& S)
2301  {
2302  return crossdifference( S.template getValue< 1, 0, 1>(), S.template getValue<-1, 0, 1>(),
2303  S.template getValue< 1, 0,-1>(), S.template getValue<-1, 0,-1>(),
2304  S.template getValue< 2, 0, 1>(), S.template getValue<-2, 0, 1>(),
2305  S.template getValue< 2, 0,-1>(), S.template getValue<-2, 0,-1>(),
2306  S.template getValue< 3, 0, 1>(), S.template getValue<-3, 0, 1>(),
2307  S.template getValue< 3, 0,-1>(), S.template getValue<-3, 0,-1>(),
2308  S.template getValue< 1, 0, 2>(), S.template getValue<-1, 0, 2>(),
2309  S.template getValue< 1, 0,-2>(), S.template getValue<-1, 0,-2>(),
2310  S.template getValue< 2, 0, 2>(), S.template getValue<-2, 0, 2>(),
2311  S.template getValue< 2, 0,-2>(), S.template getValue<-2, 0,-2>(),
2312  S.template getValue< 3, 0, 2>(), S.template getValue<-3, 0, 2>(),
2313  S.template getValue< 3, 0,-2>(), S.template getValue<-3, 0,-2>(),
2314  S.template getValue< 1, 0, 3>(), S.template getValue<-1, 0, 3>(),
2315  S.template getValue< 1, 0,-3>(), S.template getValue<-1, 0,-3>(),
2316  S.template getValue< 2, 0, 3>(), S.template getValue<-2, 0, 3>(),
2317  S.template getValue< 2, 0,-3>(), S.template getValue<-2, 0,-3>(),
2318  S.template getValue< 3, 0, 3>(), S.template getValue<-3, 0, 3>(),
2319  S.template getValue< 3, 0,-3>(), S.template getValue<-3, 0,-3>() );
2320  }
2321 
2322  template<typename Stencil>
2323  static typename Stencil::ValueType inYandZ(const Stencil& S)
2324  {
2325  return crossdifference( S.template getValue< 0, 1, 1>(), S.template getValue< 0,-1, 1>(),
2326  S.template getValue< 0, 1,-1>(), S.template getValue< 0,-1,-1>(),
2327  S.template getValue< 0, 2, 1>(), S.template getValue< 0,-2, 1>(),
2328  S.template getValue< 0, 2,-1>(), S.template getValue< 0,-2,-1>(),
2329  S.template getValue< 0, 3, 1>(), S.template getValue< 0,-3, 1>(),
2330  S.template getValue< 0, 3,-1>(), S.template getValue< 0,-3,-1>(),
2331  S.template getValue< 0, 1, 2>(), S.template getValue< 0,-1, 2>(),
2332  S.template getValue< 0, 1,-2>(), S.template getValue< 0,-1,-2>(),
2333  S.template getValue< 0, 2, 2>(), S.template getValue< 0,-2, 2>(),
2334  S.template getValue< 0, 2,-2>(), S.template getValue< 0,-2,-2>(),
2335  S.template getValue< 0, 3, 2>(), S.template getValue< 0,-3, 2>(),
2336  S.template getValue< 0, 3,-2>(), S.template getValue< 0,-3,-2>(),
2337  S.template getValue< 0, 1, 3>(), S.template getValue< 0,-1, 3>(),
2338  S.template getValue< 0, 1,-3>(), S.template getValue< 0,-1,-3>(),
2339  S.template getValue< 0, 2, 3>(), S.template getValue< 0,-2, 3>(),
2340  S.template getValue< 0, 2,-3>(), S.template getValue< 0,-2,-3>(),
2341  S.template getValue< 0, 3, 3>(), S.template getValue< 0,-3, 3>(),
2342  S.template getValue< 0, 3,-3>(), S.template getValue< 0,-3,-3>() );
2343  }
2344 
2345 };
2346 
2347 } // end math namespace
2348 } // namespace OPENVDB_VERSION_NAME
2349 } // end openvdb namespace
2350 
2351 #endif // OPENVDB_MATH_FINITEDIFFERENCE_HAS_BEEN_INCLUDED
2352 
2353 // Copyright (c) 2012-2018 DreamWorks Animation LLC
2354 // All rights reserved. This software is distributed under the
2355 // 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 x2.
Definition: Math.h:502
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)
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)
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:189
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:51
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)
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:610
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)
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)
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:135
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:549
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)