OpenVDB  3.0.0
Stencils.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2014 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 //
33 
34 #ifndef OPENVDB_MATH_STENCILS_HAS_BEEN_INCLUDED
35 #define OPENVDB_MATH_STENCILS_HAS_BEEN_INCLUDED
36 
37 #include <algorithm>
38 #include <vector>
39 #include <openvdb/math/Math.h> // for Pow2, needed by WENO and Gudonov
40 #include <openvdb/Types.h> // for Real
41 #include <openvdb/math/Coord.h> // for Coord
42 #include <openvdb/math/FiniteDifference.h> // for WENO5 and GudonovsNormSqrd
43 
44 namespace openvdb {
46 namespace OPENVDB_VERSION_NAME {
47 namespace math {
48 
49 
51 
52 
53 template<typename _GridType, typename StencilType>
55 {
56 public:
57  typedef _GridType GridType;
58  typedef typename GridType::TreeType TreeType;
59  typedef typename GridType::ValueType ValueType;
60  typedef std::vector<ValueType> BufferType;
61  typedef typename BufferType::iterator IterType;
62  typedef typename GridType::ConstAccessor AccessorType;
63 
67  inline void moveTo(const Coord& ijk)
68  {
69  mCenter = ijk;
70  mStencil[0] = mCache.getValue(ijk);
71  static_cast<StencilType&>(*this).init(mCenter);
72  }
73 
79  inline void moveTo(const Coord& ijk, const ValueType& centerValue)
80  {
81  mCenter = ijk;
82  mStencil[0] = centerValue;
83  static_cast<StencilType&>(*this).init(mCenter);
84  }
85 
91  template<typename IterType>
92  inline void moveTo(const IterType& iter)
93  {
94  mCenter = iter.getCoord();
95  mStencil[0] = *iter;
96  static_cast<StencilType&>(*this).init(mCenter);
97  }
98 
105  inline void moveTo(const Vec3R& xyz)
106  {
107  Coord ijk = openvdb::Coord::floor(xyz);
108  if (ijk != mCenter) this->moveTo(ijk);
109  }
110 
116  inline const ValueType& getValue(unsigned int pos = 0) const
117  {
118  assert(pos < mStencil.size());
119  return mStencil[pos];
120  }
121 
123  template<int i, int j, int k>
124  inline const ValueType& getValue() const
125  {
126  return mStencil[static_cast<const StencilType&>(*this).template pos<i,j,k>()];
127  }
128 
130  template<int i, int j, int k>
131  inline void setValue(const ValueType& value)
132  {
133  mStencil[static_cast<const StencilType&>(*this).template pos<i,j,k>()] = value;
134  }
135 
137  inline int size() { return mStencil.size(); }
138 
140  inline ValueType median() const
141  {
142  std::vector<ValueType> tmp(mStencil);//local copy
143  assert(!tmp.empty());
144  size_t midpoint = (tmp.size() - 1) >> 1;
145  // Partially sort the vector until the median value is at the midpoint.
146  std::nth_element(tmp.begin(), tmp.begin() + midpoint, tmp.end());
147  return tmp[midpoint];
148  }
149 
151  inline ValueType mean() const
152  {
153  ValueType sum = 0.0;
154  for (int n = 0, s = int(mStencil.size()); n < s; ++n) sum += mStencil[n];
155  return sum / mStencil.size();
156  }
157 
159  inline ValueType min() const
160  {
161  IterType iter = std::min_element(mStencil.begin(), mStencil.end());
162  return *iter;
163  }
164 
166  inline ValueType max() const
167  {
168  IterType iter = std::max_element(mStencil.begin(), mStencil.end());
169  return *iter;
170  }
171 
173  inline const Coord& getCenterCoord() const { return mCenter; }
174 
176  inline const ValueType& getCenterValue() const { return mStencil[0]; }
177 
180  inline bool intersects(const ValueType &isoValue = zeroVal<ValueType>()) const
181  {
182  const bool less = this->getValue< 0, 0, 0>() < isoValue;
183  return (less ^ (this->getValue<-1, 0, 0>() < isoValue)) ||
184  (less ^ (this->getValue< 1, 0, 0>() < isoValue)) ||
185  (less ^ (this->getValue< 0,-1, 0>() < isoValue)) ||
186  (less ^ (this->getValue< 0, 1, 0>() < isoValue)) ||
187  (less ^ (this->getValue< 0, 0,-1>() < isoValue)) ||
188  (less ^ (this->getValue< 0, 0, 1>() < isoValue)) ;
189  }
190 
193  inline const GridType& grid() const { return *mGrid; }
194 
197  inline const AccessorType& accessor() const { return mCache; }
198 
199 protected:
200  // Constructor is protected to prevent direct instantiation.
201  BaseStencil(const GridType& grid, int size):
202  mGrid(&grid), mCache(grid.getConstAccessor()),
203  mStencil(size), mCenter(Coord::max())
204  {
205  }
206 
207  const GridType* mGrid;
208  AccessorType mCache;
209  BufferType mStencil;
211 
212 }; // class BaseStencil
213 
214 
216 
217 
218 namespace { // anonymous namespace for stencil-layout map
219 
220  // the seven point stencil
221  template<int i, int j, int k> struct SevenPt {};
222  template<> struct SevenPt< 0, 0, 0> { enum { idx = 0 }; };
223  template<> struct SevenPt< 1, 0, 0> { enum { idx = 1 }; };
224  template<> struct SevenPt< 0, 1, 0> { enum { idx = 2 }; };
225  template<> struct SevenPt< 0, 0, 1> { enum { idx = 3 }; };
226  template<> struct SevenPt<-1, 0, 0> { enum { idx = 4 }; };
227  template<> struct SevenPt< 0,-1, 0> { enum { idx = 5 }; };
228  template<> struct SevenPt< 0, 0,-1> { enum { idx = 6 }; };
229 
230 }
231 
232 
233 template<typename GridType>
234 class SevenPointStencil: public BaseStencil<GridType, SevenPointStencil<GridType> >
235 {
236 public:
239  typedef typename GridType::ValueType ValueType;
241  static const int SIZE = 7;
242 
243  SevenPointStencil(const GridType& grid): BaseType(grid, SIZE) {}
244 
246  template<int i, int j, int k>
247  unsigned int pos() const { return SevenPt<i,j,k>::idx; }
248 
249 private:
250  inline void init(const Coord& ijk)
251  {
252  BaseType::template setValue<-1, 0, 0>(mCache.getValue(ijk.offsetBy(-1, 0, 0)));
253  BaseType::template setValue< 1, 0, 0>(mCache.getValue(ijk.offsetBy( 1, 0, 0)));
254 
255  BaseType::template setValue< 0,-1, 0>(mCache.getValue(ijk.offsetBy( 0,-1, 0)));
256  BaseType::template setValue< 0, 1, 0>(mCache.getValue(ijk.offsetBy( 0, 1, 0)));
257 
258  BaseType::template setValue< 0, 0,-1>(mCache.getValue(ijk.offsetBy( 0, 0,-1)));
259  BaseType::template setValue< 0, 0, 1>(mCache.getValue(ijk.offsetBy( 0, 0, 1)));
260  }
261 
262  template<typename, typename> friend class BaseStencil; // allow base class to call init()
263  using BaseType::mCache;
264  using BaseType::mStencil;
265 };
266 
267 
269 
270 
271 namespace { // anonymous namespace for stencil-layout map
272 
273  // the eight point box stencil
274  template<int i, int j, int k> struct BoxPt {};
275  template<> struct BoxPt< 0, 0, 0> { enum { idx = 0 }; };
276  template<> struct BoxPt< 0, 0, 1> { enum { idx = 1 }; };
277  template<> struct BoxPt< 0, 1, 1> { enum { idx = 2 }; };
278  template<> struct BoxPt< 0, 1, 0> { enum { idx = 3 }; };
279  template<> struct BoxPt< 1, 0, 0> { enum { idx = 4 }; };
280  template<> struct BoxPt< 1, 0, 1> { enum { idx = 5 }; };
281  template<> struct BoxPt< 1, 1, 1> { enum { idx = 6 }; };
282  template<> struct BoxPt< 1, 1, 0> { enum { idx = 7 }; };
283 }
284 
285 template<typename GridType>
286 class BoxStencil: public BaseStencil<GridType, BoxStencil<GridType> >
287 {
288 public:
291  typedef typename GridType::ValueType ValueType;
293  static const int SIZE = 8;
294 
295  BoxStencil(const GridType& grid): BaseType(grid, SIZE) {}
296 
298  template<int i, int j, int k>
299  unsigned int pos() const { return BoxPt<i,j,k>::idx; }
300 
303  inline bool intersects(const ValueType &isoValue = zeroVal<ValueType>()) const
304  {
305  const bool less = mStencil[0] < isoValue;
306  return (less ^ (mStencil[1] < isoValue)) ||
307  (less ^ (mStencil[2] < isoValue)) ||
308  (less ^ (mStencil[3] < isoValue)) ||
309  (less ^ (mStencil[4] < isoValue)) ||
310  (less ^ (mStencil[5] < isoValue)) ||
311  (less ^ (mStencil[6] < isoValue)) ||
312  (less ^ (mStencil[7] < isoValue)) ;
313  }
314 
322  inline ValueType interpolation(const Vec3Type& xyz) const
323  {
324  const Real u = xyz[0] - BaseType::mCenter[0]; assert(u>=0 && u<=1);
325  const Real v = xyz[1] - BaseType::mCenter[1]; assert(v>=0 && v<=1);
326  const Real w = xyz[2] - BaseType::mCenter[2]; assert(w>=0 && w<=1);
327 
328  ValueType V = BaseType::template getValue<0,0,0>();
329  ValueType A = static_cast<ValueType>(V + (BaseType::template getValue<0,0,1>() - V) * w);
330  V = BaseType::template getValue< 0, 1, 0>();
331  ValueType B = static_cast<ValueType>(V + (BaseType::template getValue<0,1,1>() - V) * w);
332  ValueType C = static_cast<ValueType>(A + (B - A) * v);
333 
334  V = BaseType::template getValue<1,0,0>();
335  A = static_cast<ValueType>(V + (BaseType::template getValue<1,0,1>() - V) * w);
336  V = BaseType::template getValue<1,1,0>();
337  B = static_cast<ValueType>(V + (BaseType::template getValue<1,1,1>() - V) * w);
338  ValueType D = static_cast<ValueType>(A + (B - A) * v);
339 
340  return static_cast<ValueType>(C + (D - C) * u);
341  }
342 
350  inline Vec3Type gradient(const Vec3Type& xyz) const
351  {
352  const Real u = xyz[0] - BaseType::mCenter[0]; assert(u>=0 && u<=1);
353  const Real v = xyz[1] - BaseType::mCenter[1]; assert(v>=0 && v<=1);
354  const Real w = xyz[2] - BaseType::mCenter[2]; assert(w>=0 && w<=1);
355 
356  ValueType D[4]={BaseType::template getValue<0,0,1>()-BaseType::template getValue<0,0,0>(),
357  BaseType::template getValue<0,1,1>()-BaseType::template getValue<0,1,0>(),
358  BaseType::template getValue<1,0,1>()-BaseType::template getValue<1,0,0>(),
359  BaseType::template getValue<1,1,1>()-BaseType::template getValue<1,1,0>()};
360 
361  // Z component
362  ValueType A = static_cast<ValueType>(D[0] + (D[1]- D[0]) * v);
363  ValueType B = static_cast<ValueType>(D[2] + (D[3]- D[2]) * v);
364  Vec3Type grad(
365  zeroVal<ValueType>(), zeroVal<ValueType>(), static_cast<ValueType>(A + (B - A) * u));
366 
367  D[0] = static_cast<ValueType>(BaseType::template getValue<0,0,0>() + D[0] * w);
368  D[1] = static_cast<ValueType>(BaseType::template getValue<0,1,0>() + D[1] * w);
369  D[2] = static_cast<ValueType>(BaseType::template getValue<1,0,0>() + D[2] * w);
370  D[3] = static_cast<ValueType>(BaseType::template getValue<1,1,0>() + D[3] * w);
371 
372  // X component
373  A = static_cast<ValueType>(D[0] + (D[1] - D[0]) * v);
374  B = static_cast<ValueType>(D[2] + (D[3] - D[2]) * v);
375 
376  grad[0] = B - A;
377 
378  // Y component
379  A = D[1] - D[0];
380  B = D[3] - D[2];
381 
382  grad[1] = static_cast<ValueType>(A + (B - A) * u);
383 
384  return BaseType::mGrid->transform().baseMap()->applyIJT(grad, xyz);
385  }
386 
387 private:
388  inline void init(const Coord& ijk)
389  {
390  BaseType::template setValue< 0, 0, 1>(mCache.getValue(ijk.offsetBy( 0, 0, 1)));
391  BaseType::template setValue< 0, 1, 1>(mCache.getValue(ijk.offsetBy( 0, 1, 1)));
392  BaseType::template setValue< 0, 1, 0>(mCache.getValue(ijk.offsetBy( 0, 1, 0)));
393  BaseType::template setValue< 1, 0, 0>(mCache.getValue(ijk.offsetBy( 1, 0, 0)));
394  BaseType::template setValue< 1, 0, 1>(mCache.getValue(ijk.offsetBy( 1, 0, 1)));
395  BaseType::template setValue< 1, 1, 1>(mCache.getValue(ijk.offsetBy( 1, 1, 1)));
396  BaseType::template setValue< 1, 1, 0>(mCache.getValue(ijk.offsetBy( 1, 1, 0)));
397  }
398 
399  template<typename, typename> friend class BaseStencil; // allow base class to call init()
400  using BaseType::mCache;
401  using BaseType::mStencil;
402 };
403 
404 
406 
407 
408 namespace { // anonymous namespace for stencil-layout map
409 
410  // the dense point stencil
411  template<int i, int j, int k> struct DensePt {};
412  template<> struct DensePt< 0, 0, 0> { enum { idx = 0 }; };
413 
414  template<> struct DensePt< 1, 0, 0> { enum { idx = 1 }; };
415  template<> struct DensePt< 0, 1, 0> { enum { idx = 2 }; };
416  template<> struct DensePt< 0, 0, 1> { enum { idx = 3 }; };
417 
418  template<> struct DensePt<-1, 0, 0> { enum { idx = 4 }; };
419  template<> struct DensePt< 0,-1, 0> { enum { idx = 5 }; };
420  template<> struct DensePt< 0, 0,-1> { enum { idx = 6 }; };
421 
422  template<> struct DensePt<-1,-1, 0> { enum { idx = 7 }; };
423  template<> struct DensePt< 0,-1,-1> { enum { idx = 8 }; };
424  template<> struct DensePt<-1, 0,-1> { enum { idx = 9 }; };
425 
426  template<> struct DensePt< 1,-1, 0> { enum { idx = 10 }; };
427  template<> struct DensePt< 0, 1,-1> { enum { idx = 11 }; };
428  template<> struct DensePt<-1, 0, 1> { enum { idx = 12 }; };
429 
430  template<> struct DensePt<-1, 1, 0> { enum { idx = 13 }; };
431  template<> struct DensePt< 0,-1, 1> { enum { idx = 14 }; };
432  template<> struct DensePt< 1, 0,-1> { enum { idx = 15 }; };
433 
434  template<> struct DensePt< 1, 1, 0> { enum { idx = 16 }; };
435  template<> struct DensePt< 0, 1, 1> { enum { idx = 17 }; };
436  template<> struct DensePt< 1, 0, 1> { enum { idx = 18 }; };
437 
438 }
439 
440 
441 template<typename GridType>
442 class SecondOrderDenseStencil: public BaseStencil<GridType, SecondOrderDenseStencil<GridType> >
443 {
444 public:
447  typedef typename GridType::ValueType ValueType;
449 
450  static const int SIZE = 19;
451 
452  SecondOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
453 
455  template<int i, int j, int k>
456  unsigned int pos() const { return DensePt<i,j,k>::idx; }
457 
458 private:
459  inline void init(const Coord& ijk)
460  {
461  mStencil[DensePt< 1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
462  mStencil[DensePt< 0, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
463  mStencil[DensePt< 0, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
464 
465  mStencil[DensePt<-1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
466  mStencil[DensePt< 0,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
467  mStencil[DensePt< 0, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
468 
469  mStencil[DensePt<-1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, -1, 0));
470  mStencil[DensePt< 1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, -1, 0));
471  mStencil[DensePt<-1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 1, 0));
472  mStencil[DensePt< 1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 1, 0));
473 
474  mStencil[DensePt<-1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, -1));
475  mStencil[DensePt< 1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, -1));
476  mStencil[DensePt<-1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 1));
477  mStencil[DensePt< 1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 1));
478 
479  mStencil[DensePt< 0,-1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, -1, -1));
480  mStencil[DensePt< 0, 1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, -1));
481  mStencil[DensePt< 0,-1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, -1, 1));
482  mStencil[DensePt< 0, 1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 1));
483  }
484 
485  template<typename, typename> friend class BaseStencil; // allow base class to call init()
486  using BaseType::mCache;
487  using BaseType::mStencil;
488 };
489 
490 
492 
493 
494 namespace { // anonymous namespace for stencil-layout map
495 
496  // the dense point stencil
497  template<int i, int j, int k> struct ThirteenPt {};
498  template<> struct ThirteenPt< 0, 0, 0> { enum { idx = 0 }; };
499 
500  template<> struct ThirteenPt< 1, 0, 0> { enum { idx = 1 }; };
501  template<> struct ThirteenPt< 0, 1, 0> { enum { idx = 2 }; };
502  template<> struct ThirteenPt< 0, 0, 1> { enum { idx = 3 }; };
503 
504  template<> struct ThirteenPt<-1, 0, 0> { enum { idx = 4 }; };
505  template<> struct ThirteenPt< 0,-1, 0> { enum { idx = 5 }; };
506  template<> struct ThirteenPt< 0, 0,-1> { enum { idx = 6 }; };
507 
508  template<> struct ThirteenPt< 2, 0, 0> { enum { idx = 7 }; };
509  template<> struct ThirteenPt< 0, 2, 0> { enum { idx = 8 }; };
510  template<> struct ThirteenPt< 0, 0, 2> { enum { idx = 9 }; };
511 
512  template<> struct ThirteenPt<-2, 0, 0> { enum { idx = 10 }; };
513  template<> struct ThirteenPt< 0,-2, 0> { enum { idx = 11 }; };
514  template<> struct ThirteenPt< 0, 0,-2> { enum { idx = 12 }; };
515 
516 }
517 
518 
519 template<typename GridType>
520 class ThirteenPointStencil: public BaseStencil<GridType, ThirteenPointStencil<GridType> >
521 {
522 public:
525  typedef typename GridType::ValueType ValueType;
527 
528  static const int SIZE = 13;
529 
530  ThirteenPointStencil(const GridType& grid): BaseType(grid, SIZE) {}
531 
533  template<int i, int j, int k>
534  unsigned int pos() const { return ThirteenPt<i,j,k>::idx; }
535 
536 private:
537  inline void init(const Coord& ijk)
538  {
539  mStencil[ThirteenPt< 2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 0));
540  mStencil[ThirteenPt< 1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
541  mStencil[ThirteenPt<-1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
542  mStencil[ThirteenPt<-2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 0));
543 
544  mStencil[ThirteenPt< 0, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 0));
545  mStencil[ThirteenPt< 0, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
546  mStencil[ThirteenPt< 0,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
547  mStencil[ThirteenPt< 0,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -2, 0));
548 
549  mStencil[ThirteenPt< 0, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 2));
550  mStencil[ThirteenPt< 0, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
551  mStencil[ThirteenPt< 0, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
552  mStencil[ThirteenPt< 0, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -2));
553  }
554 
555  template<typename, typename> friend class BaseStencil; // allow base class to call init()
556  using BaseType::mCache;
557  using BaseType::mStencil;
558 };
559 
560 
562 
563 
564 namespace { // anonymous namespace for stencil-layout map
565 
566  // the 4th-order dense point stencil
567  template<int i, int j, int k> struct FourthDensePt {};
568  template<> struct FourthDensePt< 0, 0, 0> { enum { idx = 0 }; };
569 
570  template<> struct FourthDensePt<-2, 2, 0> { enum { idx = 1 }; };
571  template<> struct FourthDensePt<-1, 2, 0> { enum { idx = 2 }; };
572  template<> struct FourthDensePt< 0, 2, 0> { enum { idx = 3 }; };
573  template<> struct FourthDensePt< 1, 2, 0> { enum { idx = 4 }; };
574  template<> struct FourthDensePt< 2, 2, 0> { enum { idx = 5 }; };
575 
576  template<> struct FourthDensePt<-2, 1, 0> { enum { idx = 6 }; };
577  template<> struct FourthDensePt<-1, 1, 0> { enum { idx = 7 }; };
578  template<> struct FourthDensePt< 0, 1, 0> { enum { idx = 8 }; };
579  template<> struct FourthDensePt< 1, 1, 0> { enum { idx = 9 }; };
580  template<> struct FourthDensePt< 2, 1, 0> { enum { idx = 10 }; };
581 
582  template<> struct FourthDensePt<-2, 0, 0> { enum { idx = 11 }; };
583  template<> struct FourthDensePt<-1, 0, 0> { enum { idx = 12 }; };
584  template<> struct FourthDensePt< 1, 0, 0> { enum { idx = 13 }; };
585  template<> struct FourthDensePt< 2, 0, 0> { enum { idx = 14 }; };
586 
587  template<> struct FourthDensePt<-2,-1, 0> { enum { idx = 15 }; };
588  template<> struct FourthDensePt<-1,-1, 0> { enum { idx = 16 }; };
589  template<> struct FourthDensePt< 0,-1, 0> { enum { idx = 17 }; };
590  template<> struct FourthDensePt< 1,-1, 0> { enum { idx = 18 }; };
591  template<> struct FourthDensePt< 2,-1, 0> { enum { idx = 19 }; };
592 
593  template<> struct FourthDensePt<-2,-2, 0> { enum { idx = 20 }; };
594  template<> struct FourthDensePt<-1,-2, 0> { enum { idx = 21 }; };
595  template<> struct FourthDensePt< 0,-2, 0> { enum { idx = 22 }; };
596  template<> struct FourthDensePt< 1,-2, 0> { enum { idx = 23 }; };
597  template<> struct FourthDensePt< 2,-2, 0> { enum { idx = 24 }; };
598 
599 
600  template<> struct FourthDensePt<-2, 0, 2> { enum { idx = 25 }; };
601  template<> struct FourthDensePt<-1, 0, 2> { enum { idx = 26 }; };
602  template<> struct FourthDensePt< 0, 0, 2> { enum { idx = 27 }; };
603  template<> struct FourthDensePt< 1, 0, 2> { enum { idx = 28 }; };
604  template<> struct FourthDensePt< 2, 0, 2> { enum { idx = 29 }; };
605 
606  template<> struct FourthDensePt<-2, 0, 1> { enum { idx = 30 }; };
607  template<> struct FourthDensePt<-1, 0, 1> { enum { idx = 31 }; };
608  template<> struct FourthDensePt< 0, 0, 1> { enum { idx = 32 }; };
609  template<> struct FourthDensePt< 1, 0, 1> { enum { idx = 33 }; };
610  template<> struct FourthDensePt< 2, 0, 1> { enum { idx = 34 }; };
611 
612  template<> struct FourthDensePt<-2, 0,-1> { enum { idx = 35 }; };
613  template<> struct FourthDensePt<-1, 0,-1> { enum { idx = 36 }; };
614  template<> struct FourthDensePt< 0, 0,-1> { enum { idx = 37 }; };
615  template<> struct FourthDensePt< 1, 0,-1> { enum { idx = 38 }; };
616  template<> struct FourthDensePt< 2, 0,-1> { enum { idx = 39 }; };
617 
618  template<> struct FourthDensePt<-2, 0,-2> { enum { idx = 40 }; };
619  template<> struct FourthDensePt<-1, 0,-2> { enum { idx = 41 }; };
620  template<> struct FourthDensePt< 0, 0,-2> { enum { idx = 42 }; };
621  template<> struct FourthDensePt< 1, 0,-2> { enum { idx = 43 }; };
622  template<> struct FourthDensePt< 2, 0,-2> { enum { idx = 44 }; };
623 
624 
625  template<> struct FourthDensePt< 0,-2, 2> { enum { idx = 45 }; };
626  template<> struct FourthDensePt< 0,-1, 2> { enum { idx = 46 }; };
627  template<> struct FourthDensePt< 0, 1, 2> { enum { idx = 47 }; };
628  template<> struct FourthDensePt< 0, 2, 2> { enum { idx = 48 }; };
629 
630  template<> struct FourthDensePt< 0,-2, 1> { enum { idx = 49 }; };
631  template<> struct FourthDensePt< 0,-1, 1> { enum { idx = 50 }; };
632  template<> struct FourthDensePt< 0, 1, 1> { enum { idx = 51 }; };
633  template<> struct FourthDensePt< 0, 2, 1> { enum { idx = 52 }; };
634 
635  template<> struct FourthDensePt< 0,-2,-1> { enum { idx = 53 }; };
636  template<> struct FourthDensePt< 0,-1,-1> { enum { idx = 54 }; };
637  template<> struct FourthDensePt< 0, 1,-1> { enum { idx = 55 }; };
638  template<> struct FourthDensePt< 0, 2,-1> { enum { idx = 56 }; };
639 
640  template<> struct FourthDensePt< 0,-2,-2> { enum { idx = 57 }; };
641  template<> struct FourthDensePt< 0,-1,-2> { enum { idx = 58 }; };
642  template<> struct FourthDensePt< 0, 1,-2> { enum { idx = 59 }; };
643  template<> struct FourthDensePt< 0, 2,-2> { enum { idx = 60 }; };
644 
645 }
646 
647 
648 template<typename GridType>
649 class FourthOrderDenseStencil: public BaseStencil<GridType, FourthOrderDenseStencil<GridType> >
650 {
651 public:
654  typedef typename GridType::ValueType ValueType;
656 
657  static const int SIZE = 61;
658 
659  FourthOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
660 
662  template<int i, int j, int k>
663  unsigned int pos() const { return FourthDensePt<i,j,k>::idx; }
664 
665 private:
666  inline void init(const Coord& ijk)
667  {
668  mStencil[FourthDensePt<-2, 2, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 2, 0));
669  mStencil[FourthDensePt<-1, 2, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 2, 0));
670  mStencil[FourthDensePt< 0, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 0));
671  mStencil[FourthDensePt< 1, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 2, 0));
672  mStencil[FourthDensePt< 2, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 2, 0));
673 
674  mStencil[FourthDensePt<-2, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 1, 0));
675  mStencil[FourthDensePt<-1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 1, 0));
676  mStencil[FourthDensePt< 0, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
677  mStencil[FourthDensePt< 1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 1, 0));
678  mStencil[FourthDensePt< 2, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 1, 0));
679 
680  mStencil[FourthDensePt<-2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 0));
681  mStencil[FourthDensePt<-1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
682  mStencil[FourthDensePt< 1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
683  mStencil[FourthDensePt< 2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 0));
684 
685  mStencil[FourthDensePt<-2,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-2,-1, 0));
686  mStencil[FourthDensePt<-1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1,-1, 0));
687  mStencil[FourthDensePt< 0,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 0));
688  mStencil[FourthDensePt< 1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1,-1, 0));
689  mStencil[FourthDensePt< 2,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 2,-1, 0));
690 
691  mStencil[FourthDensePt<-2,-2, 0>::idx] = mCache.getValue(ijk.offsetBy(-2,-2, 0));
692  mStencil[FourthDensePt<-1,-2, 0>::idx] = mCache.getValue(ijk.offsetBy(-1,-2, 0));
693  mStencil[FourthDensePt< 0,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 0));
694  mStencil[FourthDensePt< 1,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 1,-2, 0));
695  mStencil[FourthDensePt< 2,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 2,-2, 0));
696 
697  mStencil[FourthDensePt<-2, 0, 2>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 2));
698  mStencil[FourthDensePt<-1, 0, 2>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 2));
699  mStencil[FourthDensePt< 0, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 2));
700  mStencil[FourthDensePt< 1, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 2));
701  mStencil[FourthDensePt< 2, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 2));
702 
703  mStencil[FourthDensePt<-2, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 1));
704  mStencil[FourthDensePt<-1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 1));
705  mStencil[FourthDensePt< 0, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
706  mStencil[FourthDensePt< 1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 1));
707  mStencil[FourthDensePt< 2, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 1));
708 
709  mStencil[FourthDensePt<-2, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-2, 0,-1));
710  mStencil[FourthDensePt<-1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0,-1));
711  mStencil[FourthDensePt< 0, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0,-1));
712  mStencil[FourthDensePt< 1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0,-1));
713  mStencil[FourthDensePt< 2, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 2, 0,-1));
714 
715  mStencil[FourthDensePt<-2, 0,-2>::idx] = mCache.getValue(ijk.offsetBy(-2, 0,-2));
716  mStencil[FourthDensePt<-1, 0,-2>::idx] = mCache.getValue(ijk.offsetBy(-1, 0,-2));
717  mStencil[FourthDensePt< 0, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0,-2));
718  mStencil[FourthDensePt< 1, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 1, 0,-2));
719  mStencil[FourthDensePt< 2, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 2, 0,-2));
720 
721 
722  mStencil[FourthDensePt< 0,-2, 2>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 2));
723  mStencil[FourthDensePt< 0,-1, 2>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 2));
724  mStencil[FourthDensePt< 0, 1, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 2));
725  mStencil[FourthDensePt< 0, 2, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 2));
726 
727  mStencil[FourthDensePt< 0,-2, 1>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 1));
728  mStencil[FourthDensePt< 0,-1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 1));
729  mStencil[FourthDensePt< 0, 1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 1));
730  mStencil[FourthDensePt< 0, 2, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 1));
731 
732  mStencil[FourthDensePt< 0,-2,-1>::idx] = mCache.getValue(ijk.offsetBy( 0,-2,-1));
733  mStencil[FourthDensePt< 0,-1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0,-1,-1));
734  mStencil[FourthDensePt< 0, 1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1,-1));
735  mStencil[FourthDensePt< 0, 2,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 2,-1));
736 
737  mStencil[FourthDensePt< 0,-2,-2>::idx] = mCache.getValue(ijk.offsetBy( 0,-2,-2));
738  mStencil[FourthDensePt< 0,-1,-2>::idx] = mCache.getValue(ijk.offsetBy( 0,-1,-2));
739  mStencil[FourthDensePt< 0, 1,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 1,-2));
740  mStencil[FourthDensePt< 0, 2,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 2,-2));
741  }
742 
743  template<typename, typename> friend class BaseStencil; // allow base class to call init()
744  using BaseType::mCache;
745  using BaseType::mStencil;
746 };
747 
748 
750 
751 
752 namespace { // anonymous namespace for stencil-layout map
753 
754  // the dense point stencil
755  template<int i, int j, int k> struct NineteenPt {};
756  template<> struct NineteenPt< 0, 0, 0> { enum { idx = 0 }; };
757 
758  template<> struct NineteenPt< 1, 0, 0> { enum { idx = 1 }; };
759  template<> struct NineteenPt< 0, 1, 0> { enum { idx = 2 }; };
760  template<> struct NineteenPt< 0, 0, 1> { enum { idx = 3 }; };
761 
762  template<> struct NineteenPt<-1, 0, 0> { enum { idx = 4 }; };
763  template<> struct NineteenPt< 0,-1, 0> { enum { idx = 5 }; };
764  template<> struct NineteenPt< 0, 0,-1> { enum { idx = 6 }; };
765 
766  template<> struct NineteenPt< 2, 0, 0> { enum { idx = 7 }; };
767  template<> struct NineteenPt< 0, 2, 0> { enum { idx = 8 }; };
768  template<> struct NineteenPt< 0, 0, 2> { enum { idx = 9 }; };
769 
770  template<> struct NineteenPt<-2, 0, 0> { enum { idx = 10 }; };
771  template<> struct NineteenPt< 0,-2, 0> { enum { idx = 11 }; };
772  template<> struct NineteenPt< 0, 0,-2> { enum { idx = 12 }; };
773 
774  template<> struct NineteenPt< 3, 0, 0> { enum { idx = 13 }; };
775  template<> struct NineteenPt< 0, 3, 0> { enum { idx = 14 }; };
776  template<> struct NineteenPt< 0, 0, 3> { enum { idx = 15 }; };
777 
778  template<> struct NineteenPt<-3, 0, 0> { enum { idx = 16 }; };
779  template<> struct NineteenPt< 0,-3, 0> { enum { idx = 17 }; };
780  template<> struct NineteenPt< 0, 0,-3> { enum { idx = 18 }; };
781 
782 }
783 
784 
785 template<typename GridType>
786 class NineteenPointStencil: public BaseStencil<GridType, NineteenPointStencil<GridType> >
787 {
788 public:
791  typedef typename GridType::ValueType ValueType;
793 
794  static const int SIZE = 19;
795 
796  NineteenPointStencil(const GridType& grid): BaseType(grid, SIZE) {}
797 
799  template<int i, int j, int k>
800  unsigned int pos() const { return NineteenPt<i,j,k>::idx; }
801 
802 private:
803  inline void init(const Coord& ijk)
804  {
805  mStencil[NineteenPt< 3, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 3, 0, 0));
806  mStencil[NineteenPt< 2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 0));
807  mStencil[NineteenPt< 1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
808  mStencil[NineteenPt<-1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
809  mStencil[NineteenPt<-2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 0));
810  mStencil[NineteenPt<-3, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-3, 0, 0));
811 
812  mStencil[NineteenPt< 0, 3, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 3, 0));
813  mStencil[NineteenPt< 0, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 0));
814  mStencil[NineteenPt< 0, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
815  mStencil[NineteenPt< 0,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
816  mStencil[NineteenPt< 0,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -2, 0));
817  mStencil[NineteenPt< 0,-3, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -3, 0));
818 
819  mStencil[NineteenPt< 0, 0, 3>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 3));
820  mStencil[NineteenPt< 0, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 2));
821  mStencil[NineteenPt< 0, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
822  mStencil[NineteenPt< 0, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
823  mStencil[NineteenPt< 0, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -2));
824  mStencil[NineteenPt< 0, 0,-3>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -3));
825  }
826 
827  template<typename, typename> friend class BaseStencil; // allow base class to call init()
828  using BaseType::mCache;
829  using BaseType::mStencil;
830 };
831 
832 
834 
835 
836 namespace { // anonymous namespace for stencil-layout map
837 
838  // the 4th-order dense point stencil
839  template<int i, int j, int k> struct SixthDensePt { };
840  template<> struct SixthDensePt< 0, 0, 0> { enum { idx = 0 }; };
841 
842  template<> struct SixthDensePt<-3, 3, 0> { enum { idx = 1 }; };
843  template<> struct SixthDensePt<-2, 3, 0> { enum { idx = 2 }; };
844  template<> struct SixthDensePt<-1, 3, 0> { enum { idx = 3 }; };
845  template<> struct SixthDensePt< 0, 3, 0> { enum { idx = 4 }; };
846  template<> struct SixthDensePt< 1, 3, 0> { enum { idx = 5 }; };
847  template<> struct SixthDensePt< 2, 3, 0> { enum { idx = 6 }; };
848  template<> struct SixthDensePt< 3, 3, 0> { enum { idx = 7 }; };
849 
850  template<> struct SixthDensePt<-3, 2, 0> { enum { idx = 8 }; };
851  template<> struct SixthDensePt<-2, 2, 0> { enum { idx = 9 }; };
852  template<> struct SixthDensePt<-1, 2, 0> { enum { idx = 10 }; };
853  template<> struct SixthDensePt< 0, 2, 0> { enum { idx = 11 }; };
854  template<> struct SixthDensePt< 1, 2, 0> { enum { idx = 12 }; };
855  template<> struct SixthDensePt< 2, 2, 0> { enum { idx = 13 }; };
856  template<> struct SixthDensePt< 3, 2, 0> { enum { idx = 14 }; };
857 
858  template<> struct SixthDensePt<-3, 1, 0> { enum { idx = 15 }; };
859  template<> struct SixthDensePt<-2, 1, 0> { enum { idx = 16 }; };
860  template<> struct SixthDensePt<-1, 1, 0> { enum { idx = 17 }; };
861  template<> struct SixthDensePt< 0, 1, 0> { enum { idx = 18 }; };
862  template<> struct SixthDensePt< 1, 1, 0> { enum { idx = 19 }; };
863  template<> struct SixthDensePt< 2, 1, 0> { enum { idx = 20 }; };
864  template<> struct SixthDensePt< 3, 1, 0> { enum { idx = 21 }; };
865 
866  template<> struct SixthDensePt<-3, 0, 0> { enum { idx = 22 }; };
867  template<> struct SixthDensePt<-2, 0, 0> { enum { idx = 23 }; };
868  template<> struct SixthDensePt<-1, 0, 0> { enum { idx = 24 }; };
869  template<> struct SixthDensePt< 1, 0, 0> { enum { idx = 25 }; };
870  template<> struct SixthDensePt< 2, 0, 0> { enum { idx = 26 }; };
871  template<> struct SixthDensePt< 3, 0, 0> { enum { idx = 27 }; };
872 
873 
874  template<> struct SixthDensePt<-3,-1, 0> { enum { idx = 28 }; };
875  template<> struct SixthDensePt<-2,-1, 0> { enum { idx = 29 }; };
876  template<> struct SixthDensePt<-1,-1, 0> { enum { idx = 30 }; };
877  template<> struct SixthDensePt< 0,-1, 0> { enum { idx = 31 }; };
878  template<> struct SixthDensePt< 1,-1, 0> { enum { idx = 32 }; };
879  template<> struct SixthDensePt< 2,-1, 0> { enum { idx = 33 }; };
880  template<> struct SixthDensePt< 3,-1, 0> { enum { idx = 34 }; };
881 
882 
883  template<> struct SixthDensePt<-3,-2, 0> { enum { idx = 35 }; };
884  template<> struct SixthDensePt<-2,-2, 0> { enum { idx = 36 }; };
885  template<> struct SixthDensePt<-1,-2, 0> { enum { idx = 37 }; };
886  template<> struct SixthDensePt< 0,-2, 0> { enum { idx = 38 }; };
887  template<> struct SixthDensePt< 1,-2, 0> { enum { idx = 39 }; };
888  template<> struct SixthDensePt< 2,-2, 0> { enum { idx = 40 }; };
889  template<> struct SixthDensePt< 3,-2, 0> { enum { idx = 41 }; };
890 
891 
892  template<> struct SixthDensePt<-3,-3, 0> { enum { idx = 42 }; };
893  template<> struct SixthDensePt<-2,-3, 0> { enum { idx = 43 }; };
894  template<> struct SixthDensePt<-1,-3, 0> { enum { idx = 44 }; };
895  template<> struct SixthDensePt< 0,-3, 0> { enum { idx = 45 }; };
896  template<> struct SixthDensePt< 1,-3, 0> { enum { idx = 46 }; };
897  template<> struct SixthDensePt< 2,-3, 0> { enum { idx = 47 }; };
898  template<> struct SixthDensePt< 3,-3, 0> { enum { idx = 48 }; };
899 
900 
901  template<> struct SixthDensePt<-3, 0, 3> { enum { idx = 49 }; };
902  template<> struct SixthDensePt<-2, 0, 3> { enum { idx = 50 }; };
903  template<> struct SixthDensePt<-1, 0, 3> { enum { idx = 51 }; };
904  template<> struct SixthDensePt< 0, 0, 3> { enum { idx = 52 }; };
905  template<> struct SixthDensePt< 1, 0, 3> { enum { idx = 53 }; };
906  template<> struct SixthDensePt< 2, 0, 3> { enum { idx = 54 }; };
907  template<> struct SixthDensePt< 3, 0, 3> { enum { idx = 55 }; };
908 
909 
910  template<> struct SixthDensePt<-3, 0, 2> { enum { idx = 56 }; };
911  template<> struct SixthDensePt<-2, 0, 2> { enum { idx = 57 }; };
912  template<> struct SixthDensePt<-1, 0, 2> { enum { idx = 58 }; };
913  template<> struct SixthDensePt< 0, 0, 2> { enum { idx = 59 }; };
914  template<> struct SixthDensePt< 1, 0, 2> { enum { idx = 60 }; };
915  template<> struct SixthDensePt< 2, 0, 2> { enum { idx = 61 }; };
916  template<> struct SixthDensePt< 3, 0, 2> { enum { idx = 62 }; };
917 
918  template<> struct SixthDensePt<-3, 0, 1> { enum { idx = 63 }; };
919  template<> struct SixthDensePt<-2, 0, 1> { enum { idx = 64 }; };
920  template<> struct SixthDensePt<-1, 0, 1> { enum { idx = 65 }; };
921  template<> struct SixthDensePt< 0, 0, 1> { enum { idx = 66 }; };
922  template<> struct SixthDensePt< 1, 0, 1> { enum { idx = 67 }; };
923  template<> struct SixthDensePt< 2, 0, 1> { enum { idx = 68 }; };
924  template<> struct SixthDensePt< 3, 0, 1> { enum { idx = 69 }; };
925 
926 
927  template<> struct SixthDensePt<-3, 0,-1> { enum { idx = 70 }; };
928  template<> struct SixthDensePt<-2, 0,-1> { enum { idx = 71 }; };
929  template<> struct SixthDensePt<-1, 0,-1> { enum { idx = 72 }; };
930  template<> struct SixthDensePt< 0, 0,-1> { enum { idx = 73 }; };
931  template<> struct SixthDensePt< 1, 0,-1> { enum { idx = 74 }; };
932  template<> struct SixthDensePt< 2, 0,-1> { enum { idx = 75 }; };
933  template<> struct SixthDensePt< 3, 0,-1> { enum { idx = 76 }; };
934 
935 
936  template<> struct SixthDensePt<-3, 0,-2> { enum { idx = 77 }; };
937  template<> struct SixthDensePt<-2, 0,-2> { enum { idx = 78 }; };
938  template<> struct SixthDensePt<-1, 0,-2> { enum { idx = 79 }; };
939  template<> struct SixthDensePt< 0, 0,-2> { enum { idx = 80 }; };
940  template<> struct SixthDensePt< 1, 0,-2> { enum { idx = 81 }; };
941  template<> struct SixthDensePt< 2, 0,-2> { enum { idx = 82 }; };
942  template<> struct SixthDensePt< 3, 0,-2> { enum { idx = 83 }; };
943 
944 
945  template<> struct SixthDensePt<-3, 0,-3> { enum { idx = 84 }; };
946  template<> struct SixthDensePt<-2, 0,-3> { enum { idx = 85 }; };
947  template<> struct SixthDensePt<-1, 0,-3> { enum { idx = 86 }; };
948  template<> struct SixthDensePt< 0, 0,-3> { enum { idx = 87 }; };
949  template<> struct SixthDensePt< 1, 0,-3> { enum { idx = 88 }; };
950  template<> struct SixthDensePt< 2, 0,-3> { enum { idx = 89 }; };
951  template<> struct SixthDensePt< 3, 0,-3> { enum { idx = 90 }; };
952 
953 
954  template<> struct SixthDensePt< 0,-3, 3> { enum { idx = 91 }; };
955  template<> struct SixthDensePt< 0,-2, 3> { enum { idx = 92 }; };
956  template<> struct SixthDensePt< 0,-1, 3> { enum { idx = 93 }; };
957  template<> struct SixthDensePt< 0, 1, 3> { enum { idx = 94 }; };
958  template<> struct SixthDensePt< 0, 2, 3> { enum { idx = 95 }; };
959  template<> struct SixthDensePt< 0, 3, 3> { enum { idx = 96 }; };
960 
961  template<> struct SixthDensePt< 0,-3, 2> { enum { idx = 97 }; };
962  template<> struct SixthDensePt< 0,-2, 2> { enum { idx = 98 }; };
963  template<> struct SixthDensePt< 0,-1, 2> { enum { idx = 99 }; };
964  template<> struct SixthDensePt< 0, 1, 2> { enum { idx = 100 }; };
965  template<> struct SixthDensePt< 0, 2, 2> { enum { idx = 101 }; };
966  template<> struct SixthDensePt< 0, 3, 2> { enum { idx = 102 }; };
967 
968  template<> struct SixthDensePt< 0,-3, 1> { enum { idx = 103 }; };
969  template<> struct SixthDensePt< 0,-2, 1> { enum { idx = 104 }; };
970  template<> struct SixthDensePt< 0,-1, 1> { enum { idx = 105 }; };
971  template<> struct SixthDensePt< 0, 1, 1> { enum { idx = 106 }; };
972  template<> struct SixthDensePt< 0, 2, 1> { enum { idx = 107 }; };
973  template<> struct SixthDensePt< 0, 3, 1> { enum { idx = 108 }; };
974 
975  template<> struct SixthDensePt< 0,-3,-1> { enum { idx = 109 }; };
976  template<> struct SixthDensePt< 0,-2,-1> { enum { idx = 110 }; };
977  template<> struct SixthDensePt< 0,-1,-1> { enum { idx = 111 }; };
978  template<> struct SixthDensePt< 0, 1,-1> { enum { idx = 112 }; };
979  template<> struct SixthDensePt< 0, 2,-1> { enum { idx = 113 }; };
980  template<> struct SixthDensePt< 0, 3,-1> { enum { idx = 114 }; };
981 
982  template<> struct SixthDensePt< 0,-3,-2> { enum { idx = 115 }; };
983  template<> struct SixthDensePt< 0,-2,-2> { enum { idx = 116 }; };
984  template<> struct SixthDensePt< 0,-1,-2> { enum { idx = 117 }; };
985  template<> struct SixthDensePt< 0, 1,-2> { enum { idx = 118 }; };
986  template<> struct SixthDensePt< 0, 2,-2> { enum { idx = 119 }; };
987  template<> struct SixthDensePt< 0, 3,-2> { enum { idx = 120 }; };
988 
989  template<> struct SixthDensePt< 0,-3,-3> { enum { idx = 121 }; };
990  template<> struct SixthDensePt< 0,-2,-3> { enum { idx = 122 }; };
991  template<> struct SixthDensePt< 0,-1,-3> { enum { idx = 123 }; };
992  template<> struct SixthDensePt< 0, 1,-3> { enum { idx = 124 }; };
993  template<> struct SixthDensePt< 0, 2,-3> { enum { idx = 125 }; };
994  template<> struct SixthDensePt< 0, 3,-3> { enum { idx = 126 }; };
995 
996 }
997 
998 
999 template<typename GridType>
1000 class SixthOrderDenseStencil: public BaseStencil<GridType, SixthOrderDenseStencil<GridType> >
1001 {
1002 public:
1005  typedef typename GridType::ValueType ValueType;
1007 
1008  static const int SIZE = 127;
1009 
1010  SixthOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
1011 
1013  template<int i, int j, int k>
1014  unsigned int pos() const { return SixthDensePt<i,j,k>::idx; }
1015 
1016 private:
1017  inline void init(const Coord& ijk)
1018  {
1019  mStencil[SixthDensePt<-3, 3, 0>::idx] = mCache.getValue(ijk.offsetBy(-3, 3, 0));
1020  mStencil[SixthDensePt<-2, 3, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 3, 0));
1021  mStencil[SixthDensePt<-1, 3, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 3, 0));
1022  mStencil[SixthDensePt< 0, 3, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 3, 0));
1023  mStencil[SixthDensePt< 1, 3, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 3, 0));
1024  mStencil[SixthDensePt< 2, 3, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 3, 0));
1025  mStencil[SixthDensePt< 3, 3, 0>::idx] = mCache.getValue(ijk.offsetBy( 3, 3, 0));
1026 
1027  mStencil[SixthDensePt<-3, 2, 0>::idx] = mCache.getValue(ijk.offsetBy(-3, 2, 0));
1028  mStencil[SixthDensePt<-2, 2, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 2, 0));
1029  mStencil[SixthDensePt<-1, 2, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 2, 0));
1030  mStencil[SixthDensePt< 0, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 0));
1031  mStencil[SixthDensePt< 1, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 2, 0));
1032  mStencil[SixthDensePt< 2, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 2, 0));
1033  mStencil[SixthDensePt< 3, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 3, 2, 0));
1034 
1035  mStencil[SixthDensePt<-3, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-3, 1, 0));
1036  mStencil[SixthDensePt<-2, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 1, 0));
1037  mStencil[SixthDensePt<-1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 1, 0));
1038  mStencil[SixthDensePt< 0, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
1039  mStencil[SixthDensePt< 1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 1, 0));
1040  mStencil[SixthDensePt< 2, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 1, 0));
1041  mStencil[SixthDensePt< 3, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 3, 1, 0));
1042 
1043  mStencil[SixthDensePt<-3, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-3, 0, 0));
1044  mStencil[SixthDensePt<-2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 0));
1045  mStencil[SixthDensePt<-1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
1046  mStencil[SixthDensePt< 1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
1047  mStencil[SixthDensePt< 2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 0));
1048  mStencil[SixthDensePt< 3, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 3, 0, 0));
1049 
1050  mStencil[SixthDensePt<-3,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-3,-1, 0));
1051  mStencil[SixthDensePt<-2,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-2,-1, 0));
1052  mStencil[SixthDensePt<-1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1,-1, 0));
1053  mStencil[SixthDensePt< 0,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 0));
1054  mStencil[SixthDensePt< 1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1,-1, 0));
1055  mStencil[SixthDensePt< 2,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 2,-1, 0));
1056  mStencil[SixthDensePt< 3,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 3,-1, 0));
1057 
1058  mStencil[SixthDensePt<-3,-2, 0>::idx] = mCache.getValue(ijk.offsetBy(-3,-2, 0));
1059  mStencil[SixthDensePt<-2,-2, 0>::idx] = mCache.getValue(ijk.offsetBy(-2,-2, 0));
1060  mStencil[SixthDensePt<-1,-2, 0>::idx] = mCache.getValue(ijk.offsetBy(-1,-2, 0));
1061  mStencil[SixthDensePt< 0,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 0));
1062  mStencil[SixthDensePt< 1,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 1,-2, 0));
1063  mStencil[SixthDensePt< 2,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 2,-2, 0));
1064  mStencil[SixthDensePt< 3,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 3,-2, 0));
1065 
1066  mStencil[SixthDensePt<-3,-3, 0>::idx] = mCache.getValue(ijk.offsetBy(-3,-3, 0));
1067  mStencil[SixthDensePt<-2,-3, 0>::idx] = mCache.getValue(ijk.offsetBy(-2,-3, 0));
1068  mStencil[SixthDensePt<-1,-3, 0>::idx] = mCache.getValue(ijk.offsetBy(-1,-3, 0));
1069  mStencil[SixthDensePt< 0,-3, 0>::idx] = mCache.getValue(ijk.offsetBy( 0,-3, 0));
1070  mStencil[SixthDensePt< 1,-3, 0>::idx] = mCache.getValue(ijk.offsetBy( 1,-3, 0));
1071  mStencil[SixthDensePt< 2,-3, 0>::idx] = mCache.getValue(ijk.offsetBy( 2,-3, 0));
1072  mStencil[SixthDensePt< 3,-3, 0>::idx] = mCache.getValue(ijk.offsetBy( 3,-3, 0));
1073 
1074  mStencil[SixthDensePt<-3, 0, 3>::idx] = mCache.getValue(ijk.offsetBy(-3, 0, 3));
1075  mStencil[SixthDensePt<-2, 0, 3>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 3));
1076  mStencil[SixthDensePt<-1, 0, 3>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 3));
1077  mStencil[SixthDensePt< 0, 0, 3>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 3));
1078  mStencil[SixthDensePt< 1, 0, 3>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 3));
1079  mStencil[SixthDensePt< 2, 0, 3>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 3));
1080  mStencil[SixthDensePt< 3, 0, 3>::idx] = mCache.getValue(ijk.offsetBy( 3, 0, 3));
1081 
1082  mStencil[SixthDensePt<-3, 0, 2>::idx] = mCache.getValue(ijk.offsetBy(-3, 0, 2));
1083  mStencil[SixthDensePt<-2, 0, 2>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 2));
1084  mStencil[SixthDensePt<-1, 0, 2>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 2));
1085  mStencil[SixthDensePt< 0, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 2));
1086  mStencil[SixthDensePt< 1, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 2));
1087  mStencil[SixthDensePt< 2, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 2));
1088  mStencil[SixthDensePt< 3, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 3, 0, 2));
1089 
1090  mStencil[SixthDensePt<-3, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-3, 0, 1));
1091  mStencil[SixthDensePt<-2, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 1));
1092  mStencil[SixthDensePt<-1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 1));
1093  mStencil[SixthDensePt< 0, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
1094  mStencil[SixthDensePt< 1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 1));
1095  mStencil[SixthDensePt< 2, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 1));
1096  mStencil[SixthDensePt< 3, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 3, 0, 1));
1097 
1098  mStencil[SixthDensePt<-3, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-3, 0,-1));
1099  mStencil[SixthDensePt<-2, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-2, 0,-1));
1100  mStencil[SixthDensePt<-1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0,-1));
1101  mStencil[SixthDensePt< 0, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0,-1));
1102  mStencil[SixthDensePt< 1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0,-1));
1103  mStencil[SixthDensePt< 2, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 2, 0,-1));
1104  mStencil[SixthDensePt< 3, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 3, 0,-1));
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,-3>::idx] = mCache.getValue(ijk.offsetBy(-3, 0,-3));
1115  mStencil[SixthDensePt<-2, 0,-3>::idx] = mCache.getValue(ijk.offsetBy(-2, 0,-3));
1116  mStencil[SixthDensePt<-1, 0,-3>::idx] = mCache.getValue(ijk.offsetBy(-1, 0,-3));
1117  mStencil[SixthDensePt< 0, 0,-3>::idx] = mCache.getValue(ijk.offsetBy( 0, 0,-3));
1118  mStencil[SixthDensePt< 1, 0,-3>::idx] = mCache.getValue(ijk.offsetBy( 1, 0,-3));
1119  mStencil[SixthDensePt< 2, 0,-3>::idx] = mCache.getValue(ijk.offsetBy( 2, 0,-3));
1120  mStencil[SixthDensePt< 3, 0,-3>::idx] = mCache.getValue(ijk.offsetBy( 3, 0,-3));
1121 
1122  mStencil[SixthDensePt< 0,-3, 3>::idx] = mCache.getValue(ijk.offsetBy( 0,-3, 3));
1123  mStencil[SixthDensePt< 0,-2, 3>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 3));
1124  mStencil[SixthDensePt< 0,-1, 3>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 3));
1125  mStencil[SixthDensePt< 0, 1, 3>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 3));
1126  mStencil[SixthDensePt< 0, 2, 3>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 3));
1127  mStencil[SixthDensePt< 0, 3, 3>::idx] = mCache.getValue(ijk.offsetBy( 0, 3, 3));
1128 
1129  mStencil[SixthDensePt< 0,-3, 2>::idx] = mCache.getValue(ijk.offsetBy( 0,-3, 2));
1130  mStencil[SixthDensePt< 0,-2, 2>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 2));
1131  mStencil[SixthDensePt< 0,-1, 2>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 2));
1132  mStencil[SixthDensePt< 0, 1, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 2));
1133  mStencil[SixthDensePt< 0, 2, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 2));
1134  mStencil[SixthDensePt< 0, 3, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 3, 2));
1135 
1136  mStencil[SixthDensePt< 0,-3, 1>::idx] = mCache.getValue(ijk.offsetBy( 0,-3, 1));
1137  mStencil[SixthDensePt< 0,-2, 1>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 1));
1138  mStencil[SixthDensePt< 0,-1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 1));
1139  mStencil[SixthDensePt< 0, 1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 1));
1140  mStencil[SixthDensePt< 0, 2, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 1));
1141  mStencil[SixthDensePt< 0, 3, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 3, 1));
1142 
1143  mStencil[SixthDensePt< 0,-3,-1>::idx] = mCache.getValue(ijk.offsetBy( 0,-3,-1));
1144  mStencil[SixthDensePt< 0,-2,-1>::idx] = mCache.getValue(ijk.offsetBy( 0,-2,-1));
1145  mStencil[SixthDensePt< 0,-1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0,-1,-1));
1146  mStencil[SixthDensePt< 0, 1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1,-1));
1147  mStencil[SixthDensePt< 0, 2,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 2,-1));
1148  mStencil[SixthDensePt< 0, 3,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 3,-1));
1149 
1150  mStencil[SixthDensePt< 0,-3,-2>::idx] = mCache.getValue(ijk.offsetBy( 0,-3,-2));
1151  mStencil[SixthDensePt< 0,-2,-2>::idx] = mCache.getValue(ijk.offsetBy( 0,-2,-2));
1152  mStencil[SixthDensePt< 0,-1,-2>::idx] = mCache.getValue(ijk.offsetBy( 0,-1,-2));
1153  mStencil[SixthDensePt< 0, 1,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 1,-2));
1154  mStencil[SixthDensePt< 0, 2,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 2,-2));
1155  mStencil[SixthDensePt< 0, 3,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 3,-2));
1156 
1157  mStencil[SixthDensePt< 0,-3,-3>::idx] = mCache.getValue(ijk.offsetBy( 0,-3,-3));
1158  mStencil[SixthDensePt< 0,-2,-3>::idx] = mCache.getValue(ijk.offsetBy( 0,-2,-3));
1159  mStencil[SixthDensePt< 0,-1,-3>::idx] = mCache.getValue(ijk.offsetBy( 0,-1,-3));
1160  mStencil[SixthDensePt< 0, 1,-3>::idx] = mCache.getValue(ijk.offsetBy( 0, 1,-3));
1161  mStencil[SixthDensePt< 0, 2,-3>::idx] = mCache.getValue(ijk.offsetBy( 0, 2,-3));
1162  mStencil[SixthDensePt< 0, 3,-3>::idx] = mCache.getValue(ijk.offsetBy( 0, 3,-3));
1163  }
1164 
1165  template<typename, typename> friend class BaseStencil; // allow base class to call init()
1166  using BaseType::mCache;
1167  using BaseType::mStencil;
1168 };
1169 
1170 
1172 
1173 
1180 template<typename GridType>
1181 class GradStencil: public BaseStencil<GridType, GradStencil<GridType> >
1182 {
1183 public:
1186  typedef typename GridType::ValueType ValueType;
1188 
1189  static const int SIZE = 7;
1190 
1191  GradStencil(const GridType& grid):
1192  BaseType(grid, SIZE),
1193  mInv2Dx(ValueType(0.5 / grid.voxelSize()[0])),
1194  mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1195  {
1196  }
1197 
1198  GradStencil(const GridType& grid, Real dx):
1199  BaseType(grid, SIZE),
1200  mInv2Dx(ValueType(0.5 / dx)),
1201  mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1202  {
1203  }
1204 
1210  inline ValueType normSqGrad() const
1211  {
1212  return mInvDx2 * math::GudonovsNormSqrd(mStencil[0] > 0,
1213  mStencil[0] - mStencil[1],
1214  mStencil[2] - mStencil[0],
1215  mStencil[0] - mStencil[3],
1216  mStencil[4] - mStencil[0],
1217  mStencil[0] - mStencil[5],
1218  mStencil[6] - mStencil[0]);
1219  }
1220 
1226  inline Vec3Type gradient() const
1227  {
1228  return Vec3Type(mStencil[2] - mStencil[1],
1229  mStencil[4] - mStencil[3],
1230  mStencil[6] - mStencil[5])*mInv2Dx;
1231  }
1236  inline Vec3Type gradient(const Vec3Type& V) const
1237  {
1238  return Vec3Type(V[0]>0 ? mStencil[0] - mStencil[1] : mStencil[2] - mStencil[0],
1239  V[1]>0 ? mStencil[0] - mStencil[3] : mStencil[4] - mStencil[0],
1240  V[2]>0 ? mStencil[0] - mStencil[5] : mStencil[6] - mStencil[0])*2*mInv2Dx;
1241  }
1242 
1245  inline ValueType laplacian() const
1246  {
1247  return mInvDx2 * (mStencil[1] + mStencil[2] +
1248  mStencil[3] + mStencil[4] +
1249  mStencil[5] + mStencil[6] - 6*mStencil[0]);
1250  }
1251 
1254  inline bool zeroCrossing() const
1255  {
1256  const BufferType& v = mStencil;
1257  return (v[0]>0 ? (v[1]<0 || v[2]<0 || v[3]<0 || v[4]<0 || v[5]<0 || v[6]<0)
1258  : (v[1]>0 || v[2]>0 || v[3]>0 || v[4]>0 || v[5]>0 || v[6]>0));
1259  }
1260 
1268  inline Vec3Type cpt()
1269  {
1270  const Coord& ijk = BaseType::getCenterCoord();
1271  const ValueType d = ValueType(mStencil[0] * 0.5 * mInvDx2); // distance in voxels / (2dx^2)
1272  return Vec3Type(ijk[0] - d*(mStencil[2] - mStencil[1]),
1273  ijk[1] - d*(mStencil[4] - mStencil[3]),
1274  ijk[2] - d*(mStencil[6] - mStencil[5]));
1275  }
1276 
1277 private:
1278 
1279  inline void init(const Coord& ijk)
1280  {
1281  mStencil[1] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
1282  mStencil[2] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
1283 
1284  mStencil[3] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
1285  mStencil[4] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
1286 
1287  mStencil[5] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
1288  mStencil[6] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
1289  }
1290 
1291  template<typename, typename> friend class BaseStencil; // allow base class to call init()
1292  using BaseType::mCache;
1293  using BaseType::mStencil;
1294  const ValueType mInv2Dx, mInvDx2;
1295 }; // class GradStencil
1296 
1297 
1299 
1300 
1306 template<typename GridType>
1307 class WenoStencil: public BaseStencil<GridType, WenoStencil<GridType> >
1308 {
1309 public:
1312  typedef typename GridType::ValueType ValueType;
1314 
1315  static const int SIZE = 19;
1316 
1317  WenoStencil(const GridType& grid):
1318  BaseType(grid, SIZE),
1319  mDx2(ValueType(math::Pow2(grid.voxelSize()[0]))),
1320  mInv2Dx(ValueType(0.5 / grid.voxelSize()[0])),
1321  mInvDx2(ValueType(1.0 / mDx2))
1322  {
1323  }
1324 
1325  WenoStencil(const GridType& grid, Real dx):
1326  BaseType(grid, SIZE),
1327  mDx2(ValueType(dx * dx)),
1328  mInv2Dx(ValueType(0.5 / dx)),
1329  mInvDx2(ValueType(1.0 / mDx2))
1330  {
1331  }
1332 
1338  inline ValueType normSqGrad() const
1339  {
1340  const BufferType& v = mStencil;
1341 #ifdef DWA_OPENVDB
1342  // SSE optimized
1343  const simd::Float4
1344  v1(v[2]-v[1], v[ 8]-v[ 7], v[14]-v[13], 0),
1345  v2(v[3]-v[2], v[ 9]-v[ 8], v[15]-v[14], 0),
1346  v3(v[0]-v[3], v[ 0]-v[ 9], v[ 0]-v[15], 0),
1347  v4(v[4]-v[0], v[10]-v[ 0], v[16]-v[ 0], 0),
1348  v5(v[5]-v[4], v[11]-v[10], v[17]-v[16], 0),
1349  v6(v[6]-v[5], v[12]-v[11], v[18]-v[17], 0),
1350  dP_m = math::WENO5(v1, v2, v3, v4, v5, mDx2),
1351  dP_p = math::WENO5(v6, v5, v4, v3, v2, mDx2);
1352 
1353  return mInvDx2 * math::GudonovsNormSqrd(mStencil[0] > 0, dP_m, dP_p);
1354 #else
1355  const Real
1356  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),
1357  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),
1358  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),
1359  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),
1360  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),
1361  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);
1362  return static_cast<ValueType>(
1363  mInvDx2*math::GudonovsNormSqrd(v[0]>0,dP_xm,dP_xp,dP_ym,dP_yp,dP_zm,dP_zp));
1364 #endif
1365  }
1366 
1372  inline Vec3Type gradient(const Vec3Type& V) const
1373  {
1374  const BufferType& v = mStencil;
1375  return 2*mInv2Dx * Vec3Type(
1376  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)
1377  : math::WENO5(v[ 6]-v[ 5],v[ 5]-v[ 4],v[ 4]-v[ 0], v[ 0]-v[ 3],v[ 3]-v[ 2],mDx2),
1378  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)
1379  : math::WENO5(v[12]-v[11],v[11]-v[10],v[10]-v[ 0], v[ 0]-v[ 9],v[ 9]-v[ 8],mDx2),
1380  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)
1381  : math::WENO5(v[18]-v[17],v[17]-v[16],v[16]-v[ 0], v[ 0]-v[15],v[15]-v[14],mDx2));
1382  }
1388  inline Vec3Type gradient() const
1389  {
1390  return mInv2Dx * Vec3Type(
1391  mStencil[ 4] - mStencil[ 3],
1392  mStencil[10] - mStencil[ 9],
1393  mStencil[16] - mStencil[15]);
1394  }
1395 
1401  inline ValueType laplacian() const
1402  {
1403  return mInvDx2 * (
1404  mStencil[ 3] + mStencil[ 4] +
1405  mStencil[ 9] + mStencil[10] +
1406  mStencil[15] + mStencil[16] - 6*mStencil[0]);
1407  }
1408 
1411  inline bool zeroCrossing() const
1412  {
1413  const BufferType& v = mStencil;
1414  return (v[ 0]>0 ? (v[ 3]<0 || v[ 4]<0 || v[ 9]<0 || v[10]<0 || v[15]<0 || v[16]<0)
1415  : (v[ 3]>0 || v[ 4]>0 || v[ 9]>0 || v[10]>0 || v[15]>0 || v[16]>0));
1416  }
1417 
1418 private:
1419  inline void init(const Coord& ijk)
1420  {
1421  mStencil[ 1] = mCache.getValue(ijk.offsetBy(-3, 0, 0));
1422  mStencil[ 2] = mCache.getValue(ijk.offsetBy(-2, 0, 0));
1423  mStencil[ 3] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
1424  mStencil[ 4] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
1425  mStencil[ 5] = mCache.getValue(ijk.offsetBy( 2, 0, 0));
1426  mStencil[ 6] = mCache.getValue(ijk.offsetBy( 3, 0, 0));
1427 
1428  mStencil[ 7] = mCache.getValue(ijk.offsetBy( 0, -3, 0));
1429  mStencil[ 8] = mCache.getValue(ijk.offsetBy( 0, -2, 0));
1430  mStencil[ 9] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
1431  mStencil[10] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
1432  mStencil[11] = mCache.getValue(ijk.offsetBy( 0, 2, 0));
1433  mStencil[12] = mCache.getValue(ijk.offsetBy( 0, 3, 0));
1434 
1435  mStencil[13] = mCache.getValue(ijk.offsetBy( 0, 0, -3));
1436  mStencil[14] = mCache.getValue(ijk.offsetBy( 0, 0, -2));
1437  mStencil[15] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
1438  mStencil[16] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
1439  mStencil[17] = mCache.getValue(ijk.offsetBy( 0, 0, 2));
1440  mStencil[18] = mCache.getValue(ijk.offsetBy( 0, 0, 3));
1441  }
1442 
1443  template<typename, typename> friend class BaseStencil; // allow base class to call init()
1444  using BaseType::mCache;
1445  using BaseType::mStencil;
1446  const ValueType mDx2, mInv2Dx, mInvDx2;
1447 }; // class WenoStencil
1448 
1449 
1451 
1452 
1453 template<typename GridType>
1454 class CurvatureStencil: public BaseStencil<GridType, CurvatureStencil<GridType> >
1455 {
1456 public:
1458  typedef typename GridType::ValueType ValueType;
1459 
1460  static const int SIZE = 19;
1461 
1462  CurvatureStencil(const GridType& grid):
1463  BaseType(grid, SIZE),
1464  mInv2Dx(ValueType(0.5 / grid.voxelSize()[0])),
1465  mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1466  {
1467  }
1468 
1469  CurvatureStencil(const GridType& grid, Real dx):
1470  BaseType(grid, SIZE),
1471  mInv2Dx(ValueType(0.5 / dx)),
1472  mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1473  {
1474  }
1475 
1480  inline ValueType meanCurvature()
1481  {
1482  Real alpha, beta;
1483  return this->meanCurvature(alpha, beta) ? ValueType(alpha*mInv2Dx/math::Pow3(beta)) : 0;
1484  }
1485 
1492  inline ValueType meanCurvatureNormGrad()
1493  {
1494  Real alpha, beta;
1495  return this->meanCurvature(alpha, beta) ? ValueType(alpha*mInvDx2/(2*math::Pow2(beta))) : 0;
1496  }
1497 
1503  inline ValueType laplacian() const
1504  {
1505  return mInvDx2 * (
1506  mStencil[1] + mStencil[2] +
1507  mStencil[3] + mStencil[4] +
1508  mStencil[5] + mStencil[6] - 6*mStencil[0]);
1509  }
1510 
1517  {
1518  return math::Vec3<ValueType>(
1519  mStencil[2] - mStencil[1],
1520  mStencil[4] - mStencil[3],
1521  mStencil[6] - mStencil[5])*mInv2Dx;
1522  }
1523 
1524 private:
1525  inline void init(const Coord &ijk)
1526  {
1527  mStencil[ 1] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
1528  mStencil[ 2] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
1529 
1530  mStencil[ 3] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
1531  mStencil[ 4] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
1532 
1533  mStencil[ 5] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
1534  mStencil[ 6] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
1535 
1536  mStencil[ 7] = mCache.getValue(ijk.offsetBy(-1, -1, 0));
1537  mStencil[ 8] = mCache.getValue(ijk.offsetBy( 1, -1, 0));
1538  mStencil[ 9] = mCache.getValue(ijk.offsetBy(-1, 1, 0));
1539  mStencil[10] = mCache.getValue(ijk.offsetBy( 1, 1, 0));
1540 
1541  mStencil[11] = mCache.getValue(ijk.offsetBy(-1, 0, -1));
1542  mStencil[12] = mCache.getValue(ijk.offsetBy( 1, 0, -1));
1543  mStencil[13] = mCache.getValue(ijk.offsetBy(-1, 0, 1));
1544  mStencil[14] = mCache.getValue(ijk.offsetBy( 1, 0, 1));
1545 
1546  mStencil[15] = mCache.getValue(ijk.offsetBy( 0, -1, -1));
1547  mStencil[16] = mCache.getValue(ijk.offsetBy( 0, 1, -1));
1548  mStencil[17] = mCache.getValue(ijk.offsetBy( 0, -1, 1));
1549  mStencil[18] = mCache.getValue(ijk.offsetBy( 0, 1, 1));
1550  }
1551 
1552  inline bool meanCurvature(Real& alpha, Real& beta) const
1553  {
1554  // For performance all finite differences are unscaled wrt dx
1555  const Real
1556  Half(0.5), Quarter(0.25),
1557  Dx = Half * (mStencil[2] - mStencil[1]), Dx2 = Dx * Dx, // * 1/dx
1558  Dy = Half * (mStencil[4] - mStencil[3]), Dy2 = Dy * Dy, // * 1/dx
1559  Dz = Half * (mStencil[6] - mStencil[5]), Dz2 = Dz * Dz, // * 1/dx
1560  normGrad = Dx2 + Dy2 + Dz2;
1561  if (normGrad <= math::Tolerance<Real>::value()) {
1562  alpha = beta = 0;
1563  return false;
1564  }
1565  const Real
1566  Dxx = mStencil[2] - 2 * mStencil[0] + mStencil[1], // * 1/dx2
1567  Dyy = mStencil[4] - 2 * mStencil[0] + mStencil[3], // * 1/dx2
1568  Dzz = mStencil[6] - 2 * mStencil[0] + mStencil[5], // * 1/dx2
1569  Dxy = Quarter * (mStencil[10] - mStencil[ 8] + mStencil[7] - mStencil[ 9]), // * 1/dx2
1570  Dxz = Quarter * (mStencil[14] - mStencil[12] + mStencil[11] - mStencil[13]), // * 1/dx2
1571  Dyz = Quarter * (mStencil[18] - mStencil[16] + mStencil[15] - mStencil[17]); // * 1/dx2
1572  alpha = (Dx2*(Dyy+Dzz)+Dy2*(Dxx+Dzz)+Dz2*(Dxx+Dyy)-2*(Dx*(Dy*Dxy+Dz*Dxz)+Dy*Dz*Dyz));
1573  beta = std::sqrt(normGrad); // * 1/dx
1574  return true;
1575  }
1576 
1577  template<typename, typename> friend class BaseStencil; // allow base class to call init()
1578  using BaseType::mCache;
1579  using BaseType::mStencil;
1580  const ValueType mInv2Dx, mInvDx2;
1581 }; // class CurvatureStencil
1582 
1583 
1585 
1586 
1588 template<typename GridType>
1589 class DenseStencil: public BaseStencil<GridType, DenseStencil<GridType> >
1590 {
1591 public:
1593  typedef typename GridType::ValueType ValueType;
1594 
1595  DenseStencil(const GridType& grid, int halfWidth) :
1596  BaseType(grid, /*size=*/math::Pow3(2 * halfWidth + 1)),
1597  mHalfWidth(halfWidth)
1598  {
1599  assert(halfWidth>0);
1600  }
1601 
1602  inline const ValueType& getCenterValue() const { return mStencil[(mStencil.size()-1)>>1]; }
1603 
1606  inline void moveTo(const Coord& ijk)
1607  {
1608  BaseType::mCenter = ijk;
1609  this->init(ijk);
1610  }
1613  template<typename IterType>
1614  inline void moveTo(const IterType& iter)
1615  {
1616  BaseType::mCenter = iter.getCoord();
1617  this->init(BaseType::mCenter);
1618  }
1619 
1620 private:
1623  inline void init(const Coord& ijk)
1624  {
1625  int n = 0;
1626  for (Coord p=ijk.offsetBy(-mHalfWidth), q=ijk.offsetBy(mHalfWidth); p[0] <= q[0]; ++p[0]) {
1627  for (p[1] = ijk[1]-mHalfWidth; p[1] <= q[1]; ++p[1]) {
1628  for (p[2] = ijk[2]-mHalfWidth; p[2] <= q[2]; ++p[2]) {
1629  mStencil[n++] = mCache.getValue(p);
1630  }
1631  }
1632  }
1633  }
1634 
1635  template<typename, typename> friend class BaseStencil; // allow base class to call init()
1636  using BaseType::mCache;
1637  using BaseType::mStencil;
1638  const int mHalfWidth;
1639 };
1640 
1641 
1642 } // end math namespace
1643 } // namespace OPENVDB_VERSION_NAME
1644 } // end openvdb namespace
1645 
1646 #endif // OPENVDB_MATH_STENCILS_HAS_BEEN_INCLUDED
1647 
1648 // Copyright (c) 2012-2014 DreamWorks Animation LLC
1649 // All rights reserved. This software is distributed under the
1650 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
const AccessorType & accessor() const
Return a const reference to the ValueAccessor associated with this Stencil.
Definition: Stencils.h:197
GridType::ValueType ValueType
Definition: Stencils.h:1593
ValueType min() const
Return the smallest value in the stencil buffer.
Definition: Stencils.h:159
Definition: Stencils.h:54
GridType::ValueType ValueType
Definition: Stencils.h:654
void setValue(const ValueType &value)
Set the value at the specified location relative to the center of the stencil.
Definition: Stencils.h:131
math::Vec3< ValueType > Vec3Type
Definition: Stencils.h:292
ValueType mean() const
Return the mean value of the current stencil.
Definition: Stencils.h:151
Definition: Stencils.h:286
math::Vec3< ValueType > Vec3Type
Definition: Stencils.h:792
FourthOrderDenseStencil(const GridType &grid)
Definition: Stencils.h:659
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:79
GridType::ValueType ValueType
Definition: Stencils.h:239
Definition: Stencils.h:1181
std::vector< ValueType > BufferType
Definition: Stencils.h:60
BaseType::BufferType BufferType
Definition: Stencils.h:524
const ValueType & getValue() const
Return the value at the specified location relative to the center of the stencil. ...
Definition: Stencils.h:124
ValueType normSqGrad() const
Return the norm square of the single-sided upwind gradient (computed via Gudonov's scheme) at the pre...
Definition: Stencils.h:1210
WenoStencil(const GridType &grid, Real dx)
Definition: Stencils.h:1325
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:1307
bool zeroCrossing() const
Definition: Stencils.h:1254
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:663
ValueType laplacian() const
Definition: Stencils.h:1245
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
GridType::ValueType ValueType
Definition: Stencils.h:791
Type Pow2(Type x)
Return .
Definition: Math.h:498
BaseStencil< GridType, CurvatureStencil< GridType > > BaseType
Definition: Stencils.h:1457
BaseStencil< GridType, ThirteenPointStencil< GridType > > BaseType
Definition: Stencils.h:523
math::Vec3< ValueType > Vec3Type
Definition: Stencils.h:655
math::Vec3< ValueType > Vec3Type
Definition: Stencils.h:240
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:47
SixthOrderDenseStencil(const GridType &grid)
Definition: Stencils.h:1010
Vec3Type gradient(const Vec3Type &xyz) const
Return the gradient in world space of the trilinear interpolation kernel.
Definition: Stencils.h:350
CurvatureStencil(const GridType &grid, Real dx)
Definition: Stencils.h:1469
BaseStencil< GridType, BoxStencil< GridType > > BaseType
Definition: Stencils.h:289
math::Vec3< ValueType > Vec3Type
Definition: Stencils.h:1313
double Real
Definition: Types.h:65
Vec3Type cpt()
Compute the closest-point transform to a level set.
Definition: Stencils.h:1268
Vec3Type gradient(const Vec3Type &V) const
Return the first-order upwind gradient corresponding to the direction V.
Definition: Stencils.h:1236
ThirteenPointStencil(const GridType &grid)
Definition: Stencils.h:530
GridType::ValueType ValueType
Definition: Stencils.h:59
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:299
WenoStencil(const GridType &grid)
Definition: Stencils.h:1317
const GridType * mGrid
Definition: Stencils.h:207
BaseStencil< GridType, NineteenPointStencil< GridType > > BaseType
Definition: Stencils.h:789
ValueType meanCurvatureNormGrad()
Definition: Stencils.h:1492
NineteenPointStencil(const GridType &grid)
Definition: Stencils.h:796
BaseType::BufferType BufferType
Definition: Stencils.h:1185
GridType::ValueType ValueType
Definition: Stencils.h:1005
GridType::ValueType ValueType
Definition: Stencils.h:1458
_GridType GridType
Definition: Stencils.h:57
GradStencil(const GridType &grid)
Definition: Stencils.h:1191
BaseStencil< GridType, SecondOrderDenseStencil< GridType > > BaseType
Definition: Stencils.h:445
#define OPENVDB_VERSION_NAME
Definition: version.h:43
DenseStencil(const GridType &grid, int halfWidth)
Definition: Stencils.h:1595
Vec3Type gradient(const Vec3Type &V) const
Definition: Stencils.h:1372
void moveTo(const Vec3R &xyz)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:105
BaseStencil< GridType, DenseStencil< GridType > > BaseType
Definition: Stencils.h:1592
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:456
ValueType meanCurvature()
Return the mean curvature at the previously buffered location.
Definition: Stencils.h:1480
const ValueType & getCenterValue() const
Return the value at the center of the stencil.
Definition: Stencils.h:176
math::Vec3< ValueType > Vec3Type
Definition: Stencils.h:448
AccessorType mCache
Definition: Stencils.h:208
BaseType::BufferType BufferType
Definition: Stencils.h:1004
ValueType laplacian() const
Definition: Stencils.h:1401
Vec3Type gradient() const
Definition: Stencils.h:1388
Definition: Exceptions.h:39
void moveTo(const IterType &iter)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:92
Type Pow3(Type x)
Return .
Definition: Math.h:502
BaseType::BufferType BufferType
Definition: Stencils.h:446
BaseType::BufferType BufferType
Definition: Stencils.h:653
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:534
math::Vec3< ValueType > gradient()
Definition: Stencils.h:1516
GridType::ValueType ValueType
Definition: Stencils.h:1186
ValueType median() const
Return the median value of the current stencil.
Definition: Stencils.h:140
GridType::ValueType ValueType
Definition: Stencils.h:1312
BaseType::BufferType BufferType
Definition: Stencils.h:290
math::Vec3< ValueType > Vec3Type
Definition: Stencils.h:1006
GridType::Ptr meanCurvature(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the mean curvature of the given grid.
Definition: GridOperators.h:1033
math::Vec3< ValueType > Vec3Type
Definition: Stencils.h:526
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:247
bool zeroCrossing() const
Definition: Stencils.h:1411
GridType::ConstAccessor AccessorType
Definition: Stencils.h:62
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:800
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:67
GridType::ValueType ValueType
Definition: Stencils.h:447
const Coord & getCenterCoord() const
Return the coordinates of the center point of the stencil.
Definition: Stencils.h:173
ValueType interpolation(const Vec3Type &xyz) const
Return the trilinear interpolation at the normalized position.
Definition: Stencils.h:322
BaseStencil< GridType, GradStencil< GridType > > BaseType
Definition: Stencils.h:1184
Definition: Stencils.h:1454
Coord mCenter
Definition: Stencils.h:210
int size()
Return the size of the stencil buffer.
Definition: Stencils.h:137
BaseStencil< GridType, SevenPointStencil< GridType > > BaseType
Definition: Stencils.h:237
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:116
bool intersects(const ValueType &isoValue=zeroVal< ValueType >()) const
Return true if the center of the stencil intersects the.
Definition: Stencils.h:303
void moveTo(const IterType &iter)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:1614
const GridType & grid() const
Return a const reference to the grid from which this stencil was constructed.
Definition: Stencils.h:193
ValueType normSqGrad() const
Return the norm-square of the WENO upwind gradient (computed via WENO upwinding and Gudonov's scheme)...
Definition: Stencils.h:1338
BoxStencil(const GridType &grid)
Definition: Stencils.h:295
BaseType::BufferType BufferType
Definition: Stencils.h:790
math::Vec3< ValueType > Vec3Type
Definition: Stencils.h:1187
BaseStencil< GridType, WenoStencil< GridType > > BaseType
Definition: Stencils.h:1310
Coord offsetBy(Int32 dx, Int32 dy, Int32 dz) const
Definition: Coord.h:102
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:1014
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
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:180
BufferType::iterator IterType
Definition: Stencils.h:61
SecondOrderDenseStencil(const GridType &grid)
Definition: Stencils.h:452
const ValueType & getCenterValue() const
Definition: Stencils.h:1602
GridType::TreeType TreeType
Definition: Stencils.h:58
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Dense stencil of a given width.
Definition: Stencils.h:1589
GradStencil(const GridType &grid, Real dx)
Definition: Stencils.h:1198
Vec3Type gradient() const
Return the gradient computed at the previously buffered location by second order central differencing...
Definition: Stencils.h:1226
ValueType max() const
Return the largest value in the stencil buffer.
Definition: Stencils.h:166
BaseType::BufferType BufferType
Definition: Stencils.h:238
CurvatureStencil(const GridType &grid)
Definition: Stencils.h:1462
GridType::ValueType ValueType
Definition: Stencils.h:291
BaseStencil< GridType, SixthOrderDenseStencil< GridType > > BaseType
Definition: Stencils.h:1003
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:1606
SevenPointStencil(const GridType &grid)
Definition: Stencils.h:243
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
BaseStencil(const GridType &grid, int size)
Definition: Stencils.h:201
BaseType::BufferType BufferType
Definition: Stencils.h:1311
BufferType mStencil
Definition: Stencils.h:209
GridType::ValueType ValueType
Definition: Stencils.h:525
Tolerance for floating-point comparison.
Definition: Math.h:116
BaseStencil< GridType, FourthOrderDenseStencil< GridType > > BaseType
Definition: Stencils.h:652
ValueType laplacian() const
Definition: Stencils.h:1503