OpenVDB  1.2.0
FiniteDifference.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 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, const ValueType& v4, const ValueType& v5, float scale2 = 0.01)
332 {
333  static const double C=13.0/12.0;
334  // Weno is formulated for non-dimensional equations, here the optional scale2
335  // is a reference value (squared) for the function being interpolated. For
336  // example if 'v' is of order 1000, then scale2 = 10^6 is ok. But in practice
337  // leave scale2 = 1.
338  const double eps=1e-6*scale2;
339  // {\tilde \omega_k} = \gamma_k / ( \beta_k + \epsilon)^2 in Shu's ICASE report)
340  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),
341  A2=0.6/math::Pow2(C*math::Pow2(v2-2*v3+v4)+0.25*math::Pow2(v2-v4)+eps),
342  A3=0.3/math::Pow2(C*math::Pow2(v3-2*v4+v5)+0.25*math::Pow2(3.0*v3-4*v4+v5)+eps);
343 
344  return ValueType(A1*(2.0*v1 - 7.0*v2 + 11.0*v3) +
345  A2*(5.0*v3 - v2 + 2.0*v4) +
346  A3*(2.0*v3 + 5.0*v4 - v5))/(6.0*(A1+A2+A3));
347 }
348 
349 
350 template <typename Real>
351 inline Real GudonovsNormSqrd(bool isOutside,
352  Real dP_xm, Real dP_xp,
353  Real dP_ym, Real dP_yp,
354  Real dP_zm, Real dP_zp)
355 {
356  using math::Max;
357  using math::Min;
358  using math::Pow2;
359 
360  const Real zero(0);
361  Real dPLen2;
362  if (isOutside) { // outside
363  dPLen2 = Max(Pow2(Max(dP_xm, zero)), Pow2(Min(dP_xp,zero))); // (dP/dx)2
364  dPLen2 += Max(Pow2(Max(dP_ym, zero)), Pow2(Min(dP_yp,zero))); // (dP/dy)2
365  dPLen2 += Max(Pow2(Max(dP_zm, zero)), Pow2(Min(dP_zp,zero))); // (dP/dz)2
366  } else { // inside
367  dPLen2 = Max(Pow2(Min(dP_xm, zero)), Pow2(Max(dP_xp,zero))); // (dP/dx)2
368  dPLen2 += Max(Pow2(Min(dP_ym, zero)), Pow2(Max(dP_yp,zero))); // (dP/dy)2
369  dPLen2 += Max(Pow2(Min(dP_zm, zero)), Pow2(Max(dP_zp,zero))); // (dP/dz)2
370  }
371  return dPLen2; // |\nabla\phi|^2
372 };
373 
374 
375 template <typename Real>
376 inline Real GudonovsNormSqrd(bool isOutside, const Vec3<Real>& gradient_m, const Vec3<Real>& gradient_p)
377 {
378  return GudonovsNormSqrd<Real>(isOutside,
379  gradient_m[0], gradient_p[0],
380  gradient_m[1], gradient_p[1],
381  gradient_m[2], gradient_p[2]);
382 }
383 
384 
385 #ifdef DWA_OPENVDB
386 inline simd::Float4 simdMin(const simd::Float4& a, const simd::Float4& b) {
387  return simd::Float4(_mm_min_ps(a.base(), b.base()));
388 }
389 inline simd::Float4 simdMax(const simd::Float4& a, const simd::Float4& b) {
390  return simd::Float4(_mm_max_ps(a.base(), b.base()));
391 }
392 
393 inline float simdSum(const simd::Float4& v);
394 
395 inline simd::Float4 Pow2(const simd::Float4& v) { return v * v; }
396 
397 template<>
398 inline simd::Float4
399 WENO5<simd::Float4>(const simd::Float4& v1, const simd::Float4& v2, const simd::Float4& v3,
400  const simd::Float4& v4, const simd::Float4& v5, float scale2)
401 {
402  using math::Pow2;
403  typedef simd::Float4 F4;
404  const F4
405  C(13.0 / 12.0),
406  eps(1e-6 * scale2),
407  two(2.0), three(3.0), four(4.0), five(5.0), fourth(0.25),
408  A1 = F4(0.1) / Pow2(C*Pow2(v1-two*v2+v3) + fourth*Pow2(v1-four*v2+three*v3) + eps),
409  A2 = F4(0.6) / Pow2(C*Pow2(v2-two*v3+v4) + fourth*Pow2(v2-v4) + eps),
410  A3 = F4(0.3) / Pow2(C*Pow2(v3-two*v4+v5) + fourth*Pow2(three*v3-four*v4+v5) + eps);
411  return (A1 * (two * v1 - F4(7.0) * v2 + F4(11.0) * v3) +
412  A2 * (five * v3 - v2 + two * v4) +
413  A3 * (two * v3 + five * v4 - v5)) / (F4(6.0) * (A1 + A2 + A3));
414 }
415 
416 
417 inline float
418 simdSum(const simd::Float4& v)
419 {
420  // temp = { v3+v3, v2+v2, v1+v3, v0+v2 }
421  __m128 temp = _mm_add_ps(v.base(), _mm_movehl_ps(v.base(), v.base()));
422  // temp = { v3+v3, v2+v2, v1+v3, (v0+v2)+(v1+v3) }
423  temp = _mm_add_ss(temp, _mm_shuffle_ps(temp, temp, 1));
424  return _mm_cvtss_f32(temp);
425 }
426 
427 inline float
428 GudonovsNormSqrd(bool isOutside, const simd::Float4& dP_m, const simd::Float4& dP_p)
429 {
430  const simd::Float4 zero(0.0);
431  simd::Float4 v = isOutside
432  ? simdMax(math::Pow2(simdMax(dP_m, zero)), math::Pow2(simdMin(dP_p, zero)))
433  : simdMax(math::Pow2(simdMin(dP_m, zero)), math::Pow2(simdMax(dP_p, zero)));
434  return simdSum(v);//should be v[0]+v[1]+v[2]
435 }
436 #endif
437 
438 template<DScheme DiffScheme>
439 struct D1
440 {
441  // random access version
442  template<typename Accessor>
443  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk);
444 
445  template<typename Accessor>
446  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk);
447 
448  template<typename Accessor>
449  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk);
450 
451  // stencil access version
452  template<typename Stencil>
453  static typename Stencil::ValueType inX(const Stencil& S);
454 
455  template<typename Stencil>
456  static typename Stencil::ValueType inY(const Stencil& S);
457 
458  template<typename Stencil>
459  static typename Stencil::ValueType inZ(const Stencil& S);
460 };
461 
462 template<>
463 struct D1<CD_2NDT>
464 {
465  // the difference opperator
466  template <typename ValueType>
467  static ValueType difference(const ValueType& xp1, const ValueType& xm1) {
468  return xp1 - xm1;
469  }
470 
471  // random access version
472  template<typename Accessor>
473  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
474  {
475  return difference(grid.getValue(ijk.offsetBy(1, 0, 0)), grid.getValue(ijk.offsetBy(-1, 0, 0)));
476  }
477 
478  template<typename Accessor>
479  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
480  {
481  return difference(grid.getValue(ijk.offsetBy(0, 1, 0)), grid.getValue(ijk.offsetBy( 0, -1, 0)));
482  }
483 
484  template<typename Accessor>
485  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
486  {
487  return difference(grid.getValue(ijk.offsetBy(0, 0, 1)), grid.getValue(ijk.offsetBy( 0, 0, -1)));
488  }
489 
490  // stencil access version
491  template<typename Stencil>
492  static typename Stencil::ValueType inX(const Stencil& S)
493  {
494  return difference( S.template getValue< 1, 0, 0>(), S.template getValue<-1, 0, 0>());
495  }
496 
497  template<typename Stencil>
498  static typename Stencil::ValueType inY(const Stencil& S)
499  {
500  return difference( S.template getValue< 0, 1, 0>(), S.template getValue< 0,-1, 0>());
501  }
502 
503  template<typename Stencil>
504  static typename Stencil::ValueType inZ(const Stencil& S)
505  {
506  return difference( S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0,-1>());
507  }
508 };
509 
510 template<>
511 struct D1<CD_2ND>
512 {
513 
514  // the difference opperator
515  template <typename ValueType>
516  static ValueType difference(const ValueType& xp1, const ValueType& xm1) {
517  return (xp1 - xm1)*ValueType(0.5);
518  }
519 
520 
521  // random access
522  template<typename Accessor>
523  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
524  {
525  return difference(grid.getValue(ijk.offsetBy(1, 0, 0)), grid.getValue(ijk.offsetBy(-1, 0, 0)));
526  }
527 
528  template<typename Accessor>
529  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
530  {
531  return difference(grid.getValue(ijk.offsetBy(0, 1, 0)), grid.getValue(ijk.offsetBy( 0, -1, 0)));
532  }
533 
534  template<typename Accessor>
535  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
536  {
537  return difference(grid.getValue(ijk.offsetBy(0, 0, 1)), grid.getValue(ijk.offsetBy( 0, 0, -1)));
538  }
539 
540 
541  // stencil access version
542  template<typename Stencil>
543  static typename Stencil::ValueType inX(const Stencil& S)
544  {
545  typedef typename Stencil::ValueType ValueType;
546  return difference(S.template getValue< 1, 0, 0>(), S.template getValue<-1, 0, 0>());
547  }
548  template<typename Stencil>
549  static typename Stencil::ValueType inY(const Stencil& S)
550  {
551  typedef typename Stencil::ValueType ValueType;
552  return difference(S.template getValue< 0, 1, 0>(), S.template getValue< 0,-1, 0>());
553  }
554 
555  template<typename Stencil>
556  static typename Stencil::ValueType inZ(const Stencil& S)
557  {
558  typedef typename Stencil::ValueType ValueType;
559  return difference(S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0,-1>());
560  }
561 
562 };
563 
564 template<>
565 struct D1<CD_4TH>
566 {
567 
568  // the difference opperator
569  template <typename ValueType>
570  static ValueType difference( const ValueType& xp2, const ValueType& xp1,
571  const ValueType& xm1, const ValueType& xm2 ) {
572  return ValueType(2./3.)*(xp1 - xm1) + ValueType(1./12.)*(xm2 - xp2) ;
573  }
574 
575 
576  // random access version
577  template<typename Accessor>
578  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
579  {
580  return difference( grid.getValue(ijk.offsetBy( 2,0,0)), grid.getValue(ijk.offsetBy( 1,0,0)),
581  grid.getValue(ijk.offsetBy(-1,0,0)), grid.getValue(ijk.offsetBy(-2,0,0)) );
582 
583  }
584 
585  template<typename Accessor>
586  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
587  {
588 
589  return difference( grid.getValue(ijk.offsetBy( 0, 2, 0)), grid.getValue(ijk.offsetBy( 0, 1, 0)),
590  grid.getValue(ijk.offsetBy( 0,-1, 0)), grid.getValue(ijk.offsetBy( 0,-2, 0)) );
591 
592  }
593 
594  template<typename Accessor>
595  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
596  {
597 
598  return difference( grid.getValue(ijk.offsetBy( 0, 0, 2)), grid.getValue(ijk.offsetBy( 0, 0, 1)),
599  grid.getValue(ijk.offsetBy( 0, 0,-1)), grid.getValue(ijk.offsetBy( 0, 0,-2)) );
600  }
601 
602 
603  // stencil access version
604  template<typename Stencil>
605  static typename Stencil::ValueType inX(const Stencil& S)
606  {
607  return difference( S.template getValue< 2, 0, 0>(),
608  S.template getValue< 1, 0, 0>(),
609  S.template getValue<-1, 0, 0>(),
610  S.template getValue<-2, 0, 0>() );
611  }
612 
613  template<typename Stencil>
614  static typename Stencil::ValueType inY(const Stencil& S)
615  {
616  return difference( S.template getValue< 0, 2, 0>(),
617  S.template getValue< 0, 1, 0>(),
618  S.template getValue< 0,-1, 0>(),
619  S.template getValue< 0,-2, 0>() );
620  }
621 
622  template<typename Stencil>
623  static typename Stencil::ValueType inZ(const Stencil& S)
624  {
625  return difference( S.template getValue< 0, 0, 2>(),
626  S.template getValue< 0, 0, 1>(),
627  S.template getValue< 0, 0,-1>(),
628  S.template getValue< 0, 0,-2>() );
629  }
630 };
631 
632 template<>
633 struct D1<CD_6TH>
634 {
635 
636  // the difference opperator
637  template <typename ValueType>
638  static ValueType difference( const ValueType& xp3, const ValueType& xp2, const ValueType& xp1,
639  const ValueType& xm1, const ValueType& xm2, const ValueType& xm3 ) {
640  return ValueType(3./4.)*(xp1 - xm1) - ValueType(0.15)*(xp2 - xm2) + ValueType(1./60.)*(xp3-xm3);
641  }
642 
643 
644  // random access version
645  template<typename Accessor>
646  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
647  {
648  return difference( grid.getValue(ijk.offsetBy( 3,0,0)), grid.getValue(ijk.offsetBy( 2,0,0)),
649  grid.getValue(ijk.offsetBy( 1,0,0)), grid.getValue(ijk.offsetBy(-1,0,0)),
650  grid.getValue(ijk.offsetBy(-2,0,0)), grid.getValue(ijk.offsetBy(-3,0,0)));
651  }
652 
653  template<typename Accessor>
654  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
655  {
656  return difference( grid.getValue(ijk.offsetBy( 0, 3,0)), grid.getValue(ijk.offsetBy( 0, 2,0)),
657  grid.getValue(ijk.offsetBy( 0, 1,0)), grid.getValue(ijk.offsetBy( 0,-1,0)),
658  grid.getValue(ijk.offsetBy( 0,-2,0)), grid.getValue(ijk.offsetBy( 0,-3,0)));
659  }
660 
661  template<typename Accessor>
662  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
663  {
664  return difference( grid.getValue(ijk.offsetBy( 0, 0, 3)), grid.getValue(ijk.offsetBy( 0, 0, 2)),
665  grid.getValue(ijk.offsetBy( 0, 0, 1)), grid.getValue(ijk.offsetBy( 0, 0,-1)),
666  grid.getValue(ijk.offsetBy( 0, 0,-2)), grid.getValue(ijk.offsetBy( 0, 0,-3)));
667  }
668 
669  // stencil access version
670  template<typename Stencil>
671  static typename Stencil::ValueType inX(const Stencil& S)
672  {
673  return difference(S.template getValue< 3, 0, 0>(),
674  S.template getValue< 2, 0, 0>(),
675  S.template getValue< 1, 0, 0>(),
676  S.template getValue<-1, 0, 0>(),
677  S.template getValue<-2, 0, 0>(),
678  S.template getValue<-3, 0, 0>());
679  }
680 
681  template<typename Stencil>
682  static typename Stencil::ValueType inY(const Stencil& S)
683  {
684 
685  return difference( S.template getValue< 0, 3, 0>(),
686  S.template getValue< 0, 2, 0>(),
687  S.template getValue< 0, 1, 0>(),
688  S.template getValue< 0,-1, 0>(),
689  S.template getValue< 0,-2, 0>(),
690  S.template getValue< 0,-3, 0>());
691  }
692 
693  template<typename Stencil>
694  static typename Stencil::ValueType inZ(const Stencil& S)
695  {
696 
697  return difference( S.template getValue< 0, 0, 3>(),
698  S.template getValue< 0, 0, 2>(),
699  S.template getValue< 0, 0, 1>(),
700  S.template getValue< 0, 0,-1>(),
701  S.template getValue< 0, 0,-2>(),
702  S.template getValue< 0, 0,-3>());
703  }
704 };
705 
706 
707 template<>
708 struct D1<FD_1ST>
709 {
710 
711  // the difference opperator
712  template <typename ValueType>
713  static ValueType difference(const ValueType& xp1, const ValueType& xp0) {
714  return xp1 - xp0;
715  }
716 
717 
718  // random access version
719  template<typename Accessor>
720  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
721  {
722  return difference(grid.getValue(ijk.offsetBy(1, 0, 0)), grid.getValue(ijk));
723  }
724 
725  template<typename Accessor>
726  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
727  {
728  return difference(grid.getValue(ijk.offsetBy(0, 1, 0)), grid.getValue(ijk));
729  }
730 
731  template<typename Accessor>
732  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
733  {
734  return difference(grid.getValue(ijk.offsetBy(0, 0, 1)), grid.getValue(ijk));
735  }
736 
737  // stencil access version
738  template<typename Stencil>
739  static typename Stencil::ValueType inX(const Stencil& S)
740  {
741  return difference(S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>());
742  }
743 
744  template<typename Stencil>
745  static typename Stencil::ValueType inY(const Stencil& S)
746  {
747  return difference(S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>());
748  }
749 
750  template<typename Stencil>
751  static typename Stencil::ValueType inZ(const Stencil& S)
752  {
753  return difference(S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>());
754  }
755 };
756 
757 
758 template<>
759 struct D1<FD_2ND>
760 {
761 
762  // the difference opperator
763  template <typename ValueType>
764  static ValueType difference(const ValueType& xp2, const ValueType& xp1, const ValueType& xp0) {
765  return ValueType(2)*xp1 -(ValueType(0.5)*xp2 + ValueType(3./2.)*xp0);
766  }
767 
768 
769  // random access version
770  template<typename Accessor>
771  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
772  {
773  return difference(grid.getValue(ijk.offsetBy(2,0,0)), grid.getValue(ijk.offsetBy(1,0,0)), grid.getValue(ijk));
774  }
775 
776  template<typename Accessor>
777  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
778  {
779  return difference(grid.getValue(ijk.offsetBy(0,2,0)), grid.getValue(ijk.offsetBy(0,1,0)), grid.getValue(ijk));
780  }
781 
782  template<typename Accessor>
783  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
784  {
785  return difference(grid.getValue(ijk.offsetBy(0,0,2)), grid.getValue(ijk.offsetBy(0,0,1)), grid.getValue(ijk));
786  }
787 
788 
789  // stencil access version
790  template<typename Stencil>
791  static typename Stencil::ValueType inX(const Stencil& S)
792  {
793  return difference( S.template getValue< 2, 0, 0>(),
794  S.template getValue< 1, 0, 0>(),
795  S.template getValue< 0, 0, 0>() );
796  }
797 
798  template<typename Stencil>
799  static typename Stencil::ValueType inY(const Stencil& S)
800  {
801  return difference( S.template getValue< 0, 2, 0>(),
802  S.template getValue< 0, 1, 0>(),
803  S.template getValue< 0, 0, 0>() );
804  }
805 
806  template<typename Stencil>
807  static typename Stencil::ValueType inZ(const Stencil& S)
808  {
809  return difference( S.template getValue< 0, 0, 2>(),
810  S.template getValue< 0, 0, 1>(),
811  S.template getValue< 0, 0, 0>() );
812  }
813 
814 };
815 
816 
817 template<>
818 struct D1<FD_3RD>
819 {
820 
821  // the difference opperator
822  template <typename ValueType>
823  static ValueType difference(const ValueType& xp3, const ValueType& xp2, const ValueType& xp1, const ValueType& xp0) {
824  return ValueType(1./3.)*xp3 - ValueType(1.5)*xp2 + ValueType(3.)*xp1 - ValueType(11./6.)*xp0;
825  }
826 
827 
828  // random access version
829  template<typename Accessor>
830  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
831  {
832  return difference( grid.getValue(ijk.offsetBy(3,0,0)),
833  grid.getValue(ijk.offsetBy(2,0,0)),
834  grid.getValue(ijk.offsetBy(1,0,0)),
835  grid.getValue(ijk) );
836  }
837 
838  template<typename Accessor>
839  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
840  {
841  return difference( grid.getValue(ijk.offsetBy(0,3,0)),
842  grid.getValue(ijk.offsetBy(0,2,0)),
843  grid.getValue(ijk.offsetBy(0,1,0)),
844  grid.getValue(ijk) );
845  }
846 
847  template<typename Accessor>
848  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
849  {
850  return difference( grid.getValue(ijk.offsetBy(0,0,3)),
851  grid.getValue(ijk.offsetBy(0,0,2)),
852  grid.getValue(ijk.offsetBy(0,0,1)),
853  grid.getValue(ijk) );
854  }
855 
856 
857  // stencil access version
858  template<typename Stencil>
859  static typename Stencil::ValueType inX(const Stencil& S)
860  {
861  return difference(S.template getValue< 3, 0, 0>(),
862  S.template getValue< 2, 0, 0>(),
863  S.template getValue< 1, 0, 0>(),
864  S.template getValue< 0, 0, 0>() );
865  }
866 
867  template<typename Stencil>
868  static typename Stencil::ValueType inY(const Stencil& S)
869  {
870  return difference(S.template getValue< 0, 3, 0>(),
871  S.template getValue< 0, 2, 0>(),
872  S.template getValue< 0, 1, 0>(),
873  S.template getValue< 0, 0, 0>() );
874  }
875 
876  template<typename Stencil>
877  static typename Stencil::ValueType inZ(const Stencil& S)
878  {
879  return difference( S.template getValue< 0, 0, 3>(),
880  S.template getValue< 0, 0, 2>(),
881  S.template getValue< 0, 0, 1>(),
882  S.template getValue< 0, 0, 0>() );
883  }
884 };
885 
886 
887 template<>
888 struct D1<BD_1ST>
889 {
890 
891  // the difference opperator
892  template <typename ValueType>
893  static ValueType difference(const ValueType& xm1, const ValueType& xm0) {
894  return -D1<FD_1ST>::difference(xm1, xm0);
895  }
896 
897 
898  // random access version
899  template<typename Accessor>
900  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
901  {
902  return difference(grid.getValue(ijk.offsetBy(-1,0,0)), grid.getValue(ijk));
903  }
904 
905  template<typename Accessor>
906  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
907  {
908  return difference(grid.getValue(ijk.offsetBy(0,-1,0)), grid.getValue(ijk));
909  }
910 
911  template<typename Accessor>
912  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
913  {
914  return difference(grid.getValue(ijk.offsetBy(0, 0,-1)), grid.getValue(ijk));
915  }
916 
917 
918  // stencil access version
919  template<typename Stencil>
920  static typename Stencil::ValueType inX(const Stencil& S)
921  {
922  return difference(S.template getValue<-1, 0, 0>(), S.template getValue< 0, 0, 0>());
923  }
924 
925  template<typename Stencil>
926  static typename Stencil::ValueType inY(const Stencil& S)
927  {
928  return difference(S.template getValue< 0,-1, 0>(), S.template getValue< 0, 0, 0>());
929  }
930 
931  template<typename Stencil>
932  static typename Stencil::ValueType inZ(const Stencil& S)
933  {
934  return difference(S.template getValue< 0, 0,-1>(), S.template getValue< 0, 0, 0>());
935  }
936 };
937 
938 
939 template<>
940 struct D1<BD_2ND>
941 {
942 
943  // the difference opperator
944  template <typename ValueType>
945  static ValueType difference(const ValueType& xm2, const ValueType& xm1, const ValueType& xm0) {
946  return -D1<FD_2ND>::difference(xm2, xm1, xm0);
947  }
948 
949 
950  // random access version
951  template<typename Accessor>
952  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
953  {
954  return difference( grid.getValue(ijk.offsetBy(-2,0,0)),
955  grid.getValue(ijk.offsetBy(-1,0,0)),
956  grid.getValue(ijk) );
957  }
958 
959  template<typename Accessor>
960  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
961  {
962  return difference( grid.getValue(ijk.offsetBy(0,-2,0)),
963  grid.getValue(ijk.offsetBy(0,-1,0)),
964  grid.getValue(ijk) );
965  }
966 
967  template<typename Accessor>
968  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
969  {
970  return difference( grid.getValue(ijk.offsetBy(0,0,-2)),
971  grid.getValue(ijk.offsetBy(0,0,-1)),
972  grid.getValue(ijk) );
973  }
974 
975  // stencil access version
976  template<typename Stencil>
977  static typename Stencil::ValueType inX(const Stencil& S)
978  {
979  return difference( S.template getValue<-2, 0, 0>(),
980  S.template getValue<-1, 0, 0>(),
981  S.template getValue< 0, 0, 0>() );
982  }
983 
984  template<typename Stencil>
985  static typename Stencil::ValueType inY(const Stencil& S)
986  {
987  return difference( S.template getValue< 0,-2, 0>(),
988  S.template getValue< 0,-1, 0>(),
989  S.template getValue< 0, 0, 0>() );
990  }
991 
992  template<typename Stencil>
993  static typename Stencil::ValueType inZ(const Stencil& S)
994  {
995  return difference( S.template getValue< 0, 0,-2>(),
996  S.template getValue< 0, 0,-1>(),
997  S.template getValue< 0, 0, 0>() );
998  }
999 };
1000 
1001 
1002 template<>
1003 struct D1<BD_3RD>
1004 {
1005 
1006  // the difference opperator
1007  template <typename ValueType>
1008  static ValueType difference(const ValueType& xm3, const ValueType& xm2, const ValueType& xm1, const ValueType& xm0){
1009  return -D1<FD_3RD>::difference(xm3, xm2, xm1, xm0);
1010  }
1011 
1012  // random access version
1013  template<typename Accessor>
1014  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1015  {
1016  return difference( grid.getValue(ijk.offsetBy(-3,0,0)),
1017  grid.getValue(ijk.offsetBy(-2,0,0)),
1018  grid.getValue(ijk.offsetBy(-1,0,0)),
1019  grid.getValue(ijk) );
1020  }
1021 
1022  template<typename Accessor>
1023  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1024  {
1025  return difference( grid.getValue(ijk.offsetBy( 0,-3,0)),
1026  grid.getValue(ijk.offsetBy( 0,-2,0)),
1027  grid.getValue(ijk.offsetBy( 0,-1,0)),
1028  grid.getValue(ijk) );
1029  }
1030 
1031  template<typename Accessor>
1032  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1033  {
1034  return difference( grid.getValue(ijk.offsetBy( 0, 0,-3)),
1035  grid.getValue(ijk.offsetBy( 0, 0,-2)),
1036  grid.getValue(ijk.offsetBy( 0, 0,-1)),
1037  grid.getValue(ijk) );
1038  }
1039 
1040  // stencil access version
1041  template<typename Stencil>
1042  static typename Stencil::ValueType inX(const Stencil& S)
1043  {
1044  return difference( S.template getValue<-3, 0, 0>(),
1045  S.template getValue<-2, 0, 0>(),
1046  S.template getValue<-1, 0, 0>(),
1047  S.template getValue< 0, 0, 0>() );
1048  }
1049 
1050  template<typename Stencil>
1051  static typename Stencil::ValueType inY(const Stencil& S)
1052  {
1053  return difference( S.template getValue< 0,-3, 0>(),
1054  S.template getValue< 0,-2, 0>(),
1055  S.template getValue< 0,-1, 0>(),
1056  S.template getValue< 0, 0, 0>() );
1057  }
1058 
1059  template<typename Stencil>
1060  static typename Stencil::ValueType inZ(const Stencil& S)
1061  {
1062  return difference( S.template getValue< 0, 0,-3>(),
1063  S.template getValue< 0, 0,-2>(),
1064  S.template getValue< 0, 0,-1>(),
1065  S.template getValue< 0, 0, 0>() );
1066  }
1067 
1068 };
1069 
1070 template<>
1071 struct D1<FD_WENO5>
1072 {
1073  // the difference opperator
1074  template <typename ValueType>
1075  static ValueType difference(const ValueType& xp3, const ValueType& xp2,
1076  const ValueType& xp1, const ValueType& xp0,
1077  const ValueType& xm1, const ValueType& xm2) {
1078  return WENO5<ValueType>(xp3, xp2, xp1, xp0, xm1)
1079  - WENO5<ValueType>(xp2, xp1, xp0, xm1, xm2);
1080  }
1081 
1082 
1083  // random access version
1084  template<typename Accessor>
1085  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1086  {
1087  typedef typename Accessor::ValueType ValueType;
1088  ValueType V[6];
1089  V[0] = grid.getValue(ijk.offsetBy(3,0,0));
1090  V[1] = grid.getValue(ijk.offsetBy(2,0,0));
1091  V[2] = grid.getValue(ijk.offsetBy(1,0,0));
1092  V[3] = grid.getValue(ijk);
1093  V[4] = grid.getValue(ijk.offsetBy(-1,0,0));
1094  V[5] = grid.getValue(ijk.offsetBy(-2,0,0));
1095 
1096  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1097  }
1098 
1099  template<typename Accessor>
1100  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1101  {
1102  typedef typename Accessor::ValueType ValueType;
1103  ValueType V[6];
1104  V[0] = grid.getValue(ijk.offsetBy(0,3,0));
1105  V[1] = grid.getValue(ijk.offsetBy(0,2,0));
1106  V[2] = grid.getValue(ijk.offsetBy(0,1,0));
1107  V[3] = grid.getValue(ijk);
1108  V[4] = grid.getValue(ijk.offsetBy(0,-1,0));
1109  V[5] = grid.getValue(ijk.offsetBy(0,-2,0));
1110 
1111  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1112  }
1113 
1114  template<typename Accessor>
1115  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1116  {
1117  typedef typename Accessor::ValueType ValueType;
1118  ValueType V[6];
1119  V[0] = grid.getValue(ijk.offsetBy(0,0,3));
1120  V[1] = grid.getValue(ijk.offsetBy(0,0,2));
1121  V[2] = grid.getValue(ijk.offsetBy(0,0,1));
1122  V[3] = grid.getValue(ijk);
1123  V[4] = grid.getValue(ijk.offsetBy(0,0,-1));
1124  V[5] = grid.getValue(ijk.offsetBy(0,0,-2));
1125 
1126  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1127  }
1128 
1129  // stencil access version
1130  template<typename Stencil>
1131  static typename Stencil::ValueType inX(const Stencil& S)
1132  {
1133 
1134  return difference( S.template getValue< 3, 0, 0>(),
1135  S.template getValue< 2, 0, 0>(),
1136  S.template getValue< 1, 0, 0>(),
1137  S.template getValue< 0, 0, 0>(),
1138  S.template getValue<-1, 0, 0>(),
1139  S.template getValue<-2, 0, 0>() );
1140 
1141  }
1142 
1143  template<typename Stencil>
1144  static typename Stencil::ValueType inY(const Stencil& S)
1145  {
1146  return difference( S.template getValue< 0, 3, 0>(),
1147  S.template getValue< 0, 2, 0>(),
1148  S.template getValue< 0, 1, 0>(),
1149  S.template getValue< 0, 0, 0>(),
1150  S.template getValue< 0,-1, 0>(),
1151  S.template getValue< 0,-2, 0>() );
1152  }
1153 
1154  template<typename Stencil>
1155  static typename Stencil::ValueType inZ(const Stencil& S)
1156  {
1157 
1158  return difference( S.template getValue< 0, 0, 3>(),
1159  S.template getValue< 0, 0, 2>(),
1160  S.template getValue< 0, 0, 1>(),
1161  S.template getValue< 0, 0, 0>(),
1162  S.template getValue< 0, 0,-1>(),
1163  S.template getValue< 0, 0,-2>() );
1164  }
1165 };
1166 
1167 template<>
1169 {
1170 
1171  // the difference opperator
1172  template <typename ValueType>
1173  static ValueType difference(const ValueType& xp3, const ValueType& xp2,
1174  const ValueType& xp1, const ValueType& xp0,
1175  const ValueType& xm1, const ValueType& xm2) {
1176  return WENO5<ValueType>(xp3 - xp2, xp2 - xp1, xp1 - xp0, xp0-xm1, xm1-xm2);
1177  }
1178 
1179  // random access version
1180  template<typename Accessor>
1181  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1182  {
1183  typedef typename Accessor::ValueType ValueType;
1184  ValueType V[6];
1185  V[0] = grid.getValue(ijk.offsetBy(3,0,0));
1186  V[1] = grid.getValue(ijk.offsetBy(2,0,0));
1187  V[2] = grid.getValue(ijk.offsetBy(1,0,0));
1188  V[3] = grid.getValue(ijk);
1189  V[4] = grid.getValue(ijk.offsetBy(-1,0,0));
1190  V[5] = grid.getValue(ijk.offsetBy(-2,0,0));
1191 
1192  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1193 
1194  }
1195 
1196  template<typename Accessor>
1197  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1198  {
1199  typedef typename Accessor::ValueType ValueType;
1200  ValueType V[6];
1201  V[0] = grid.getValue(ijk.offsetBy(0,3,0));
1202  V[1] = grid.getValue(ijk.offsetBy(0,2,0));
1203  V[2] = grid.getValue(ijk.offsetBy(0,1,0));
1204  V[3] = grid.getValue(ijk);
1205  V[4] = grid.getValue(ijk.offsetBy(0,-1,0));
1206  V[5] = grid.getValue(ijk.offsetBy(0,-2,0));
1207 
1208  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1209  }
1210 
1211  template<typename Accessor>
1212  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1213  {
1214  typedef typename Accessor::ValueType ValueType;
1215  ValueType V[6];
1216  V[0] = grid.getValue(ijk.offsetBy(0,0,3));
1217  V[1] = grid.getValue(ijk.offsetBy(0,0,2));
1218  V[2] = grid.getValue(ijk.offsetBy(0,0,1));
1219  V[3] = grid.getValue(ijk);
1220  V[4] = grid.getValue(ijk.offsetBy(0,0,-1));
1221  V[5] = grid.getValue(ijk.offsetBy(0,0,-2));
1222 
1223  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1224  }
1225 
1226  // stencil access version
1227  template<typename Stencil>
1228  static typename Stencil::ValueType inX(const Stencil& S)
1229  {
1230 
1231  return difference( S.template getValue< 3, 0, 0>(),
1232  S.template getValue< 2, 0, 0>(),
1233  S.template getValue< 1, 0, 0>(),
1234  S.template getValue< 0, 0, 0>(),
1235  S.template getValue<-1, 0, 0>(),
1236  S.template getValue<-2, 0, 0>() );
1237 
1238  }
1239 
1240  template<typename Stencil>
1241  static typename Stencil::ValueType inY(const Stencil& S)
1242  {
1243  return difference( S.template getValue< 0, 3, 0>(),
1244  S.template getValue< 0, 2, 0>(),
1245  S.template getValue< 0, 1, 0>(),
1246  S.template getValue< 0, 0, 0>(),
1247  S.template getValue< 0,-1, 0>(),
1248  S.template getValue< 0,-2, 0>() );
1249  }
1250 
1251  template<typename Stencil>
1252  static typename Stencil::ValueType inZ(const Stencil& S)
1253  {
1254 
1255  return difference( S.template getValue< 0, 0, 3>(),
1256  S.template getValue< 0, 0, 2>(),
1257  S.template getValue< 0, 0, 1>(),
1258  S.template getValue< 0, 0, 0>(),
1259  S.template getValue< 0, 0,-1>(),
1260  S.template getValue< 0, 0,-2>() );
1261  }
1262 
1263 };
1264 
1265 template<>
1266 struct D1<BD_WENO5>
1267 {
1268 
1269  template<typename ValueType>
1270  static ValueType difference(const ValueType& xm3, const ValueType& xm2, const ValueType& xm1,
1271  const ValueType& xm0, const ValueType& xp1, const ValueType& xp2) {
1272  return -D1<FD_WENO5>::difference(xm3, xm2, xm1, xm0, xp1, xp2);
1273  }
1274 
1275 
1276  // random access version
1277  template<typename Accessor>
1278  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1279  {
1280  typedef typename Accessor::ValueType ValueType;
1281  ValueType V[6];
1282  V[0] = grid.getValue(ijk.offsetBy(-3,0,0));
1283  V[1] = grid.getValue(ijk.offsetBy(-2,0,0));
1284  V[2] = grid.getValue(ijk.offsetBy(-1,0,0));
1285  V[3] = grid.getValue(ijk);
1286  V[4] = grid.getValue(ijk.offsetBy(1,0,0));
1287  V[5] = grid.getValue(ijk.offsetBy(2,0,0));
1288 
1289  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1290  }
1291 
1292  template<typename Accessor>
1293  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1294  {
1295  typedef typename Accessor::ValueType ValueType;
1296  ValueType V[6];
1297  V[0] = grid.getValue(ijk.offsetBy(0,-3,0));
1298  V[1] = grid.getValue(ijk.offsetBy(0,-2,0));
1299  V[2] = grid.getValue(ijk.offsetBy(0,-1,0));
1300  V[3] = grid.getValue(ijk);
1301  V[4] = grid.getValue(ijk.offsetBy(0,1,0));
1302  V[5] = grid.getValue(ijk.offsetBy(0,2,0));
1303 
1304  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1305  }
1306 
1307  template<typename Accessor>
1308  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1309  {
1310  typedef typename Accessor::ValueType ValueType;
1311  ValueType V[6];
1312  V[0] = grid.getValue(ijk.offsetBy(0,0,-3));
1313  V[1] = grid.getValue(ijk.offsetBy(0,0,-2));
1314  V[2] = grid.getValue(ijk.offsetBy(0,0,-1));
1315  V[3] = grid.getValue(ijk);
1316  V[4] = grid.getValue(ijk.offsetBy(0,0,1));
1317  V[5] = grid.getValue(ijk.offsetBy(0,0,2));
1318 
1319  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1320  }
1321 
1322  // stencil access version
1323  template<typename Stencil>
1324  static typename Stencil::ValueType inX(const Stencil& S)
1325  {
1326  typedef typename Stencil::ValueType ValueType;
1327  ValueType V[6];
1328  V[0] = S.template getValue<-3, 0, 0>();
1329  V[1] = S.template getValue<-2, 0, 0>();
1330  V[2] = S.template getValue<-1, 0, 0>();
1331  V[3] = S.template getValue< 0, 0, 0>();
1332  V[4] = S.template getValue< 1, 0, 0>();
1333  V[5] = S.template getValue< 2, 0, 0>();
1334 
1335  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1336  }
1337 
1338  template<typename Stencil>
1339  static typename Stencil::ValueType inY(const Stencil& S)
1340  {
1341  typedef typename Stencil::ValueType ValueType;
1342  ValueType V[6];
1343  V[0] = S.template getValue< 0,-3, 0>();
1344  V[1] = S.template getValue< 0,-2, 0>();
1345  V[2] = S.template getValue< 0,-1, 0>();
1346  V[3] = S.template getValue< 0, 0, 0>();
1347  V[4] = S.template getValue< 0, 1, 0>();
1348  V[5] = S.template getValue< 0, 2, 0>();
1349 
1350  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1351  }
1352 
1353  template<typename Stencil>
1354  static typename Stencil::ValueType inZ(const Stencil& S)
1355  {
1356  typedef typename Stencil::ValueType ValueType;
1357  ValueType V[6];
1358  V[0] = S.template getValue< 0, 0,-3>();
1359  V[1] = S.template getValue< 0, 0,-2>();
1360  V[2] = S.template getValue< 0, 0,-1>();
1361  V[3] = S.template getValue< 0, 0, 0>();
1362  V[4] = S.template getValue< 0, 0, 1>();
1363  V[5] = S.template getValue< 0, 0, 2>();
1364 
1365  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1366  }
1367 };
1368 
1369 
1370 template<>
1372 {
1373  template<typename ValueType>
1374  static ValueType difference(const ValueType& xm3, const ValueType& xm2, const ValueType& xm1,
1375  const ValueType& xm0, const ValueType& xp1, const ValueType& xp2) {
1376  return -D1<FD_HJWENO5>::difference(xm3, xm2, xm1, xm0, xp1, xp2);
1377  }
1378 
1379  // random access version
1380  template<typename Accessor>
1381  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1382  {
1383  typedef typename Accessor::ValueType ValueType;
1384  ValueType V[6];
1385  V[0] = grid.getValue(ijk.offsetBy(-3,0,0));
1386  V[1] = grid.getValue(ijk.offsetBy(-2,0,0));
1387  V[2] = grid.getValue(ijk.offsetBy(-1,0,0));
1388  V[3] = grid.getValue(ijk);
1389  V[4] = grid.getValue(ijk.offsetBy(1,0,0));
1390  V[5] = grid.getValue(ijk.offsetBy(2,0,0));
1391 
1392  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1393  }
1394 
1395  template<typename Accessor>
1396  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1397  {
1398  typedef typename Accessor::ValueType ValueType;
1399  ValueType V[6];
1400  V[0] = grid.getValue(ijk.offsetBy(0,-3,0));
1401  V[1] = grid.getValue(ijk.offsetBy(0,-2,0));
1402  V[2] = grid.getValue(ijk.offsetBy(0,-1,0));
1403  V[3] = grid.getValue(ijk);
1404  V[4] = grid.getValue(ijk.offsetBy(0,1,0));
1405  V[5] = grid.getValue(ijk.offsetBy(0,2,0));
1406 
1407  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1408  }
1409 
1410  template<typename Accessor>
1411  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1412  {
1413  typedef typename Accessor::ValueType ValueType;
1414  ValueType V[6];
1415  V[0] = grid.getValue(ijk.offsetBy(0,0,-3));
1416  V[1] = grid.getValue(ijk.offsetBy(0,0,-2));
1417  V[2] = grid.getValue(ijk.offsetBy(0,0,-1));
1418  V[3] = grid.getValue(ijk);
1419  V[4] = grid.getValue(ijk.offsetBy(0,0,1));
1420  V[5] = grid.getValue(ijk.offsetBy(0,0,2));
1421 
1422  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1423  }
1424 
1425  // stencil access version
1426  template<typename Stencil>
1427  static typename Stencil::ValueType inX(const Stencil& S)
1428  {
1429  typedef typename Stencil::ValueType ValueType;
1430  ValueType V[6];
1431  V[0] = S.template getValue<-3, 0, 0>();
1432  V[1] = S.template getValue<-2, 0, 0>();
1433  V[2] = S.template getValue<-1, 0, 0>();
1434  V[3] = S.template getValue< 0, 0, 0>();
1435  V[4] = S.template getValue< 1, 0, 0>();
1436  V[5] = S.template getValue< 2, 0, 0>();
1437 
1438  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1439  }
1440 
1441  template<typename Stencil>
1442  static typename Stencil::ValueType inY(const Stencil& S)
1443  {
1444  typedef typename Stencil::ValueType ValueType;
1445  ValueType V[6];
1446  V[0] = S.template getValue< 0,-3, 0>();
1447  V[1] = S.template getValue< 0,-2, 0>();
1448  V[2] = S.template getValue< 0,-1, 0>();
1449  V[3] = S.template getValue< 0, 0, 0>();
1450  V[4] = S.template getValue< 0, 1, 0>();
1451  V[5] = S.template getValue< 0, 2, 0>();
1452 
1453  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1454  }
1455 
1456  template<typename Stencil>
1457  static typename Stencil::ValueType inZ(const Stencil& S)
1458  {
1459  typedef typename Stencil::ValueType ValueType;
1460  ValueType V[6];
1461  V[0] = S.template getValue< 0, 0,-3>();
1462  V[1] = S.template getValue< 0, 0,-2>();
1463  V[2] = S.template getValue< 0, 0,-1>();
1464  V[3] = S.template getValue< 0, 0, 0>();
1465  V[4] = S.template getValue< 0, 0, 1>();
1466  V[5] = S.template getValue< 0, 0, 2>();
1467 
1468  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1469  }
1470 };
1471 
1472 
1473 template<DScheme DiffScheme>
1474 struct D1Vec
1475 {
1476  // random access version
1477  template<typename Accessor>
1478  static typename Accessor::ValueType::value_type inX(const Accessor& grid, const Coord& ijk, int n)
1479  {
1480  return D1<DiffScheme>::inX(grid, ijk)[n];
1481  }
1482 
1483  template<typename Accessor>
1484  static typename Accessor::ValueType::value_type inY(const Accessor& grid, const Coord& ijk, int n)
1485  {
1486  return D1<DiffScheme>::inY(grid, ijk)[n];
1487  }
1488  template<typename Accessor>
1489  static typename Accessor::ValueType::value_type inZ(const Accessor& grid, const Coord& ijk, int n)
1490  {
1491  return D1<DiffScheme>::inZ(grid, ijk)[n];
1492  }
1493 
1494 
1495  // stencil access version
1496  template<typename Stencil>
1497  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1498  {
1499  return D1<DiffScheme>::inX(S)[n];
1500  }
1501 
1502  template<typename Stencil>
1503  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1504  {
1505  return D1<DiffScheme>::inY(S)[n];
1506  }
1507 
1508  template<typename Stencil>
1509  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1510  {
1511  return D1<DiffScheme>::inZ(S)[n];
1512  }
1513 };
1514 
1515 
1516 template<>
1518 {
1519 
1520  // random access version
1521  template<typename Accessor>
1522  static typename Accessor::ValueType::value_type inX(const Accessor& grid, const Coord& ijk, int n)
1523  {
1524  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy( 1, 0, 0))[n],
1525  grid.getValue(ijk.offsetBy(-1, 0, 0))[n] );
1526  }
1527 
1528  template<typename Accessor>
1529  static typename Accessor::ValueType::value_type inY(const Accessor& grid, const Coord& ijk, int n)
1530  {
1531  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy(0, 1, 0))[n],
1532  grid.getValue(ijk.offsetBy(0,-1, 0))[n] );
1533  }
1534 
1535  template<typename Accessor>
1536  static typename Accessor::ValueType::value_type inZ(const Accessor& grid, const Coord& ijk, int n)
1537  {
1538  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy(0, 0, 1))[n],
1539  grid.getValue(ijk.offsetBy(0, 0,-1))[n] );
1540  }
1541 
1542  // stencil access version
1543  template<typename Stencil>
1544  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1545  {
1546  return D1<CD_2NDT>::difference( S.template getValue< 1, 0, 0>()[n],
1547  S.template getValue<-1, 0, 0>()[n] );
1548  }
1549 
1550  template<typename Stencil>
1551  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1552  {
1553  return D1<CD_2NDT>::difference( S.template getValue< 0, 1, 0>()[n],
1554  S.template getValue< 0,-1, 0>()[n] );
1555  }
1556 
1557  template<typename Stencil>
1558  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1559  {
1560  return D1<CD_2NDT>::difference( S.template getValue< 0, 0, 1>()[n],
1561  S.template getValue< 0, 0,-1>()[n] );
1562  }
1563 };
1564 
1565 template<>
1566 struct D1Vec<CD_2ND>
1567 {
1568 
1569  // random access version
1570  template<typename Accessor>
1571  static typename Accessor::ValueType::value_type inX(const Accessor& grid, const Coord& ijk, int n)
1572  {
1573  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy( 1, 0, 0))[n] ,
1574  grid.getValue(ijk.offsetBy(-1, 0, 0))[n] );
1575  }
1576 
1577  template<typename Accessor>
1578  static typename Accessor::ValueType::value_type inY(const Accessor& grid, const Coord& ijk, int n)
1579  {
1580  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy(0, 1, 0))[n] ,
1581  grid.getValue(ijk.offsetBy(0,-1, 0))[n] );
1582  }
1583 
1584  template<typename Accessor>
1585  static typename Accessor::ValueType::value_type inZ(const Accessor& grid, const Coord& ijk, int n)
1586  {
1587  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy(0, 0, 1))[n] ,
1588  grid.getValue(ijk.offsetBy(0, 0,-1))[n] );
1589  }
1590 
1591 
1592  // stencil access version
1593  template<typename Stencil>
1594  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1595  {
1596  return D1<CD_2ND>::difference( S.template getValue< 1, 0, 0>()[n],
1597  S.template getValue<-1, 0, 0>()[n] );
1598  }
1599 
1600  template<typename Stencil>
1601  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1602  {
1603  return D1<CD_2ND>::difference( S.template getValue< 0, 1, 0>()[n],
1604  S.template getValue< 0,-1, 0>()[n] );
1605  }
1606 
1607  template<typename Stencil>
1608  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1609  {
1610  return D1<CD_2ND>::difference( S.template getValue< 0, 0, 1>()[n],
1611  S.template getValue< 0, 0,-1>()[n] );
1612  }
1613 };
1614 
1615 
1616 template<>
1617 struct D1Vec<CD_4TH> {
1618  // typedef typename Accessor::ValueType::value_type value_type;
1619 
1620 
1621  // random access version
1622  template<typename Accessor>
1623  static typename Accessor::ValueType::value_type inX(const Accessor& grid, const Coord& ijk, int n) {
1624  return D1<CD_4TH>::difference(grid.getValue(ijk.offsetBy(2, 0, 0))[n], grid.getValue(ijk.offsetBy( 1, 0, 0))[n],
1625  grid.getValue(ijk.offsetBy(-1,0, 0))[n], grid.getValue(ijk.offsetBy(-2, 0, 0))[n]);
1626  }
1627 
1628  template<typename Accessor>
1629  static typename Accessor::ValueType::value_type inY(const Accessor& grid, const Coord& ijk, int n) {
1630  return D1<CD_4TH>::difference(grid.getValue(ijk.offsetBy( 0, 2, 0))[n], grid.getValue(ijk.offsetBy( 0, 1, 0))[n],
1631  grid.getValue(ijk.offsetBy( 0,-1, 0))[n], grid.getValue(ijk.offsetBy( 0,-2, 0))[n]);
1632  }
1633 
1634  template<typename Accessor>
1635  static typename Accessor::ValueType::value_type inZ(const Accessor& grid, const Coord& ijk, int n) {
1636  return D1<CD_4TH>::difference(grid.getValue(ijk.offsetBy(0,0, 2))[n], grid.getValue(ijk.offsetBy( 0, 0, 1))[n],
1637  grid.getValue(ijk.offsetBy(0,0,-1))[n], grid.getValue(ijk.offsetBy( 0, 0,-2))[n]);
1638  }
1639 
1640  // stencil access version
1641  template<typename Stencil>
1642  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n) {
1643  return D1<CD_4TH>::difference(S.template getValue< 2, 0, 0>()[n], S.template getValue< 1, 0, 0>()[n],
1644  S.template getValue<-1, 0, 0>()[n], S.template getValue<-2, 0, 0>()[n] );
1645  }
1646 
1647  template<typename Stencil>
1648  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n) {
1649  return D1<CD_4TH>::difference(S.template getValue< 0, 2, 0>()[n], S.template getValue< 0, 1, 0>()[n],
1650  S.template getValue< 0,-1, 0>()[n], S.template getValue< 0,-2, 0>()[n]);
1651  }
1652 
1653  template<typename Stencil>
1654  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n) {
1655  return D1<CD_4TH>::difference(S.template getValue< 0, 0, 2>()[n], S.template getValue< 0, 0, 1>()[n],
1656  S.template getValue< 0, 0,-1>()[n], S.template getValue< 0, 0,-2>()[n]);
1657  }
1658 };
1659 
1660 template<>
1661 struct D1Vec<CD_6TH>
1662 {
1663  //typedef typename Accessor::ValueType::value_type::value_type ValueType;
1664 
1665  // random access version
1666  template<typename Accessor>
1667  static typename Accessor::ValueType::value_type inX(const Accessor& grid, const Coord& ijk, int n)
1668  {
1669  return D1<CD_6TH>::difference( grid.getValue(ijk.offsetBy( 3, 0, 0))[n], grid.getValue(ijk.offsetBy( 2, 0, 0))[n],
1670  grid.getValue(ijk.offsetBy( 1, 0, 0))[n], grid.getValue(ijk.offsetBy(-1, 0, 0))[n],
1671  grid.getValue(ijk.offsetBy(-2, 0, 0))[n], grid.getValue(ijk.offsetBy(-3, 0, 0))[n] );
1672  }
1673 
1674  template<typename Accessor>
1675  static typename Accessor::ValueType::value_type inY(const Accessor& grid, const Coord& ijk, int n)
1676  {
1677 
1678  return D1<CD_6TH>::difference( grid.getValue(ijk.offsetBy( 0, 3, 0))[n], grid.getValue(ijk.offsetBy( 0, 2, 0))[n],
1679  grid.getValue(ijk.offsetBy( 0, 1, 0))[n], grid.getValue(ijk.offsetBy( 0,-1, 0))[n],
1680  grid.getValue(ijk.offsetBy( 0,-2, 0))[n], grid.getValue(ijk.offsetBy( 0,-3, 0))[n] );
1681  }
1682 
1683  template<typename Accessor>
1684  static typename Accessor::ValueType::value_type inZ(const Accessor& grid, const Coord& ijk, int n)
1685  {
1686 
1687  return D1<CD_6TH>::difference( grid.getValue(ijk.offsetBy( 0, 0, 3))[n], grid.getValue(ijk.offsetBy( 0, 0, 2))[n],
1688  grid.getValue(ijk.offsetBy( 0, 0, 1))[n], grid.getValue(ijk.offsetBy( 0, 0,-1))[n],
1689  grid.getValue(ijk.offsetBy( 0, 0,-2))[n], grid.getValue(ijk.offsetBy( 0, 0,-3))[n] );
1690  }
1691 
1692 
1693  // stencil access version
1694  template<typename Stencil>
1695  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1696  {
1697  return D1<CD_6TH>::difference( S.template getValue< 3, 0, 0>()[n], S.template getValue< 2, 0, 0>()[n],
1698  S.template getValue< 1, 0, 0>()[n], S.template getValue<-1, 0, 0>()[n],
1699  S.template getValue<-2, 0, 0>()[n], S.template getValue<-3, 0, 0>()[n] );
1700  }
1701 
1702  template<typename Stencil>
1703  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1704  {
1705  return D1<CD_6TH>::difference( S.template getValue< 0, 3, 0>()[n], S.template getValue< 0, 2, 0>()[n],
1706  S.template getValue< 0, 1, 0>()[n], S.template getValue< 0,-1, 0>()[n],
1707  S.template getValue< 0,-2, 0>()[n], S.template getValue< 0,-3, 0>()[n] );
1708  }
1709 
1710  template<typename Stencil>
1711  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1712  {
1713  return D1<CD_6TH>::difference( S.template getValue< 0, 0, 3>()[n], S.template getValue< 0, 0, 2>()[n],
1714  S.template getValue< 0, 0, 1>()[n], S.template getValue< 0, 0,-1>()[n],
1715  S.template getValue< 0, 0,-2>()[n], S.template getValue< 0, 0,-3>()[n] );
1716  }
1717 };
1718 
1719 template<DDScheme DiffScheme>
1720 struct D2
1721 {
1722 
1723  template<typename Accessor>
1724  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk);
1725  template<typename Accessor>
1726  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk);
1727  template<typename Accessor>
1728  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk);
1729 
1730  // cross derivatives
1731  template<typename Accessor>
1732  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk);
1733 
1734  template<typename Accessor>
1735  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk);
1736 
1737  template<typename Accessor>
1738  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk);
1739 
1740 
1741  // stencil access version
1742  template<typename Stencil>
1743  static typename Stencil::ValueType inX(const Stencil& S);
1744  template<typename Stencil>
1745  static typename Stencil::ValueType inY(const Stencil& S);
1746  template<typename Stencil>
1747  static typename Stencil::ValueType inZ(const Stencil& S);
1748 
1749  // cross derivatives
1750  template<typename Stencil>
1751  static typename Stencil::ValueType inXandY(const Stencil& S);
1752 
1753  template<typename Stencil>
1754  static typename Stencil::ValueType inXandZ(const Stencil& S);
1755 
1756  template<typename Stencil>
1757  static typename Stencil::ValueType inYandZ(const Stencil& S);
1758 };
1759 
1760 template<>
1761 struct D2<CD_SECOND>
1762 {
1763 
1764  // the difference opperator
1765  template <typename ValueType>
1766  static ValueType difference(const ValueType& xp1, const ValueType& xp0, const ValueType& xm1) {
1767  return xp1 + xm1 - ValueType(2)*xp0;
1768  }
1769 
1770  template <typename ValueType>
1771  static ValueType crossdifference(const ValueType& xpyp, const ValueType& xpym,
1772  const ValueType& xmyp, const ValueType& xmym) {
1773  return ValueType(0.25)*(xpyp + xmym - xpym - xmyp);
1774  }
1775 
1776  // random access version
1777  template<typename Accessor>
1778  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1779  {
1780  return difference( grid.getValue(ijk.offsetBy( 1,0,0)), grid.getValue(ijk),
1781  grid.getValue(ijk.offsetBy(-1,0,0)) );
1782  }
1783 
1784  template<typename Accessor>
1785  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1786  {
1787 
1788  return difference( grid.getValue(ijk.offsetBy(0, 1,0)), grid.getValue(ijk),
1789  grid.getValue(ijk.offsetBy(0,-1,0)) );
1790  }
1791 
1792  template<typename Accessor>
1793  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1794  {
1795  return difference( grid.getValue(ijk.offsetBy( 0,0, 1)), grid.getValue(ijk),
1796  grid.getValue(ijk.offsetBy( 0,0,-1)) );
1797  }
1798 
1799  // cross derivatives
1800  template<typename Accessor>
1801  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
1802  {
1803  return crossdifference(grid.getValue(ijk.offsetBy(1, 1,0)), grid.getValue(ijk.offsetBy( 1,-1,0)),
1804  grid.getValue(ijk.offsetBy(-1,1,0)), grid.getValue(ijk.offsetBy(-1,-1,0)));
1805 
1806  }
1807 
1808  template<typename Accessor>
1809  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
1810  {
1811  return crossdifference(grid.getValue(ijk.offsetBy(1,0, 1)), grid.getValue(ijk.offsetBy(1, 0,-1)),
1812  grid.getValue(ijk.offsetBy(-1,0,1)), grid.getValue(ijk.offsetBy(-1,0,-1)) );
1813  }
1814 
1815  template<typename Accessor>
1816  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
1817  {
1818  return crossdifference(grid.getValue(ijk.offsetBy(0, 1,1)), grid.getValue(ijk.offsetBy(0, 1,-1)),
1819  grid.getValue(ijk.offsetBy(0,-1,1)), grid.getValue(ijk.offsetBy(0,-1,-1)) );
1820  }
1821 
1822 
1823  // stencil access version
1824  template<typename Stencil>
1825  static typename Stencil::ValueType inX(const Stencil& S)
1826  {
1827  return difference( S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>(),
1828  S.template getValue<-1, 0, 0>() );
1829  }
1830 
1831  template<typename Stencil>
1832  static typename Stencil::ValueType inY(const Stencil& S)
1833  {
1834  return difference( S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>(),
1835  S.template getValue< 0,-1, 0>() );
1836  }
1837 
1838  template<typename Stencil>
1839  static typename Stencil::ValueType inZ(const Stencil& S)
1840  {
1841  return difference( S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>(),
1842  S.template getValue< 0, 0,-1>() );
1843  }
1844 
1845  // cross derivatives
1846  template<typename Stencil>
1847  static typename Stencil::ValueType inXandY(const Stencil& S)
1848  {
1849  return crossdifference(S.template getValue< 1, 1, 0>(), S.template getValue< 1,-1, 0>(),
1850  S.template getValue<-1, 1, 0>(), S.template getValue<-1,-1, 0>() );
1851  }
1852 
1853  template<typename Stencil>
1854  static typename Stencil::ValueType inXandZ(const Stencil& S)
1855  {
1856  return crossdifference(S.template getValue< 1, 0, 1>(), S.template getValue< 1, 0,-1>(),
1857  S.template getValue<-1, 0, 1>(), S.template getValue<-1, 0,-1>() );
1858  }
1859 
1860  template<typename Stencil>
1861  static typename Stencil::ValueType inYandZ(const Stencil& S)
1862  {
1863  return crossdifference(S.template getValue< 0, 1, 1>(), S.template getValue< 0, 1,-1>(),
1864  S.template getValue< 0,-1, 1>(), S.template getValue< 0,-1,-1>() );
1865  }
1866 };
1867 
1868 
1869 template<>
1870 struct D2<CD_FOURTH>
1871 {
1872 
1873  // the difference opperator
1874  template <typename ValueType>
1875  static ValueType difference(const ValueType& xp2, const ValueType& xp1, const ValueType& xp0,
1876  const ValueType& xm1, const ValueType& xm2) {
1877  return ValueType(-1./12.)*(xp2 + xm2) + ValueType(4./3.)*(xp1 + xm1) -ValueType(2.5)*xp0;
1878  }
1879 
1880  template <typename ValueType>
1881  static ValueType crossdifference(const ValueType& xp2yp2, const ValueType& xp2yp1,
1882  const ValueType& xp2ym1, const ValueType& xp2ym2,
1883  const ValueType& xp1yp2, const ValueType& xp1yp1,
1884  const ValueType& xp1ym1, const ValueType& xp1ym2,
1885  const ValueType& xm2yp2, const ValueType& xm2yp1,
1886  const ValueType& xm2ym1, const ValueType& xm2ym2,
1887  const ValueType& xm1yp2, const ValueType& xm1yp1,
1888  const ValueType& xm1ym1, const ValueType& xm1ym2 ) {
1889  ValueType tmp1 =
1890  ValueType(2./3.0)*(xp1yp1 - xm1yp1 - xp1ym1 + xm1ym1)-
1891  ValueType(1./12.)*(xp2yp1 - xm2yp1 - xp2ym1 + xm2ym1);
1892  ValueType tmp2 =
1893  ValueType(2./3.0)*(xp1yp2 - xm1yp2 - xp1ym2 + xm1ym2)-
1894  ValueType(1./12.)*(xp2yp2 - xm2yp2 - xp2ym2 + xm2ym2);
1895 
1896  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
1897  }
1898 
1899 
1900 
1901  // random access version
1902  template<typename Accessor>
1903  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1904  {
1905  return difference(grid.getValue(ijk.offsetBy(2,0,0)), grid.getValue(ijk.offsetBy( 1,0,0)),
1906  grid.getValue(ijk),
1907  grid.getValue(ijk.offsetBy(-1,0,0)), grid.getValue(ijk.offsetBy(-2, 0, 0)));
1908  }
1909 
1910  template<typename Accessor>
1911  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1912  {
1913  return difference(grid.getValue(ijk.offsetBy(0, 2,0)), grid.getValue(ijk.offsetBy(0, 1,0)),
1914  grid.getValue(ijk),
1915  grid.getValue(ijk.offsetBy(0,-1,0)), grid.getValue(ijk.offsetBy(0,-2, 0)));
1916  }
1917 
1918  template<typename Accessor>
1919  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1920  {
1921  return difference(grid.getValue(ijk.offsetBy(0,0, 2)), grid.getValue(ijk.offsetBy(0, 0,1)),
1922  grid.getValue(ijk),
1923  grid.getValue(ijk.offsetBy(0,0,-1)), grid.getValue(ijk.offsetBy(0,0,-2)));
1924  }
1925 
1926  // cross derivatives
1927  template<typename Accessor>
1928  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
1929  {
1930  typedef typename Accessor::ValueType ValueType;
1931  typename Accessor::ValueType tmp1 =
1932  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 1, 0)) -
1933  D1<CD_4TH>::inX(grid, ijk.offsetBy(0,-1, 0));
1934  typename Accessor::ValueType tmp2 =
1935  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 2, 0)) -
1936  D1<CD_4TH>::inX(grid, ijk.offsetBy(0,-2, 0));
1937  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
1938  }
1939 
1940  template<typename Accessor>
1941  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
1942  {
1943  typedef typename Accessor::ValueType ValueType;
1944  typename Accessor::ValueType tmp1 =
1945  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0, 1)) -
1946  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0,-1));
1947  typename Accessor::ValueType tmp2 =
1948  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0, 2)) -
1949  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0,-2));
1950  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
1951  }
1952 
1953  template<typename Accessor>
1954  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
1955  {
1956  typedef typename Accessor::ValueType ValueType;
1957  typename Accessor::ValueType tmp1 =
1958  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0, 1)) -
1959  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0,-1));
1960  typename Accessor::ValueType tmp2 =
1961  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0, 2)) -
1962  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0,-2));
1963  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
1964  }
1965 
1966 
1967  // stencil access version
1968  template<typename Stencil>
1969  static typename Stencil::ValueType inX(const Stencil& S)
1970  {
1971  return difference(S.template getValue< 2, 0, 0>(), S.template getValue< 1, 0, 0>(),
1972  S.template getValue< 0, 0, 0>(),
1973  S.template getValue<-1, 0, 0>(), S.template getValue<-2, 0, 0>() );
1974  }
1975 
1976  template<typename Stencil>
1977  static typename Stencil::ValueType inY(const Stencil& S)
1978  {
1979  return difference(S.template getValue< 0, 2, 0>(), S.template getValue< 0, 1, 0>(),
1980  S.template getValue< 0, 0, 0>(),
1981  S.template getValue< 0,-1, 0>(), S.template getValue< 0,-2, 0>() );
1982  }
1983 
1984  template<typename Stencil>
1985  static typename Stencil::ValueType inZ(const Stencil& S)
1986  {
1987  return difference(S.template getValue< 0, 0, 2>(), S.template getValue< 0, 0, 1>(),
1988  S.template getValue< 0, 0, 0>(),
1989  S.template getValue< 0, 0,-1>(), S.template getValue< 0, 0,-2>() );
1990  }
1991 
1992  // cross derivatives
1993  template<typename Stencil>
1994  static typename Stencil::ValueType inXandY(const Stencil& S)
1995  {
1996  return crossdifference( S.template getValue< 2, 2, 0>(), S.template getValue< 2, 1, 0>(),
1997  S.template getValue< 2,-1, 0>(), S.template getValue< 2,-2, 0>(),
1998  S.template getValue< 1, 2, 0>(), S.template getValue< 1, 1, 0>(),
1999  S.template getValue< 1,-1, 0>(), S.template getValue< 1,-2, 0>(),
2000  S.template getValue<-2, 2, 0>(), S.template getValue<-2, 1, 0>(),
2001  S.template getValue<-2,-1, 0>(), S.template getValue<-2,-2, 0>(),
2002  S.template getValue<-1, 2, 0>(), S.template getValue<-1, 1, 0>(),
2003  S.template getValue<-1,-1, 0>(), S.template getValue<-1,-2, 0>() );
2004  }
2005 
2006  template<typename Stencil>
2007  static typename Stencil::ValueType inXandZ(const Stencil& S)
2008  {
2009  return crossdifference( S.template getValue< 2, 0, 2>(), S.template getValue< 2, 0, 1>(),
2010  S.template getValue< 2, 0,-1>(), S.template getValue< 2, 0,-2>(),
2011  S.template getValue< 1, 0, 2>(), S.template getValue< 1, 0, 1>(),
2012  S.template getValue< 1, 0,-1>(), S.template getValue< 1, 0,-2>(),
2013  S.template getValue<-2, 0, 2>(), S.template getValue<-2, 0, 1>(),
2014  S.template getValue<-2, 0,-1>(), S.template getValue<-2, 0,-2>(),
2015  S.template getValue<-1, 0, 2>(), S.template getValue<-1, 0, 1>(),
2016  S.template getValue<-1, 0,-1>(), S.template getValue<-1, 0,-2>() );
2017  }
2018 
2019  template<typename Stencil>
2020  static typename Stencil::ValueType inYandZ(const Stencil& S)
2021  {
2022  return crossdifference( S.template getValue< 0, 2, 2>(), S.template getValue< 0, 2, 1>(),
2023  S.template getValue< 0, 2,-1>(), S.template getValue< 0, 2,-2>(),
2024  S.template getValue< 0, 1, 2>(), S.template getValue< 0, 1, 1>(),
2025  S.template getValue< 0, 1,-1>(), S.template getValue< 0, 1,-2>(),
2026  S.template getValue< 0,-2, 2>(), S.template getValue< 0,-2, 1>(),
2027  S.template getValue< 0,-2,-1>(), S.template getValue< 0,-2,-2>(),
2028  S.template getValue< 0,-1, 2>(), S.template getValue< 0,-1, 1>(),
2029  S.template getValue< 0,-1,-1>(), S.template getValue< 0,-1,-2>() );
2030  }
2031 
2032 
2033 
2034 };
2035 
2036 template<>
2037 struct D2<CD_SIXTH>
2038 {
2039 
2040  // the difference opperator
2041  template <typename ValueType>
2042  static ValueType difference(const ValueType& xp3, const ValueType& xp2, const ValueType& xp1,
2043  const ValueType& xp0,
2044  const ValueType& xm1, const ValueType& xm2, const ValueType& xm3) {
2045  return ValueType(1./90.)*(xp3 + xm3) - ValueType(3./20.)*(xp2 + xm2)
2046  + ValueType(1.5)*(xp1 + xm1) - ValueType(49./18.)*xp0;
2047  }
2048 
2049  template <typename ValueType>
2050  static ValueType crossdifference( const ValueType& xp1yp1,const ValueType& xm1yp1,
2051  const ValueType& xp1ym1,const ValueType& xm1ym1,
2052  const ValueType& xp2yp1,const ValueType& xm2yp1,
2053  const ValueType& xp2ym1,const ValueType& xm2ym1,
2054  const ValueType& xp3yp1,const ValueType& xm3yp1,
2055  const ValueType& xp3ym1,const ValueType& xm3ym1,
2056  const ValueType& xp1yp2,const ValueType& xm1yp2,
2057  const ValueType& xp1ym2,const ValueType& xm1ym2,
2058  const ValueType& xp2yp2,const ValueType& xm2yp2,
2059  const ValueType& xp2ym2,const ValueType& xm2ym2,
2060  const ValueType& xp3yp2,const ValueType& xm3yp2,
2061  const ValueType& xp3ym2,const ValueType& xm3ym2,
2062  const ValueType& xp1yp3,const ValueType& xm1yp3,
2063  const ValueType& xp1ym3,const ValueType& xm1ym3,
2064  const ValueType& xp2yp3,const ValueType& xm2yp3,
2065  const ValueType& xp2ym3,const ValueType& xm2ym3,
2066  const ValueType& xp3yp3,const ValueType& xm3yp3,
2067  const ValueType& xp3ym3,const ValueType& xm3ym3 )
2068  {
2069  ValueType tmp1 =
2070  ValueType(0.7500)*(xp1yp1 - xm1yp1 - xp1ym1 + xm1ym1) -
2071  ValueType(0.1500)*(xp2yp1 - xm2yp1 - xp2ym1 + xm2ym1) +
2072  ValueType(1./60.)*(xp3yp1 - xm3yp1 - xp3ym1 + xm3ym1);
2073 
2074  ValueType tmp2 =
2075  ValueType(0.7500)*(xp1yp2 - xm1yp2 - xp1ym2 + xm1ym2) -
2076  ValueType(0.1500)*(xp2yp2 - xm2yp2 - xp2ym2 + xm2ym2) +
2077  ValueType(1./60.)*(xp3yp2 - xm3yp2 - xp3ym2 + xm3ym2);
2078 
2079  ValueType tmp3 =
2080  ValueType(0.7500)*(xp1yp3 - xm1yp3 - xp1ym3 + xm1ym3) -
2081  ValueType(0.1500)*(xp2yp3 - xm2yp3 - xp2ym3 + xm2ym3) +
2082  ValueType(1./60.)*(xp3yp3 - xm3yp3 - xp3ym3 + xm3ym3);
2083 
2084  return ValueType(0.75)*tmp1 - ValueType(0.15)*tmp2 + ValueType(1./60)*tmp3;
2085  }
2086 
2087  // random access version
2088 
2089  template<typename Accessor>
2090  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
2091  {
2092  return difference(grid.getValue(ijk.offsetBy( 3, 0, 0)), grid.getValue(ijk.offsetBy( 2, 0, 0)),
2093  grid.getValue(ijk.offsetBy( 1, 0, 0)), grid.getValue(ijk),
2094  grid.getValue(ijk.offsetBy(-1, 0, 0)), grid.getValue(ijk.offsetBy(-2, 0, 0)),
2095  grid.getValue(ijk.offsetBy(-3, 0, 0)) );
2096 
2097  }
2098 
2099  template<typename Accessor>
2100  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
2101  {
2102 
2103  return difference(grid.getValue(ijk.offsetBy( 0, 3, 0)), grid.getValue(ijk.offsetBy( 0, 2, 0)),
2104  grid.getValue(ijk.offsetBy( 0, 1, 0)), grid.getValue(ijk),
2105  grid.getValue(ijk.offsetBy( 0,-1, 0)), grid.getValue(ijk.offsetBy( 0,-2, 0)),
2106  grid.getValue(ijk.offsetBy( 0,-3, 0)) );
2107 
2108  }
2109 
2110  template<typename Accessor>
2111  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
2112  {
2113 
2114  return difference(grid.getValue(ijk.offsetBy( 0, 0, 3)), grid.getValue(ijk.offsetBy( 0, 0, 2)),
2115  grid.getValue(ijk.offsetBy( 0, 0, 1)), grid.getValue(ijk),
2116  grid.getValue(ijk.offsetBy( 0, 0,-1)), grid.getValue(ijk.offsetBy( 0, 0,-2)),
2117  grid.getValue(ijk.offsetBy( 0, 0,-3)) );
2118 
2119  }
2120 
2121  template<typename Accessor>
2122  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
2123  {
2124  typename Accessor::ValueType tmp1 =
2125  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 1, 0)) -
2126  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-1, 0));
2127  typename Accessor::ValueType tmp2 =
2128  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 2, 0)) -
2129  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-2, 0));
2130  typename Accessor::ValueType tmp3 =
2131  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 3, 0)) -
2132  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-3, 0));
2133  return 0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3;
2134  }
2135 
2136  template<typename Accessor>
2137  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
2138  {
2139  typename Accessor::ValueType tmp1 =
2140  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 1)) -
2141  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-1));
2142  typename Accessor::ValueType tmp2 =
2143  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 2)) -
2144  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-2));
2145  typename Accessor::ValueType tmp3 =
2146  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 3)) -
2147  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-3));
2148  return 0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3;
2149  }
2150 
2151  template<typename Accessor>
2152  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
2153  {
2154  typename Accessor::ValueType tmp1 =
2155  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 1)) -
2156  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-1));
2157  typename Accessor::ValueType tmp2 =
2158  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 2)) -
2159  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-2));
2160  typename Accessor::ValueType tmp3 =
2161  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 3)) -
2162  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-3));
2163  return 0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3;
2164  }
2165 
2166 
2167  // stencil access version
2168  template<typename Stencil>
2169  static typename Stencil::ValueType inX(const Stencil& S)
2170  {
2171  return difference( S.template getValue< 3, 0, 0>(), S.template getValue< 2, 0, 0>(),
2172  S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>(),
2173  S.template getValue<-1, 0, 0>(), S.template getValue<-2, 0, 0>(),
2174  S.template getValue<-3, 0, 0>() );
2175  }
2176 
2177  template<typename Stencil>
2178  static typename Stencil::ValueType inY(const Stencil& S)
2179  {
2180  return difference( S.template getValue< 0, 3, 0>(), S.template getValue< 0, 2, 0>(),
2181  S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>(),
2182  S.template getValue< 0,-1, 0>(), S.template getValue< 0,-2, 0>(),
2183  S.template getValue< 0,-3, 0>() );
2184 
2185  }
2186 
2187  template<typename Stencil>
2188  static typename Stencil::ValueType inZ(const Stencil& S)
2189  {
2190  return difference( S.template getValue< 0, 0, 3>(), S.template getValue< 0, 0, 2>(),
2191  S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>(),
2192  S.template getValue< 0, 0,-1>(), S.template getValue< 0, 0,-2>(),
2193  S.template getValue< 0, 0,-3>() );
2194  }
2195 
2196  template<typename Stencil>
2197  static typename Stencil::ValueType inXandY(const Stencil& S)
2198  {
2199  return crossdifference( S.template getValue< 1, 1, 0>(), S.template getValue<-1, 1, 0>(),
2200  S.template getValue< 1,-1, 0>(), S.template getValue<-1,-1, 0>(),
2201  S.template getValue< 2, 1, 0>(), S.template getValue<-2, 1, 0>(),
2202  S.template getValue< 2,-1, 0>(), S.template getValue<-2,-1, 0>(),
2203  S.template getValue< 3, 1, 0>(), S.template getValue<-3, 1, 0>(),
2204  S.template getValue< 3,-1, 0>(), S.template getValue<-3,-1, 0>(),
2205  S.template getValue< 1, 2, 0>(), S.template getValue<-1, 2, 0>(),
2206  S.template getValue< 1,-2, 0>(), S.template getValue<-1,-2, 0>(),
2207  S.template getValue< 2, 2, 0>(), S.template getValue<-2, 2, 0>(),
2208  S.template getValue< 2,-2, 0>(), S.template getValue<-2,-2, 0>(),
2209  S.template getValue< 3, 2, 0>(), S.template getValue<-3, 2, 0>(),
2210  S.template getValue< 3,-2, 0>(), S.template getValue<-3,-2, 0>(),
2211  S.template getValue< 1, 3, 0>(), S.template getValue<-1, 3, 0>(),
2212  S.template getValue< 1,-3, 0>(), S.template getValue<-1,-3, 0>(),
2213  S.template getValue< 2, 3, 0>(), S.template getValue<-2, 3, 0>(),
2214  S.template getValue< 2,-3, 0>(), S.template getValue<-2,-3, 0>(),
2215  S.template getValue< 3, 3, 0>(), S.template getValue<-3, 3, 0>(),
2216  S.template getValue< 3,-3, 0>(), S.template getValue<-3,-3, 0>() );
2217 
2218  }
2219 
2220  template<typename Stencil>
2221  static typename Stencil::ValueType inXandZ(const Stencil& S)
2222  {
2223  return crossdifference( S.template getValue< 1, 0, 1>(), S.template getValue<-1, 0, 1>(),
2224  S.template getValue< 1, 0,-1>(), S.template getValue<-1, 0,-1>(),
2225  S.template getValue< 2, 0, 1>(), S.template getValue<-2, 0, 1>(),
2226  S.template getValue< 2, 0,-1>(), S.template getValue<-2, 0,-1>(),
2227  S.template getValue< 3, 0, 1>(), S.template getValue<-3, 0, 1>(),
2228  S.template getValue< 3, 0,-1>(), S.template getValue<-3, 0,-1>(),
2229  S.template getValue< 1, 0, 2>(), S.template getValue<-1, 0, 2>(),
2230  S.template getValue< 1, 0,-2>(), S.template getValue<-1, 0,-2>(),
2231  S.template getValue< 2, 0, 2>(), S.template getValue<-2, 0, 2>(),
2232  S.template getValue< 2, 0,-2>(), S.template getValue<-2, 0,-2>(),
2233  S.template getValue< 3, 0, 2>(), S.template getValue<-3, 0, 2>(),
2234  S.template getValue< 3, 0,-2>(), S.template getValue<-3, 0,-2>(),
2235  S.template getValue< 1, 0, 3>(), S.template getValue<-1, 0, 3>(),
2236  S.template getValue< 1, 0,-3>(), S.template getValue<-1, 0,-3>(),
2237  S.template getValue< 2, 0, 3>(), S.template getValue<-2, 0, 3>(),
2238  S.template getValue< 2, 0,-3>(), S.template getValue<-2, 0,-3>(),
2239  S.template getValue< 3, 0, 3>(), S.template getValue<-3, 0, 3>(),
2240  S.template getValue< 3, 0,-3>(), S.template getValue<-3, 0,-3>() );
2241 
2242  }
2243 
2244  template<typename Stencil>
2245  static typename Stencil::ValueType inYandZ(const Stencil& S)
2246  {
2247  return crossdifference( S.template getValue< 0, 1, 1>(), S.template getValue< 0,-1, 1>(),
2248  S.template getValue< 0, 1,-1>(), S.template getValue< 0,-1,-1>(),
2249  S.template getValue< 0, 2, 1>(), S.template getValue< 0,-2, 1>(),
2250  S.template getValue< 0, 2,-1>(), S.template getValue< 0,-2,-1>(),
2251  S.template getValue< 0, 3, 1>(), S.template getValue< 0,-3, 1>(),
2252  S.template getValue< 0, 3,-1>(), S.template getValue< 0,-3,-1>(),
2253  S.template getValue< 0, 1, 2>(), S.template getValue< 0,-1, 2>(),
2254  S.template getValue< 0, 1,-2>(), S.template getValue< 0,-1,-2>(),
2255  S.template getValue< 0, 2, 2>(), S.template getValue< 0,-2, 2>(),
2256  S.template getValue< 0, 2,-2>(), S.template getValue< 0,-2,-2>(),
2257  S.template getValue< 0, 3, 2>(), S.template getValue< 0,-3, 2>(),
2258  S.template getValue< 0, 3,-2>(), S.template getValue< 0,-3,-2>(),
2259  S.template getValue< 0, 1, 3>(), S.template getValue< 0,-1, 3>(),
2260  S.template getValue< 0, 1,-3>(), S.template getValue< 0,-1,-3>(),
2261  S.template getValue< 0, 2, 3>(), S.template getValue< 0,-2, 3>(),
2262  S.template getValue< 0, 2,-3>(), S.template getValue< 0,-2,-3>(),
2263  S.template getValue< 0, 3, 3>(), S.template getValue< 0,-3, 3>(),
2264  S.template getValue< 0, 3,-3>(), S.template getValue< 0,-3,-3>() );
2265  }
2266 
2267 };
2268 
2269 
2270 
2271 } // end math namespace
2272 } // namespace OPENVDB_VERSION_NAME
2273 } // end openvdb namespace
2274 
2275 #endif // OPENVDB_MATH_FINITEDIFFERENCE_HAS_BEEN_INCLUDED
2276 
2277 // Copyright (c) 2012-2013 DreamWorks Animation LLC
2278 // All rights reserved. This software is distributed under the
2279 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
2280