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