OpenVDB 11.0.0
Loading...
Searching...
No Matches
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
20namespace openvdb {
22namespace OPENVDB_VERSION_NAME {
23namespace 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.
31enum DScheme {
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
50
51
52inline 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
76inline DScheme
77stringToDScheme(const std::string& s)
78{
79 DScheme ret = UNKNOWN_DS;
80
81 std::string str = s;
82 openvdb::string::trim(str);
83 openvdb::string::to_lower(str);
84
85 if (str == dsSchemeToString(CD_2NDT)) {
86 ret = CD_2NDT;
87 } else if (str == dsSchemeToString(CD_2ND)) {
88 ret = CD_2ND;
89 } else if (str == dsSchemeToString(CD_4TH)) {
90 ret = CD_4TH;
91 } else if (str == dsSchemeToString(CD_6TH)) {
92 ret = CD_6TH;
93 } else if (str == dsSchemeToString(FD_1ST)) {
94 ret = FD_1ST;
95 } else if (str == dsSchemeToString(FD_2ND)) {
96 ret = FD_2ND;
97 } else if (str == dsSchemeToString(FD_3RD)) {
98 ret = FD_3RD;
99 } else if (str == dsSchemeToString(BD_1ST)) {
100 ret = BD_1ST;
101 } else if (str == dsSchemeToString(BD_2ND)) {
102 ret = BD_2ND;
103 } else if (str == dsSchemeToString(BD_3RD)) {
104 ret = BD_3RD;
105 } else if (str == dsSchemeToString(FD_WENO5)) {
106 ret = FD_WENO5;
107 } else if (str == dsSchemeToString(BD_WENO5)) {
108 ret = BD_WENO5;
109 } else if (str == dsSchemeToString(FD_HJWENO5)) {
110 ret = FD_HJWENO5;
111 } else if (str == dsSchemeToString(BD_HJWENO5)) {
112 ret = BD_HJWENO5;
113 }
114
115 return ret;
116}
117
118inline 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.
151 CD_SECOND = 0, // center difference, 2nd order
152 CD_FOURTH, // center difference, 4th order
153 CD_SIXTH // center difference, 6th order
155
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
172
174
175inline 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
190inline BiasedGradientScheme
191stringToBiasedGradientScheme(const std::string& s)
192{
194
195 std::string str = s;
196 openvdb::string::trim(str);
197 openvdb::string::to_lower(str);
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
213inline 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
237 TVD_RK3
239
241
242inline 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
255inline TemporalIntegrationScheme
257{
259
260 std::string str = s;
261 openvdb::string::trim(str);
262 openvdb::string::to_lower(str);
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
275inline 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)
301template<typename ValueType>
302inline ValueType
303WENO5(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
324template <typename Real>
325inline 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
349template<typename Real>
350inline Real
351GodunovsNormSqrd(bool isOutside, const Vec3<Real>& gradient_m, const Vec3<Real>& gradient_p)
352{
353 return GodunovsNormSqrd<Real>(isOutside,
354 gradient_m[0], gradient_p[0],
355 gradient_m[1], gradient_p[1],
356 gradient_m[2], gradient_p[2]);
357}
358
359
360#ifdef DWA_OPENVDB
361inline simd::Float4 simdMin(const simd::Float4& a, const simd::Float4& b) {
362 return simd::Float4(_mm_min_ps(a.base(), b.base()));
363}
364inline simd::Float4 simdMax(const simd::Float4& a, const simd::Float4& b) {
365 return simd::Float4(_mm_max_ps(a.base(), b.base()));
366}
367
368inline float simdSum(const simd::Float4& v);
369
370inline simd::Float4 Pow2(const simd::Float4& v) { return v * v; }
371
372template<>
373inline simd::Float4
374WENO5<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
392inline float
393simdSum(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
402inline float
403GodunovsNormSqrd(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
414template<DScheme DiffScheme>
415struct 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
438template<>
439struct 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
492template<>
493struct 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
552template<>
553struct 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
621template<>
622struct 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
701template<>
702struct 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
752template<>
753struct 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
820template<>
821struct 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
892template<>
893struct 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
944template<>
945struct 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
1008template<>
1009struct 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
1078template<>
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
1177template<>
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
1275template<>
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
1381template<>
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
1485template<DScheme DiffScheme>
1486struct 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
1531template<>
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
1583template<>
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
1637template<>
1638struct 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 {
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 {
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 {
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 {
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 {
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 {
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
1697template<>
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 {
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 {
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 {
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 {
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 {
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 {
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
1763template<DDScheme DiffScheme>
1764struct 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
1804template<>
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
1918template<>
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
2089template<>
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
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Signed (x, y, z) 32-bit integer coordinates.
Definition Coord.h:25
Coord offsetBy(Int32 dx, Int32 dy, Int32 dz) const
Definition Coord.h:91
Definition Vec3.h:24
TemporalIntegrationScheme
Temporal integration schemes.
Definition FiniteDifference.h:233
@ TVD_RK1
Definition FiniteDifference.h:235
@ UNKNOWN_TIS
Definition FiniteDifference.h:234
@ TVD_RK2
Definition FiniteDifference.h:236
@ TVD_RK3
Definition FiniteDifference.h:237
ValueType WENO5(const ValueType &v1, const ValueType &v2, const ValueType &v3, const ValueType &v4, const ValueType &v5, float scale2=0.01f)
Implementation of nominally fifth-order finite-difference WENO.
Definition FiniteDifference.h:303
std::string dsSchemeToMenuName(DScheme dss)
Definition FiniteDifference.h:119
DScheme stringToDScheme(const std::string &s)
Definition FiniteDifference.h:77
Real GodunovsNormSqrd(bool isOutside, Real dP_xm, Real dP_xp, Real dP_ym, Real dP_yp, Real dP_zm, Real dP_zp)
Definition FiniteDifference.h:325
DScheme
Different discrete schemes used in the first derivatives.
Definition FiniteDifference.h:31
@ FD_2ND
Definition FiniteDifference.h:38
@ BD_3RD
Definition FiniteDifference.h:42
@ UNKNOWN_DS
Definition FiniteDifference.h:32
@ CD_6TH
Definition FiniteDifference.h:36
@ CD_2NDT
Definition FiniteDifference.h:33
@ FD_HJWENO5
Definition FiniteDifference.h:45
@ FD_WENO5
Definition FiniteDifference.h:43
@ BD_2ND
Definition FiniteDifference.h:41
@ FD_1ST
Definition FiniteDifference.h:37
@ BD_HJWENO5
Definition FiniteDifference.h:46
@ CD_2ND
Definition FiniteDifference.h:34
@ BD_WENO5
Definition FiniteDifference.h:44
@ BD_1ST
Definition FiniteDifference.h:40
@ FD_3RD
Definition FiniteDifference.h:39
@ CD_4TH
Definition FiniteDifference.h:35
DDScheme
Different discrete schemes used in the second derivatives.
Definition FiniteDifference.h:149
@ UNKNOWN_DD
Definition FiniteDifference.h:150
@ CD_SIXTH
Definition FiniteDifference.h:153
@ CD_FOURTH
Definition FiniteDifference.h:152
@ CD_SECOND
Definition FiniteDifference.h:151
@ NUM_BIAS_SCHEMES
Definition FiniteDifference.h:173
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition Math.h:595
std::string biasedGradientSchemeToMenuName(BiasedGradientScheme bgs)
Definition FiniteDifference.h:214
@ NUM_DS_SCHEMES
Definition FiniteDifference.h:49
std::string dsSchemeToString(DScheme dss)
Definition FiniteDifference.h:53
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition FiniteDifference.h:164
@ FIRST_BIAS
Definition FiniteDifference.h:166
@ THIRD_BIAS
Definition FiniteDifference.h:168
@ WENO5_BIAS
Definition FiniteDifference.h:169
@ UNKNOWN_BIAS
Definition FiniteDifference.h:165
@ SECOND_BIAS
Definition FiniteDifference.h:167
@ HJWENO5_BIAS
Definition FiniteDifference.h:170
@ NUM_DD_SCHEMES
Definition FiniteDifference.h:156
@ NUM_TEMPORAL_SCHEMES
Definition FiniteDifference.h:240
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition Math.h:656
std::string temporalIntegrationSchemeToString(TemporalIntegrationScheme tis)
Definition FiniteDifference.h:243
BiasedGradientScheme stringToBiasedGradientScheme(const std::string &s)
Definition FiniteDifference.h:191
std::string temporalIntegrationSchemeToMenuName(TemporalIntegrationScheme tis)
Definition FiniteDifference.h:276
TemporalIntegrationScheme stringToTemporalIntegrationScheme(const std::string &s)
Definition FiniteDifference.h:256
std::string biasedGradientSchemeToString(BiasedGradientScheme bgs)
Definition FiniteDifference.h:176
Type Pow2(Type x)
Return x2.
Definition Math.h:548
double Real
Definition Types.h:60
Definition Exceptions.h:13
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
Definition FiniteDifference.h:1576
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
Definition FiniteDifference.h:1562
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1554
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1538
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1546
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
Definition FiniteDifference.h:1569
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
Definition FiniteDifference.h:1629
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
Definition FiniteDifference.h:1615
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1606
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1590
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1598
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
Definition FiniteDifference.h:1622
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
Definition FiniteDifference.h:1688
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
Definition FiniteDifference.h:1672
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1663
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1645
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1654
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
Definition FiniteDifference.h:1680
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
Definition FiniteDifference.h:1754
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
Definition FiniteDifference.h:1736
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1725
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1705
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1715
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
Definition FiniteDifference.h:1745
Definition FiniteDifference.h:1487
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
Definition FiniteDifference.h:1524
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
Definition FiniteDifference.h:1512
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1504
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1491
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
Definition FiniteDifference.h:1498
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
Definition FiniteDifference.h:1518
static ValueType difference(const ValueType &xm1, const ValueType &xm0)
Definition FiniteDifference.h:898
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:905
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:917
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:937
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:911
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:931
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:925
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:958
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:974
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:999
static ValueType difference(const ValueType &xm2, const ValueType &xm1, const ValueType &xm0)
Definition FiniteDifference.h:950
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:966
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:991
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:983
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1022
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1040
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:1068
static ValueType difference(const ValueType &xm3, const ValueType &xm2, const ValueType &xm1, const ValueType &xm0)
Definition FiniteDifference.h:1014
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1031
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:1059
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:1050
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1393
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1423
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:1469
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1408
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:1454
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:1439
static ValueType difference(const ValueType &xm3, const ValueType &xm2, const ValueType &xm1, const ValueType &xm0, const ValueType &xp1, const ValueType &xp2)
Definition FiniteDifference.h:1385
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1289
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1319
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:1365
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1304
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:1350
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:1335
static ValueType difference(const ValueType &xm3, const ValueType &xm2, const ValueType &xm1, const ValueType &xm0, const ValueType &xp1, const ValueType &xp2)
Definition FiniteDifference.h:1280
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:449
static ValueType difference(const ValueType &xp1, const ValueType &xm1)
Definition FiniteDifference.h:443
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:465
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:486
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:457
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:480
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:474
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:508
static ValueType difference(const ValueType &xp1, const ValueType &xm1)
Definition FiniteDifference.h:498
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:524
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:545
static bool difference(const bool &xp1, const bool &)
Definition FiniteDifference.h:501
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:516
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:539
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:534
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:566
static ValueType difference(const ValueType &xp2, const ValueType &xp1, const ValueType &xm1, const ValueType &xm2)
Definition FiniteDifference.h:558
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:583
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:612
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:574
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:603
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:594
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:637
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:655
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:688
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xm1, const ValueType &xm2, const ValueType &xm3)
Definition FiniteDifference.h:627
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:646
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:676
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:665
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:714
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:726
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:745
static ValueType difference(const ValueType &xp1, const ValueType &xp0)
Definition FiniteDifference.h:707
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:720
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:739
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:733
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:765
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:783
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:810
static ValueType difference(const ValueType &xp2, const ValueType &xp1, const ValueType &xp0)
Definition FiniteDifference.h:757
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:774
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:802
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:794
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:835
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:853
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xp0)
Definition FiniteDifference.h:826
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:882
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:844
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:873
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:864
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1191
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1222
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:1262
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xp0, const ValueType &xm1, const ValueType &xm2)
Definition FiniteDifference.h:1183
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1207
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:1251
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:1238
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1093
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1123
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:1165
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xp0, const ValueType &xm1, const ValueType &xm2)
Definition FiniteDifference.h:1083
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1108
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:1153
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:1139
Definition FiniteDifference.h:416
static Accessor::ValueType inX(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 inY(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inY(const Stencil &S)
static Stencil::ValueType inX(const Stencil &S)
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1952
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1970
static ValueType crossdifference(const ValueType &xp2yp2, const ValueType &xp2yp1, const ValueType &xp2ym1, const ValueType &xp2ym2, const ValueType &xp1yp2, const ValueType &xp1yp1, const ValueType &xp1ym1, const ValueType &xp1ym2, const ValueType &xm2yp2, const ValueType &xm2yp1, const ValueType &xm2ym1, const ValueType &xm2ym2, const ValueType &xm1yp2, const ValueType &xm1yp1, const ValueType &xm1ym1, const ValueType &xm1ym2)
Definition FiniteDifference.h:1930
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:2037
static Accessor::ValueType inXandZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1993
static Accessor::ValueType inYandZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:2006
static Stencil::ValueType inXandY(const Stencil &S)
Definition FiniteDifference.h:2046
static Accessor::ValueType inXandY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1980
static Stencil::ValueType inXandZ(const Stencil &S)
Definition FiniteDifference.h:2060
static ValueType difference(const ValueType &xp2, const ValueType &xp1, const ValueType &xp0, const ValueType &xm1, const ValueType &xm2)
Definition FiniteDifference.h:1924
static Stencil::ValueType inYandZ(const Stencil &S)
Definition FiniteDifference.h:2074
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1961
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:2029
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:2021
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1824
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1839
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:1888
static Accessor::ValueType inXandZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1856
static ValueType crossdifference(const ValueType &xpyp, const ValueType &xpym, const ValueType &xmyp, const ValueType &xmym)
Definition FiniteDifference.h:1816
static Accessor::ValueType inYandZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1864
static Stencil::ValueType inXandY(const Stencil &S)
Definition FiniteDifference.h:1896
static Accessor::ValueType inXandY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1847
static Stencil::ValueType inXandZ(const Stencil &S)
Definition FiniteDifference.h:1903
static ValueType difference(const ValueType &xp1, const ValueType &xp0, const ValueType &xm1)
Definition FiniteDifference.h:1810
static Stencil::ValueType inYandZ(const Stencil &S)
Definition FiniteDifference.h:1910
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:1831
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:1881
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:1874
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:2143
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:2163
static ValueType crossdifference(const ValueType &xp1yp1, const ValueType &xm1yp1, const ValueType &xp1ym1, const ValueType &xm1ym1, const ValueType &xp2yp1, const ValueType &xm2yp1, const ValueType &xp2ym1, const ValueType &xm2ym1, const ValueType &xp3yp1, const ValueType &xm3yp1, const ValueType &xp3ym1, const ValueType &xm3ym1, const ValueType &xp1yp2, const ValueType &xm1yp2, const ValueType &xp1ym2, const ValueType &xm1ym2, const ValueType &xp2yp2, const ValueType &xm2yp2, const ValueType &xp2ym2, const ValueType &xm2ym2, const ValueType &xp3yp2, const ValueType &xm3yp2, const ValueType &xp3ym2, const ValueType &xm3ym2, const ValueType &xp1yp3, const ValueType &xm1yp3, const ValueType &xp1ym3, const ValueType &xm1ym3, const ValueType &xp2yp3, const ValueType &xm2yp3, const ValueType &xp2ym3, const ValueType &xm2ym3, const ValueType &xp3yp3, const ValueType &xm3yp3, const ValueType &xp3ym3, const ValueType &xm3ym3)
Definition FiniteDifference.h:2103
static Stencil::ValueType inZ(const Stencil &S)
Definition FiniteDifference.h:2243
static Accessor::ValueType inXandZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:2190
static Accessor::ValueType inYandZ(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:2206
static Stencil::ValueType inXandY(const Stencil &S)
Definition FiniteDifference.h:2252
static Accessor::ValueType inXandY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:2174
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xp0, const ValueType &xm1, const ValueType &xm2, const ValueType &xm3)
Definition FiniteDifference.h:2094
static Stencil::ValueType inXandZ(const Stencil &S)
Definition FiniteDifference.h:2275
static Stencil::ValueType inYandZ(const Stencil &S)
Definition FiniteDifference.h:2298
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition FiniteDifference.h:2153
static Stencil::ValueType inY(const Stencil &S)
Definition FiniteDifference.h:2233
static Stencil::ValueType inX(const Stencil &S)
Definition FiniteDifference.h:2224
Definition FiniteDifference.h:1765
static Accessor::ValueType inX(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 inXandZ(const Accessor &grid, const Coord &ijk)
static Accessor::ValueType inYandZ(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inXandY(const Stencil &S)
static Accessor::ValueType inXandY(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inXandZ(const Stencil &S)
static Stencil::ValueType inYandZ(const Stencil &S)
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
static Stencil::ValueType inY(const Stencil &S)
static Stencil::ValueType inX(const Stencil &S)
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212