OpenVDB  2.3.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,
332  const ValueType& v4, const ValueType& v5, float scale2 = 0.01)
333 {
334  static 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 ValueType(A1*(2.0*v1 - 7.0*v2 + 11.0*v3) +
346  A2*(5.0*v3 - v2 + 2.0*v4) +
347  A3*(2.0*v3 + 5.0*v4 - v5))/(6.0*(A1+A2+A3));
348 }
349 
350 
351 template <typename Real>
352 inline Real GudonovsNormSqrd(bool isOutside,
353  Real dP_xm, Real dP_xp,
354  Real dP_ym, Real dP_yp,
355  Real dP_zm, Real dP_zp)
356 {
357  using math::Max;
358  using math::Min;
359  using math::Pow2;
360 
361  const Real zero(0);
362  Real dPLen2;
363  if (isOutside) { // outside
364  dPLen2 = Max(Pow2(Max(dP_xm, zero)), Pow2(Min(dP_xp,zero))); // (dP/dx)2
365  dPLen2 += Max(Pow2(Max(dP_ym, zero)), Pow2(Min(dP_yp,zero))); // (dP/dy)2
366  dPLen2 += Max(Pow2(Max(dP_zm, zero)), Pow2(Min(dP_zp,zero))); // (dP/dz)2
367  } else { // inside
368  dPLen2 = Max(Pow2(Min(dP_xm, zero)), Pow2(Max(dP_xp,zero))); // (dP/dx)2
369  dPLen2 += Max(Pow2(Min(dP_ym, zero)), Pow2(Max(dP_yp,zero))); // (dP/dy)2
370  dPLen2 += Max(Pow2(Min(dP_zm, zero)), Pow2(Max(dP_zp,zero))); // (dP/dz)2
371  }
372  return dPLen2; // |\nabla\phi|^2
373 }
374 
375 
376 template<typename Real>
377 inline Real
378 GudonovsNormSqrd(bool isOutside, const Vec3<Real>& gradient_m, const Vec3<Real>& gradient_p)
379 {
380  return GudonovsNormSqrd<Real>(isOutside,
381  gradient_m[0], gradient_p[0],
382  gradient_m[1], gradient_p[1],
383  gradient_m[2], gradient_p[2]);
384 }
385 
386 
387 #ifdef DWA_OPENVDB
388 inline simd::Float4 simdMin(const simd::Float4& a, const simd::Float4& b) {
389  return simd::Float4(_mm_min_ps(a.base(), b.base()));
390 }
391 inline simd::Float4 simdMax(const simd::Float4& a, const simd::Float4& b) {
392  return simd::Float4(_mm_max_ps(a.base(), b.base()));
393 }
394 
395 inline float simdSum(const simd::Float4& v);
396 
397 inline simd::Float4 Pow2(const simd::Float4& v) { return v * v; }
398 
399 template<>
400 inline simd::Float4
401 WENO5<simd::Float4>(const simd::Float4& v1, const simd::Float4& v2, const simd::Float4& v3,
402  const simd::Float4& v4, const simd::Float4& v5, float scale2)
403 {
404  using math::Pow2;
405  typedef simd::Float4 F4;
406  const F4
407  C(13.0 / 12.0),
408  eps(1e-6 * scale2),
409  two(2.0), three(3.0), four(4.0), five(5.0), fourth(0.25),
410  A1 = F4(0.1) / Pow2(C*Pow2(v1-two*v2+v3) + fourth*Pow2(v1-four*v2+three*v3) + eps),
411  A2 = F4(0.6) / Pow2(C*Pow2(v2-two*v3+v4) + fourth*Pow2(v2-v4) + eps),
412  A3 = F4(0.3) / Pow2(C*Pow2(v3-two*v4+v5) + fourth*Pow2(three*v3-four*v4+v5) + eps);
413  return (A1 * (two * v1 - F4(7.0) * v2 + F4(11.0) * v3) +
414  A2 * (five * v3 - v2 + two * v4) +
415  A3 * (two * v3 + five * v4 - v5)) / (F4(6.0) * (A1 + A2 + A3));
416 }
417 
418 
419 inline float
420 simdSum(const simd::Float4& v)
421 {
422  // temp = { v3+v3, v2+v2, v1+v3, v0+v2 }
423  __m128 temp = _mm_add_ps(v.base(), _mm_movehl_ps(v.base(), v.base()));
424  // temp = { v3+v3, v2+v2, v1+v3, (v0+v2)+(v1+v3) }
425  temp = _mm_add_ss(temp, _mm_shuffle_ps(temp, temp, 1));
426  return _mm_cvtss_f32(temp);
427 }
428 
429 inline float
430 GudonovsNormSqrd(bool isOutside, const simd::Float4& dP_m, const simd::Float4& dP_p)
431 {
432  const simd::Float4 zero(0.0);
433  simd::Float4 v = isOutside
434  ? simdMax(math::Pow2(simdMax(dP_m, zero)), math::Pow2(simdMin(dP_p, zero)))
435  : simdMax(math::Pow2(simdMin(dP_m, zero)), math::Pow2(simdMax(dP_p, zero)));
436  return simdSum(v);//should be v[0]+v[1]+v[2]
437 }
438 #endif
439 
440 template<DScheme DiffScheme>
441 struct D1
442 {
443  // random access version
444  template<typename Accessor>
445  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk);
446 
447  template<typename Accessor>
448  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk);
449 
450  template<typename Accessor>
451  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk);
452 
453  // stencil access version
454  template<typename Stencil>
455  static typename Stencil::ValueType inX(const Stencil& S);
456 
457  template<typename Stencil>
458  static typename Stencil::ValueType inY(const Stencil& S);
459 
460  template<typename Stencil>
461  static typename Stencil::ValueType inZ(const Stencil& S);
462 };
463 
464 template<>
465 struct D1<CD_2NDT>
466 {
467  // the difference opperator
468  template <typename ValueType>
469  static ValueType difference(const ValueType& xp1, const ValueType& xm1) {
470  return xp1 - xm1;
471  }
472 
473  // random access version
474  template<typename Accessor>
475  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
476  {
477  return difference(
478  grid.getValue(ijk.offsetBy(1, 0, 0)),
479  grid.getValue(ijk.offsetBy(-1, 0, 0)));
480  }
481 
482  template<typename Accessor>
483  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
484  {
485  return difference(
486  grid.getValue(ijk.offsetBy(0, 1, 0)),
487  grid.getValue(ijk.offsetBy( 0, -1, 0)));
488  }
489 
490  template<typename Accessor>
491  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
492  {
493  return difference(
494  grid.getValue(ijk.offsetBy(0, 0, 1)),
495  grid.getValue(ijk.offsetBy( 0, 0, -1)));
496  }
497 
498  // stencil access version
499  template<typename Stencil>
500  static typename Stencil::ValueType inX(const Stencil& S)
501  {
502  return difference( S.template getValue< 1, 0, 0>(), S.template getValue<-1, 0, 0>());
503  }
504 
505  template<typename Stencil>
506  static typename Stencil::ValueType inY(const Stencil& S)
507  {
508  return difference( S.template getValue< 0, 1, 0>(), S.template getValue< 0,-1, 0>());
509  }
510 
511  template<typename Stencil>
512  static typename Stencil::ValueType inZ(const Stencil& S)
513  {
514  return difference( S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0,-1>());
515  }
516 };
517 
518 template<>
519 struct D1<CD_2ND>
520 {
521 
522  // the difference opperator
523  template <typename ValueType>
524  static ValueType difference(const ValueType& xp1, const ValueType& xm1) {
525  return (xp1 - xm1)*ValueType(0.5);
526  }
527 
528 
529  // random access
530  template<typename Accessor>
531  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
532  {
533  return difference(
534  grid.getValue(ijk.offsetBy(1, 0, 0)),
535  grid.getValue(ijk.offsetBy(-1, 0, 0)));
536  }
537 
538  template<typename Accessor>
539  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
540  {
541  return difference(
542  grid.getValue(ijk.offsetBy(0, 1, 0)),
543  grid.getValue(ijk.offsetBy( 0, -1, 0)));
544  }
545 
546  template<typename Accessor>
547  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
548  {
549  return difference(
550  grid.getValue(ijk.offsetBy(0, 0, 1)),
551  grid.getValue(ijk.offsetBy( 0, 0, -1)));
552  }
553 
554 
555  // stencil access version
556  template<typename Stencil>
557  static typename Stencil::ValueType inX(const Stencil& S)
558  {
559  return difference(S.template getValue< 1, 0, 0>(), S.template getValue<-1, 0, 0>());
560  }
561  template<typename Stencil>
562  static typename Stencil::ValueType inY(const Stencil& S)
563  {
564  return difference(S.template getValue< 0, 1, 0>(), S.template getValue< 0,-1, 0>());
565  }
566 
567  template<typename Stencil>
568  static typename Stencil::ValueType inZ(const Stencil& S)
569  {
570  return difference(S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0,-1>());
571  }
572 
573 };
574 
575 template<>
576 struct D1<CD_4TH>
577 {
578 
579  // the difference opperator
580  template <typename ValueType>
581  static ValueType difference( const ValueType& xp2, const ValueType& xp1,
582  const ValueType& xm1, const ValueType& xm2 ) {
583  return ValueType(2./3.)*(xp1 - xm1) + ValueType(1./12.)*(xm2 - xp2) ;
584  }
585 
586 
587  // random access version
588  template<typename Accessor>
589  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
590  {
591  return difference(
592  grid.getValue(ijk.offsetBy( 2,0,0)), grid.getValue(ijk.offsetBy( 1,0,0)),
593  grid.getValue(ijk.offsetBy(-1,0,0)), grid.getValue(ijk.offsetBy(-2,0,0)) );
594  }
595 
596  template<typename Accessor>
597  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
598  {
599 
600  return difference(
601  grid.getValue(ijk.offsetBy( 0, 2, 0)), grid.getValue(ijk.offsetBy( 0, 1, 0)),
602  grid.getValue(ijk.offsetBy( 0,-1, 0)), grid.getValue(ijk.offsetBy( 0,-2, 0)) );
603  }
604 
605  template<typename Accessor>
606  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
607  {
608 
609  return difference(
610  grid.getValue(ijk.offsetBy( 0, 0, 2)), grid.getValue(ijk.offsetBy( 0, 0, 1)),
611  grid.getValue(ijk.offsetBy( 0, 0,-1)), grid.getValue(ijk.offsetBy( 0, 0,-2)) );
612  }
613 
614 
615  // stencil access version
616  template<typename Stencil>
617  static typename Stencil::ValueType inX(const Stencil& S)
618  {
619  return difference( S.template getValue< 2, 0, 0>(),
620  S.template getValue< 1, 0, 0>(),
621  S.template getValue<-1, 0, 0>(),
622  S.template getValue<-2, 0, 0>() );
623  }
624 
625  template<typename Stencil>
626  static typename Stencil::ValueType inY(const Stencil& S)
627  {
628  return difference( S.template getValue< 0, 2, 0>(),
629  S.template getValue< 0, 1, 0>(),
630  S.template getValue< 0,-1, 0>(),
631  S.template getValue< 0,-2, 0>() );
632  }
633 
634  template<typename Stencil>
635  static typename Stencil::ValueType inZ(const Stencil& S)
636  {
637  return difference( S.template getValue< 0, 0, 2>(),
638  S.template getValue< 0, 0, 1>(),
639  S.template getValue< 0, 0,-1>(),
640  S.template getValue< 0, 0,-2>() );
641  }
642 };
643 
644 template<>
645 struct D1<CD_6TH>
646 {
647 
648  // the difference opperator
649  template <typename ValueType>
650  static ValueType difference( const ValueType& xp3, const ValueType& xp2, const ValueType& xp1,
651  const ValueType& xm1, const ValueType& xm2, const ValueType& xm3 )
652  {
653  return ValueType(3./4.)*(xp1 - xm1) - ValueType(0.15)*(xp2 - xm2)
654  + ValueType(1./60.)*(xp3-xm3);
655  }
656 
657 
658  // random access version
659  template<typename Accessor>
660  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
661  {
662  return difference(
663  grid.getValue(ijk.offsetBy( 3,0,0)), grid.getValue(ijk.offsetBy( 2,0,0)),
664  grid.getValue(ijk.offsetBy( 1,0,0)), grid.getValue(ijk.offsetBy(-1,0,0)),
665  grid.getValue(ijk.offsetBy(-2,0,0)), grid.getValue(ijk.offsetBy(-3,0,0)));
666  }
667 
668  template<typename Accessor>
669  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
670  {
671  return difference(
672  grid.getValue(ijk.offsetBy( 0, 3, 0)), grid.getValue(ijk.offsetBy( 0, 2, 0)),
673  grid.getValue(ijk.offsetBy( 0, 1, 0)), grid.getValue(ijk.offsetBy( 0,-1, 0)),
674  grid.getValue(ijk.offsetBy( 0,-2, 0)), grid.getValue(ijk.offsetBy( 0,-3, 0)));
675  }
676 
677  template<typename Accessor>
678  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
679  {
680  return difference(
681  grid.getValue(ijk.offsetBy( 0, 0, 3)), grid.getValue(ijk.offsetBy( 0, 0, 2)),
682  grid.getValue(ijk.offsetBy( 0, 0, 1)), grid.getValue(ijk.offsetBy( 0, 0,-1)),
683  grid.getValue(ijk.offsetBy( 0, 0,-2)), grid.getValue(ijk.offsetBy( 0, 0,-3)));
684  }
685 
686  // stencil access version
687  template<typename Stencil>
688  static typename Stencil::ValueType inX(const Stencil& S)
689  {
690  return difference(S.template getValue< 3, 0, 0>(),
691  S.template getValue< 2, 0, 0>(),
692  S.template getValue< 1, 0, 0>(),
693  S.template getValue<-1, 0, 0>(),
694  S.template getValue<-2, 0, 0>(),
695  S.template getValue<-3, 0, 0>());
696  }
697 
698  template<typename Stencil>
699  static typename Stencil::ValueType inY(const Stencil& S)
700  {
701 
702  return difference( S.template getValue< 0, 3, 0>(),
703  S.template getValue< 0, 2, 0>(),
704  S.template getValue< 0, 1, 0>(),
705  S.template getValue< 0,-1, 0>(),
706  S.template getValue< 0,-2, 0>(),
707  S.template getValue< 0,-3, 0>());
708  }
709 
710  template<typename Stencil>
711  static typename Stencil::ValueType inZ(const Stencil& S)
712  {
713 
714  return difference( S.template getValue< 0, 0, 3>(),
715  S.template getValue< 0, 0, 2>(),
716  S.template getValue< 0, 0, 1>(),
717  S.template getValue< 0, 0,-1>(),
718  S.template getValue< 0, 0,-2>(),
719  S.template getValue< 0, 0,-3>());
720  }
721 };
722 
723 
724 template<>
725 struct D1<FD_1ST>
726 {
727 
728  // the difference opperator
729  template <typename ValueType>
730  static ValueType difference(const ValueType& xp1, const ValueType& xp0) {
731  return xp1 - xp0;
732  }
733 
734 
735  // random access version
736  template<typename Accessor>
737  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
738  {
739  return difference(grid.getValue(ijk.offsetBy(1, 0, 0)), grid.getValue(ijk));
740  }
741 
742  template<typename Accessor>
743  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
744  {
745  return difference(grid.getValue(ijk.offsetBy(0, 1, 0)), grid.getValue(ijk));
746  }
747 
748  template<typename Accessor>
749  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
750  {
751  return difference(grid.getValue(ijk.offsetBy(0, 0, 1)), grid.getValue(ijk));
752  }
753 
754  // stencil access version
755  template<typename Stencil>
756  static typename Stencil::ValueType inX(const Stencil& S)
757  {
758  return difference(S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>());
759  }
760 
761  template<typename Stencil>
762  static typename Stencil::ValueType inY(const Stencil& S)
763  {
764  return difference(S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>());
765  }
766 
767  template<typename Stencil>
768  static typename Stencil::ValueType inZ(const Stencil& S)
769  {
770  return difference(S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>());
771  }
772 };
773 
774 
775 template<>
776 struct D1<FD_2ND>
777 {
778  // the difference opperator
779  template <typename ValueType>
780  static ValueType difference(const ValueType& xp2, const ValueType& xp1, const ValueType& xp0)
781  {
782  return ValueType(2)*xp1 -(ValueType(0.5)*xp2 + ValueType(3./2.)*xp0);
783  }
784 
785 
786  // random access version
787  template<typename Accessor>
788  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
789  {
790  return difference(
791  grid.getValue(ijk.offsetBy(2,0,0)),
792  grid.getValue(ijk.offsetBy(1,0,0)),
793  grid.getValue(ijk));
794  }
795 
796  template<typename Accessor>
797  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
798  {
799  return difference(
800  grid.getValue(ijk.offsetBy(0,2,0)),
801  grid.getValue(ijk.offsetBy(0,1,0)),
802  grid.getValue(ijk));
803  }
804 
805  template<typename Accessor>
806  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
807  {
808  return difference(
809  grid.getValue(ijk.offsetBy(0,0,2)),
810  grid.getValue(ijk.offsetBy(0,0,1)),
811  grid.getValue(ijk));
812  }
813 
814 
815  // stencil access version
816  template<typename Stencil>
817  static typename Stencil::ValueType inX(const Stencil& S)
818  {
819  return difference( S.template getValue< 2, 0, 0>(),
820  S.template getValue< 1, 0, 0>(),
821  S.template getValue< 0, 0, 0>() );
822  }
823 
824  template<typename Stencil>
825  static typename Stencil::ValueType inY(const Stencil& S)
826  {
827  return difference( S.template getValue< 0, 2, 0>(),
828  S.template getValue< 0, 1, 0>(),
829  S.template getValue< 0, 0, 0>() );
830  }
831 
832  template<typename Stencil>
833  static typename Stencil::ValueType inZ(const Stencil& S)
834  {
835  return difference( S.template getValue< 0, 0, 2>(),
836  S.template getValue< 0, 0, 1>(),
837  S.template getValue< 0, 0, 0>() );
838  }
839 
840 };
841 
842 
843 template<>
844 struct D1<FD_3RD>
845 {
846 
847  // the difference opperator
848  template <typename ValueType>
849  static ValueType difference(const ValueType& xp3, const ValueType& xp2,
850  const ValueType& xp1, const ValueType& xp0)
851  {
852  return ValueType(1./3.)*xp3 - ValueType(1.5)*xp2
853  + ValueType(3.)*xp1 - ValueType(11./6.)*xp0;
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 opperator
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 difference( S.template getValue< 3, 0, 0>(),
1167  S.template getValue< 2, 0, 0>(),
1168  S.template getValue< 1, 0, 0>(),
1169  S.template getValue< 0, 0, 0>(),
1170  S.template getValue<-1, 0, 0>(),
1171  S.template getValue<-2, 0, 0>() );
1172 
1173  }
1174 
1175  template<typename Stencil>
1176  static typename Stencil::ValueType inY(const Stencil& S)
1177  {
1178  return difference( S.template getValue< 0, 3, 0>(),
1179  S.template getValue< 0, 2, 0>(),
1180  S.template getValue< 0, 1, 0>(),
1181  S.template getValue< 0, 0, 0>(),
1182  S.template getValue< 0,-1, 0>(),
1183  S.template getValue< 0,-2, 0>() );
1184  }
1185 
1186  template<typename Stencil>
1187  static typename Stencil::ValueType inZ(const Stencil& S)
1188  {
1189 
1190  return difference( S.template getValue< 0, 0, 3>(),
1191  S.template getValue< 0, 0, 2>(),
1192  S.template getValue< 0, 0, 1>(),
1193  S.template getValue< 0, 0, 0>(),
1194  S.template getValue< 0, 0,-1>(),
1195  S.template getValue< 0, 0,-2>() );
1196  }
1197 };
1198 
1199 template<>
1201 {
1202 
1203  // the difference opperator
1204  template <typename ValueType>
1205  static ValueType difference(const ValueType& xp3, const ValueType& xp2,
1206  const ValueType& xp1, const ValueType& xp0,
1207  const ValueType& xm1, const ValueType& xm2) {
1208  return WENO5<ValueType>(xp3 - xp2, xp2 - xp1, xp1 - xp0, xp0-xm1, xm1-xm2);
1209  }
1210 
1211  // random access version
1212  template<typename Accessor>
1213  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1214  {
1215  typedef typename Accessor::ValueType ValueType;
1216  ValueType V[6];
1217  V[0] = grid.getValue(ijk.offsetBy(3,0,0));
1218  V[1] = grid.getValue(ijk.offsetBy(2,0,0));
1219  V[2] = grid.getValue(ijk.offsetBy(1,0,0));
1220  V[3] = grid.getValue(ijk);
1221  V[4] = grid.getValue(ijk.offsetBy(-1,0,0));
1222  V[5] = grid.getValue(ijk.offsetBy(-2,0,0));
1223 
1224  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1225 
1226  }
1227 
1228  template<typename Accessor>
1229  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1230  {
1231  typedef typename Accessor::ValueType ValueType;
1232  ValueType V[6];
1233  V[0] = grid.getValue(ijk.offsetBy(0,3,0));
1234  V[1] = grid.getValue(ijk.offsetBy(0,2,0));
1235  V[2] = grid.getValue(ijk.offsetBy(0,1,0));
1236  V[3] = grid.getValue(ijk);
1237  V[4] = grid.getValue(ijk.offsetBy(0,-1,0));
1238  V[5] = grid.getValue(ijk.offsetBy(0,-2,0));
1239 
1240  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1241  }
1242 
1243  template<typename Accessor>
1244  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1245  {
1246  typedef typename Accessor::ValueType ValueType;
1247  ValueType V[6];
1248  V[0] = grid.getValue(ijk.offsetBy(0,0,3));
1249  V[1] = grid.getValue(ijk.offsetBy(0,0,2));
1250  V[2] = grid.getValue(ijk.offsetBy(0,0,1));
1251  V[3] = grid.getValue(ijk);
1252  V[4] = grid.getValue(ijk.offsetBy(0,0,-1));
1253  V[5] = grid.getValue(ijk.offsetBy(0,0,-2));
1254 
1255  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1256  }
1257 
1258  // stencil access version
1259  template<typename Stencil>
1260  static typename Stencil::ValueType inX(const Stencil& S)
1261  {
1262 
1263  return difference( S.template getValue< 3, 0, 0>(),
1264  S.template getValue< 2, 0, 0>(),
1265  S.template getValue< 1, 0, 0>(),
1266  S.template getValue< 0, 0, 0>(),
1267  S.template getValue<-1, 0, 0>(),
1268  S.template getValue<-2, 0, 0>() );
1269 
1270  }
1271 
1272  template<typename Stencil>
1273  static typename Stencil::ValueType inY(const Stencil& S)
1274  {
1275  return difference( S.template getValue< 0, 3, 0>(),
1276  S.template getValue< 0, 2, 0>(),
1277  S.template getValue< 0, 1, 0>(),
1278  S.template getValue< 0, 0, 0>(),
1279  S.template getValue< 0,-1, 0>(),
1280  S.template getValue< 0,-2, 0>() );
1281  }
1282 
1283  template<typename Stencil>
1284  static typename Stencil::ValueType inZ(const Stencil& S)
1285  {
1286 
1287  return difference( S.template getValue< 0, 0, 3>(),
1288  S.template getValue< 0, 0, 2>(),
1289  S.template getValue< 0, 0, 1>(),
1290  S.template getValue< 0, 0, 0>(),
1291  S.template getValue< 0, 0,-1>(),
1292  S.template getValue< 0, 0,-2>() );
1293  }
1294 
1295 };
1296 
1297 template<>
1298 struct D1<BD_WENO5>
1299 {
1300 
1301  template<typename ValueType>
1302  static ValueType difference(const ValueType& xm3, const ValueType& xm2, const ValueType& xm1,
1303  const ValueType& xm0, const ValueType& xp1, const ValueType& xp2)
1304  {
1305  return -D1<FD_WENO5>::difference(xm3, xm2, xm1, xm0, xp1, xp2);
1306  }
1307 
1308 
1309  // random access version
1310  template<typename Accessor>
1311  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1312  {
1313  typedef typename Accessor::ValueType ValueType;
1314  ValueType V[6];
1315  V[0] = grid.getValue(ijk.offsetBy(-3,0,0));
1316  V[1] = grid.getValue(ijk.offsetBy(-2,0,0));
1317  V[2] = grid.getValue(ijk.offsetBy(-1,0,0));
1318  V[3] = grid.getValue(ijk);
1319  V[4] = grid.getValue(ijk.offsetBy(1,0,0));
1320  V[5] = grid.getValue(ijk.offsetBy(2,0,0));
1321 
1322  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1323  }
1324 
1325  template<typename Accessor>
1326  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1327  {
1328  typedef typename Accessor::ValueType ValueType;
1329  ValueType V[6];
1330  V[0] = grid.getValue(ijk.offsetBy(0,-3,0));
1331  V[1] = grid.getValue(ijk.offsetBy(0,-2,0));
1332  V[2] = grid.getValue(ijk.offsetBy(0,-1,0));
1333  V[3] = grid.getValue(ijk);
1334  V[4] = grid.getValue(ijk.offsetBy(0,1,0));
1335  V[5] = grid.getValue(ijk.offsetBy(0,2,0));
1336 
1337  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1338  }
1339 
1340  template<typename Accessor>
1341  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1342  {
1343  typedef typename Accessor::ValueType ValueType;
1344  ValueType V[6];
1345  V[0] = grid.getValue(ijk.offsetBy(0,0,-3));
1346  V[1] = grid.getValue(ijk.offsetBy(0,0,-2));
1347  V[2] = grid.getValue(ijk.offsetBy(0,0,-1));
1348  V[3] = grid.getValue(ijk);
1349  V[4] = grid.getValue(ijk.offsetBy(0,0,1));
1350  V[5] = grid.getValue(ijk.offsetBy(0,0,2));
1351 
1352  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1353  }
1354 
1355  // stencil access version
1356  template<typename Stencil>
1357  static typename Stencil::ValueType inX(const Stencil& S)
1358  {
1359  typedef typename Stencil::ValueType ValueType;
1360  ValueType V[6];
1361  V[0] = S.template getValue<-3, 0, 0>();
1362  V[1] = S.template getValue<-2, 0, 0>();
1363  V[2] = S.template getValue<-1, 0, 0>();
1364  V[3] = S.template getValue< 0, 0, 0>();
1365  V[4] = S.template getValue< 1, 0, 0>();
1366  V[5] = S.template getValue< 2, 0, 0>();
1367 
1368  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1369  }
1370 
1371  template<typename Stencil>
1372  static typename Stencil::ValueType inY(const Stencil& S)
1373  {
1374  typedef typename Stencil::ValueType ValueType;
1375  ValueType V[6];
1376  V[0] = S.template getValue< 0,-3, 0>();
1377  V[1] = S.template getValue< 0,-2, 0>();
1378  V[2] = S.template getValue< 0,-1, 0>();
1379  V[3] = S.template getValue< 0, 0, 0>();
1380  V[4] = S.template getValue< 0, 1, 0>();
1381  V[5] = S.template getValue< 0, 2, 0>();
1382 
1383  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1384  }
1385 
1386  template<typename Stencil>
1387  static typename Stencil::ValueType inZ(const Stencil& S)
1388  {
1389  typedef typename Stencil::ValueType ValueType;
1390  ValueType V[6];
1391  V[0] = S.template getValue< 0, 0,-3>();
1392  V[1] = S.template getValue< 0, 0,-2>();
1393  V[2] = S.template getValue< 0, 0,-1>();
1394  V[3] = S.template getValue< 0, 0, 0>();
1395  V[4] = S.template getValue< 0, 0, 1>();
1396  V[5] = S.template getValue< 0, 0, 2>();
1397 
1398  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1399  }
1400 };
1401 
1402 
1403 template<>
1405 {
1406  template<typename ValueType>
1407  static ValueType difference(const ValueType& xm3, const ValueType& xm2, const ValueType& xm1,
1408  const ValueType& xm0, const ValueType& xp1, const ValueType& xp2)
1409  {
1410  return -D1<FD_HJWENO5>::difference(xm3, xm2, xm1, xm0, xp1, xp2);
1411  }
1412 
1413  // random access version
1414  template<typename Accessor>
1415  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1416  {
1417  typedef typename Accessor::ValueType ValueType;
1418  ValueType V[6];
1419  V[0] = grid.getValue(ijk.offsetBy(-3,0,0));
1420  V[1] = grid.getValue(ijk.offsetBy(-2,0,0));
1421  V[2] = grid.getValue(ijk.offsetBy(-1,0,0));
1422  V[3] = grid.getValue(ijk);
1423  V[4] = grid.getValue(ijk.offsetBy(1,0,0));
1424  V[5] = grid.getValue(ijk.offsetBy(2,0,0));
1425 
1426  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1427  }
1428 
1429  template<typename Accessor>
1430  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1431  {
1432  typedef typename Accessor::ValueType ValueType;
1433  ValueType V[6];
1434  V[0] = grid.getValue(ijk.offsetBy(0,-3,0));
1435  V[1] = grid.getValue(ijk.offsetBy(0,-2,0));
1436  V[2] = grid.getValue(ijk.offsetBy(0,-1,0));
1437  V[3] = grid.getValue(ijk);
1438  V[4] = grid.getValue(ijk.offsetBy(0,1,0));
1439  V[5] = grid.getValue(ijk.offsetBy(0,2,0));
1440 
1441  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1442  }
1443 
1444  template<typename Accessor>
1445  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1446  {
1447  typedef typename Accessor::ValueType ValueType;
1448  ValueType V[6];
1449  V[0] = grid.getValue(ijk.offsetBy(0,0,-3));
1450  V[1] = grid.getValue(ijk.offsetBy(0,0,-2));
1451  V[2] = grid.getValue(ijk.offsetBy(0,0,-1));
1452  V[3] = grid.getValue(ijk);
1453  V[4] = grid.getValue(ijk.offsetBy(0,0,1));
1454  V[5] = grid.getValue(ijk.offsetBy(0,0,2));
1455 
1456  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1457  }
1458 
1459  // stencil access version
1460  template<typename Stencil>
1461  static typename Stencil::ValueType inX(const Stencil& S)
1462  {
1463  typedef typename Stencil::ValueType ValueType;
1464  ValueType V[6];
1465  V[0] = S.template getValue<-3, 0, 0>();
1466  V[1] = S.template getValue<-2, 0, 0>();
1467  V[2] = S.template getValue<-1, 0, 0>();
1468  V[3] = S.template getValue< 0, 0, 0>();
1469  V[4] = S.template getValue< 1, 0, 0>();
1470  V[5] = S.template getValue< 2, 0, 0>();
1471 
1472  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1473  }
1474 
1475  template<typename Stencil>
1476  static typename Stencil::ValueType inY(const Stencil& S)
1477  {
1478  typedef typename Stencil::ValueType ValueType;
1479  ValueType V[6];
1480  V[0] = S.template getValue< 0,-3, 0>();
1481  V[1] = S.template getValue< 0,-2, 0>();
1482  V[2] = S.template getValue< 0,-1, 0>();
1483  V[3] = S.template getValue< 0, 0, 0>();
1484  V[4] = S.template getValue< 0, 1, 0>();
1485  V[5] = S.template getValue< 0, 2, 0>();
1486 
1487  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1488  }
1489 
1490  template<typename Stencil>
1491  static typename Stencil::ValueType inZ(const Stencil& S)
1492  {
1493  typedef typename Stencil::ValueType ValueType;
1494  ValueType V[6];
1495  V[0] = S.template getValue< 0, 0,-3>();
1496  V[1] = S.template getValue< 0, 0,-2>();
1497  V[2] = S.template getValue< 0, 0,-1>();
1498  V[3] = S.template getValue< 0, 0, 0>();
1499  V[4] = S.template getValue< 0, 0, 1>();
1500  V[5] = S.template getValue< 0, 0, 2>();
1501 
1502  return difference(V[0], V[1], V[2], V[3], V[4], V[5]);
1503  }
1504 };
1505 
1506 
1507 template<DScheme DiffScheme>
1508 struct D1Vec
1509 {
1510  // random access version
1511  template<typename Accessor>
1512  static typename Accessor::ValueType::value_type
1513  inX(const Accessor& grid, const Coord& ijk, int n)
1514  {
1515  return D1<DiffScheme>::inX(grid, ijk)[n];
1516  }
1517 
1518  template<typename Accessor>
1519  static typename Accessor::ValueType::value_type
1520  inY(const Accessor& grid, const Coord& ijk, int n)
1521  {
1522  return D1<DiffScheme>::inY(grid, ijk)[n];
1523  }
1524  template<typename Accessor>
1525  static typename Accessor::ValueType::value_type
1526  inZ(const Accessor& grid, const Coord& ijk, int n)
1527  {
1528  return D1<DiffScheme>::inZ(grid, ijk)[n];
1529  }
1530 
1531 
1532  // stencil access version
1533  template<typename Stencil>
1534  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1535  {
1536  return D1<DiffScheme>::inX(S)[n];
1537  }
1538 
1539  template<typename Stencil>
1540  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1541  {
1542  return D1<DiffScheme>::inY(S)[n];
1543  }
1544 
1545  template<typename Stencil>
1546  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1547  {
1548  return D1<DiffScheme>::inZ(S)[n];
1549  }
1550 };
1551 
1552 
1553 template<>
1555 {
1556 
1557  // random access version
1558  template<typename Accessor>
1559  static typename Accessor::ValueType::value_type
1560  inX(const Accessor& grid, const Coord& ijk, int n)
1561  {
1562  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy( 1, 0, 0))[n],
1563  grid.getValue(ijk.offsetBy(-1, 0, 0))[n] );
1564  }
1565 
1566  template<typename Accessor>
1567  static typename Accessor::ValueType::value_type
1568  inY(const Accessor& grid, const Coord& ijk, int n)
1569  {
1570  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy(0, 1, 0))[n],
1571  grid.getValue(ijk.offsetBy(0,-1, 0))[n] );
1572  }
1573 
1574  template<typename Accessor>
1575  static typename Accessor::ValueType::value_type
1576  inZ(const Accessor& grid, const Coord& ijk, int n)
1577  {
1578  return D1<CD_2NDT>::difference( grid.getValue(ijk.offsetBy(0, 0, 1))[n],
1579  grid.getValue(ijk.offsetBy(0, 0,-1))[n] );
1580  }
1581 
1582  // stencil access version
1583  template<typename Stencil>
1584  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1585  {
1586  return D1<CD_2NDT>::difference( S.template getValue< 1, 0, 0>()[n],
1587  S.template getValue<-1, 0, 0>()[n] );
1588  }
1589 
1590  template<typename Stencil>
1591  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1592  {
1593  return D1<CD_2NDT>::difference( S.template getValue< 0, 1, 0>()[n],
1594  S.template getValue< 0,-1, 0>()[n] );
1595  }
1596 
1597  template<typename Stencil>
1598  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1599  {
1600  return D1<CD_2NDT>::difference( S.template getValue< 0, 0, 1>()[n],
1601  S.template getValue< 0, 0,-1>()[n] );
1602  }
1603 };
1604 
1605 template<>
1606 struct D1Vec<CD_2ND>
1607 {
1608 
1609  // random access version
1610  template<typename Accessor>
1611  static typename Accessor::ValueType::value_type
1612  inX(const Accessor& grid, const Coord& ijk, int n)
1613  {
1614  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy( 1, 0, 0))[n] ,
1615  grid.getValue(ijk.offsetBy(-1, 0, 0))[n] );
1616  }
1617 
1618  template<typename Accessor>
1619  static typename Accessor::ValueType::value_type
1620  inY(const Accessor& grid, const Coord& ijk, int n)
1621  {
1622  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy(0, 1, 0))[n] ,
1623  grid.getValue(ijk.offsetBy(0,-1, 0))[n] );
1624  }
1625 
1626  template<typename Accessor>
1627  static typename Accessor::ValueType::value_type
1628  inZ(const Accessor& grid, const Coord& ijk, int n)
1629  {
1630  return D1<CD_2ND>::difference( grid.getValue(ijk.offsetBy(0, 0, 1))[n] ,
1631  grid.getValue(ijk.offsetBy(0, 0,-1))[n] );
1632  }
1633 
1634 
1635  // stencil access version
1636  template<typename Stencil>
1637  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1638  {
1639  return D1<CD_2ND>::difference( S.template getValue< 1, 0, 0>()[n],
1640  S.template getValue<-1, 0, 0>()[n] );
1641  }
1642 
1643  template<typename Stencil>
1644  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1645  {
1646  return D1<CD_2ND>::difference( S.template getValue< 0, 1, 0>()[n],
1647  S.template getValue< 0,-1, 0>()[n] );
1648  }
1649 
1650  template<typename Stencil>
1651  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1652  {
1653  return D1<CD_2ND>::difference( S.template getValue< 0, 0, 1>()[n],
1654  S.template getValue< 0, 0,-1>()[n] );
1655  }
1656 };
1657 
1658 
1659 template<>
1660 struct D1Vec<CD_4TH> {
1661  // typedef typename Accessor::ValueType::value_type value_type;
1662 
1663 
1664  // random access version
1665  template<typename Accessor>
1666  static typename Accessor::ValueType::value_type
1667  inX(const Accessor& grid, const Coord& ijk, int n)
1668  {
1669  return D1<CD_4TH>::difference(
1670  grid.getValue(ijk.offsetBy(2, 0, 0))[n], grid.getValue(ijk.offsetBy( 1, 0, 0))[n],
1671  grid.getValue(ijk.offsetBy(-1,0, 0))[n], grid.getValue(ijk.offsetBy(-2, 0, 0))[n]);
1672  }
1673 
1674  template<typename Accessor>
1675  static typename Accessor::ValueType::value_type
1676  inY(const Accessor& grid, const Coord& ijk, int n)
1677  {
1678  return D1<CD_4TH>::difference(
1679  grid.getValue(ijk.offsetBy( 0, 2, 0))[n], grid.getValue(ijk.offsetBy( 0, 1, 0))[n],
1680  grid.getValue(ijk.offsetBy( 0,-1, 0))[n], grid.getValue(ijk.offsetBy( 0,-2, 0))[n]);
1681  }
1682 
1683  template<typename Accessor>
1684  static typename Accessor::ValueType::value_type
1685  inZ(const Accessor& grid, const Coord& ijk, int n)
1686  {
1687  return D1<CD_4TH>::difference(
1688  grid.getValue(ijk.offsetBy(0,0, 2))[n], grid.getValue(ijk.offsetBy( 0, 0, 1))[n],
1689  grid.getValue(ijk.offsetBy(0,0,-1))[n], grid.getValue(ijk.offsetBy( 0, 0,-2))[n]);
1690  }
1691 
1692  // stencil access version
1693  template<typename Stencil>
1694  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1695  {
1696  return D1<CD_4TH>::difference(
1697  S.template getValue< 2, 0, 0>()[n], S.template getValue< 1, 0, 0>()[n],
1698  S.template getValue<-1, 0, 0>()[n], S.template getValue<-2, 0, 0>()[n] );
1699  }
1700 
1701  template<typename Stencil>
1702  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1703  {
1704  return D1<CD_4TH>::difference(
1705  S.template getValue< 0, 2, 0>()[n], S.template getValue< 0, 1, 0>()[n],
1706  S.template getValue< 0,-1, 0>()[n], S.template getValue< 0,-2, 0>()[n]);
1707  }
1708 
1709  template<typename Stencil>
1710  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1711  {
1712  return D1<CD_4TH>::difference(
1713  S.template getValue< 0, 0, 2>()[n], S.template getValue< 0, 0, 1>()[n],
1714  S.template getValue< 0, 0,-1>()[n], S.template getValue< 0, 0,-2>()[n]);
1715  }
1716 };
1717 
1718 
1719 template<>
1720 struct D1Vec<CD_6TH>
1721 {
1722  //typedef typename Accessor::ValueType::value_type::value_type ValueType;
1723 
1724  // random access version
1725  template<typename Accessor>
1726  static typename Accessor::ValueType::value_type
1727  inX(const Accessor& grid, const Coord& ijk, int n)
1728  {
1729  return D1<CD_6TH>::difference(
1730  grid.getValue(ijk.offsetBy( 3, 0, 0))[n], grid.getValue(ijk.offsetBy( 2, 0, 0))[n],
1731  grid.getValue(ijk.offsetBy( 1, 0, 0))[n], grid.getValue(ijk.offsetBy(-1, 0, 0))[n],
1732  grid.getValue(ijk.offsetBy(-2, 0, 0))[n], grid.getValue(ijk.offsetBy(-3, 0, 0))[n] );
1733  }
1734 
1735  template<typename Accessor>
1736  static typename Accessor::ValueType::value_type
1737  inY(const Accessor& grid, const Coord& ijk, int n)
1738  {
1739  return D1<CD_6TH>::difference(
1740  grid.getValue(ijk.offsetBy( 0, 3, 0))[n], grid.getValue(ijk.offsetBy( 0, 2, 0))[n],
1741  grid.getValue(ijk.offsetBy( 0, 1, 0))[n], grid.getValue(ijk.offsetBy( 0,-1, 0))[n],
1742  grid.getValue(ijk.offsetBy( 0,-2, 0))[n], grid.getValue(ijk.offsetBy( 0,-3, 0))[n] );
1743  }
1744 
1745  template<typename Accessor>
1746  static typename Accessor::ValueType::value_type
1747  inZ(const Accessor& grid, const Coord& ijk, int n)
1748  {
1749  return D1<CD_6TH>::difference(
1750  grid.getValue(ijk.offsetBy( 0, 0, 3))[n], grid.getValue(ijk.offsetBy( 0, 0, 2))[n],
1751  grid.getValue(ijk.offsetBy( 0, 0, 1))[n], grid.getValue(ijk.offsetBy( 0, 0,-1))[n],
1752  grid.getValue(ijk.offsetBy( 0, 0,-2))[n], grid.getValue(ijk.offsetBy( 0, 0,-3))[n] );
1753  }
1754 
1755 
1756  // stencil access version
1757  template<typename Stencil>
1758  static typename Stencil::ValueType::value_type inX(const Stencil& S, int n)
1759  {
1760  return D1<CD_6TH>::difference(
1761  S.template getValue< 3, 0, 0>()[n], S.template getValue< 2, 0, 0>()[n],
1762  S.template getValue< 1, 0, 0>()[n], S.template getValue<-1, 0, 0>()[n],
1763  S.template getValue<-2, 0, 0>()[n], S.template getValue<-3, 0, 0>()[n] );
1764  }
1765 
1766  template<typename Stencil>
1767  static typename Stencil::ValueType::value_type inY(const Stencil& S, int n)
1768  {
1769  return D1<CD_6TH>::difference(
1770  S.template getValue< 0, 3, 0>()[n], S.template getValue< 0, 2, 0>()[n],
1771  S.template getValue< 0, 1, 0>()[n], S.template getValue< 0,-1, 0>()[n],
1772  S.template getValue< 0,-2, 0>()[n], S.template getValue< 0,-3, 0>()[n] );
1773  }
1774 
1775  template<typename Stencil>
1776  static typename Stencil::ValueType::value_type inZ(const Stencil& S, int n)
1777  {
1778  return D1<CD_6TH>::difference(
1779  S.template getValue< 0, 0, 3>()[n], S.template getValue< 0, 0, 2>()[n],
1780  S.template getValue< 0, 0, 1>()[n], S.template getValue< 0, 0,-1>()[n],
1781  S.template getValue< 0, 0,-2>()[n], S.template getValue< 0, 0,-3>()[n] );
1782  }
1783 };
1784 
1785 template<DDScheme DiffScheme>
1786 struct D2
1787 {
1788 
1789  template<typename Accessor>
1790  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk);
1791  template<typename Accessor>
1792  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk);
1793  template<typename Accessor>
1794  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk);
1795 
1796  // cross derivatives
1797  template<typename Accessor>
1798  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk);
1799 
1800  template<typename Accessor>
1801  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk);
1802 
1803  template<typename Accessor>
1804  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk);
1805 
1806 
1807  // stencil access version
1808  template<typename Stencil>
1809  static typename Stencil::ValueType inX(const Stencil& S);
1810  template<typename Stencil>
1811  static typename Stencil::ValueType inY(const Stencil& S);
1812  template<typename Stencil>
1813  static typename Stencil::ValueType inZ(const Stencil& S);
1814 
1815  // cross derivatives
1816  template<typename Stencil>
1817  static typename Stencil::ValueType inXandY(const Stencil& S);
1818 
1819  template<typename Stencil>
1820  static typename Stencil::ValueType inXandZ(const Stencil& S);
1821 
1822  template<typename Stencil>
1823  static typename Stencil::ValueType inYandZ(const Stencil& S);
1824 };
1825 
1826 template<>
1827 struct D2<CD_SECOND>
1828 {
1829 
1830  // the difference opperator
1831  template <typename ValueType>
1832  static ValueType difference(const ValueType& xp1, const ValueType& xp0, const ValueType& xm1)
1833  {
1834  return xp1 + xm1 - ValueType(2)*xp0;
1835  }
1836 
1837  template <typename ValueType>
1838  static ValueType crossdifference(const ValueType& xpyp, const ValueType& xpym,
1839  const ValueType& xmyp, const ValueType& xmym)
1840  {
1841  return ValueType(0.25)*(xpyp + xmym - xpym - xmyp);
1842  }
1843 
1844  // random access version
1845  template<typename Accessor>
1846  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1847  {
1848  return difference( grid.getValue(ijk.offsetBy( 1,0,0)), grid.getValue(ijk),
1849  grid.getValue(ijk.offsetBy(-1,0,0)) );
1850  }
1851 
1852  template<typename Accessor>
1853  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1854  {
1855 
1856  return difference( grid.getValue(ijk.offsetBy(0, 1,0)), grid.getValue(ijk),
1857  grid.getValue(ijk.offsetBy(0,-1,0)) );
1858  }
1859 
1860  template<typename Accessor>
1861  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1862  {
1863  return difference( grid.getValue(ijk.offsetBy( 0,0, 1)), grid.getValue(ijk),
1864  grid.getValue(ijk.offsetBy( 0,0,-1)) );
1865  }
1866 
1867  // cross derivatives
1868  template<typename Accessor>
1869  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
1870  {
1871  return crossdifference(
1872  grid.getValue(ijk.offsetBy(1, 1,0)), grid.getValue(ijk.offsetBy( 1,-1,0)),
1873  grid.getValue(ijk.offsetBy(-1,1,0)), grid.getValue(ijk.offsetBy(-1,-1,0)));
1874 
1875  }
1876 
1877  template<typename Accessor>
1878  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
1879  {
1880  return crossdifference(
1881  grid.getValue(ijk.offsetBy(1,0, 1)), grid.getValue(ijk.offsetBy(1, 0,-1)),
1882  grid.getValue(ijk.offsetBy(-1,0,1)), grid.getValue(ijk.offsetBy(-1,0,-1)) );
1883  }
1884 
1885  template<typename Accessor>
1886  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
1887  {
1888  return crossdifference(
1889  grid.getValue(ijk.offsetBy(0, 1,1)), grid.getValue(ijk.offsetBy(0, 1,-1)),
1890  grid.getValue(ijk.offsetBy(0,-1,1)), grid.getValue(ijk.offsetBy(0,-1,-1)) );
1891  }
1892 
1893 
1894  // stencil access version
1895  template<typename Stencil>
1896  static typename Stencil::ValueType inX(const Stencil& S)
1897  {
1898  return difference( S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>(),
1899  S.template getValue<-1, 0, 0>() );
1900  }
1901 
1902  template<typename Stencil>
1903  static typename Stencil::ValueType inY(const Stencil& S)
1904  {
1905  return difference( S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>(),
1906  S.template getValue< 0,-1, 0>() );
1907  }
1908 
1909  template<typename Stencil>
1910  static typename Stencil::ValueType inZ(const Stencil& S)
1911  {
1912  return difference( S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>(),
1913  S.template getValue< 0, 0,-1>() );
1914  }
1915 
1916  // cross derivatives
1917  template<typename Stencil>
1918  static typename Stencil::ValueType inXandY(const Stencil& S)
1919  {
1920  return crossdifference(S.template getValue< 1, 1, 0>(), S.template getValue< 1,-1, 0>(),
1921  S.template getValue<-1, 1, 0>(), S.template getValue<-1,-1, 0>() );
1922  }
1923 
1924  template<typename Stencil>
1925  static typename Stencil::ValueType inXandZ(const Stencil& S)
1926  {
1927  return crossdifference(S.template getValue< 1, 0, 1>(), S.template getValue< 1, 0,-1>(),
1928  S.template getValue<-1, 0, 1>(), S.template getValue<-1, 0,-1>() );
1929  }
1930 
1931  template<typename Stencil>
1932  static typename Stencil::ValueType inYandZ(const Stencil& S)
1933  {
1934  return crossdifference(S.template getValue< 0, 1, 1>(), S.template getValue< 0, 1,-1>(),
1935  S.template getValue< 0,-1, 1>(), S.template getValue< 0,-1,-1>() );
1936  }
1937 };
1938 
1939 
1940 template<>
1941 struct D2<CD_FOURTH>
1942 {
1943 
1944  // the difference opperator
1945  template <typename ValueType>
1946  static ValueType difference(const ValueType& xp2, const ValueType& xp1, const ValueType& xp0,
1947  const ValueType& xm1, const ValueType& xm2) {
1948  return ValueType(-1./12.)*(xp2 + xm2) + ValueType(4./3.)*(xp1 + xm1) -ValueType(2.5)*xp0;
1949  }
1950 
1951  template <typename ValueType>
1952  static ValueType crossdifference(const ValueType& xp2yp2, const ValueType& xp2yp1,
1953  const ValueType& xp2ym1, const ValueType& xp2ym2,
1954  const ValueType& xp1yp2, const ValueType& xp1yp1,
1955  const ValueType& xp1ym1, const ValueType& xp1ym2,
1956  const ValueType& xm2yp2, const ValueType& xm2yp1,
1957  const ValueType& xm2ym1, const ValueType& xm2ym2,
1958  const ValueType& xm1yp2, const ValueType& xm1yp1,
1959  const ValueType& xm1ym1, const ValueType& xm1ym2 ) {
1960  ValueType tmp1 =
1961  ValueType(2./3.0)*(xp1yp1 - xm1yp1 - xp1ym1 + xm1ym1)-
1962  ValueType(1./12.)*(xp2yp1 - xm2yp1 - xp2ym1 + xm2ym1);
1963  ValueType tmp2 =
1964  ValueType(2./3.0)*(xp1yp2 - xm1yp2 - xp1ym2 + xm1ym2)-
1965  ValueType(1./12.)*(xp2yp2 - xm2yp2 - xp2ym2 + xm2ym2);
1966 
1967  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
1968  }
1969 
1970 
1971 
1972  // random access version
1973  template<typename Accessor>
1974  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
1975  {
1976  return difference(
1977  grid.getValue(ijk.offsetBy(2,0,0)), grid.getValue(ijk.offsetBy( 1,0,0)),
1978  grid.getValue(ijk),
1979  grid.getValue(ijk.offsetBy(-1,0,0)), grid.getValue(ijk.offsetBy(-2, 0, 0)));
1980  }
1981 
1982  template<typename Accessor>
1983  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
1984  {
1985  return difference(
1986  grid.getValue(ijk.offsetBy(0, 2,0)), grid.getValue(ijk.offsetBy(0, 1,0)),
1987  grid.getValue(ijk),
1988  grid.getValue(ijk.offsetBy(0,-1,0)), grid.getValue(ijk.offsetBy(0,-2, 0)));
1989  }
1990 
1991  template<typename Accessor>
1992  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
1993  {
1994  return difference(
1995  grid.getValue(ijk.offsetBy(0,0, 2)), grid.getValue(ijk.offsetBy(0, 0,1)),
1996  grid.getValue(ijk),
1997  grid.getValue(ijk.offsetBy(0,0,-1)), grid.getValue(ijk.offsetBy(0,0,-2)));
1998  }
1999 
2000  // cross derivatives
2001  template<typename Accessor>
2002  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
2003  {
2004  typedef typename Accessor::ValueType ValueType;
2005  typename Accessor::ValueType tmp1 =
2006  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 1, 0)) -
2007  D1<CD_4TH>::inX(grid, ijk.offsetBy(0,-1, 0));
2008  typename Accessor::ValueType tmp2 =
2009  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 2, 0)) -
2010  D1<CD_4TH>::inX(grid, ijk.offsetBy(0,-2, 0));
2011  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
2012  }
2013 
2014  template<typename Accessor>
2015  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
2016  {
2017  typedef typename Accessor::ValueType ValueType;
2018  typename Accessor::ValueType tmp1 =
2019  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0, 1)) -
2020  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0,-1));
2021  typename Accessor::ValueType tmp2 =
2022  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0, 2)) -
2023  D1<CD_4TH>::inX(grid, ijk.offsetBy(0, 0,-2));
2024  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
2025  }
2026 
2027  template<typename Accessor>
2028  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
2029  {
2030  typedef typename Accessor::ValueType ValueType;
2031  typename Accessor::ValueType tmp1 =
2032  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0, 1)) -
2033  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0,-1));
2034  typename Accessor::ValueType tmp2 =
2035  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0, 2)) -
2036  D1<CD_4TH>::inY(grid, ijk.offsetBy(0, 0,-2));
2037  return ValueType(2./3.)*tmp1 - ValueType(1./12.)*tmp2;
2038  }
2039 
2040 
2041  // stencil access version
2042  template<typename Stencil>
2043  static typename Stencil::ValueType inX(const Stencil& S)
2044  {
2045  return difference(S.template getValue< 2, 0, 0>(), S.template getValue< 1, 0, 0>(),
2046  S.template getValue< 0, 0, 0>(),
2047  S.template getValue<-1, 0, 0>(), S.template getValue<-2, 0, 0>() );
2048  }
2049 
2050  template<typename Stencil>
2051  static typename Stencil::ValueType inY(const Stencil& S)
2052  {
2053  return difference(S.template getValue< 0, 2, 0>(), S.template getValue< 0, 1, 0>(),
2054  S.template getValue< 0, 0, 0>(),
2055  S.template getValue< 0,-1, 0>(), S.template getValue< 0,-2, 0>() );
2056  }
2057 
2058  template<typename Stencil>
2059  static typename Stencil::ValueType inZ(const Stencil& S)
2060  {
2061  return difference(S.template getValue< 0, 0, 2>(), S.template getValue< 0, 0, 1>(),
2062  S.template getValue< 0, 0, 0>(),
2063  S.template getValue< 0, 0,-1>(), S.template getValue< 0, 0,-2>() );
2064  }
2065 
2066  // cross derivatives
2067  template<typename Stencil>
2068  static typename Stencil::ValueType inXandY(const Stencil& S)
2069  {
2070  return crossdifference(
2071  S.template getValue< 2, 2, 0>(), S.template getValue< 2, 1, 0>(),
2072  S.template getValue< 2,-1, 0>(), S.template getValue< 2,-2, 0>(),
2073  S.template getValue< 1, 2, 0>(), S.template getValue< 1, 1, 0>(),
2074  S.template getValue< 1,-1, 0>(), S.template getValue< 1,-2, 0>(),
2075  S.template getValue<-2, 2, 0>(), S.template getValue<-2, 1, 0>(),
2076  S.template getValue<-2,-1, 0>(), S.template getValue<-2,-2, 0>(),
2077  S.template getValue<-1, 2, 0>(), S.template getValue<-1, 1, 0>(),
2078  S.template getValue<-1,-1, 0>(), S.template getValue<-1,-2, 0>() );
2079  }
2080 
2081  template<typename Stencil>
2082  static typename Stencil::ValueType inXandZ(const Stencil& S)
2083  {
2084  return crossdifference(
2085  S.template getValue< 2, 0, 2>(), S.template getValue< 2, 0, 1>(),
2086  S.template getValue< 2, 0,-1>(), S.template getValue< 2, 0,-2>(),
2087  S.template getValue< 1, 0, 2>(), S.template getValue< 1, 0, 1>(),
2088  S.template getValue< 1, 0,-1>(), S.template getValue< 1, 0,-2>(),
2089  S.template getValue<-2, 0, 2>(), S.template getValue<-2, 0, 1>(),
2090  S.template getValue<-2, 0,-1>(), S.template getValue<-2, 0,-2>(),
2091  S.template getValue<-1, 0, 2>(), S.template getValue<-1, 0, 1>(),
2092  S.template getValue<-1, 0,-1>(), S.template getValue<-1, 0,-2>() );
2093  }
2094 
2095  template<typename Stencil>
2096  static typename Stencil::ValueType inYandZ(const Stencil& S)
2097  {
2098  return crossdifference(
2099  S.template getValue< 0, 2, 2>(), S.template getValue< 0, 2, 1>(),
2100  S.template getValue< 0, 2,-1>(), S.template getValue< 0, 2,-2>(),
2101  S.template getValue< 0, 1, 2>(), S.template getValue< 0, 1, 1>(),
2102  S.template getValue< 0, 1,-1>(), S.template getValue< 0, 1,-2>(),
2103  S.template getValue< 0,-2, 2>(), S.template getValue< 0,-2, 1>(),
2104  S.template getValue< 0,-2,-1>(), S.template getValue< 0,-2,-2>(),
2105  S.template getValue< 0,-1, 2>(), S.template getValue< 0,-1, 1>(),
2106  S.template getValue< 0,-1,-1>(), S.template getValue< 0,-1,-2>() );
2107  }
2108 };
2109 
2110 
2111 template<>
2112 struct D2<CD_SIXTH>
2113 {
2114  // the difference opperator
2115  template <typename ValueType>
2116  static ValueType difference(const ValueType& xp3, const ValueType& xp2, const ValueType& xp1,
2117  const ValueType& xp0,
2118  const ValueType& xm1, const ValueType& xm2, const ValueType& xm3)
2119  {
2120  return ValueType(1./90.)*(xp3 + xm3) - ValueType(3./20.)*(xp2 + xm2)
2121  + ValueType(1.5)*(xp1 + xm1) - ValueType(49./18.)*xp0;
2122  }
2123 
2124  template <typename ValueType>
2125  static ValueType crossdifference( const ValueType& xp1yp1,const ValueType& xm1yp1,
2126  const ValueType& xp1ym1,const ValueType& xm1ym1,
2127  const ValueType& xp2yp1,const ValueType& xm2yp1,
2128  const ValueType& xp2ym1,const ValueType& xm2ym1,
2129  const ValueType& xp3yp1,const ValueType& xm3yp1,
2130  const ValueType& xp3ym1,const ValueType& xm3ym1,
2131  const ValueType& xp1yp2,const ValueType& xm1yp2,
2132  const ValueType& xp1ym2,const ValueType& xm1ym2,
2133  const ValueType& xp2yp2,const ValueType& xm2yp2,
2134  const ValueType& xp2ym2,const ValueType& xm2ym2,
2135  const ValueType& xp3yp2,const ValueType& xm3yp2,
2136  const ValueType& xp3ym2,const ValueType& xm3ym2,
2137  const ValueType& xp1yp3,const ValueType& xm1yp3,
2138  const ValueType& xp1ym3,const ValueType& xm1ym3,
2139  const ValueType& xp2yp3,const ValueType& xm2yp3,
2140  const ValueType& xp2ym3,const ValueType& xm2ym3,
2141  const ValueType& xp3yp3,const ValueType& xm3yp3,
2142  const ValueType& xp3ym3,const ValueType& xm3ym3 )
2143  {
2144  ValueType tmp1 =
2145  ValueType(0.7500)*(xp1yp1 - xm1yp1 - xp1ym1 + xm1ym1) -
2146  ValueType(0.1500)*(xp2yp1 - xm2yp1 - xp2ym1 + xm2ym1) +
2147  ValueType(1./60.)*(xp3yp1 - xm3yp1 - xp3ym1 + xm3ym1);
2148 
2149  ValueType tmp2 =
2150  ValueType(0.7500)*(xp1yp2 - xm1yp2 - xp1ym2 + xm1ym2) -
2151  ValueType(0.1500)*(xp2yp2 - xm2yp2 - xp2ym2 + xm2ym2) +
2152  ValueType(1./60.)*(xp3yp2 - xm3yp2 - xp3ym2 + xm3ym2);
2153 
2154  ValueType tmp3 =
2155  ValueType(0.7500)*(xp1yp3 - xm1yp3 - xp1ym3 + xm1ym3) -
2156  ValueType(0.1500)*(xp2yp3 - xm2yp3 - xp2ym3 + xm2ym3) +
2157  ValueType(1./60.)*(xp3yp3 - xm3yp3 - xp3ym3 + xm3ym3);
2158 
2159  return ValueType(0.75)*tmp1 - ValueType(0.15)*tmp2 + ValueType(1./60)*tmp3;
2160  }
2161 
2162  // random access version
2163 
2164  template<typename Accessor>
2165  static typename Accessor::ValueType inX(const Accessor& grid, const Coord& ijk)
2166  {
2167  return difference(
2168  grid.getValue(ijk.offsetBy( 3, 0, 0)), grid.getValue(ijk.offsetBy( 2, 0, 0)),
2169  grid.getValue(ijk.offsetBy( 1, 0, 0)), grid.getValue(ijk),
2170  grid.getValue(ijk.offsetBy(-1, 0, 0)), grid.getValue(ijk.offsetBy(-2, 0, 0)),
2171  grid.getValue(ijk.offsetBy(-3, 0, 0)) );
2172  }
2173 
2174  template<typename Accessor>
2175  static typename Accessor::ValueType inY(const Accessor& grid, const Coord& ijk)
2176  {
2177  return difference(
2178  grid.getValue(ijk.offsetBy( 0, 3, 0)), grid.getValue(ijk.offsetBy( 0, 2, 0)),
2179  grid.getValue(ijk.offsetBy( 0, 1, 0)), grid.getValue(ijk),
2180  grid.getValue(ijk.offsetBy( 0,-1, 0)), grid.getValue(ijk.offsetBy( 0,-2, 0)),
2181  grid.getValue(ijk.offsetBy( 0,-3, 0)) );
2182  }
2183 
2184  template<typename Accessor>
2185  static typename Accessor::ValueType inZ(const Accessor& grid, const Coord& ijk)
2186  {
2187 
2188  return difference(
2189  grid.getValue(ijk.offsetBy( 0, 0, 3)), grid.getValue(ijk.offsetBy( 0, 0, 2)),
2190  grid.getValue(ijk.offsetBy( 0, 0, 1)), grid.getValue(ijk),
2191  grid.getValue(ijk.offsetBy( 0, 0,-1)), grid.getValue(ijk.offsetBy( 0, 0,-2)),
2192  grid.getValue(ijk.offsetBy( 0, 0,-3)) );
2193  }
2194 
2195  template<typename Accessor>
2196  static typename Accessor::ValueType inXandY(const Accessor& grid, const Coord& ijk)
2197  {
2198  typename Accessor::ValueType tmp1 =
2199  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 1, 0)) -
2200  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-1, 0));
2201  typename Accessor::ValueType tmp2 =
2202  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 2, 0)) -
2203  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-2, 0));
2204  typename Accessor::ValueType tmp3 =
2205  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 3, 0)) -
2206  D1<CD_6TH>::inX(grid, ijk.offsetBy(0,-3, 0));
2207  return 0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3;
2208  }
2209 
2210  template<typename Accessor>
2211  static typename Accessor::ValueType inXandZ(const Accessor& grid, const Coord& ijk)
2212  {
2213  typename Accessor::ValueType tmp1 =
2214  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 1)) -
2215  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-1));
2216  typename Accessor::ValueType tmp2 =
2217  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 2)) -
2218  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-2));
2219  typename Accessor::ValueType tmp3 =
2220  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0, 3)) -
2221  D1<CD_6TH>::inX(grid, ijk.offsetBy(0, 0,-3));
2222  return 0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3;
2223  }
2224 
2225  template<typename Accessor>
2226  static typename Accessor::ValueType inYandZ(const Accessor& grid, const Coord& ijk)
2227  {
2228  typename Accessor::ValueType tmp1 =
2229  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 1)) -
2230  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-1));
2231  typename Accessor::ValueType tmp2 =
2232  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 2)) -
2233  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-2));
2234  typename Accessor::ValueType tmp3 =
2235  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0, 3)) -
2236  D1<CD_6TH>::inY(grid, ijk.offsetBy(0, 0,-3));
2237  return 0.75*tmp1 - 0.15*tmp2 + 1./60*tmp3;
2238  }
2239 
2240 
2241  // stencil access version
2242  template<typename Stencil>
2243  static typename Stencil::ValueType inX(const Stencil& S)
2244  {
2245  return difference( S.template getValue< 3, 0, 0>(), S.template getValue< 2, 0, 0>(),
2246  S.template getValue< 1, 0, 0>(), S.template getValue< 0, 0, 0>(),
2247  S.template getValue<-1, 0, 0>(), S.template getValue<-2, 0, 0>(),
2248  S.template getValue<-3, 0, 0>() );
2249  }
2250 
2251  template<typename Stencil>
2252  static typename Stencil::ValueType inY(const Stencil& S)
2253  {
2254  return difference( S.template getValue< 0, 3, 0>(), S.template getValue< 0, 2, 0>(),
2255  S.template getValue< 0, 1, 0>(), S.template getValue< 0, 0, 0>(),
2256  S.template getValue< 0,-1, 0>(), S.template getValue< 0,-2, 0>(),
2257  S.template getValue< 0,-3, 0>() );
2258 
2259  }
2260 
2261  template<typename Stencil>
2262  static typename Stencil::ValueType inZ(const Stencil& S)
2263  {
2264  return difference( S.template getValue< 0, 0, 3>(), S.template getValue< 0, 0, 2>(),
2265  S.template getValue< 0, 0, 1>(), S.template getValue< 0, 0, 0>(),
2266  S.template getValue< 0, 0,-1>(), S.template getValue< 0, 0,-2>(),
2267  S.template getValue< 0, 0,-3>() );
2268  }
2269 
2270  template<typename Stencil>
2271  static typename Stencil::ValueType inXandY(const Stencil& S)
2272  {
2273  return crossdifference( S.template getValue< 1, 1, 0>(), S.template getValue<-1, 1, 0>(),
2274  S.template getValue< 1,-1, 0>(), S.template getValue<-1,-1, 0>(),
2275  S.template getValue< 2, 1, 0>(), S.template getValue<-2, 1, 0>(),
2276  S.template getValue< 2,-1, 0>(), S.template getValue<-2,-1, 0>(),
2277  S.template getValue< 3, 1, 0>(), S.template getValue<-3, 1, 0>(),
2278  S.template getValue< 3,-1, 0>(), S.template getValue<-3,-1, 0>(),
2279  S.template getValue< 1, 2, 0>(), S.template getValue<-1, 2, 0>(),
2280  S.template getValue< 1,-2, 0>(), S.template getValue<-1,-2, 0>(),
2281  S.template getValue< 2, 2, 0>(), S.template getValue<-2, 2, 0>(),
2282  S.template getValue< 2,-2, 0>(), S.template getValue<-2,-2, 0>(),
2283  S.template getValue< 3, 2, 0>(), S.template getValue<-3, 2, 0>(),
2284  S.template getValue< 3,-2, 0>(), S.template getValue<-3,-2, 0>(),
2285  S.template getValue< 1, 3, 0>(), S.template getValue<-1, 3, 0>(),
2286  S.template getValue< 1,-3, 0>(), S.template getValue<-1,-3, 0>(),
2287  S.template getValue< 2, 3, 0>(), S.template getValue<-2, 3, 0>(),
2288  S.template getValue< 2,-3, 0>(), S.template getValue<-2,-3, 0>(),
2289  S.template getValue< 3, 3, 0>(), S.template getValue<-3, 3, 0>(),
2290  S.template getValue< 3,-3, 0>(), S.template getValue<-3,-3, 0>() );
2291  }
2292 
2293  template<typename Stencil>
2294  static typename Stencil::ValueType inXandZ(const Stencil& S)
2295  {
2296  return crossdifference( S.template getValue< 1, 0, 1>(), S.template getValue<-1, 0, 1>(),
2297  S.template getValue< 1, 0,-1>(), S.template getValue<-1, 0,-1>(),
2298  S.template getValue< 2, 0, 1>(), S.template getValue<-2, 0, 1>(),
2299  S.template getValue< 2, 0,-1>(), S.template getValue<-2, 0,-1>(),
2300  S.template getValue< 3, 0, 1>(), S.template getValue<-3, 0, 1>(),
2301  S.template getValue< 3, 0,-1>(), S.template getValue<-3, 0,-1>(),
2302  S.template getValue< 1, 0, 2>(), S.template getValue<-1, 0, 2>(),
2303  S.template getValue< 1, 0,-2>(), S.template getValue<-1, 0,-2>(),
2304  S.template getValue< 2, 0, 2>(), S.template getValue<-2, 0, 2>(),
2305  S.template getValue< 2, 0,-2>(), S.template getValue<-2, 0,-2>(),
2306  S.template getValue< 3, 0, 2>(), S.template getValue<-3, 0, 2>(),
2307  S.template getValue< 3, 0,-2>(), S.template getValue<-3, 0,-2>(),
2308  S.template getValue< 1, 0, 3>(), S.template getValue<-1, 0, 3>(),
2309  S.template getValue< 1, 0,-3>(), S.template getValue<-1, 0,-3>(),
2310  S.template getValue< 2, 0, 3>(), S.template getValue<-2, 0, 3>(),
2311  S.template getValue< 2, 0,-3>(), S.template getValue<-2, 0,-3>(),
2312  S.template getValue< 3, 0, 3>(), S.template getValue<-3, 0, 3>(),
2313  S.template getValue< 3, 0,-3>(), S.template getValue<-3, 0,-3>() );
2314  }
2315 
2316  template<typename Stencil>
2317  static typename Stencil::ValueType inYandZ(const Stencil& S)
2318  {
2319  return crossdifference( S.template getValue< 0, 1, 1>(), S.template getValue< 0,-1, 1>(),
2320  S.template getValue< 0, 1,-1>(), S.template getValue< 0,-1,-1>(),
2321  S.template getValue< 0, 2, 1>(), S.template getValue< 0,-2, 1>(),
2322  S.template getValue< 0, 2,-1>(), S.template getValue< 0,-2,-1>(),
2323  S.template getValue< 0, 3, 1>(), S.template getValue< 0,-3, 1>(),
2324  S.template getValue< 0, 3,-1>(), S.template getValue< 0,-3,-1>(),
2325  S.template getValue< 0, 1, 2>(), S.template getValue< 0,-1, 2>(),
2326  S.template getValue< 0, 1,-2>(), S.template getValue< 0,-1,-2>(),
2327  S.template getValue< 0, 2, 2>(), S.template getValue< 0,-2, 2>(),
2328  S.template getValue< 0, 2,-2>(), S.template getValue< 0,-2,-2>(),
2329  S.template getValue< 0, 3, 2>(), S.template getValue< 0,-3, 2>(),
2330  S.template getValue< 0, 3,-2>(), S.template getValue< 0,-3,-2>(),
2331  S.template getValue< 0, 1, 3>(), S.template getValue< 0,-1, 3>(),
2332  S.template getValue< 0, 1,-3>(), S.template getValue< 0,-1,-3>(),
2333  S.template getValue< 0, 2, 3>(), S.template getValue< 0,-2, 3>(),
2334  S.template getValue< 0, 2,-3>(), S.template getValue< 0,-2,-3>(),
2335  S.template getValue< 0, 3, 3>(), S.template getValue< 0,-3, 3>(),
2336  S.template getValue< 0, 3,-3>(), S.template getValue< 0,-3,-3>() );
2337  }
2338 
2339 };
2340 
2341 } // end math namespace
2342 } // namespace OPENVDB_VERSION_NAME
2343 } // end openvdb namespace
2344 
2345 #endif // OPENVDB_MATH_FINITEDIFFERENCE_HAS_BEEN_INCLUDED
2346 
2347 // Copyright (c) 2012-2013 DreamWorks Animation LLC
2348 // All rights reserved. This software is distributed under the
2349 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
static ValueType difference(const ValueType &xm2, const ValueType &xm1, const ValueType &xm0)
Definition: FiniteDifference.h:974
static Stencil::ValueType inYandZ(const Stencil &S)
Definition: FiniteDifference.h:2096
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:888
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:877
Definition: FiniteDifference.h:268
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:859
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:500
std::string dsSchemeToString(DScheme dss)
Definition: FiniteDifference.h:81
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:635
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1229
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:678
static Stencil::ValueType inXandY(const Stencil &S)
Definition: FiniteDifference.h:2271
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1846
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:806
std::string biasedGradientSchemeToMenuName(BiasedGradientScheme bgs)
Definition: FiniteDifference.h:242
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:788
static ValueType difference(const ValueType &xp1, const ValueType &xp0, const ValueType &xm1)
Definition: FiniteDifference.h:1832
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:941
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1747
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1147
Definition: FiniteDifference.h:201
Definition: FiniteDifference.h:1786
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:547
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1727
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1520
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
Definition: FiniteDifference.h:1767
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1974
Definition: FiniteDifference.h:64
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:1387
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1311
static Accessor::ValueType inXandZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:2211
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:762
Definition: FiniteDifference.h:70
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:825
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1064
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:491
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
Definition: FiniteDifference.h:1598
static Stencil::ValueType inYandZ(const Stencil &S)
Definition: FiniteDifference.h:1932
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:711
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1445
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:743
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:2043
std::string temporalIntegrationSchemeToString(TemporalIntegrationScheme tis)
Definition: FiniteDifference.h:271
DScheme stringToDScheme(const std::string &s)
Definition: FiniteDifference.h:105
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:998
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:688
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1685
static ValueType crossdifference(const ValueType &xpyp, const ValueType &xpym, const ValueType &xmyp, const ValueType &xmym)
Definition: FiniteDifference.h:1838
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1415
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:483
static Stencil::ValueType inXandZ(const Stencil &S)
Definition: FiniteDifference.h:2294
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:1372
static Stencil::ValueType inXandY(const Stencil &S)
Definition: FiniteDifference.h:1918
static ValueType difference(const ValueType &xm3, const ValueType &xm2, const ValueType &xm1, const ValueType &xm0, const ValueType &xp1, const ValueType &xp2)
Definition: FiniteDifference.h:1302
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:868
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:1326
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:557
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:2125
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:935
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:1007
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:660
Definition: FiniteDifference.h:264
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:617
Definition: FiniteDifference.h:69
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:1083
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:961
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:1491
Type Pow2(Type x)
Return .
Definition: Math.h:458
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:512
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1983
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:626
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:568
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1861
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1055
static ValueType difference(const ValueType &xp2, const ValueType &xp1, const ValueType &xm1, const ValueType &xm2)
Definition: FiniteDifference.h:581
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
Definition: FiniteDifference.h:1540
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1046
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:531
Definition: FiniteDifference.h:196
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:817
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:949
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:589
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:1461
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:562
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:982
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:1903
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:1023
static Accessor::ValueType inXandZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:2015
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:699
static Accessor::ValueType inXandY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:2002
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:1092
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1132
Real GudonovsNormSqrd(bool isOutside, const Vec3< Real > &gradient_m, const Vec3< Real > &gradient_p)
Definition: FiniteDifference.h:378
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:1896
std::string biasedGradientSchemeToString(BiasedGradientScheme bgs)
Definition: FiniteDifference.h:204
DScheme
Different discrete schemes used in the first derivatives.
Definition: FiniteDifference.h:59
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:955
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:1273
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:1015
static ValueType difference(const ValueType &xp2, const ValueType &xp1, const ValueType &xp0, const ValueType &xm1, const ValueType &xm2)
Definition: FiniteDifference.h:1946
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:505
TemporalIntegrationScheme
Temporal integration schemes.
Definition: FiniteDifference.h:261
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1992
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:756
Definition: FiniteDifference.h:265
#define OPENVDB_VERSION_NAME
Definition: version.h:45
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:47
Definition: FiniteDifference.h:1508
static ValueType difference(const ValueType &xp1, const ValueType &xp0)
Definition: FiniteDifference.h:730
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1341
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
Definition: FiniteDifference.h:1651
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
static ValueType difference(const ValueType &xp2, const ValueType &xp1, const ValueType &xp0)
Definition: FiniteDifference.h:780
Definition: FiniteDifference.h:193
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:797
static Accessor::ValueType inYandZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:2226
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:737
Definition: FiniteDifference.h:67
Definition: FiniteDifference.h:61
Definition: FiniteDifference.h:198
Definition: FiniteDifference.h:77
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:2116
Definition: FiniteDifference.h:60
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:1357
Definition: FiniteDifference.h:68
static Accessor::ValueType inYandZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:2028
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:768
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:475
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:1952
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:1176
std::string temporalIntegrationSchemeToMenuName(TemporalIntegrationScheme tis)
Definition: FiniteDifference.h:304
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1560
Definition: FiniteDifference.h:262
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:990
static Stencil::ValueType inXandY(const Stencil &S)
Definition: FiniteDifference.h:2068
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
Definition: FiniteDifference.h:1758
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:597
Definition: FiniteDifference.h:72
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:539
static ValueType difference(const ValueType &xp1, const ValueType &xm1)
Definition: FiniteDifference.h:469
static Stencil::ValueType inYandZ(const Stencil &S)
Definition: FiniteDifference.h:2317
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1737
Definition: FiniteDifference.h:66
Definition: FiniteDifference.h:195
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1628
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:2262
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:833
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1612
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:2165
Definition: FiniteDifference.h:63
Coord offsetBy(Int32 dx, Int32 dy, Int32 dz) const
Definition: Coord.h:102
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:906
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:1910
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1244
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1526
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
Definition: FiniteDifference.h:1694
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:669
Definition: FiniteDifference.h:178
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xp0, const ValueType &xm1, const ValueType &xm2)
Definition: FiniteDifference.h:1205
Definition: FiniteDifference.h:180
Definition: FiniteDifference.h:71
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:2051
std::string dsSchemeToMenuName(DScheme dss)
Definition: FiniteDifference.h:147
Definition: FiniteDifference.h:74
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1213
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1676
static Accessor::ValueType::value_type inZ(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1576
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:506
TemporalIntegrationScheme stringToTemporalIntegrationScheme(const std::string &s)
Definition: FiniteDifference.h:284
static Stencil::ValueType inXandZ(const Stencil &S)
Definition: FiniteDifference.h:2082
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
Definition: FiniteDifference.h:1534
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
Definition: FiniteDifference.h:1584
double Real
Definition: Types.h:63
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:2243
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
Definition: FiniteDifference.h:1546
Definition: FiniteDifference.h:263
static Accessor::ValueType inXandY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1869
BiasedGradientScheme stringToBiasedGradientScheme(const std::string &s)
Definition: FiniteDifference.h:219
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:929
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:2059
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:1284
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1430
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
Definition: FiniteDifference.h:1644
static Stencil::ValueType::value_type inX(const Stencil &S, int n)
Definition: FiniteDifference.h:1637
Definition: FiniteDifference.h:441
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:2175
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:2252
Definition: FiniteDifference.h:179
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
Definition: FiniteDifference.h:1591
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
Definition: FiniteDifference.h:1776
static Accessor::ValueType inXandY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:2196
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
static Accessor::ValueType inYandZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1886
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1513
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:2185
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1620
Definition: FiniteDifference.h:184
static Accessor::ValueType::value_type inY(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1568
static Stencil::ValueType inZ(const Stencil &S)
Definition: FiniteDifference.h:1187
static Accessor::ValueType inXandZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1878
static Accessor::ValueType inX(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1117
static Stencil::ValueType inXandZ(const Stencil &S)
Definition: FiniteDifference.h:1925
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xp0)
Definition: FiniteDifference.h:849
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:1074
Definition: FiniteDifference.h:62
Definition: FiniteDifference.h:181
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:1260
static Accessor::ValueType inY(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:1853
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:606
static ValueType difference(const ValueType &xp3, const ValueType &xp2, const ValueType &xp1, const ValueType &xm1, const ValueType &xm2, const ValueType &xm3)
Definition: FiniteDifference.h:650
static ValueType difference(const ValueType &xm3, const ValueType &xm2, const ValueType &xm1, const ValueType &xm0, const ValueType &xp1, const ValueType &xp2)
Definition: FiniteDifference.h:1407
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:897
static Accessor::ValueType inZ(const Accessor &grid, const Coord &ijk)
Definition: FiniteDifference.h:749
DDScheme
Different discrete schemes used in the second derivatives.
Definition: FiniteDifference.h:177
static ValueType difference(const ValueType &xm3, const ValueType &xm2, const ValueType &xm1, const ValueType &xm0)
Definition: FiniteDifference.h:1038
static Stencil::ValueType inX(const Stencil &S)
Definition: FiniteDifference.h:1163
static Stencil::ValueType::value_type inZ(const Stencil &S, int n)
Definition: FiniteDifference.h:1710
Definition: FiniteDifference.h:65
static Accessor::ValueType::value_type inX(const Accessor &grid, const Coord &ijk, int n)
Definition: FiniteDifference.h:1667
Definition: FiniteDifference.h:197
Definition: FiniteDifference.h:73
static Stencil::ValueType::value_type inY(const Stencil &S, int n)
Definition: FiniteDifference.h:1702
Definition: FiniteDifference.h:194
static ValueType difference(const ValueType &xp1, const ValueType &xm1)
Definition: FiniteDifference.h:524
static Stencil::ValueType inY(const Stencil &S)
Definition: FiniteDifference.h:1476
ValueType WENO5(const ValueType &v1, const ValueType &v2, const ValueType &v3, const ValueType &v4, const ValueType &v5, float scale2=0.01)
Implementation of nominally fifth-order finite-difference WENO.
Definition: FiniteDifference.h:331
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:566