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