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