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