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