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