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