OpenVDB  8.1.1
Stencils.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
12 
13 #ifndef OPENVDB_MATH_STENCILS_HAS_BEEN_INCLUDED
14 #define OPENVDB_MATH_STENCILS_HAS_BEEN_INCLUDED
15 
16 #include <algorithm>
17 #include <vector> // for std::vector
18 #include <bitset> // for std::bitset
19 #include <openvdb/math/Math.h> // for Pow2, needed by WENO and Godunov
20 #include <openvdb/Types.h> // for Real
21 #include <openvdb/math/Coord.h> // for Coord
22 #include <openvdb/math/FiniteDifference.h> // for WENO5 and GodunovsNormSqrd
24 
25 namespace openvdb {
27 namespace OPENVDB_VERSION_NAME {
28 namespace math {
29 
30 
32 
33 template<typename DerivedType, typename GridT, bool IsSafe>
35 {
36 public:
37  typedef GridT GridType;
38  typedef typename GridT::TreeType TreeType;
39  typedef typename GridT::ValueType ValueType;
41  typedef std::vector<ValueType> BufferType;
42  typedef typename BufferType::iterator IterType;
43 
47  inline void moveTo(const Coord& ijk)
48  {
49  mCenter = ijk;
50  mValues[0] = mAcc.getValue(ijk);
51  static_cast<DerivedType&>(*this).init(mCenter);
52  }
53 
59  inline void moveTo(const Coord& ijk, const ValueType& centerValue)
60  {
61  mCenter = ijk;
62  mValues[0] = centerValue;
63  static_cast<DerivedType&>(*this).init(mCenter);
64  }
65 
71  template<typename IterType>
72  inline void moveTo(const IterType& iter)
73  {
74  mCenter = iter.getCoord();
75  mValues[0] = *iter;
76  static_cast<DerivedType&>(*this).init(mCenter);
77  }
78 
85  template<typename RealType>
86  inline void moveTo(const Vec3<RealType>& xyz)
87  {
88  Coord ijk = Coord::floor(xyz);
89  if (ijk != mCenter) this->moveTo(ijk);
90  }
91 
97  inline const ValueType& getValue(unsigned int pos = 0) const
98  {
99  assert(pos < mValues.size());
100  return mValues[pos];
101  }
102 
104  template<int i, int j, int k>
105  inline const ValueType& getValue() const
106  {
107  return mValues[static_cast<const DerivedType&>(*this).template pos<i,j,k>()];
108  }
109 
111  template<int i, int j, int k>
112  inline void setValue(const ValueType& value)
113  {
114  mValues[static_cast<const DerivedType&>(*this).template pos<i,j,k>()] = value;
115  }
116 
118  inline int size() { return mValues.size(); }
119 
121  inline ValueType median() const
122  {
123  BufferType tmp(mValues);//local copy
124  assert(!tmp.empty());
125  size_t midpoint = (tmp.size() - 1) >> 1;
126  // Partially sort the vector until the median value is at the midpoint.
127 #if !defined(_MSC_VER) || _MSC_VER < 1924
128  std::nth_element(tmp.begin(), tmp.begin() + midpoint, tmp.end());
129 #else
130  // Workaround MSVC bool warning C4804 unsafe use of type 'bool'
131  std::nth_element(tmp.begin(), tmp.begin() + midpoint, tmp.end(),
132  std::less<ValueType>());
133 #endif
134  return tmp[midpoint];
135  }
136 
138  inline ValueType mean() const
139  {
140  ValueType sum = 0.0;
141  for (int n = 0, s = int(mValues.size()); n < s; ++n) sum += mValues[n];
142  return sum / ValueType(mValues.size());
143  }
144 
146  inline ValueType min() const
147  {
148  IterType iter = std::min_element(mValues.begin(), mValues.end());
149  return *iter;
150  }
151 
153  inline ValueType max() const
154  {
155  IterType iter = std::max_element(mValues.begin(), mValues.end());
156  return *iter;
157  }
158 
160  inline const Coord& getCenterCoord() const { return mCenter; }
161 
163  inline const ValueType& getCenterValue() const { return mValues[0]; }
164 
167  inline bool intersects(const ValueType &isoValue = zeroVal<ValueType>()) const
168  {
169  const bool less = this->getValue< 0, 0, 0>() < isoValue;
170  return (less ^ (this->getValue<-1, 0, 0>() < isoValue)) ||
171  (less ^ (this->getValue< 1, 0, 0>() < isoValue)) ||
172  (less ^ (this->getValue< 0,-1, 0>() < isoValue)) ||
173  (less ^ (this->getValue< 0, 1, 0>() < isoValue)) ||
174  (less ^ (this->getValue< 0, 0,-1>() < isoValue)) ||
175  (less ^ (this->getValue< 0, 0, 1>() < isoValue)) ;
176  }
177 
188  inline std::bitset<6> intersectionMask(const ValueType &isoValue = zeroVal<ValueType>()) const
189  {
190  std::bitset<6> mask;
191  const bool less = this->getValue< 0, 0, 0>() < isoValue;
192  mask[0] = less ^ (this->getValue<-1, 0, 0>() < isoValue);
193  mask[1] = less ^ (this->getValue< 1, 0, 0>() < isoValue);
194  mask[2] = less ^ (this->getValue< 0,-1, 0>() < isoValue);
195  mask[3] = less ^ (this->getValue< 0, 1, 0>() < isoValue);
196  mask[4] = less ^ (this->getValue< 0, 0,-1>() < isoValue);
197  mask[5] = less ^ (this->getValue< 0, 0, 1>() < isoValue);
198  return mask;
199  }
200 
203  inline const GridType& grid() const { return *mGrid; }
204 
207  inline const AccessorType& accessor() const { return mAcc; }
208 
209 protected:
210  // Constructor is protected to prevent direct instantiation.
211  BaseStencil(const GridType& grid, int size)
212  : mGrid(&grid)
213  , mAcc(grid.tree())
214  , mValues(size)
215  , mCenter(Coord::max())
216  {
217  }
218 
219  const GridType* mGrid;
220  AccessorType mAcc;
221  BufferType mValues;
223 
224 }; // BaseStencil class
225 
226 
228 
229 
230 namespace { // anonymous namespace for stencil-layout map
231 
232  // the seven point stencil
233  template<int i, int j, int k> struct SevenPt {};
234  template<> struct SevenPt< 0, 0, 0> { enum { idx = 0 }; };
235  template<> struct SevenPt< 1, 0, 0> { enum { idx = 1 }; };
236  template<> struct SevenPt< 0, 1, 0> { enum { idx = 2 }; };
237  template<> struct SevenPt< 0, 0, 1> { enum { idx = 3 }; };
238  template<> struct SevenPt<-1, 0, 0> { enum { idx = 4 }; };
239  template<> struct SevenPt< 0,-1, 0> { enum { idx = 5 }; };
240  template<> struct SevenPt< 0, 0,-1> { enum { idx = 6 }; };
241 
242 }
243 
244 
245 template<typename GridT, bool IsSafe = true>
246 class SevenPointStencil: public BaseStencil<SevenPointStencil<GridT, IsSafe>, GridT, IsSafe>
247 {
250 public:
251  typedef GridT GridType;
252  typedef typename GridT::TreeType TreeType;
253  typedef typename GridT::ValueType ValueType;
254 
255  static const int SIZE = 7;
256 
257  SevenPointStencil(const GridT& grid): BaseType(grid, SIZE) {}
258 
260  template<int i, int j, int k>
261  unsigned int pos() const { return SevenPt<i,j,k>::idx; }
262 
263 private:
264  inline void init(const Coord& ijk)
265  {
266  BaseType::template setValue<-1, 0, 0>(mAcc.getValue(ijk.offsetBy(-1, 0, 0)));
267  BaseType::template setValue< 1, 0, 0>(mAcc.getValue(ijk.offsetBy( 1, 0, 0)));
268 
269  BaseType::template setValue< 0,-1, 0>(mAcc.getValue(ijk.offsetBy( 0,-1, 0)));
270  BaseType::template setValue< 0, 1, 0>(mAcc.getValue(ijk.offsetBy( 0, 1, 0)));
271 
272  BaseType::template setValue< 0, 0,-1>(mAcc.getValue(ijk.offsetBy( 0, 0,-1)));
273  BaseType::template setValue< 0, 0, 1>(mAcc.getValue(ijk.offsetBy( 0, 0, 1)));
274  }
275 
276  template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
277  using BaseType::mAcc;
278  using BaseType::mValues;
279 };// SevenPointStencil class
280 
281 
283 
284 
285 namespace { // anonymous namespace for stencil-layout map
286 
287  // the eight point box stencil
288  template<int i, int j, int k> struct BoxPt {};
289  template<> struct BoxPt< 0, 0, 0> { enum { idx = 0 }; };
290  template<> struct BoxPt< 0, 0, 1> { enum { idx = 1 }; };
291  template<> struct BoxPt< 0, 1, 1> { enum { idx = 2 }; };
292  template<> struct BoxPt< 0, 1, 0> { enum { idx = 3 }; };
293  template<> struct BoxPt< 1, 0, 0> { enum { idx = 4 }; };
294  template<> struct BoxPt< 1, 0, 1> { enum { idx = 5 }; };
295  template<> struct BoxPt< 1, 1, 1> { enum { idx = 6 }; };
296  template<> struct BoxPt< 1, 1, 0> { enum { idx = 7 }; };
297 }
298 
299 template<typename GridT, bool IsSafe = true>
300 class BoxStencil: public BaseStencil<BoxStencil<GridT, IsSafe>, GridT, IsSafe>
301 {
304 public:
305  typedef GridT GridType;
306  typedef typename GridT::TreeType TreeType;
307  typedef typename GridT::ValueType ValueType;
308 
309  static const int SIZE = 8;
310 
311  BoxStencil(const GridType& grid): BaseType(grid, SIZE) {}
312 
314  template<int i, int j, int k>
315  unsigned int pos() const { return BoxPt<i,j,k>::idx; }
316 
319  inline bool intersects(const ValueType &isoValue = zeroVal<ValueType>()) const
320  {
321  const bool less = mValues[0] < isoValue;
322  return (less ^ (mValues[1] < isoValue)) ||
323  (less ^ (mValues[2] < isoValue)) ||
324  (less ^ (mValues[3] < isoValue)) ||
325  (less ^ (mValues[4] < isoValue)) ||
326  (less ^ (mValues[5] < isoValue)) ||
327  (less ^ (mValues[6] < isoValue)) ||
328  (less ^ (mValues[7] < isoValue)) ;
329  }
330 
338  inline ValueType interpolation(const math::Vec3<ValueType>& xyz) const
339  {
341  const ValueType u = xyz[0] - BaseType::mCenter[0];
342  const ValueType v = xyz[1] - BaseType::mCenter[1];
343  const ValueType w = xyz[2] - BaseType::mCenter[2];
345 
346  assert(u>=0 && u<=1);
347  assert(v>=0 && v<=1);
348  assert(w>=0 && w<=1);
349 
350  ValueType V = BaseType::template getValue<0,0,0>();
351  ValueType A = static_cast<ValueType>(V + (BaseType::template getValue<0,0,1>() - V) * w);
352  V = BaseType::template getValue< 0, 1, 0>();
353  ValueType B = static_cast<ValueType>(V + (BaseType::template getValue<0,1,1>() - V) * w);
354  ValueType C = static_cast<ValueType>(A + (B - A) * v);
355 
356  V = BaseType::template getValue<1,0,0>();
357  A = static_cast<ValueType>(V + (BaseType::template getValue<1,0,1>() - V) * w);
358  V = BaseType::template getValue<1,1,0>();
359  B = static_cast<ValueType>(V + (BaseType::template getValue<1,1,1>() - V) * w);
360  ValueType D = static_cast<ValueType>(A + (B - A) * v);
361 
362  return static_cast<ValueType>(C + (D - C) * u);
363  }
364 
373  {
375  const ValueType u = xyz[0] - BaseType::mCenter[0];
376  const ValueType v = xyz[1] - BaseType::mCenter[1];
377  const ValueType w = xyz[2] - BaseType::mCenter[2];
379 
380  assert(u>=0 && u<=1);
381  assert(v>=0 && v<=1);
382  assert(w>=0 && w<=1);
383 
384  ValueType D[4]={BaseType::template getValue<0,0,1>()-BaseType::template getValue<0,0,0>(),
385  BaseType::template getValue<0,1,1>()-BaseType::template getValue<0,1,0>(),
386  BaseType::template getValue<1,0,1>()-BaseType::template getValue<1,0,0>(),
387  BaseType::template getValue<1,1,1>()-BaseType::template getValue<1,1,0>()};
388 
389  // Z component
390  ValueType A = static_cast<ValueType>(D[0] + (D[1]- D[0]) * v);
391  ValueType B = static_cast<ValueType>(D[2] + (D[3]- D[2]) * v);
392  math::Vec3<ValueType> grad(zeroVal<ValueType>(),
393  zeroVal<ValueType>(),
394  static_cast<ValueType>(A + (B - A) * u));
395 
396  D[0] = static_cast<ValueType>(BaseType::template getValue<0,0,0>() + D[0] * w);
397  D[1] = static_cast<ValueType>(BaseType::template getValue<0,1,0>() + D[1] * w);
398  D[2] = static_cast<ValueType>(BaseType::template getValue<1,0,0>() + D[2] * w);
399  D[3] = static_cast<ValueType>(BaseType::template getValue<1,1,0>() + D[3] * w);
400 
401  // X component
402  A = static_cast<ValueType>(D[0] + (D[1] - D[0]) * v);
403  B = static_cast<ValueType>(D[2] + (D[3] - D[2]) * v);
404 
405  grad[0] = B - A;
406 
407  // Y component
408  A = D[1] - D[0];
409  B = D[3] - D[2];
410 
411  grad[1] = static_cast<ValueType>(A + (B - A) * u);
412 
413  return BaseType::mGrid->transform().baseMap()->applyIJT(grad, xyz);
414  }
415 
416 private:
417  inline void init(const Coord& ijk)
418  {
419  BaseType::template setValue< 0, 0, 1>(mAcc.getValue(ijk.offsetBy( 0, 0, 1)));
420  BaseType::template setValue< 0, 1, 1>(mAcc.getValue(ijk.offsetBy( 0, 1, 1)));
421  BaseType::template setValue< 0, 1, 0>(mAcc.getValue(ijk.offsetBy( 0, 1, 0)));
422  BaseType::template setValue< 1, 0, 0>(mAcc.getValue(ijk.offsetBy( 1, 0, 0)));
423  BaseType::template setValue< 1, 0, 1>(mAcc.getValue(ijk.offsetBy( 1, 0, 1)));
424  BaseType::template setValue< 1, 1, 1>(mAcc.getValue(ijk.offsetBy( 1, 1, 1)));
425  BaseType::template setValue< 1, 1, 0>(mAcc.getValue(ijk.offsetBy( 1, 1, 0)));
426  }
427 
428  template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
429  using BaseType::mAcc;
430  using BaseType::mValues;
431 };// BoxStencil class
432 
433 
435 
436 
437 namespace { // anonymous namespace for stencil-layout map
438 
439  // the dense point stencil
440  template<int i, int j, int k> struct DensePt {};
441  template<> struct DensePt< 0, 0, 0> { enum { idx = 0 }; };
442 
443  template<> struct DensePt< 1, 0, 0> { enum { idx = 1 }; };
444  template<> struct DensePt< 0, 1, 0> { enum { idx = 2 }; };
445  template<> struct DensePt< 0, 0, 1> { enum { idx = 3 }; };
446 
447  template<> struct DensePt<-1, 0, 0> { enum { idx = 4 }; };
448  template<> struct DensePt< 0,-1, 0> { enum { idx = 5 }; };
449  template<> struct DensePt< 0, 0,-1> { enum { idx = 6 }; };
450 
451  template<> struct DensePt<-1,-1, 0> { enum { idx = 7 }; };
452  template<> struct DensePt< 0,-1,-1> { enum { idx = 8 }; };
453  template<> struct DensePt<-1, 0,-1> { enum { idx = 9 }; };
454 
455  template<> struct DensePt< 1,-1, 0> { enum { idx = 10 }; };
456  template<> struct DensePt< 0, 1,-1> { enum { idx = 11 }; };
457  template<> struct DensePt<-1, 0, 1> { enum { idx = 12 }; };
458 
459  template<> struct DensePt<-1, 1, 0> { enum { idx = 13 }; };
460  template<> struct DensePt< 0,-1, 1> { enum { idx = 14 }; };
461  template<> struct DensePt< 1, 0,-1> { enum { idx = 15 }; };
462 
463  template<> struct DensePt< 1, 1, 0> { enum { idx = 16 }; };
464  template<> struct DensePt< 0, 1, 1> { enum { idx = 17 }; };
465  template<> struct DensePt< 1, 0, 1> { enum { idx = 18 }; };
466 }
467 
468 
469 template<typename GridT, bool IsSafe = true>
471  : public BaseStencil<SecondOrderDenseStencil<GridT, IsSafe>, GridT, IsSafe >
472 {
474  typedef BaseStencil<SelfT, GridT, IsSafe > BaseType;
475 public:
476  typedef GridT GridType;
477  typedef typename GridT::TreeType TreeType;
478  typedef typename GridType::ValueType ValueType;
479 
480  static const int SIZE = 19;
481 
482  SecondOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
483 
485  template<int i, int j, int k>
486  unsigned int pos() const { return DensePt<i,j,k>::idx; }
487 
488 private:
489  inline void init(const Coord& ijk)
490  {
491  mValues[DensePt< 1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
492  mValues[DensePt< 0, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
493  mValues[DensePt< 0, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
494 
495  mValues[DensePt<-1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
496  mValues[DensePt< 0,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
497  mValues[DensePt< 0, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
498 
499  mValues[DensePt<-1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, -1, 0));
500  mValues[DensePt< 1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, -1, 0));
501  mValues[DensePt<-1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 1, 0));
502  mValues[DensePt< 1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 1, 0));
503 
504  mValues[DensePt<-1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, -1));
505  mValues[DensePt< 1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, -1));
506  mValues[DensePt<-1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 1));
507  mValues[DensePt< 1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 1));
508 
509  mValues[DensePt< 0,-1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, -1, -1));
510  mValues[DensePt< 0, 1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, -1));
511  mValues[DensePt< 0,-1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, -1, 1));
512  mValues[DensePt< 0, 1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 1));
513  }
514 
515  template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
516  using BaseType::mAcc;
517  using BaseType::mValues;
518 };// SecondOrderDenseStencil class
519 
520 
522 
523 
524 namespace { // anonymous namespace for stencil-layout map
525 
526  // the dense point stencil
527  template<int i, int j, int k> struct ThirteenPt {};
528  template<> struct ThirteenPt< 0, 0, 0> { enum { idx = 0 }; };
529 
530  template<> struct ThirteenPt< 1, 0, 0> { enum { idx = 1 }; };
531  template<> struct ThirteenPt< 0, 1, 0> { enum { idx = 2 }; };
532  template<> struct ThirteenPt< 0, 0, 1> { enum { idx = 3 }; };
533 
534  template<> struct ThirteenPt<-1, 0, 0> { enum { idx = 4 }; };
535  template<> struct ThirteenPt< 0,-1, 0> { enum { idx = 5 }; };
536  template<> struct ThirteenPt< 0, 0,-1> { enum { idx = 6 }; };
537 
538  template<> struct ThirteenPt< 2, 0, 0> { enum { idx = 7 }; };
539  template<> struct ThirteenPt< 0, 2, 0> { enum { idx = 8 }; };
540  template<> struct ThirteenPt< 0, 0, 2> { enum { idx = 9 }; };
541 
542  template<> struct ThirteenPt<-2, 0, 0> { enum { idx = 10 }; };
543  template<> struct ThirteenPt< 0,-2, 0> { enum { idx = 11 }; };
544  template<> struct ThirteenPt< 0, 0,-2> { enum { idx = 12 }; };
545 
546 }
547 
548 
549 template<typename GridT, bool IsSafe = true>
551  : public BaseStencil<ThirteenPointStencil<GridT, IsSafe>, GridT, IsSafe>
552 {
554  typedef BaseStencil<SelfT, GridT, IsSafe > BaseType;
555 public:
556  typedef GridT GridType;
557  typedef typename GridT::TreeType TreeType;
558  typedef typename GridType::ValueType ValueType;
559 
560  static const int SIZE = 13;
561 
562  ThirteenPointStencil(const GridType& grid): BaseType(grid, SIZE) {}
563 
565  template<int i, int j, int k>
566  unsigned int pos() const { return ThirteenPt<i,j,k>::idx; }
567 
568 private:
569  inline void init(const Coord& ijk)
570  {
571  mValues[ThirteenPt< 2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 0));
572  mValues[ThirteenPt< 1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
573  mValues[ThirteenPt<-1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
574  mValues[ThirteenPt<-2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 0));
575 
576  mValues[ThirteenPt< 0, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 0));
577  mValues[ThirteenPt< 0, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
578  mValues[ThirteenPt< 0,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
579  mValues[ThirteenPt< 0,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -2, 0));
580 
581  mValues[ThirteenPt< 0, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 2));
582  mValues[ThirteenPt< 0, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
583  mValues[ThirteenPt< 0, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
584  mValues[ThirteenPt< 0, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -2));
585  }
586 
587  template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
588  using BaseType::mAcc;
589  using BaseType::mValues;
590 };// ThirteenPointStencil class
591 
592 
594 
595 
596 namespace { // anonymous namespace for stencil-layout map
597 
598  // the 4th-order dense point stencil
599  template<int i, int j, int k> struct FourthDensePt {};
600  template<> struct FourthDensePt< 0, 0, 0> { enum { idx = 0 }; };
601 
602  template<> struct FourthDensePt<-2, 2, 0> { enum { idx = 1 }; };
603  template<> struct FourthDensePt<-1, 2, 0> { enum { idx = 2 }; };
604  template<> struct FourthDensePt< 0, 2, 0> { enum { idx = 3 }; };
605  template<> struct FourthDensePt< 1, 2, 0> { enum { idx = 4 }; };
606  template<> struct FourthDensePt< 2, 2, 0> { enum { idx = 5 }; };
607 
608  template<> struct FourthDensePt<-2, 1, 0> { enum { idx = 6 }; };
609  template<> struct FourthDensePt<-1, 1, 0> { enum { idx = 7 }; };
610  template<> struct FourthDensePt< 0, 1, 0> { enum { idx = 8 }; };
611  template<> struct FourthDensePt< 1, 1, 0> { enum { idx = 9 }; };
612  template<> struct FourthDensePt< 2, 1, 0> { enum { idx = 10 }; };
613 
614  template<> struct FourthDensePt<-2, 0, 0> { enum { idx = 11 }; };
615  template<> struct FourthDensePt<-1, 0, 0> { enum { idx = 12 }; };
616  template<> struct FourthDensePt< 1, 0, 0> { enum { idx = 13 }; };
617  template<> struct FourthDensePt< 2, 0, 0> { enum { idx = 14 }; };
618 
619  template<> struct FourthDensePt<-2,-1, 0> { enum { idx = 15 }; };
620  template<> struct FourthDensePt<-1,-1, 0> { enum { idx = 16 }; };
621  template<> struct FourthDensePt< 0,-1, 0> { enum { idx = 17 }; };
622  template<> struct FourthDensePt< 1,-1, 0> { enum { idx = 18 }; };
623  template<> struct FourthDensePt< 2,-1, 0> { enum { idx = 19 }; };
624 
625  template<> struct FourthDensePt<-2,-2, 0> { enum { idx = 20 }; };
626  template<> struct FourthDensePt<-1,-2, 0> { enum { idx = 21 }; };
627  template<> struct FourthDensePt< 0,-2, 0> { enum { idx = 22 }; };
628  template<> struct FourthDensePt< 1,-2, 0> { enum { idx = 23 }; };
629  template<> struct FourthDensePt< 2,-2, 0> { enum { idx = 24 }; };
630 
631 
632  template<> struct FourthDensePt<-2, 0, 2> { enum { idx = 25 }; };
633  template<> struct FourthDensePt<-1, 0, 2> { enum { idx = 26 }; };
634  template<> struct FourthDensePt< 0, 0, 2> { enum { idx = 27 }; };
635  template<> struct FourthDensePt< 1, 0, 2> { enum { idx = 28 }; };
636  template<> struct FourthDensePt< 2, 0, 2> { enum { idx = 29 }; };
637 
638  template<> struct FourthDensePt<-2, 0, 1> { enum { idx = 30 }; };
639  template<> struct FourthDensePt<-1, 0, 1> { enum { idx = 31 }; };
640  template<> struct FourthDensePt< 0, 0, 1> { enum { idx = 32 }; };
641  template<> struct FourthDensePt< 1, 0, 1> { enum { idx = 33 }; };
642  template<> struct FourthDensePt< 2, 0, 1> { enum { idx = 34 }; };
643 
644  template<> struct FourthDensePt<-2, 0,-1> { enum { idx = 35 }; };
645  template<> struct FourthDensePt<-1, 0,-1> { enum { idx = 36 }; };
646  template<> struct FourthDensePt< 0, 0,-1> { enum { idx = 37 }; };
647  template<> struct FourthDensePt< 1, 0,-1> { enum { idx = 38 }; };
648  template<> struct FourthDensePt< 2, 0,-1> { enum { idx = 39 }; };
649 
650  template<> struct FourthDensePt<-2, 0,-2> { enum { idx = 40 }; };
651  template<> struct FourthDensePt<-1, 0,-2> { enum { idx = 41 }; };
652  template<> struct FourthDensePt< 0, 0,-2> { enum { idx = 42 }; };
653  template<> struct FourthDensePt< 1, 0,-2> { enum { idx = 43 }; };
654  template<> struct FourthDensePt< 2, 0,-2> { enum { idx = 44 }; };
655 
656 
657  template<> struct FourthDensePt< 0,-2, 2> { enum { idx = 45 }; };
658  template<> struct FourthDensePt< 0,-1, 2> { enum { idx = 46 }; };
659  template<> struct FourthDensePt< 0, 1, 2> { enum { idx = 47 }; };
660  template<> struct FourthDensePt< 0, 2, 2> { enum { idx = 48 }; };
661 
662  template<> struct FourthDensePt< 0,-2, 1> { enum { idx = 49 }; };
663  template<> struct FourthDensePt< 0,-1, 1> { enum { idx = 50 }; };
664  template<> struct FourthDensePt< 0, 1, 1> { enum { idx = 51 }; };
665  template<> struct FourthDensePt< 0, 2, 1> { enum { idx = 52 }; };
666 
667  template<> struct FourthDensePt< 0,-2,-1> { enum { idx = 53 }; };
668  template<> struct FourthDensePt< 0,-1,-1> { enum { idx = 54 }; };
669  template<> struct FourthDensePt< 0, 1,-1> { enum { idx = 55 }; };
670  template<> struct FourthDensePt< 0, 2,-1> { enum { idx = 56 }; };
671 
672  template<> struct FourthDensePt< 0,-2,-2> { enum { idx = 57 }; };
673  template<> struct FourthDensePt< 0,-1,-2> { enum { idx = 58 }; };
674  template<> struct FourthDensePt< 0, 1,-2> { enum { idx = 59 }; };
675  template<> struct FourthDensePt< 0, 2,-2> { enum { idx = 60 }; };
676 
677 }
678 
679 
680 template<typename GridT, bool IsSafe = true>
682  : public BaseStencil<FourthOrderDenseStencil<GridT, IsSafe>, GridT, IsSafe>
683 {
685  typedef BaseStencil<SelfT, GridT, IsSafe > BaseType;
686 public:
687  typedef GridT GridType;
688  typedef typename GridT::TreeType TreeType;
689  typedef typename GridType::ValueType ValueType;
690 
691  static const int SIZE = 61;
692 
693  FourthOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
694 
696  template<int i, int j, int k>
697  unsigned int pos() const { return FourthDensePt<i,j,k>::idx; }
698 
699 private:
700  inline void init(const Coord& ijk)
701  {
702  mValues[FourthDensePt<-2, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 2, 0));
703  mValues[FourthDensePt<-1, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 2, 0));
704  mValues[FourthDensePt< 0, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 0));
705  mValues[FourthDensePt< 1, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 2, 0));
706  mValues[FourthDensePt< 2, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 2, 0));
707 
708  mValues[FourthDensePt<-2, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 1, 0));
709  mValues[FourthDensePt<-1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 1, 0));
710  mValues[FourthDensePt< 0, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
711  mValues[FourthDensePt< 1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 1, 0));
712  mValues[FourthDensePt< 2, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 1, 0));
713 
714  mValues[FourthDensePt<-2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 0));
715  mValues[FourthDensePt<-1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
716  mValues[FourthDensePt< 1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
717  mValues[FourthDensePt< 2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 0));
718 
719  mValues[FourthDensePt<-2,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2,-1, 0));
720  mValues[FourthDensePt<-1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1,-1, 0));
721  mValues[FourthDensePt< 0,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 0));
722  mValues[FourthDensePt< 1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1,-1, 0));
723  mValues[FourthDensePt< 2,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2,-1, 0));
724 
725  mValues[FourthDensePt<-2,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2,-2, 0));
726  mValues[FourthDensePt<-1,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1,-2, 0));
727  mValues[FourthDensePt< 0,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 0));
728  mValues[FourthDensePt< 1,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1,-2, 0));
729  mValues[FourthDensePt< 2,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2,-2, 0));
730 
731  mValues[FourthDensePt<-2, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 2));
732  mValues[FourthDensePt<-1, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 2));
733  mValues[FourthDensePt< 0, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 2));
734  mValues[FourthDensePt< 1, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 2));
735  mValues[FourthDensePt< 2, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 2));
736 
737  mValues[FourthDensePt<-2, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 1));
738  mValues[FourthDensePt<-1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 1));
739  mValues[FourthDensePt< 0, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
740  mValues[FourthDensePt< 1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 1));
741  mValues[FourthDensePt< 2, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 1));
742 
743  mValues[FourthDensePt<-2, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0,-1));
744  mValues[FourthDensePt<-1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0,-1));
745  mValues[FourthDensePt< 0, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0,-1));
746  mValues[FourthDensePt< 1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0,-1));
747  mValues[FourthDensePt< 2, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0,-1));
748 
749  mValues[FourthDensePt<-2, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0,-2));
750  mValues[FourthDensePt<-1, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0,-2));
751  mValues[FourthDensePt< 0, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0,-2));
752  mValues[FourthDensePt< 1, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0,-2));
753  mValues[FourthDensePt< 2, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0,-2));
754 
755 
756  mValues[FourthDensePt< 0,-2, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 2));
757  mValues[FourthDensePt< 0,-1, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 2));
758  mValues[FourthDensePt< 0, 1, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 2));
759  mValues[FourthDensePt< 0, 2, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 2));
760 
761  mValues[FourthDensePt< 0,-2, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 1));
762  mValues[FourthDensePt< 0,-1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 1));
763  mValues[FourthDensePt< 0, 1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 1));
764  mValues[FourthDensePt< 0, 2, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 1));
765 
766  mValues[FourthDensePt< 0,-2,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2,-1));
767  mValues[FourthDensePt< 0,-1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1,-1));
768  mValues[FourthDensePt< 0, 1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1,-1));
769  mValues[FourthDensePt< 0, 2,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2,-1));
770 
771  mValues[FourthDensePt< 0,-2,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2,-2));
772  mValues[FourthDensePt< 0,-1,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1,-2));
773  mValues[FourthDensePt< 0, 1,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1,-2));
774  mValues[FourthDensePt< 0, 2,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2,-2));
775  }
776 
777  template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
778  using BaseType::mAcc;
779  using BaseType::mValues;
780 };// FourthOrderDenseStencil class
781 
782 
784 
785 
786 namespace { // anonymous namespace for stencil-layout map
787 
788  // the dense point stencil
789  template<int i, int j, int k> struct NineteenPt {};
790  template<> struct NineteenPt< 0, 0, 0> { enum { idx = 0 }; };
791 
792  template<> struct NineteenPt< 1, 0, 0> { enum { idx = 1 }; };
793  template<> struct NineteenPt< 0, 1, 0> { enum { idx = 2 }; };
794  template<> struct NineteenPt< 0, 0, 1> { enum { idx = 3 }; };
795 
796  template<> struct NineteenPt<-1, 0, 0> { enum { idx = 4 }; };
797  template<> struct NineteenPt< 0,-1, 0> { enum { idx = 5 }; };
798  template<> struct NineteenPt< 0, 0,-1> { enum { idx = 6 }; };
799 
800  template<> struct NineteenPt< 2, 0, 0> { enum { idx = 7 }; };
801  template<> struct NineteenPt< 0, 2, 0> { enum { idx = 8 }; };
802  template<> struct NineteenPt< 0, 0, 2> { enum { idx = 9 }; };
803 
804  template<> struct NineteenPt<-2, 0, 0> { enum { idx = 10 }; };
805  template<> struct NineteenPt< 0,-2, 0> { enum { idx = 11 }; };
806  template<> struct NineteenPt< 0, 0,-2> { enum { idx = 12 }; };
807 
808  template<> struct NineteenPt< 3, 0, 0> { enum { idx = 13 }; };
809  template<> struct NineteenPt< 0, 3, 0> { enum { idx = 14 }; };
810  template<> struct NineteenPt< 0, 0, 3> { enum { idx = 15 }; };
811 
812  template<> struct NineteenPt<-3, 0, 0> { enum { idx = 16 }; };
813  template<> struct NineteenPt< 0,-3, 0> { enum { idx = 17 }; };
814  template<> struct NineteenPt< 0, 0,-3> { enum { idx = 18 }; };
815 
816 }
817 
818 
819 template<typename GridT, bool IsSafe = true>
821  : public BaseStencil<NineteenPointStencil<GridT, IsSafe>, GridT, IsSafe>
822 {
824  typedef BaseStencil<SelfT, GridT, IsSafe > BaseType;
825 public:
826  typedef GridT GridType;
827  typedef typename GridT::TreeType TreeType;
828  typedef typename GridType::ValueType ValueType;
829 
830  static const int SIZE = 19;
831 
832  NineteenPointStencil(const GridType& grid): BaseType(grid, SIZE) {}
833 
835  template<int i, int j, int k>
836  unsigned int pos() const { return NineteenPt<i,j,k>::idx; }
837 
838 private:
839  inline void init(const Coord& ijk)
840  {
841  mValues[NineteenPt< 3, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0, 0));
842  mValues[NineteenPt< 2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 0));
843  mValues[NineteenPt< 1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
844  mValues[NineteenPt<-1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
845  mValues[NineteenPt<-2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 0));
846  mValues[NineteenPt<-3, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0, 0));
847 
848  mValues[NineteenPt< 0, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3, 0));
849  mValues[NineteenPt< 0, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 0));
850  mValues[NineteenPt< 0, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
851  mValues[NineteenPt< 0,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
852  mValues[NineteenPt< 0,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -2, 0));
853  mValues[NineteenPt< 0,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -3, 0));
854 
855  mValues[NineteenPt< 0, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 3));
856  mValues[NineteenPt< 0, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 2));
857  mValues[NineteenPt< 0, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
858  mValues[NineteenPt< 0, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
859  mValues[NineteenPt< 0, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -2));
860  mValues[NineteenPt< 0, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -3));
861  }
862 
863  template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
864  using BaseType::mAcc;
865  using BaseType::mValues;
866 };// NineteenPointStencil class
867 
868 
870 
871 
872 namespace { // anonymous namespace for stencil-layout map
873 
874  // the 4th-order dense point stencil
875  template<int i, int j, int k> struct SixthDensePt { };
876  template<> struct SixthDensePt< 0, 0, 0> { enum { idx = 0 }; };
877 
878  template<> struct SixthDensePt<-3, 3, 0> { enum { idx = 1 }; };
879  template<> struct SixthDensePt<-2, 3, 0> { enum { idx = 2 }; };
880  template<> struct SixthDensePt<-1, 3, 0> { enum { idx = 3 }; };
881  template<> struct SixthDensePt< 0, 3, 0> { enum { idx = 4 }; };
882  template<> struct SixthDensePt< 1, 3, 0> { enum { idx = 5 }; };
883  template<> struct SixthDensePt< 2, 3, 0> { enum { idx = 6 }; };
884  template<> struct SixthDensePt< 3, 3, 0> { enum { idx = 7 }; };
885 
886  template<> struct SixthDensePt<-3, 2, 0> { enum { idx = 8 }; };
887  template<> struct SixthDensePt<-2, 2, 0> { enum { idx = 9 }; };
888  template<> struct SixthDensePt<-1, 2, 0> { enum { idx = 10 }; };
889  template<> struct SixthDensePt< 0, 2, 0> { enum { idx = 11 }; };
890  template<> struct SixthDensePt< 1, 2, 0> { enum { idx = 12 }; };
891  template<> struct SixthDensePt< 2, 2, 0> { enum { idx = 13 }; };
892  template<> struct SixthDensePt< 3, 2, 0> { enum { idx = 14 }; };
893 
894  template<> struct SixthDensePt<-3, 1, 0> { enum { idx = 15 }; };
895  template<> struct SixthDensePt<-2, 1, 0> { enum { idx = 16 }; };
896  template<> struct SixthDensePt<-1, 1, 0> { enum { idx = 17 }; };
897  template<> struct SixthDensePt< 0, 1, 0> { enum { idx = 18 }; };
898  template<> struct SixthDensePt< 1, 1, 0> { enum { idx = 19 }; };
899  template<> struct SixthDensePt< 2, 1, 0> { enum { idx = 20 }; };
900  template<> struct SixthDensePt< 3, 1, 0> { enum { idx = 21 }; };
901 
902  template<> struct SixthDensePt<-3, 0, 0> { enum { idx = 22 }; };
903  template<> struct SixthDensePt<-2, 0, 0> { enum { idx = 23 }; };
904  template<> struct SixthDensePt<-1, 0, 0> { enum { idx = 24 }; };
905  template<> struct SixthDensePt< 1, 0, 0> { enum { idx = 25 }; };
906  template<> struct SixthDensePt< 2, 0, 0> { enum { idx = 26 }; };
907  template<> struct SixthDensePt< 3, 0, 0> { enum { idx = 27 }; };
908 
909 
910  template<> struct SixthDensePt<-3,-1, 0> { enum { idx = 28 }; };
911  template<> struct SixthDensePt<-2,-1, 0> { enum { idx = 29 }; };
912  template<> struct SixthDensePt<-1,-1, 0> { enum { idx = 30 }; };
913  template<> struct SixthDensePt< 0,-1, 0> { enum { idx = 31 }; };
914  template<> struct SixthDensePt< 1,-1, 0> { enum { idx = 32 }; };
915  template<> struct SixthDensePt< 2,-1, 0> { enum { idx = 33 }; };
916  template<> struct SixthDensePt< 3,-1, 0> { enum { idx = 34 }; };
917 
918 
919  template<> struct SixthDensePt<-3,-2, 0> { enum { idx = 35 }; };
920  template<> struct SixthDensePt<-2,-2, 0> { enum { idx = 36 }; };
921  template<> struct SixthDensePt<-1,-2, 0> { enum { idx = 37 }; };
922  template<> struct SixthDensePt< 0,-2, 0> { enum { idx = 38 }; };
923  template<> struct SixthDensePt< 1,-2, 0> { enum { idx = 39 }; };
924  template<> struct SixthDensePt< 2,-2, 0> { enum { idx = 40 }; };
925  template<> struct SixthDensePt< 3,-2, 0> { enum { idx = 41 }; };
926 
927 
928  template<> struct SixthDensePt<-3,-3, 0> { enum { idx = 42 }; };
929  template<> struct SixthDensePt<-2,-3, 0> { enum { idx = 43 }; };
930  template<> struct SixthDensePt<-1,-3, 0> { enum { idx = 44 }; };
931  template<> struct SixthDensePt< 0,-3, 0> { enum { idx = 45 }; };
932  template<> struct SixthDensePt< 1,-3, 0> { enum { idx = 46 }; };
933  template<> struct SixthDensePt< 2,-3, 0> { enum { idx = 47 }; };
934  template<> struct SixthDensePt< 3,-3, 0> { enum { idx = 48 }; };
935 
936 
937  template<> struct SixthDensePt<-3, 0, 3> { enum { idx = 49 }; };
938  template<> struct SixthDensePt<-2, 0, 3> { enum { idx = 50 }; };
939  template<> struct SixthDensePt<-1, 0, 3> { enum { idx = 51 }; };
940  template<> struct SixthDensePt< 0, 0, 3> { enum { idx = 52 }; };
941  template<> struct SixthDensePt< 1, 0, 3> { enum { idx = 53 }; };
942  template<> struct SixthDensePt< 2, 0, 3> { enum { idx = 54 }; };
943  template<> struct SixthDensePt< 3, 0, 3> { enum { idx = 55 }; };
944 
945 
946  template<> struct SixthDensePt<-3, 0, 2> { enum { idx = 56 }; };
947  template<> struct SixthDensePt<-2, 0, 2> { enum { idx = 57 }; };
948  template<> struct SixthDensePt<-1, 0, 2> { enum { idx = 58 }; };
949  template<> struct SixthDensePt< 0, 0, 2> { enum { idx = 59 }; };
950  template<> struct SixthDensePt< 1, 0, 2> { enum { idx = 60 }; };
951  template<> struct SixthDensePt< 2, 0, 2> { enum { idx = 61 }; };
952  template<> struct SixthDensePt< 3, 0, 2> { enum { idx = 62 }; };
953 
954  template<> struct SixthDensePt<-3, 0, 1> { enum { idx = 63 }; };
955  template<> struct SixthDensePt<-2, 0, 1> { enum { idx = 64 }; };
956  template<> struct SixthDensePt<-1, 0, 1> { enum { idx = 65 }; };
957  template<> struct SixthDensePt< 0, 0, 1> { enum { idx = 66 }; };
958  template<> struct SixthDensePt< 1, 0, 1> { enum { idx = 67 }; };
959  template<> struct SixthDensePt< 2, 0, 1> { enum { idx = 68 }; };
960  template<> struct SixthDensePt< 3, 0, 1> { enum { idx = 69 }; };
961 
962 
963  template<> struct SixthDensePt<-3, 0,-1> { enum { idx = 70 }; };
964  template<> struct SixthDensePt<-2, 0,-1> { enum { idx = 71 }; };
965  template<> struct SixthDensePt<-1, 0,-1> { enum { idx = 72 }; };
966  template<> struct SixthDensePt< 0, 0,-1> { enum { idx = 73 }; };
967  template<> struct SixthDensePt< 1, 0,-1> { enum { idx = 74 }; };
968  template<> struct SixthDensePt< 2, 0,-1> { enum { idx = 75 }; };
969  template<> struct SixthDensePt< 3, 0,-1> { enum { idx = 76 }; };
970 
971 
972  template<> struct SixthDensePt<-3, 0,-2> { enum { idx = 77 }; };
973  template<> struct SixthDensePt<-2, 0,-2> { enum { idx = 78 }; };
974  template<> struct SixthDensePt<-1, 0,-2> { enum { idx = 79 }; };
975  template<> struct SixthDensePt< 0, 0,-2> { enum { idx = 80 }; };
976  template<> struct SixthDensePt< 1, 0,-2> { enum { idx = 81 }; };
977  template<> struct SixthDensePt< 2, 0,-2> { enum { idx = 82 }; };
978  template<> struct SixthDensePt< 3, 0,-2> { enum { idx = 83 }; };
979 
980 
981  template<> struct SixthDensePt<-3, 0,-3> { enum { idx = 84 }; };
982  template<> struct SixthDensePt<-2, 0,-3> { enum { idx = 85 }; };
983  template<> struct SixthDensePt<-1, 0,-3> { enum { idx = 86 }; };
984  template<> struct SixthDensePt< 0, 0,-3> { enum { idx = 87 }; };
985  template<> struct SixthDensePt< 1, 0,-3> { enum { idx = 88 }; };
986  template<> struct SixthDensePt< 2, 0,-3> { enum { idx = 89 }; };
987  template<> struct SixthDensePt< 3, 0,-3> { enum { idx = 90 }; };
988 
989 
990  template<> struct SixthDensePt< 0,-3, 3> { enum { idx = 91 }; };
991  template<> struct SixthDensePt< 0,-2, 3> { enum { idx = 92 }; };
992  template<> struct SixthDensePt< 0,-1, 3> { enum { idx = 93 }; };
993  template<> struct SixthDensePt< 0, 1, 3> { enum { idx = 94 }; };
994  template<> struct SixthDensePt< 0, 2, 3> { enum { idx = 95 }; };
995  template<> struct SixthDensePt< 0, 3, 3> { enum { idx = 96 }; };
996 
997  template<> struct SixthDensePt< 0,-3, 2> { enum { idx = 97 }; };
998  template<> struct SixthDensePt< 0,-2, 2> { enum { idx = 98 }; };
999  template<> struct SixthDensePt< 0,-1, 2> { enum { idx = 99 }; };
1000  template<> struct SixthDensePt< 0, 1, 2> { enum { idx = 100 }; };
1001  template<> struct SixthDensePt< 0, 2, 2> { enum { idx = 101 }; };
1002  template<> struct SixthDensePt< 0, 3, 2> { enum { idx = 102 }; };
1003 
1004  template<> struct SixthDensePt< 0,-3, 1> { enum { idx = 103 }; };
1005  template<> struct SixthDensePt< 0,-2, 1> { enum { idx = 104 }; };
1006  template<> struct SixthDensePt< 0,-1, 1> { enum { idx = 105 }; };
1007  template<> struct SixthDensePt< 0, 1, 1> { enum { idx = 106 }; };
1008  template<> struct SixthDensePt< 0, 2, 1> { enum { idx = 107 }; };
1009  template<> struct SixthDensePt< 0, 3, 1> { enum { idx = 108 }; };
1010 
1011  template<> struct SixthDensePt< 0,-3,-1> { enum { idx = 109 }; };
1012  template<> struct SixthDensePt< 0,-2,-1> { enum { idx = 110 }; };
1013  template<> struct SixthDensePt< 0,-1,-1> { enum { idx = 111 }; };
1014  template<> struct SixthDensePt< 0, 1,-1> { enum { idx = 112 }; };
1015  template<> struct SixthDensePt< 0, 2,-1> { enum { idx = 113 }; };
1016  template<> struct SixthDensePt< 0, 3,-1> { enum { idx = 114 }; };
1017 
1018  template<> struct SixthDensePt< 0,-3,-2> { enum { idx = 115 }; };
1019  template<> struct SixthDensePt< 0,-2,-2> { enum { idx = 116 }; };
1020  template<> struct SixthDensePt< 0,-1,-2> { enum { idx = 117 }; };
1021  template<> struct SixthDensePt< 0, 1,-2> { enum { idx = 118 }; };
1022  template<> struct SixthDensePt< 0, 2,-2> { enum { idx = 119 }; };
1023  template<> struct SixthDensePt< 0, 3,-2> { enum { idx = 120 }; };
1024 
1025  template<> struct SixthDensePt< 0,-3,-3> { enum { idx = 121 }; };
1026  template<> struct SixthDensePt< 0,-2,-3> { enum { idx = 122 }; };
1027  template<> struct SixthDensePt< 0,-1,-3> { enum { idx = 123 }; };
1028  template<> struct SixthDensePt< 0, 1,-3> { enum { idx = 124 }; };
1029  template<> struct SixthDensePt< 0, 2,-3> { enum { idx = 125 }; };
1030  template<> struct SixthDensePt< 0, 3,-3> { enum { idx = 126 }; };
1031 
1032 }
1033 
1034 
1035 template<typename GridT, bool IsSafe = true>
1037  : public BaseStencil<SixthOrderDenseStencil<GridT, IsSafe>, GridT, IsSafe>
1038 {
1040  typedef BaseStencil<SelfT, GridT, IsSafe > BaseType;
1041 public:
1042  typedef GridT GridType;
1043  typedef typename GridT::TreeType TreeType;
1044  typedef typename GridType::ValueType ValueType;
1045 
1046  static const int SIZE = 127;
1047 
1048  SixthOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
1049 
1051  template<int i, int j, int k>
1052  unsigned int pos() const { return SixthDensePt<i,j,k>::idx; }
1053 
1054 private:
1055  inline void init(const Coord& ijk)
1056  {
1057  mValues[SixthDensePt<-3, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3, 3, 0));
1058  mValues[SixthDensePt<-2, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 3, 0));
1059  mValues[SixthDensePt<-1, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 3, 0));
1060  mValues[SixthDensePt< 0, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3, 0));
1061  mValues[SixthDensePt< 1, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 3, 0));
1062  mValues[SixthDensePt< 2, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 3, 0));
1063  mValues[SixthDensePt< 3, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3, 3, 0));
1064 
1065  mValues[SixthDensePt<-3, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3, 2, 0));
1066  mValues[SixthDensePt<-2, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 2, 0));
1067  mValues[SixthDensePt<-1, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 2, 0));
1068  mValues[SixthDensePt< 0, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 0));
1069  mValues[SixthDensePt< 1, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 2, 0));
1070  mValues[SixthDensePt< 2, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 2, 0));
1071  mValues[SixthDensePt< 3, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3, 2, 0));
1072 
1073  mValues[SixthDensePt<-3, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3, 1, 0));
1074  mValues[SixthDensePt<-2, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 1, 0));
1075  mValues[SixthDensePt<-1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 1, 0));
1076  mValues[SixthDensePt< 0, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
1077  mValues[SixthDensePt< 1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 1, 0));
1078  mValues[SixthDensePt< 2, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 1, 0));
1079  mValues[SixthDensePt< 3, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3, 1, 0));
1080 
1081  mValues[SixthDensePt<-3, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0, 0));
1082  mValues[SixthDensePt<-2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 0));
1083  mValues[SixthDensePt<-1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
1084  mValues[SixthDensePt< 1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
1085  mValues[SixthDensePt< 2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 0));
1086  mValues[SixthDensePt< 3, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0, 0));
1087 
1088  mValues[SixthDensePt<-3,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3,-1, 0));
1089  mValues[SixthDensePt<-2,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2,-1, 0));
1090  mValues[SixthDensePt<-1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1,-1, 0));
1091  mValues[SixthDensePt< 0,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 0));
1092  mValues[SixthDensePt< 1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1,-1, 0));
1093  mValues[SixthDensePt< 2,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2,-1, 0));
1094  mValues[SixthDensePt< 3,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3,-1, 0));
1095 
1096  mValues[SixthDensePt<-3,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3,-2, 0));
1097  mValues[SixthDensePt<-2,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2,-2, 0));
1098  mValues[SixthDensePt<-1,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1,-2, 0));
1099  mValues[SixthDensePt< 0,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 0));
1100  mValues[SixthDensePt< 1,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1,-2, 0));
1101  mValues[SixthDensePt< 2,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2,-2, 0));
1102  mValues[SixthDensePt< 3,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3,-2, 0));
1103 
1104  mValues[SixthDensePt<-3,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3,-3, 0));
1105  mValues[SixthDensePt<-2,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2,-3, 0));
1106  mValues[SixthDensePt<-1,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1,-3, 0));
1107  mValues[SixthDensePt< 0,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3, 0));
1108  mValues[SixthDensePt< 1,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1,-3, 0));
1109  mValues[SixthDensePt< 2,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2,-3, 0));
1110  mValues[SixthDensePt< 3,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3,-3, 0));
1111 
1112  mValues[SixthDensePt<-3, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0, 3));
1113  mValues[SixthDensePt<-2, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 3));
1114  mValues[SixthDensePt<-1, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 3));
1115  mValues[SixthDensePt< 0, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 3));
1116  mValues[SixthDensePt< 1, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 3));
1117  mValues[SixthDensePt< 2, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 3));
1118  mValues[SixthDensePt< 3, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0, 3));
1119 
1120  mValues[SixthDensePt<-3, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0, 2));
1121  mValues[SixthDensePt<-2, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 2));
1122  mValues[SixthDensePt<-1, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 2));
1123  mValues[SixthDensePt< 0, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 2));
1124  mValues[SixthDensePt< 1, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 2));
1125  mValues[SixthDensePt< 2, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 2));
1126  mValues[SixthDensePt< 3, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0, 2));
1127 
1128  mValues[SixthDensePt<-3, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0, 1));
1129  mValues[SixthDensePt<-2, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 1));
1130  mValues[SixthDensePt<-1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 1));
1131  mValues[SixthDensePt< 0, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
1132  mValues[SixthDensePt< 1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 1));
1133  mValues[SixthDensePt< 2, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 1));
1134  mValues[SixthDensePt< 3, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0, 1));
1135 
1136  mValues[SixthDensePt<-3, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0,-1));
1137  mValues[SixthDensePt<-2, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0,-1));
1138  mValues[SixthDensePt<-1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0,-1));
1139  mValues[SixthDensePt< 0, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0,-1));
1140  mValues[SixthDensePt< 1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0,-1));
1141  mValues[SixthDensePt< 2, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0,-1));
1142  mValues[SixthDensePt< 3, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0,-1));
1143 
1144  mValues[SixthDensePt<-3, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0,-2));
1145  mValues[SixthDensePt<-2, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0,-2));
1146  mValues[SixthDensePt<-1, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0,-2));
1147  mValues[SixthDensePt< 0, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0,-2));
1148  mValues[SixthDensePt< 1, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0,-2));
1149  mValues[SixthDensePt< 2, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0,-2));
1150  mValues[SixthDensePt< 3, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0,-2));
1151 
1152  mValues[SixthDensePt<-3, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0,-3));
1153  mValues[SixthDensePt<-2, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0,-3));
1154  mValues[SixthDensePt<-1, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0,-3));
1155  mValues[SixthDensePt< 0, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0,-3));
1156  mValues[SixthDensePt< 1, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0,-3));
1157  mValues[SixthDensePt< 2, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0,-3));
1158  mValues[SixthDensePt< 3, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0,-3));
1159 
1160  mValues[SixthDensePt< 0,-3, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3, 3));
1161  mValues[SixthDensePt< 0,-2, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 3));
1162  mValues[SixthDensePt< 0,-1, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 3));
1163  mValues[SixthDensePt< 0, 1, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 3));
1164  mValues[SixthDensePt< 0, 2, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 3));
1165  mValues[SixthDensePt< 0, 3, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3, 3));
1166 
1167  mValues[SixthDensePt< 0,-3, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3, 2));
1168  mValues[SixthDensePt< 0,-2, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 2));
1169  mValues[SixthDensePt< 0,-1, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 2));
1170  mValues[SixthDensePt< 0, 1, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 2));
1171  mValues[SixthDensePt< 0, 2, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 2));
1172  mValues[SixthDensePt< 0, 3, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3, 2));
1173 
1174  mValues[SixthDensePt< 0,-3, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3, 1));
1175  mValues[SixthDensePt< 0,-2, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 1));
1176  mValues[SixthDensePt< 0,-1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 1));
1177  mValues[SixthDensePt< 0, 1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 1));
1178  mValues[SixthDensePt< 0, 2, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 1));
1179  mValues[SixthDensePt< 0, 3, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3, 1));
1180 
1181  mValues[SixthDensePt< 0,-3,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3,-1));
1182  mValues[SixthDensePt< 0,-2,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2,-1));
1183  mValues[SixthDensePt< 0,-1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1,-1));
1184  mValues[SixthDensePt< 0, 1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1,-1));
1185  mValues[SixthDensePt< 0, 2,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2,-1));
1186  mValues[SixthDensePt< 0, 3,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3,-1));
1187 
1188  mValues[SixthDensePt< 0,-3,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3,-2));
1189  mValues[SixthDensePt< 0,-2,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2,-2));
1190  mValues[SixthDensePt< 0,-1,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1,-2));
1191  mValues[SixthDensePt< 0, 1,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1,-2));
1192  mValues[SixthDensePt< 0, 2,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2,-2));
1193  mValues[SixthDensePt< 0, 3,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3,-2));
1194 
1195  mValues[SixthDensePt< 0,-3,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3,-3));
1196  mValues[SixthDensePt< 0,-2,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2,-3));
1197  mValues[SixthDensePt< 0,-1,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1,-3));
1198  mValues[SixthDensePt< 0, 1,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1,-3));
1199  mValues[SixthDensePt< 0, 2,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2,-3));
1200  mValues[SixthDensePt< 0, 3,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3,-3));
1201  }
1202 
1203  template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
1204  using BaseType::mAcc;
1205  using BaseType::mValues;
1206 };// SixthOrderDenseStencil class
1207 
1208 
1210 
1211 namespace { // anonymous namespace for stencil-layout map
1212 
1213  // the seven point stencil with a different layout from SevenPt
1214  template<int i, int j, int k> struct GradPt {};
1215  template<> struct GradPt< 0, 0, 0> { enum { idx = 0 }; };
1216  template<> struct GradPt< 1, 0, 0> { enum { idx = 2 }; };
1217  template<> struct GradPt< 0, 1, 0> { enum { idx = 4 }; };
1218  template<> struct GradPt< 0, 0, 1> { enum { idx = 6 }; };
1219  template<> struct GradPt<-1, 0, 0> { enum { idx = 1 }; };
1220  template<> struct GradPt< 0,-1, 0> { enum { idx = 3 }; };
1221  template<> struct GradPt< 0, 0,-1> { enum { idx = 5 }; };
1222 }
1223 
1230 template<typename GridT, bool IsSafe = true>
1231 class GradStencil : public BaseStencil<GradStencil<GridT, IsSafe>, GridT, IsSafe>
1232 {
1233  typedef GradStencil<GridT, IsSafe> SelfT;
1234  typedef BaseStencil<SelfT, GridT, IsSafe > BaseType;
1235 public:
1236  typedef GridT GridType;
1237  typedef typename GridT::TreeType TreeType;
1238  typedef typename GridType::ValueType ValueType;
1239 
1240  static const int SIZE = 7;
1241 
1242  GradStencil(const GridType& grid)
1243  : BaseType(grid, SIZE)
1244  , mInv2Dx(ValueType(0.5 / grid.voxelSize()[0]))
1245  , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1246  {
1247  }
1248 
1249  GradStencil(const GridType& grid, Real dx)
1250  : BaseType(grid, SIZE)
1251  , mInv2Dx(ValueType(0.5 / dx))
1252  , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1253  {
1254  }
1255 
1261  inline ValueType normSqGrad() const
1262  {
1263  return mInvDx2 * math::GodunovsNormSqrd(mValues[0] > zeroVal<ValueType>(),
1264  mValues[0] - mValues[1],
1265  mValues[2] - mValues[0],
1266  mValues[0] - mValues[3],
1267  mValues[4] - mValues[0],
1268  mValues[0] - mValues[5],
1269  mValues[6] - mValues[0]);
1270  }
1271 
1278  {
1279  return math::Vec3<ValueType>(mValues[2] - mValues[1],
1280  mValues[4] - mValues[3],
1281  mValues[6] - mValues[5])*mInv2Dx;
1282  }
1288  {
1289  return math::Vec3<ValueType>(
1290  V[0]>0 ? mValues[0] - mValues[1] : mValues[2] - mValues[0],
1291  V[1]>0 ? mValues[0] - mValues[3] : mValues[4] - mValues[0],
1292  V[2]>0 ? mValues[0] - mValues[5] : mValues[6] - mValues[0])*2*mInv2Dx;
1293  }
1294 
1297  inline ValueType laplacian() const
1298  {
1299  return mInvDx2 * (mValues[1] + mValues[2] +
1300  mValues[3] + mValues[4] +
1301  mValues[5] + mValues[6] - 6*mValues[0]);
1302  }
1303 
1306  inline bool zeroCrossing() const
1307  {
1308  const typename BaseType::BufferType& v = mValues;
1309  return (v[0]>0 ? (v[1]<0 || v[2]<0 || v[3]<0 || v[4]<0 || v[5]<0 || v[6]<0)
1310  : (v[1]>0 || v[2]>0 || v[3]>0 || v[4]>0 || v[5]>0 || v[6]>0));
1311  }
1312 
1321  {
1322  const Coord& ijk = BaseType::getCenterCoord();
1323  const ValueType d = ValueType(mValues[0] * 0.5 * mInvDx2); // distance in voxels / (2dx^2)
1325  const auto value = math::Vec3<ValueType>( ijk[0] - d*(mValues[2] - mValues[1]),
1326  ijk[1] - d*(mValues[4] - mValues[3]),
1327  ijk[2] - d*(mValues[6] - mValues[5]));
1329  return value;
1330  }
1331 
1333  template<int i, int j, int k>
1334  unsigned int pos() const { return GradPt<i,j,k>::idx; }
1335 
1336 private:
1337 
1338  inline void init(const Coord& ijk)
1339  {
1340  BaseType::template setValue<-1, 0, 0>(mAcc.getValue(ijk.offsetBy(-1, 0, 0)));
1341  BaseType::template setValue< 1, 0, 0>(mAcc.getValue(ijk.offsetBy( 1, 0, 0)));
1342 
1343  BaseType::template setValue< 0,-1, 0>(mAcc.getValue(ijk.offsetBy( 0,-1, 0)));
1344  BaseType::template setValue< 0, 1, 0>(mAcc.getValue(ijk.offsetBy( 0, 1, 0)));
1345 
1346  BaseType::template setValue< 0, 0,-1>(mAcc.getValue(ijk.offsetBy( 0, 0,-1)));
1347  BaseType::template setValue< 0, 0, 1>(mAcc.getValue(ijk.offsetBy( 0, 0, 1)));
1348  }
1349 
1350  template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
1351  using BaseType::mAcc;
1352  using BaseType::mValues;
1353  const ValueType mInv2Dx, mInvDx2;
1354 }; // GradStencil class
1355 
1357 
1358 
1364 template<typename GridT, bool IsSafe = true>
1365 class WenoStencil: public BaseStencil<WenoStencil<GridT, IsSafe>, GridT, IsSafe>
1366 {
1367  typedef WenoStencil<GridT, IsSafe> SelfT;
1368  typedef BaseStencil<SelfT, GridT, IsSafe > BaseType;
1369 public:
1370  typedef GridT GridType;
1371  typedef typename GridT::TreeType TreeType;
1372  typedef typename GridType::ValueType ValueType;
1373 
1374  static const int SIZE = 19;
1375 
1376  WenoStencil(const GridType& grid)
1377  : BaseType(grid, SIZE)
1378  , mDx2(ValueType(math::Pow2(grid.voxelSize()[0])))
1379  , mInv2Dx(ValueType(0.5 / grid.voxelSize()[0]))
1380  , mInvDx2(ValueType(1.0 / mDx2))
1381  {
1382  }
1383 
1384  WenoStencil(const GridType& grid, Real dx)
1385  : BaseType(grid, SIZE)
1386  , mDx2(ValueType(dx * dx))
1387  , mInv2Dx(ValueType(0.5 / dx))
1388  , mInvDx2(ValueType(1.0 / mDx2))
1389  {
1390  }
1391 
1397  inline ValueType normSqGrad(const ValueType &isoValue = zeroVal<ValueType>()) const
1398  {
1399  const typename BaseType::BufferType& v = mValues;
1400 #ifdef DWA_OPENVDB
1401  // SSE optimized
1402  const simd::Float4
1403  v1(v[2]-v[1], v[ 8]-v[ 7], v[14]-v[13], 0),
1404  v2(v[3]-v[2], v[ 9]-v[ 8], v[15]-v[14], 0),
1405  v3(v[0]-v[3], v[ 0]-v[ 9], v[ 0]-v[15], 0),
1406  v4(v[4]-v[0], v[10]-v[ 0], v[16]-v[ 0], 0),
1407  v5(v[5]-v[4], v[11]-v[10], v[17]-v[16], 0),
1408  v6(v[6]-v[5], v[12]-v[11], v[18]-v[17], 0),
1409  dP_m = math::WENO5(v1, v2, v3, v4, v5, mDx2),
1410  dP_p = math::WENO5(v6, v5, v4, v3, v2, mDx2);
1411 
1412  return mInvDx2 * math::GodunovsNormSqrd(mValues[0] > isoValue, dP_m, dP_p);
1413 #else
1414  const Real
1415  dP_xm = math::WENO5(v[ 2]-v[ 1],v[ 3]-v[ 2],v[ 0]-v[ 3],v[ 4]-v[ 0],v[ 5]-v[ 4],mDx2),
1416  dP_xp = math::WENO5(v[ 6]-v[ 5],v[ 5]-v[ 4],v[ 4]-v[ 0],v[ 0]-v[ 3],v[ 3]-v[ 2],mDx2),
1417  dP_ym = math::WENO5(v[ 8]-v[ 7],v[ 9]-v[ 8],v[ 0]-v[ 9],v[10]-v[ 0],v[11]-v[10],mDx2),
1418  dP_yp = math::WENO5(v[12]-v[11],v[11]-v[10],v[10]-v[ 0],v[ 0]-v[ 9],v[ 9]-v[ 8],mDx2),
1419  dP_zm = math::WENO5(v[14]-v[13],v[15]-v[14],v[ 0]-v[15],v[16]-v[ 0],v[17]-v[16],mDx2),
1420  dP_zp = math::WENO5(v[18]-v[17],v[17]-v[16],v[16]-v[ 0],v[ 0]-v[15],v[15]-v[14],mDx2);
1421  return static_cast<ValueType>(
1422  mInvDx2*math::GodunovsNormSqrd(v[0]>isoValue, dP_xm, dP_xp, dP_ym, dP_yp, dP_zm, dP_zp));
1423 #endif
1424  }
1425 
1432  {
1433  const typename BaseType::BufferType& v = mValues;
1434  return 2*mInv2Dx * math::Vec3<ValueType>(
1435  V[0]>0 ? math::WENO5(v[ 2]-v[ 1],v[ 3]-v[ 2],v[ 0]-v[ 3], v[ 4]-v[ 0],v[ 5]-v[ 4],mDx2)
1436  : math::WENO5(v[ 6]-v[ 5],v[ 5]-v[ 4],v[ 4]-v[ 0], v[ 0]-v[ 3],v[ 3]-v[ 2],mDx2),
1437  V[1]>0 ? math::WENO5(v[ 8]-v[ 7],v[ 9]-v[ 8],v[ 0]-v[ 9], v[10]-v[ 0],v[11]-v[10],mDx2)
1438  : math::WENO5(v[12]-v[11],v[11]-v[10],v[10]-v[ 0], v[ 0]-v[ 9],v[ 9]-v[ 8],mDx2),
1439  V[2]>0 ? math::WENO5(v[14]-v[13],v[15]-v[14],v[ 0]-v[15], v[16]-v[ 0],v[17]-v[16],mDx2)
1440  : math::WENO5(v[18]-v[17],v[17]-v[16],v[16]-v[ 0], v[ 0]-v[15],v[15]-v[14],mDx2));
1441  }
1448  {
1449  return mInv2Dx * math::Vec3<ValueType>(mValues[ 4] - mValues[ 3],
1450  mValues[10] - mValues[ 9],
1451  mValues[16] - mValues[15]);
1452  }
1453 
1459  inline ValueType laplacian() const
1460  {
1461  return mInvDx2 * (
1462  mValues[ 3] + mValues[ 4] +
1463  mValues[ 9] + mValues[10] +
1464  mValues[15] + mValues[16] - 6*mValues[0]);
1465  }
1466 
1469  inline bool zeroCrossing() const
1470  {
1471  const typename BaseType::BufferType& v = mValues;
1472  return (v[ 0]>0 ? (v[ 3]<0 || v[ 4]<0 || v[ 9]<0 || v[10]<0 || v[15]<0 || v[16]<0)
1473  : (v[ 3]>0 || v[ 4]>0 || v[ 9]>0 || v[10]>0 || v[15]>0 || v[16]>0));
1474  }
1475 
1476 private:
1477  inline void init(const Coord& ijk)
1478  {
1479  mValues[ 1] = mAcc.getValue(ijk.offsetBy(-3, 0, 0));
1480  mValues[ 2] = mAcc.getValue(ijk.offsetBy(-2, 0, 0));
1481  mValues[ 3] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
1482  mValues[ 4] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
1483  mValues[ 5] = mAcc.getValue(ijk.offsetBy( 2, 0, 0));
1484  mValues[ 6] = mAcc.getValue(ijk.offsetBy( 3, 0, 0));
1485 
1486  mValues[ 7] = mAcc.getValue(ijk.offsetBy( 0, -3, 0));
1487  mValues[ 8] = mAcc.getValue(ijk.offsetBy( 0, -2, 0));
1488  mValues[ 9] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
1489  mValues[10] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
1490  mValues[11] = mAcc.getValue(ijk.offsetBy( 0, 2, 0));
1491  mValues[12] = mAcc.getValue(ijk.offsetBy( 0, 3, 0));
1492 
1493  mValues[13] = mAcc.getValue(ijk.offsetBy( 0, 0, -3));
1494  mValues[14] = mAcc.getValue(ijk.offsetBy( 0, 0, -2));
1495  mValues[15] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
1496  mValues[16] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
1497  mValues[17] = mAcc.getValue(ijk.offsetBy( 0, 0, 2));
1498  mValues[18] = mAcc.getValue(ijk.offsetBy( 0, 0, 3));
1499  }
1500 
1501  template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
1502  using BaseType::mAcc;
1503  using BaseType::mValues;
1504  const ValueType mDx2, mInv2Dx, mInvDx2;
1505 }; // WenoStencil class
1506 
1507 
1509 
1510 
1511 template<typename GridT, bool IsSafe = true>
1512 class CurvatureStencil: public BaseStencil<CurvatureStencil<GridT, IsSafe>, GridT, IsSafe>
1513 {
1514  typedef CurvatureStencil<GridT, IsSafe> SelfT;
1515  typedef BaseStencil<SelfT, GridT, IsSafe> BaseType;
1516 public:
1517  typedef GridT GridType;
1518  typedef typename GridT::TreeType TreeType;
1519  typedef typename GridT::ValueType ValueType;
1520 
1521  static const int SIZE = 19;
1522 
1523  CurvatureStencil(const GridType& grid)
1524  : BaseType(grid, SIZE)
1525  , mInv2Dx(ValueType(0.5 / grid.voxelSize()[0]))
1526  , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1527  {
1528  }
1529 
1530  CurvatureStencil(const GridType& grid, Real dx)
1531  : BaseType(grid, SIZE)
1532  , mInv2Dx(ValueType(0.5 / dx))
1533  , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1534  {
1535  }
1536 
1541  inline ValueType meanCurvature() const
1542  {
1543  Real alpha, normGrad;
1544  return this->meanCurvature(alpha, normGrad) ?
1545  ValueType(alpha*mInv2Dx/math::Pow3(normGrad)) : 0;
1546  }
1547 
1552  inline ValueType gaussianCurvature() const
1553  {
1554  Real alpha, normGrad;
1555  return this->gaussianCurvature(alpha, normGrad) ?
1556  ValueType(alpha*mInvDx2/math::Pow4(normGrad)) : 0;
1557  }
1558 
1564  inline void curvatures(ValueType &mean, ValueType& gauss) const
1565  {
1566  Real alphaM, alphaG, normGrad;
1567  if (this->curvatures(alphaM, alphaG, normGrad)) {
1568  mean = ValueType(alphaM*mInv2Dx/math::Pow3(normGrad));
1569  gauss = ValueType(alphaG*mInvDx2/math::Pow4(normGrad));
1570  } else {
1571  mean = gauss = 0;
1572  }
1573  }
1574 
1581  inline ValueType meanCurvatureNormGrad() const
1582  {
1583  Real alpha, normGrad;
1584  return this->meanCurvature(alpha, normGrad) ?
1585  ValueType(alpha*mInvDx2/(2*math::Pow2(normGrad))) : 0;
1586  }
1587 
1593  inline ValueType gaussianCurvatureNormGrad() const
1594  {
1595  Real alpha, normGrad;
1596  return this->gaussianCurvature(alpha, normGrad) ?
1597  ValueType(2*alpha*mInv2Dx*mInvDx2/math::Pow3(normGrad)) : 0;
1598  }
1599 
1605  inline void curvaturesNormGrad(ValueType &mean, ValueType& gauss) const
1606  {
1607  Real alphaM, alphaG, normGrad;
1608  if (this->curvatures(alphaM, alphaG, normGrad)) {
1609  mean = ValueType(alphaM*mInvDx2/(2*math::Pow2(normGrad)));
1610  gauss = ValueType(2*alphaG*mInv2Dx*mInvDx2/math::Pow3(normGrad));
1611  } else {
1612  mean = gauss = 0;
1613  }
1614  }
1615 
1621  inline std::pair<ValueType, ValueType> principalCurvatures() const
1622  {
1623  std::pair<ValueType, ValueType> pair(0, 0);// min, max
1624  Real alphaM, alphaG, normGrad;
1625  if (this->curvatures(alphaM, alphaG, normGrad)) {
1626  const Real mean = alphaM*mInv2Dx/math::Pow3(normGrad);
1627  const Real tmp = std::sqrt(mean*mean - alphaG*mInvDx2/math::Pow4(normGrad));
1628  pair.first = ValueType(mean - tmp);
1629  pair.second = ValueType(mean + tmp);
1630  }
1631  return pair;// min, max
1632  }
1633 
1639  inline ValueType laplacian() const
1640  {
1641  return mInvDx2 * (
1642  mValues[1] + mValues[2] +
1643  mValues[3] + mValues[4] +
1644  mValues[5] + mValues[6] - 6*mValues[0]);
1645  }
1646 
1653  {
1654  return math::Vec3<ValueType>(
1655  mValues[2] - mValues[1],
1656  mValues[4] - mValues[3],
1657  mValues[6] - mValues[5])*mInv2Dx;
1658  }
1659 
1660 private:
1661  inline void init(const Coord &ijk)
1662  {
1663  mValues[ 1] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
1664  mValues[ 2] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
1665 
1666  mValues[ 3] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
1667  mValues[ 4] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
1668 
1669  mValues[ 5] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
1670  mValues[ 6] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
1671 
1672  mValues[ 7] = mAcc.getValue(ijk.offsetBy(-1, -1, 0));
1673  mValues[ 8] = mAcc.getValue(ijk.offsetBy( 1, -1, 0));
1674  mValues[ 9] = mAcc.getValue(ijk.offsetBy(-1, 1, 0));
1675  mValues[10] = mAcc.getValue(ijk.offsetBy( 1, 1, 0));
1676 
1677  mValues[11] = mAcc.getValue(ijk.offsetBy(-1, 0, -1));
1678  mValues[12] = mAcc.getValue(ijk.offsetBy( 1, 0, -1));
1679  mValues[13] = mAcc.getValue(ijk.offsetBy(-1, 0, 1));
1680  mValues[14] = mAcc.getValue(ijk.offsetBy( 1, 0, 1));
1681 
1682  mValues[15] = mAcc.getValue(ijk.offsetBy( 0, -1, -1));
1683  mValues[16] = mAcc.getValue(ijk.offsetBy( 0, 1, -1));
1684  mValues[17] = mAcc.getValue(ijk.offsetBy( 0, -1, 1));
1685  mValues[18] = mAcc.getValue(ijk.offsetBy( 0, 1, 1));
1686  }
1687 
1688  inline Real Dx() const { return 0.5*(mValues[2] - mValues[1]); }// * 1/dx
1689  inline Real Dy() const { return 0.5*(mValues[4] - mValues[3]); }// * 1/dx
1690  inline Real Dz() const { return 0.5*(mValues[6] - mValues[5]); }// * 1/dx
1691  inline Real Dxx() const { return mValues[2] - 2 * mValues[0] + mValues[1]; }// * 1/dx2
1692  inline Real Dyy() const { return mValues[4] - 2 * mValues[0] + mValues[3]; }// * 1/dx2}
1693  inline Real Dzz() const { return mValues[6] - 2 * mValues[0] + mValues[5]; }// * 1/dx2
1694  inline Real Dxy() const { return 0.25 * (mValues[10] - mValues[ 8] + mValues[ 7] - mValues[ 9]); }// * 1/dx2
1695  inline Real Dxz() const { return 0.25 * (mValues[14] - mValues[12] + mValues[11] - mValues[13]); }// * 1/dx2
1696  inline Real Dyz() const { return 0.25 * (mValues[18] - mValues[16] + mValues[15] - mValues[17]); }// * 1/dx2
1697 
1698  inline bool meanCurvature(Real& alpha, Real& normGrad) const
1699  {
1700  // For performance all finite differences are unscaled wrt dx
1701  const Real Dx = this->Dx(), Dy = this->Dy(), Dz = this->Dz(),
1702  Dx2 = Dx*Dx, Dy2 = Dy*Dy, Dz2 = Dz*Dz, normGrad2 = Dx2 + Dy2 + Dz2;
1703  if (normGrad2 <= math::Tolerance<Real>::value()) {
1704  alpha = normGrad = 0;
1705  return false;
1706  }
1707  const Real Dxx = this->Dxx(), Dyy = this->Dyy(), Dzz = this->Dzz();
1708  alpha = Dx2*(Dyy + Dzz) + Dy2*(Dxx + Dzz) + Dz2*(Dxx + Dyy) -
1709  2*(Dx*(Dy*this->Dxy() + Dz*this->Dxz()) + Dy*Dz*this->Dyz());// * 1/dx^4
1710  normGrad = std::sqrt(normGrad2); // * 1/dx
1711  return true;
1712  }
1713 
1714  inline bool gaussianCurvature(Real& alpha, Real& normGrad) const
1715  {
1716  // For performance all finite differences are unscaled wrt dx
1717  const Real Dx = this->Dx(), Dy = this->Dy(), Dz = this->Dz(),
1718  Dx2 = Dx*Dx, Dy2 = Dy*Dy, Dz2 = Dz*Dz, normGrad2 = Dx2 + Dy2 + Dz2;
1719  if (normGrad2 <= math::Tolerance<Real>::value()) {
1720  alpha = normGrad = 0;
1721  return false;
1722  }
1723  const Real Dxx = this->Dxx(), Dyy = this->Dyy(), Dzz = this->Dzz(),
1724  Dxy = this->Dxy(), Dxz = this->Dxz(), Dyz = this->Dyz();
1725  alpha = Dx2*(Dyy*Dzz - Dyz*Dyz) + Dy2*(Dxx*Dzz - Dxz*Dxz) + Dz2*(Dxx*Dyy - Dxy*Dxy) +
1726  2*( Dy*Dz*(Dxy*Dxz - Dyz*Dxx) + Dx*Dz*(Dxy*Dyz - Dxz*Dyy) + Dx*Dy*(Dxz*Dyz - Dxy*Dzz) );// * 1/dx^6
1727  normGrad = std::sqrt(normGrad2); // * 1/dx
1728  return true;
1729  }
1730  inline bool curvatures(Real& alphaM, Real& alphaG, Real& normGrad) const
1731  {
1732  // For performance all finite differences are unscaled wrt dx
1733  const Real Dx = this->Dx(), Dy = this->Dy(), Dz = this->Dz(),
1734  Dx2 = Dx*Dx, Dy2 = Dy*Dy, Dz2 = Dz*Dz, normGrad2 = Dx2 + Dy2 + Dz2;
1735  if (normGrad2 <= math::Tolerance<Real>::value()) {
1736  alphaM = alphaG =normGrad = 0;
1737  return false;
1738  }
1739  const Real Dxx = this->Dxx(), Dyy = this->Dyy(), Dzz = this->Dzz(),
1740  Dxy = this->Dxy(), Dxz = this->Dxz(), Dyz = this->Dyz();
1741  alphaM = Dx2*(Dyy + Dzz) + Dy2*(Dxx + Dzz) + Dz2*(Dxx + Dyy) -
1742  2*(Dx*(Dy*Dxy + Dz*Dxz) + Dy*Dz*Dyz);// *1/dx^4
1743  alphaG = Dx2*(Dyy*Dzz - Dyz*Dyz) + Dy2*(Dxx*Dzz - Dxz*Dxz) + Dz2*(Dxx*Dyy - Dxy*Dxy) +
1744  2*( Dy*Dz*(Dxy*Dxz - Dyz*Dxx) + Dx*Dz*(Dxy*Dyz - Dxz*Dyy) + Dx*Dy*(Dxz*Dyz - Dxy*Dzz) );// *1/dx^6
1745  normGrad = std::sqrt(normGrad2); // * 1/dx
1746  return true;
1747  }
1748 
1749  template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
1750  using BaseType::mAcc;
1751  using BaseType::mValues;
1752  const ValueType mInv2Dx, mInvDx2;
1753 }; // CurvatureStencil class
1754 
1755 
1757 
1758 
1760 template<typename GridT, bool IsSafe = true>
1761 class DenseStencil: public BaseStencil<DenseStencil<GridT, IsSafe>, GridT, IsSafe>
1762 {
1763  typedef DenseStencil<GridT, IsSafe> SelfT;
1764  typedef BaseStencil<SelfT, GridT, IsSafe> BaseType;
1765 public:
1766  typedef GridT GridType;
1767  typedef typename GridT::TreeType TreeType;
1768  typedef typename GridType::ValueType ValueType;
1769 
1770  DenseStencil(const GridType& grid, int halfWidth)
1771  : BaseType(grid, /*size=*/math::Pow3(2 * halfWidth + 1))
1772  , mHalfWidth(halfWidth)
1773  {
1774  assert(halfWidth>0);
1775  }
1776 
1777  inline const ValueType& getCenterValue() const { return mValues[(mValues.size()-1)>>1]; }
1778 
1781  inline void moveTo(const Coord& ijk)
1782  {
1783  BaseType::mCenter = ijk;
1784  this->init(ijk);
1785  }
1788  template<typename IterType>
1789  inline void moveTo(const IterType& iter)
1790  {
1791  BaseType::mCenter = iter.getCoord();
1792  this->init(BaseType::mCenter);
1793  }
1794 
1795 private:
1798  inline void init(const Coord& ijk)
1799  {
1800  int n = 0;
1801  for (Coord p=ijk.offsetBy(-mHalfWidth), q=ijk.offsetBy(mHalfWidth); p[0] <= q[0]; ++p[0]) {
1802  for (p[1] = ijk[1]-mHalfWidth; p[1] <= q[1]; ++p[1]) {
1803  for (p[2] = ijk[2]-mHalfWidth; p[2] <= q[2]; ++p[2]) {
1804  mValues[n++] = mAcc.getValue(p);
1805  }
1806  }
1807  }
1808  }
1809 
1810  template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
1811  using BaseType::mAcc;
1812  using BaseType::mValues;
1813  const int mHalfWidth;
1814 };// DenseStencil class
1815 
1816 
1817 } // end math namespace
1818 } // namespace OPENVDB_VERSION_NAME
1819 } // end openvdb namespace
1820 
1821 #endif // OPENVDB_MATH_STENCILS_HAS_BEEN_INCLUDED
bool zeroCrossing() const
Definition: Stencils.h:1306
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
Bracket code with OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN/_END, to inhibit warnings about type conve...
Definition: Platform.h:206
GridType::Ptr meanCurvature(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the mean curvature of the given grid.
Definition: GridOperators.h:1030
GridT::TreeType TreeType
Definition: Stencils.h:306
ValueType laplacian() const
Definition: Stencils.h:1297
GridT GridType
Definition: Stencils.h:1766
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:1781
ValueType normSqGrad(const ValueType &isoValue=zeroVal< ValueType >()) const
Return the norm-square of the WENO upwind gradient (computed via WENO upwinding and Godunov&#39;s scheme)...
Definition: Stencils.h:1397
Definition: Stencils.h:246
bool intersects(const ValueType &isoValue=zeroVal< ValueType >()) const
Return true if the center of the stencil intersects the.
Definition: Stencils.h:319
SevenPointStencil(const GridT &grid)
Definition: Stencils.h:257
GridT::TreeType TreeType
Definition: Stencils.h:1518
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:261
GridT::ValueType ValueType
Definition: Stencils.h:1519
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
GridType::ValueType ValueType
Definition: Stencils.h:1044
GridT GridType
Definition: Stencils.h:826
GridT GridType
Definition: Stencils.h:556
BaseStencil(const GridType &grid, int size)
Definition: Stencils.h:211
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:315
AccessorType mAcc
Definition: Stencils.h:220
const AccessorType & accessor() const
Return a const reference to the ValueAccessor associated with this Stencil.
Definition: Stencils.h:207
const ValueType & getCenterValue() const
Definition: Stencils.h:1777
GridT::ValueType ValueType
Definition: Stencils.h:253
int size()
Return the size of the stencil buffer.
Definition: Stencils.h:118
GridT::TreeType TreeType
Definition: Stencils.h:477
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:566
WenoStencil(const GridType &grid, Real dx)
Definition: Stencils.h:1384
GridType::ValueType ValueType
Definition: Stencils.h:1238
const ValueType & getValue() const
Return the value at the specified location relative to the center of the stencil. ...
Definition: Stencils.h:105
GridType::ValueType ValueType
Definition: Stencils.h:1768
GridT::ValueType ValueType
Definition: Stencils.h:39
bool intersects(const ValueType &isoValue=zeroVal< ValueType >()) const
Return true if the center of the stencil intersects the iso-contour specified by the isoValue...
Definition: Stencils.h:167
math::Vec3< ValueType > cpt()
Compute the closest-point transform to a level set.
Definition: Stencils.h:1320
void moveTo(const IterType &iter)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:1789
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:1334
ValueType interpolation(const math::Vec3< ValueType > &xyz) const
Return the trilinear interpolation at the normalized position.
Definition: Stencils.h:338
GridT::TreeType TreeType
Definition: Stencils.h:1043
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:97
SecondOrderDenseStencil(const GridType &grid)
Definition: Stencils.h:482
GridT::TreeType TreeType
Definition: Stencils.h:252
Definition: Stencils.h:1231
ThirteenPointStencil(const GridType &grid)
Definition: Stencils.h:562
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:47
const GridType & grid() const
Return a const reference to the grid from which this stencil was constructed.
Definition: Stencils.h:203
ValueType WENO5(const ValueType &v1, const ValueType &v2, const ValueType &v3, const ValueType &v4, const ValueType &v5, float scale2=0.01f)
Implementation of nominally fifth-order finite-difference WENO.
Definition: FiniteDifference.h:304
void moveTo(const IterType &iter)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:72
GridT::ValueType ValueType
Definition: Stencils.h:307
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:1052
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_END
Definition: Platform.h:207
CurvatureStencil(const GridType &grid, Real dx)
Definition: Stencils.h:1530
Tolerance for floating-point comparison.
Definition: Math.h:147
Coord mCenter
Definition: Stencils.h:222
GridT GridType
Definition: Stencils.h:1517
ValueType normSqGrad() const
Return the norm square of the single-sided upwind gradient (computed via Godunov&#39;s scheme) at the pre...
Definition: Stencils.h:1261
FourthOrderDenseStencil(const GridType &grid)
Definition: Stencils.h:693
math::Vec3< ValueType > gradient() const
Return the gradient computed at the previously buffered location by second order central differencing...
Definition: Stencils.h:1277
const GridType * mGrid
Definition: Stencils.h:219
tree::ValueAccessor< const TreeType, IsSafe > AccessorType
Definition: Stencils.h:40
WenoStencil(const GridType &grid)
Definition: Stencils.h:1376
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:486
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:697
BoxStencil(const GridType &grid)
Definition: Stencils.h:311
GridT GridType
Definition: Stencils.h:251
GridT::TreeType TreeType
Definition: Stencils.h:1767
ValueType gaussianCurvature() const
Return the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:1552
Definition: Stencils.h:300
ValueType max() const
Return the largest value in the stencil buffer.
Definition: Stencils.h:153
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:107
void curvatures(ValueType &mean, ValueType &gauss) const
Return both the mean and the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:1564
ValueType meanCurvature() const
Return the mean curvature at the previously buffered location.
Definition: Stencils.h:1541
GridT::TreeType TreeType
Definition: Stencils.h:1237
GridT::TreeType TreeType
Definition: Stencils.h:688
GridType::ValueType ValueType
Definition: Stencils.h:689
Type Pow4(Type x)
Return x4.
Definition: Math.h:559
Definition: Mat.h:187
GridT::TreeType TreeType
Definition: Stencils.h:557
std::vector< ValueType > BufferType
Definition: Stencils.h:41
ValueType mean() const
Return the mean value of the current stencil.
Definition: Stencils.h:138
GridT::TreeType TreeType
Definition: Stencils.h:827
CurvatureStencil(const GridType &grid)
Definition: Stencils.h:1523
math::Vec3< ValueType > gradient(const math::Vec3< ValueType > &xyz) const
Return the gradient in world space of the trilinear interpolation kernel.
Definition: Stencils.h:372
Definition: openvdb/Exceptions.h:13
ValueType median() const
Return the median value of the current stencil.
Definition: Stencils.h:121
NineteenPointStencil(const GridType &grid)
Definition: Stencils.h:832
const Coord & getCenterCoord() const
Return the coordinates of the center point of the stencil.
Definition: Stencils.h:160
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
GridType::ValueType ValueType
Definition: Stencils.h:478
Type Pow3(Type x)
Return x3.
Definition: Math.h:555
ValueType meanCurvatureNormGrad() const
Definition: Stencils.h:1581
GridT::TreeType TreeType
Definition: Stencils.h:38
void moveTo(const Coord &ijk, const ValueType &centerValue)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors. The method also takes a value of the center element of the stencil, assuming it is already known.
Definition: Stencils.h:59
double Real
Definition: openvdb/Types.h:56
Coord offsetBy(Int32 dx, Int32 dy, Int32 dz) const
Definition: Coord.h:92
GridType::ValueType ValueType
Definition: Stencils.h:828
Definition: Stencils.h:1512
GridT::TreeType TreeType
Definition: Stencils.h:1371
Dense stencil of a given width.
Definition: Stencils.h:1761
GridT GridType
Definition: Stencils.h:37
math::Vec3< ValueType > gradient() const
Definition: Stencils.h:1652
GradStencil(const GridType &grid)
Definition: Stencils.h:1242
std::pair< ValueType, ValueType > principalCurvatures() const
Return the pair (minimum, maximum) principal curvature at the previously buffered location...
Definition: Stencils.h:1621
This is a special 19-point stencil that supports optimal fifth-order WENO upwinding, second-order central differencing, Laplacian, and zero-crossing test.
Definition: Stencils.h:1365
void curvaturesNormGrad(ValueType &mean, ValueType &gauss) const
Return both the mean and the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:1605
SixthOrderDenseStencil(const GridType &grid)
Definition: Stencils.h:1048
math::Vec3< ValueType > gradient() const
Definition: Stencils.h:1447
void moveTo(const Vec3< RealType > &xyz)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:86
const ValueType & getCenterValue() const
Return the value at the center of the stencil.
Definition: Stencils.h:163
GridT GridType
Definition: Stencils.h:1370
GridType::ValueType ValueType
Definition: Stencils.h:558
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
GridT GridType
Definition: Stencils.h:687
BufferType mValues
Definition: Stencils.h:221
ValueType gaussianCurvatureNormGrad() const
Definition: Stencils.h:1593
ValueType laplacian() const
Definition: Stencils.h:1459
GridT GridType
Definition: Stencils.h:1236
GridT GridType
Definition: Stencils.h:1042
std::bitset< 6 > intersectionMask(const ValueType &isoValue=zeroVal< ValueType >()) const
Return true a bit-mask where the 6 bits indicates if the center of the stencil intersects the iso-con...
Definition: Stencils.h:188
bool zeroCrossing() const
Definition: Stencils.h:1469
ValueType min() const
Return the smallest value in the stencil buffer.
Definition: Stencils.h:146
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:836
GridT GridType
Definition: Stencils.h:476
math::Vec3< ValueType > gradient(const math::Vec3< ValueType > &V) const
Return the first-order upwind gradient corresponding to the direction V.
Definition: Stencils.h:1287
void setValue(const ValueType &value)
Set the value at the specified location relative to the center of the stencil.
Definition: Stencils.h:112
ValueType laplacian() const
Definition: Stencils.h:1639
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
Definition: Stencils.h:34
math::Vec3< ValueType > gradient(const math::Vec3< ValueType > &V) const
Definition: Stencils.h:1431
GradStencil(const GridType &grid, Real dx)
Definition: Stencils.h:1249
Real GodunovsNormSqrd(bool isOutside, Real dP_xm, Real dP_xp, Real dP_ym, Real dP_yp, Real dP_zm, Real dP_zp)
Definition: FiniteDifference.h:326
BufferType::iterator IterType
Definition: Stencils.h:42
DenseStencil(const GridType &grid, int halfWidth)
Definition: Stencils.h:1770
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:178
GridType::ValueType ValueType
Definition: Stencils.h:1372
GridT GridType
Definition: Stencils.h:305