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