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