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