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