OpenVDB  1.2.0
Stencils.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 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 
65  inline void moveTo(const Coord& ijk)
66  {
67  mCenter = ijk;
68  mStencil[0] = mCache.getValue(ijk);
69  static_cast<StencilType&>(*this).init(mCenter);
70  }
76  template<typename IterType>
77  inline void moveTo(const IterType& iter)
78  {
79  mCenter = iter.getCoord();
80  mStencil[0] = *iter;
81  static_cast<StencilType&>(*this).init(mCenter);
82  }
83 
89  inline ValueType getValue(unsigned int pos = 0) const
90  {
91  assert(pos < mStencil.size());
92  return mStencil[pos];
93  }
94 
96  template<int i, int j, int k>
97  const ValueType& getValue() const
98  {
99  return mStencil[static_cast<const StencilType&>(*this).template pos<i,j,k>()];
100  }
101 
103  template<int i, int j, int k>
104  void setValue(const ValueType& value)
105  {
106  mStencil[static_cast<const StencilType&>(*this).template pos<i,j,k>()] = value;
107  }
108 
110  inline int size() { return mStencil.size(); }
111 
113  inline ValueType median() const
114  {
115  std::vector<ValueType> tmp(mStencil);//local copy
116  assert(!tmp.empty());
117  size_t midpoint = (tmp.size() - 1) >> 1;
118  // Partially sort the vector until the median value is at the midpoint.
119  std::nth_element(tmp.begin(), tmp.begin() + midpoint, tmp.end());
120  return tmp[midpoint];
121  }
122 
124  inline ValueType mean() const
125  {
126  ValueType sum = 0.0;
127  for (int n=0, s=mStencil.size(); n<s; ++n) sum += mStencil[n];
128  return sum / mStencil.size();
129  }
130 
132  inline ValueType min() const
133  {
134  IterType iter = std::min_element(mStencil.begin(), mStencil.end());
135  return *iter;
136  }
137 
139  inline ValueType max() const
140  {
141  IterType iter = std::max_element(mStencil.begin(), mStencil.end());
142  return *iter;
143  }
144 
146  inline const Coord& getCenterCoord() const { return mCenter; }
147 
149  inline const ValueType& getCenterValue() const { return mStencil[0]; }
150 
153  inline bool intersects(const ValueType &isoValue = zeroVal<ValueType>()) const
154  {
155  const bool less = this->getValue< 0, 0, 0>() < isoValue;
156  return (less ^ (this->getValue<-1, 0, 0>() < isoValue)) ||
157  (less ^ (this->getValue< 1, 0, 0>() < isoValue)) ||
158  (less ^ (this->getValue< 0,-1, 0>() < isoValue)) ||
159  (less ^ (this->getValue< 0, 1, 0>() < isoValue)) ||
160  (less ^ (this->getValue< 0, 0,-1>() < isoValue)) ||
161  (less ^ (this->getValue< 0, 0, 1>() < isoValue)) ;
162  }
163 
164 protected:
165  // Constructor is protected to prevent direct instantiation.
166  BaseStencil(const GridType& grid, int size):
167  mCache(grid.getConstAccessor()),
168  mStencil(size)
169  {
170  }
171 
172  typename GridType::ConstAccessor mCache;
175 
176 }; // class BaseStencil
177 
178 
180 
181 
182 namespace { // anonymous namespace for stencil-layout map
183 
184  // the seven point stencil
185  template<int i, int j, int k> struct SevenPt {};
186  template<> struct SevenPt< 0, 0, 0> { enum { idx = 0 }; };
187  template<> struct SevenPt< 1, 0, 0> { enum { idx = 1 }; };
188  template<> struct SevenPt< 0, 1, 0> { enum { idx = 2 }; };
189  template<> struct SevenPt< 0, 0, 1> { enum { idx = 3 }; };
190  template<> struct SevenPt<-1, 0, 0> { enum { idx = 4 }; };
191  template<> struct SevenPt< 0,-1, 0> { enum { idx = 5 }; };
192  template<> struct SevenPt< 0, 0,-1> { enum { idx = 6 }; };
193 
194 }
195 
196 
197 template<typename GridType>
198 class SevenPointStencil: public BaseStencil<GridType, SevenPointStencil<GridType> >
199 {
200 public:
203  typedef typename GridType::ValueType ValueType;
205  static const int SIZE = 7;
206 
207  SevenPointStencil(const GridType& grid): BaseType(grid, SIZE) {}
208 
210  template<int i, int j, int k>
211  unsigned int pos() const { return SevenPt<i,j,k>::idx; }
212 
213 private:
214  inline void init(const Coord& ijk)
215  {
216  BaseType::template setValue<-1, 0, 0>(mCache.getValue(ijk.offsetBy(-1, 0, 0)));
217  BaseType::template setValue< 1, 0, 0>(mCache.getValue(ijk.offsetBy( 1, 0, 0)));
218 
219  BaseType::template setValue< 0,-1, 0>(mCache.getValue(ijk.offsetBy( 0,-1, 0)));
220  BaseType::template setValue< 0, 1, 0>(mCache.getValue(ijk.offsetBy( 0, 1, 0)));
221 
222  BaseType::template setValue< 0, 0,-1>(mCache.getValue(ijk.offsetBy( 0, 0,-1)));
223  BaseType::template setValue< 0, 0, 1>(mCache.getValue(ijk.offsetBy( 0, 0, 1)));
224  }
225 
226  template<typename, typename> friend class BaseStencil; // allow base class to call init()
227  using BaseType::mCache;
228  using BaseType::mStencil;
229 };
230 
231 
233 
234 
235 namespace { // anonymous namespace for stencil-layout map
236 
237  // the dense point stencil
238  template<int i, int j, int k> struct DensePt {};
239  template<> struct DensePt< 0, 0, 0> { enum { idx = 0 }; };
240 
241  template<> struct DensePt< 1, 0, 0> { enum { idx = 1 }; };
242  template<> struct DensePt< 0, 1, 0> { enum { idx = 2 }; };
243  template<> struct DensePt< 0, 0, 1> { enum { idx = 3 }; };
244 
245  template<> struct DensePt<-1, 0, 0> { enum { idx = 4 }; };
246  template<> struct DensePt< 0,-1, 0> { enum { idx = 5 }; };
247  template<> struct DensePt< 0, 0,-1> { enum { idx = 6 }; };
248 
249  template<> struct DensePt<-1,-1, 0> { enum { idx = 7 }; };
250  template<> struct DensePt< 0,-1,-1> { enum { idx = 8 }; };
251  template<> struct DensePt<-1, 0,-1> { enum { idx = 9 }; };
252 
253  template<> struct DensePt< 1,-1, 0> { enum { idx = 10 }; };
254  template<> struct DensePt< 0, 1,-1> { enum { idx = 11 }; };
255  template<> struct DensePt<-1, 0, 1> { enum { idx = 12 }; };
256 
257  template<> struct DensePt<-1, 1, 0> { enum { idx = 13 }; };
258  template<> struct DensePt< 0,-1, 1> { enum { idx = 14 }; };
259  template<> struct DensePt< 1, 0,-1> { enum { idx = 15 }; };
260 
261  template<> struct DensePt< 1, 1, 0> { enum { idx = 16 }; };
262  template<> struct DensePt< 0, 1, 1> { enum { idx = 17 }; };
263  template<> struct DensePt< 1, 0, 1> { enum { idx = 18 }; };
264 
265 }
266 
267 
268 template<typename GridType>
269 class SecondOrderDenseStencil: public BaseStencil<GridType, SecondOrderDenseStencil<GridType> >
270 {
271 public:
274  typedef typename GridType::ValueType ValueType;
276 
277  static const int SIZE = 19;
278 
279  SecondOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
280 
282  template<int i, int j, int k>
283  unsigned int pos() const { return DensePt<i,j,k>::idx; }
284 
285 private:
286  inline void init(const Coord& ijk)
287  {
288  mStencil[DensePt< 1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
289  mStencil[DensePt< 0, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
290  mStencil[DensePt< 0, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
291 
292  mStencil[DensePt<-1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
293  mStencil[DensePt< 0,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
294  mStencil[DensePt< 0, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
295 
296  mStencil[DensePt<-1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, -1, 0));
297  mStencil[DensePt< 1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, -1, 0));
298  mStencil[DensePt<-1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 1, 0));
299  mStencil[DensePt< 1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 1, 0));
300 
301  mStencil[DensePt<-1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, -1));
302  mStencil[DensePt< 1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, -1));
303  mStencil[DensePt<-1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 1));
304  mStencil[DensePt< 1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 1));
305 
306  mStencil[DensePt< 0,-1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, -1, -1));
307  mStencil[DensePt< 0, 1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, -1));
308  mStencil[DensePt< 0,-1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, -1, 1));
309  mStencil[DensePt< 0, 1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 1));
310  }
311 
312  template<typename, typename> friend class BaseStencil; // allow base class to call init()
313  using BaseType::mCache;
314  using BaseType::mStencil;
315 };
316 
317 
319 
320 
321 namespace { // anonymous namespace for stencil-layout map
322 
323  // the dense point stencil
324  template<int i, int j, int k> struct ThirteenPt {};
325  template<> struct ThirteenPt< 0, 0, 0> { enum { idx = 0 }; };
326 
327  template<> struct ThirteenPt< 1, 0, 0> { enum { idx = 1 }; };
328  template<> struct ThirteenPt< 0, 1, 0> { enum { idx = 2 }; };
329  template<> struct ThirteenPt< 0, 0, 1> { enum { idx = 3 }; };
330 
331  template<> struct ThirteenPt<-1, 0, 0> { enum { idx = 4 }; };
332  template<> struct ThirteenPt< 0,-1, 0> { enum { idx = 5 }; };
333  template<> struct ThirteenPt< 0, 0,-1> { enum { idx = 6 }; };
334 
335  template<> struct ThirteenPt< 2, 0, 0> { enum { idx = 7 }; };
336  template<> struct ThirteenPt< 0, 2, 0> { enum { idx = 8 }; };
337  template<> struct ThirteenPt< 0, 0, 2> { enum { idx = 9 }; };
338 
339  template<> struct ThirteenPt<-2, 0, 0> { enum { idx = 10 }; };
340  template<> struct ThirteenPt< 0,-2, 0> { enum { idx = 11 }; };
341  template<> struct ThirteenPt< 0, 0,-2> { enum { idx = 12 }; };
342 
343 }
344 
345 
346 template<typename GridType>
347 class ThirteenPointStencil: public BaseStencil<GridType, ThirteenPointStencil<GridType> >
348 {
349 public:
352  typedef typename GridType::ValueType ValueType;
354 
355  static const int SIZE = 13;
356 
357  ThirteenPointStencil(const GridType& grid): BaseType(grid, SIZE) {}
358 
360  template<int i, int j, int k>
361  unsigned int pos() const { return ThirteenPt<i,j,k>::idx; }
362 
363 private:
364  inline void init(const Coord& ijk)
365  {
366  mStencil[ThirteenPt< 2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 0));
367  mStencil[ThirteenPt< 1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
368  mStencil[ThirteenPt<-1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
369  mStencil[ThirteenPt<-2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 0));
370 
371  mStencil[ThirteenPt< 0, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 0));
372  mStencil[ThirteenPt< 0, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
373  mStencil[ThirteenPt< 0,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
374  mStencil[ThirteenPt< 0,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -2, 0));
375 
376  mStencil[ThirteenPt< 0, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 2));
377  mStencil[ThirteenPt< 0, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
378  mStencil[ThirteenPt< 0, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
379  mStencil[ThirteenPt< 0, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -2));
380  }
381 
382  template<typename, typename> friend class BaseStencil; // allow base class to call init()
383  using BaseType::mCache;
384  using BaseType::mStencil;
385 };
386 
387 
389 
390 
391 namespace { // anonymous namespace for stencil-layout map
392 
393  // the 4th-order dense point stencil
394  template<int i, int j, int k> struct FourthDensePt {};
395  template<> struct FourthDensePt< 0, 0, 0> { enum { idx = 0 }; };
396 
397  template<> struct FourthDensePt<-2, 2, 0> { enum { idx = 1 }; };
398  template<> struct FourthDensePt<-1, 2, 0> { enum { idx = 2 }; };
399  template<> struct FourthDensePt< 0, 2, 0> { enum { idx = 3 }; };
400  template<> struct FourthDensePt< 1, 2, 0> { enum { idx = 4 }; };
401  template<> struct FourthDensePt< 2, 2, 0> { enum { idx = 5 }; };
402 
403  template<> struct FourthDensePt<-2, 1, 0> { enum { idx = 6 }; };
404  template<> struct FourthDensePt<-1, 1, 0> { enum { idx = 7 }; };
405  template<> struct FourthDensePt< 0, 1, 0> { enum { idx = 8 }; };
406  template<> struct FourthDensePt< 1, 1, 0> { enum { idx = 9 }; };
407  template<> struct FourthDensePt< 2, 1, 0> { enum { idx = 10 }; };
408 
409  template<> struct FourthDensePt<-2, 0, 0> { enum { idx = 11 }; };
410  template<> struct FourthDensePt<-1, 0, 0> { enum { idx = 12 }; };
411  template<> struct FourthDensePt< 1, 0, 0> { enum { idx = 13 }; };
412  template<> struct FourthDensePt< 2, 0, 0> { enum { idx = 14 }; };
413 
414  template<> struct FourthDensePt<-2,-1, 0> { enum { idx = 15 }; };
415  template<> struct FourthDensePt<-1,-1, 0> { enum { idx = 16 }; };
416  template<> struct FourthDensePt< 0,-1, 0> { enum { idx = 17 }; };
417  template<> struct FourthDensePt< 1,-1, 0> { enum { idx = 18 }; };
418  template<> struct FourthDensePt< 2,-1, 0> { enum { idx = 19 }; };
419 
420  template<> struct FourthDensePt<-2,-2, 0> { enum { idx = 20 }; };
421  template<> struct FourthDensePt<-1,-2, 0> { enum { idx = 21 }; };
422  template<> struct FourthDensePt< 0,-2, 0> { enum { idx = 22 }; };
423  template<> struct FourthDensePt< 1,-2, 0> { enum { idx = 23 }; };
424  template<> struct FourthDensePt< 2,-2, 0> { enum { idx = 24 }; };
425 
426 
427  template<> struct FourthDensePt<-2, 0, 2> { enum { idx = 25 }; };
428  template<> struct FourthDensePt<-1, 0, 2> { enum { idx = 26 }; };
429  template<> struct FourthDensePt< 0, 0, 2> { enum { idx = 27 }; };
430  template<> struct FourthDensePt< 1, 0, 2> { enum { idx = 28 }; };
431  template<> struct FourthDensePt< 2, 0, 2> { enum { idx = 29 }; };
432 
433  template<> struct FourthDensePt<-2, 0, 1> { enum { idx = 30 }; };
434  template<> struct FourthDensePt<-1, 0, 1> { enum { idx = 31 }; };
435  template<> struct FourthDensePt< 0, 0, 1> { enum { idx = 32 }; };
436  template<> struct FourthDensePt< 1, 0, 1> { enum { idx = 33 }; };
437  template<> struct FourthDensePt< 2, 0, 1> { enum { idx = 34 }; };
438 
439  template<> struct FourthDensePt<-2, 0,-1> { enum { idx = 35 }; };
440  template<> struct FourthDensePt<-1, 0,-1> { enum { idx = 36 }; };
441  template<> struct FourthDensePt< 0, 0,-1> { enum { idx = 37 }; };
442  template<> struct FourthDensePt< 1, 0,-1> { enum { idx = 38 }; };
443  template<> struct FourthDensePt< 2, 0,-1> { enum { idx = 39 }; };
444 
445  template<> struct FourthDensePt<-2, 0,-2> { enum { idx = 40 }; };
446  template<> struct FourthDensePt<-1, 0,-2> { enum { idx = 41 }; };
447  template<> struct FourthDensePt< 0, 0,-2> { enum { idx = 42 }; };
448  template<> struct FourthDensePt< 1, 0,-2> { enum { idx = 43 }; };
449  template<> struct FourthDensePt< 2, 0,-2> { enum { idx = 44 }; };
450 
451 
452  template<> struct FourthDensePt< 0,-2, 2> { enum { idx = 45 }; };
453  template<> struct FourthDensePt< 0,-1, 2> { enum { idx = 46 }; };
454  template<> struct FourthDensePt< 0, 1, 2> { enum { idx = 47 }; };
455  template<> struct FourthDensePt< 0, 2, 2> { enum { idx = 48 }; };
456 
457  template<> struct FourthDensePt< 0,-2, 1> { enum { idx = 49 }; };
458  template<> struct FourthDensePt< 0,-1, 1> { enum { idx = 50 }; };
459  template<> struct FourthDensePt< 0, 1, 1> { enum { idx = 51 }; };
460  template<> struct FourthDensePt< 0, 2, 1> { enum { idx = 52 }; };
461 
462  template<> struct FourthDensePt< 0,-2,-1> { enum { idx = 53 }; };
463  template<> struct FourthDensePt< 0,-1,-1> { enum { idx = 54 }; };
464  template<> struct FourthDensePt< 0, 1,-1> { enum { idx = 55 }; };
465  template<> struct FourthDensePt< 0, 2,-1> { enum { idx = 56 }; };
466 
467  template<> struct FourthDensePt< 0,-2,-2> { enum { idx = 57 }; };
468  template<> struct FourthDensePt< 0,-1,-2> { enum { idx = 58 }; };
469  template<> struct FourthDensePt< 0, 1,-2> { enum { idx = 59 }; };
470  template<> struct FourthDensePt< 0, 2,-2> { enum { idx = 60 }; };
471 
472 }
473 
474 
475 template<typename GridType>
476 class FourthOrderDenseStencil: public BaseStencil<GridType, FourthOrderDenseStencil<GridType> >
477 {
478 public:
481  typedef typename GridType::ValueType ValueType;
483 
484  static const int SIZE = 61;
485 
486  FourthOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
487 
489  template<int i, int j, int k>
490  unsigned int pos() const { return FourthDensePt<i,j,k>::idx; }
491 
492 private:
493  inline void init(const Coord& ijk)
494  {
495  mStencil[FourthDensePt<-2, 2, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 2, 0));
496  mStencil[FourthDensePt<-1, 2, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 2, 0));
497  mStencil[FourthDensePt< 0, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 0));
498  mStencil[FourthDensePt< 1, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 2, 0));
499  mStencil[FourthDensePt< 2, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 2, 0));
500 
501  mStencil[FourthDensePt<-2, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 1, 0));
502  mStencil[FourthDensePt<-1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 1, 0));
503  mStencil[FourthDensePt< 0, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
504  mStencil[FourthDensePt< 1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 1, 0));
505  mStencil[FourthDensePt< 2, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 1, 0));
506 
507  mStencil[FourthDensePt<-2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 0));
508  mStencil[FourthDensePt<-1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
509  mStencil[FourthDensePt< 1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
510  mStencil[FourthDensePt< 2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 0));
511 
512  mStencil[FourthDensePt<-2,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-2,-1, 0));
513  mStencil[FourthDensePt<-1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1,-1, 0));
514  mStencil[FourthDensePt< 0,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 0));
515  mStencil[FourthDensePt< 1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1,-1, 0));
516  mStencil[FourthDensePt< 2,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 2,-1, 0));
517 
518  mStencil[FourthDensePt<-2,-2, 0>::idx] = mCache.getValue(ijk.offsetBy(-2,-2, 0));
519  mStencil[FourthDensePt<-1,-2, 0>::idx] = mCache.getValue(ijk.offsetBy(-1,-2, 0));
520  mStencil[FourthDensePt< 0,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 0));
521  mStencil[FourthDensePt< 1,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 1,-2, 0));
522  mStencil[FourthDensePt< 2,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 2,-2, 0));
523 
524  mStencil[FourthDensePt<-2, 0, 2>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 2));
525  mStencil[FourthDensePt<-1, 0, 2>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 2));
526  mStencil[FourthDensePt< 0, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 2));
527  mStencil[FourthDensePt< 1, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 2));
528  mStencil[FourthDensePt< 2, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 2));
529 
530  mStencil[FourthDensePt<-2, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 1));
531  mStencil[FourthDensePt<-1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 1));
532  mStencil[FourthDensePt< 0, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
533  mStencil[FourthDensePt< 1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 1));
534  mStencil[FourthDensePt< 2, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 1));
535 
536  mStencil[FourthDensePt<-2, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-2, 0,-1));
537  mStencil[FourthDensePt<-1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0,-1));
538  mStencil[FourthDensePt< 0, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0,-1));
539  mStencil[FourthDensePt< 1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0,-1));
540  mStencil[FourthDensePt< 2, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 2, 0,-1));
541 
542  mStencil[FourthDensePt<-2, 0,-2>::idx] = mCache.getValue(ijk.offsetBy(-2, 0,-2));
543  mStencil[FourthDensePt<-1, 0,-2>::idx] = mCache.getValue(ijk.offsetBy(-1, 0,-2));
544  mStencil[FourthDensePt< 0, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0,-2));
545  mStencil[FourthDensePt< 1, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 1, 0,-2));
546  mStencil[FourthDensePt< 2, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 2, 0,-2));
547 
548 
549  mStencil[FourthDensePt< 0,-2, 2>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 2));
550  mStencil[FourthDensePt< 0,-1, 2>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 2));
551  mStencil[FourthDensePt< 0, 1, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 2));
552  mStencil[FourthDensePt< 0, 2, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 2));
553 
554  mStencil[FourthDensePt< 0,-2, 1>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 1));
555  mStencil[FourthDensePt< 0,-1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 1));
556  mStencil[FourthDensePt< 0, 1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 1));
557  mStencil[FourthDensePt< 0, 2, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 1));
558 
559  mStencil[FourthDensePt< 0,-2,-1>::idx] = mCache.getValue(ijk.offsetBy( 0,-2,-1));
560  mStencil[FourthDensePt< 0,-1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0,-1,-1));
561  mStencil[FourthDensePt< 0, 1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1,-1));
562  mStencil[FourthDensePt< 0, 2,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 2,-1));
563 
564  mStencil[FourthDensePt< 0,-2,-2>::idx] = mCache.getValue(ijk.offsetBy( 0,-2,-2));
565  mStencil[FourthDensePt< 0,-1,-2>::idx] = mCache.getValue(ijk.offsetBy( 0,-1,-2));
566  mStencil[FourthDensePt< 0, 1,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 1,-2));
567  mStencil[FourthDensePt< 0, 2,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 2,-2));
568  }
569 
570  template<typename, typename> friend class BaseStencil; // allow base class to call init()
571  using BaseType::mCache;
572  using BaseType::mStencil;
573 };
574 
575 
577 
578 
579 namespace { // anonymous namespace for stencil-layout map
580 
581  // the dense point stencil
582  template<int i, int j, int k> struct NineteenPt {};
583  template<> struct NineteenPt< 0, 0, 0> { enum { idx = 0 }; };
584 
585  template<> struct NineteenPt< 1, 0, 0> { enum { idx = 1 }; };
586  template<> struct NineteenPt< 0, 1, 0> { enum { idx = 2 }; };
587  template<> struct NineteenPt< 0, 0, 1> { enum { idx = 3 }; };
588 
589  template<> struct NineteenPt<-1, 0, 0> { enum { idx = 4 }; };
590  template<> struct NineteenPt< 0,-1, 0> { enum { idx = 5 }; };
591  template<> struct NineteenPt< 0, 0,-1> { enum { idx = 6 }; };
592 
593  template<> struct NineteenPt< 2, 0, 0> { enum { idx = 7 }; };
594  template<> struct NineteenPt< 0, 2, 0> { enum { idx = 8 }; };
595  template<> struct NineteenPt< 0, 0, 2> { enum { idx = 9 }; };
596 
597  template<> struct NineteenPt<-2, 0, 0> { enum { idx = 10 }; };
598  template<> struct NineteenPt< 0,-2, 0> { enum { idx = 11 }; };
599  template<> struct NineteenPt< 0, 0,-2> { enum { idx = 12 }; };
600 
601  template<> struct NineteenPt< 3, 0, 0> { enum { idx = 13 }; };
602  template<> struct NineteenPt< 0, 3, 0> { enum { idx = 14 }; };
603  template<> struct NineteenPt< 0, 0, 3> { enum { idx = 15 }; };
604 
605  template<> struct NineteenPt<-3, 0, 0> { enum { idx = 16 }; };
606  template<> struct NineteenPt< 0,-3, 0> { enum { idx = 17 }; };
607  template<> struct NineteenPt< 0, 0,-3> { enum { idx = 18 }; };
608 
609 }
610 
611 
612 template<typename GridType>
613 class NineteenPointStencil: public BaseStencil<GridType, NineteenPointStencil<GridType> >
614 {
615 public:
618  typedef typename GridType::ValueType ValueType;
620 
621  static const int SIZE = 19;
622 
623  NineteenPointStencil(const GridType& grid): BaseType(grid, SIZE) {}
624 
626  template<int i, int j, int k>
627  unsigned int pos() const { return NineteenPt<i,j,k>::idx; }
628 
629 private:
630  inline void init(const Coord& ijk)
631  {
632  mStencil[NineteenPt< 3, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 3, 0, 0));
633  mStencil[NineteenPt< 2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 0));
634  mStencil[NineteenPt< 1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
635  mStencil[NineteenPt<-1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
636  mStencil[NineteenPt<-2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 0));
637  mStencil[NineteenPt<-3, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-3, 0, 0));
638 
639  mStencil[NineteenPt< 0, 3, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 3, 0));
640  mStencil[NineteenPt< 0, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 0));
641  mStencil[NineteenPt< 0, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
642  mStencil[NineteenPt< 0,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
643  mStencil[NineteenPt< 0,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -2, 0));
644  mStencil[NineteenPt< 0,-3, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, -3, 0));
645 
646  mStencil[NineteenPt< 0, 0, 3>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 3));
647  mStencil[NineteenPt< 0, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 2));
648  mStencil[NineteenPt< 0, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
649  mStencil[NineteenPt< 0, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
650  mStencil[NineteenPt< 0, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -2));
651  mStencil[NineteenPt< 0, 0,-3>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, -3));
652  }
653 
654  template<typename, typename> friend class BaseStencil; // allow base class to call init()
655  using BaseType::mCache;
656  using BaseType::mStencil;
657 };
658 
659 
661 
662 
663 namespace { // anonymous namespace for stencil-layout map
664 
665  // the 4th-order dense point stencil
666  template<int i, int j, int k> struct SixthDensePt { };
667  template<> struct SixthDensePt< 0, 0, 0> { enum { idx = 0 }; };
668 
669  template<> struct SixthDensePt<-3, 3, 0> { enum { idx = 1 }; };
670  template<> struct SixthDensePt<-2, 3, 0> { enum { idx = 2 }; };
671  template<> struct SixthDensePt<-1, 3, 0> { enum { idx = 3 }; };
672  template<> struct SixthDensePt< 0, 3, 0> { enum { idx = 4 }; };
673  template<> struct SixthDensePt< 1, 3, 0> { enum { idx = 5 }; };
674  template<> struct SixthDensePt< 2, 3, 0> { enum { idx = 6 }; };
675  template<> struct SixthDensePt< 3, 3, 0> { enum { idx = 7 }; };
676 
677  template<> struct SixthDensePt<-3, 2, 0> { enum { idx = 8 }; };
678  template<> struct SixthDensePt<-2, 2, 0> { enum { idx = 9 }; };
679  template<> struct SixthDensePt<-1, 2, 0> { enum { idx = 10 }; };
680  template<> struct SixthDensePt< 0, 2, 0> { enum { idx = 11 }; };
681  template<> struct SixthDensePt< 1, 2, 0> { enum { idx = 12 }; };
682  template<> struct SixthDensePt< 2, 2, 0> { enum { idx = 13 }; };
683  template<> struct SixthDensePt< 3, 2, 0> { enum { idx = 14 }; };
684 
685  template<> struct SixthDensePt<-3, 1, 0> { enum { idx = 15 }; };
686  template<> struct SixthDensePt<-2, 1, 0> { enum { idx = 16 }; };
687  template<> struct SixthDensePt<-1, 1, 0> { enum { idx = 17 }; };
688  template<> struct SixthDensePt< 0, 1, 0> { enum { idx = 18 }; };
689  template<> struct SixthDensePt< 1, 1, 0> { enum { idx = 19 }; };
690  template<> struct SixthDensePt< 2, 1, 0> { enum { idx = 20 }; };
691  template<> struct SixthDensePt< 3, 1, 0> { enum { idx = 21 }; };
692 
693  template<> struct SixthDensePt<-3, 0, 0> { enum { idx = 22 }; };
694  template<> struct SixthDensePt<-2, 0, 0> { enum { idx = 23 }; };
695  template<> struct SixthDensePt<-1, 0, 0> { enum { idx = 24 }; };
696  template<> struct SixthDensePt< 1, 0, 0> { enum { idx = 25 }; };
697  template<> struct SixthDensePt< 2, 0, 0> { enum { idx = 26 }; };
698  template<> struct SixthDensePt< 3, 0, 0> { enum { idx = 27 }; };
699 
700 
701  template<> struct SixthDensePt<-3,-1, 0> { enum { idx = 28 }; };
702  template<> struct SixthDensePt<-2,-1, 0> { enum { idx = 29 }; };
703  template<> struct SixthDensePt<-1,-1, 0> { enum { idx = 30 }; };
704  template<> struct SixthDensePt< 0,-1, 0> { enum { idx = 31 }; };
705  template<> struct SixthDensePt< 1,-1, 0> { enum { idx = 32 }; };
706  template<> struct SixthDensePt< 2,-1, 0> { enum { idx = 33 }; };
707  template<> struct SixthDensePt< 3,-1, 0> { enum { idx = 34 }; };
708 
709 
710  template<> struct SixthDensePt<-3,-2, 0> { enum { idx = 35 }; };
711  template<> struct SixthDensePt<-2,-2, 0> { enum { idx = 36 }; };
712  template<> struct SixthDensePt<-1,-2, 0> { enum { idx = 37 }; };
713  template<> struct SixthDensePt< 0,-2, 0> { enum { idx = 38 }; };
714  template<> struct SixthDensePt< 1,-2, 0> { enum { idx = 39 }; };
715  template<> struct SixthDensePt< 2,-2, 0> { enum { idx = 40 }; };
716  template<> struct SixthDensePt< 3,-2, 0> { enum { idx = 41 }; };
717 
718 
719  template<> struct SixthDensePt<-3,-3, 0> { enum { idx = 42 }; };
720  template<> struct SixthDensePt<-2,-3, 0> { enum { idx = 43 }; };
721  template<> struct SixthDensePt<-1,-3, 0> { enum { idx = 44 }; };
722  template<> struct SixthDensePt< 0,-3, 0> { enum { idx = 45 }; };
723  template<> struct SixthDensePt< 1,-3, 0> { enum { idx = 46 }; };
724  template<> struct SixthDensePt< 2,-3, 0> { enum { idx = 47 }; };
725  template<> struct SixthDensePt< 3,-3, 0> { enum { idx = 48 }; };
726 
727 
728  template<> struct SixthDensePt<-3, 0, 3> { enum { idx = 49 }; };
729  template<> struct SixthDensePt<-2, 0, 3> { enum { idx = 50 }; };
730  template<> struct SixthDensePt<-1, 0, 3> { enum { idx = 51 }; };
731  template<> struct SixthDensePt< 0, 0, 3> { enum { idx = 52 }; };
732  template<> struct SixthDensePt< 1, 0, 3> { enum { idx = 53 }; };
733  template<> struct SixthDensePt< 2, 0, 3> { enum { idx = 54 }; };
734  template<> struct SixthDensePt< 3, 0, 3> { enum { idx = 55 }; };
735 
736 
737  template<> struct SixthDensePt<-3, 0, 2> { enum { idx = 56 }; };
738  template<> struct SixthDensePt<-2, 0, 2> { enum { idx = 57 }; };
739  template<> struct SixthDensePt<-1, 0, 2> { enum { idx = 58 }; };
740  template<> struct SixthDensePt< 0, 0, 2> { enum { idx = 59 }; };
741  template<> struct SixthDensePt< 1, 0, 2> { enum { idx = 60 }; };
742  template<> struct SixthDensePt< 2, 0, 2> { enum { idx = 61 }; };
743  template<> struct SixthDensePt< 3, 0, 2> { enum { idx = 62 }; };
744 
745  template<> struct SixthDensePt<-3, 0, 1> { enum { idx = 63 }; };
746  template<> struct SixthDensePt<-2, 0, 1> { enum { idx = 64 }; };
747  template<> struct SixthDensePt<-1, 0, 1> { enum { idx = 65 }; };
748  template<> struct SixthDensePt< 0, 0, 1> { enum { idx = 66 }; };
749  template<> struct SixthDensePt< 1, 0, 1> { enum { idx = 67 }; };
750  template<> struct SixthDensePt< 2, 0, 1> { enum { idx = 68 }; };
751  template<> struct SixthDensePt< 3, 0, 1> { enum { idx = 69 }; };
752 
753 
754  template<> struct SixthDensePt<-3, 0,-1> { enum { idx = 70 }; };
755  template<> struct SixthDensePt<-2, 0,-1> { enum { idx = 71 }; };
756  template<> struct SixthDensePt<-1, 0,-1> { enum { idx = 72 }; };
757  template<> struct SixthDensePt< 0, 0,-1> { enum { idx = 73 }; };
758  template<> struct SixthDensePt< 1, 0,-1> { enum { idx = 74 }; };
759  template<> struct SixthDensePt< 2, 0,-1> { enum { idx = 75 }; };
760  template<> struct SixthDensePt< 3, 0,-1> { enum { idx = 76 }; };
761 
762 
763  template<> struct SixthDensePt<-3, 0,-2> { enum { idx = 77 }; };
764  template<> struct SixthDensePt<-2, 0,-2> { enum { idx = 78 }; };
765  template<> struct SixthDensePt<-1, 0,-2> { enum { idx = 79 }; };
766  template<> struct SixthDensePt< 0, 0,-2> { enum { idx = 80 }; };
767  template<> struct SixthDensePt< 1, 0,-2> { enum { idx = 81 }; };
768  template<> struct SixthDensePt< 2, 0,-2> { enum { idx = 82 }; };
769  template<> struct SixthDensePt< 3, 0,-2> { enum { idx = 83 }; };
770 
771 
772  template<> struct SixthDensePt<-3, 0,-3> { enum { idx = 84 }; };
773  template<> struct SixthDensePt<-2, 0,-3> { enum { idx = 85 }; };
774  template<> struct SixthDensePt<-1, 0,-3> { enum { idx = 86 }; };
775  template<> struct SixthDensePt< 0, 0,-3> { enum { idx = 87 }; };
776  template<> struct SixthDensePt< 1, 0,-3> { enum { idx = 88 }; };
777  template<> struct SixthDensePt< 2, 0,-3> { enum { idx = 89 }; };
778  template<> struct SixthDensePt< 3, 0,-3> { enum { idx = 90 }; };
779 
780 
781  template<> struct SixthDensePt< 0,-3, 3> { enum { idx = 91 }; };
782  template<> struct SixthDensePt< 0,-2, 3> { enum { idx = 92 }; };
783  template<> struct SixthDensePt< 0,-1, 3> { enum { idx = 93 }; };
784  template<> struct SixthDensePt< 0, 1, 3> { enum { idx = 94 }; };
785  template<> struct SixthDensePt< 0, 2, 3> { enum { idx = 95 }; };
786  template<> struct SixthDensePt< 0, 3, 3> { enum { idx = 96 }; };
787 
788  template<> struct SixthDensePt< 0,-3, 2> { enum { idx = 97 }; };
789  template<> struct SixthDensePt< 0,-2, 2> { enum { idx = 98 }; };
790  template<> struct SixthDensePt< 0,-1, 2> { enum { idx = 99 }; };
791  template<> struct SixthDensePt< 0, 1, 2> { enum { idx = 100 }; };
792  template<> struct SixthDensePt< 0, 2, 2> { enum { idx = 101 }; };
793  template<> struct SixthDensePt< 0, 3, 2> { enum { idx = 102 }; };
794 
795  template<> struct SixthDensePt< 0,-3, 1> { enum { idx = 103 }; };
796  template<> struct SixthDensePt< 0,-2, 1> { enum { idx = 104 }; };
797  template<> struct SixthDensePt< 0,-1, 1> { enum { idx = 105 }; };
798  template<> struct SixthDensePt< 0, 1, 1> { enum { idx = 106 }; };
799  template<> struct SixthDensePt< 0, 2, 1> { enum { idx = 107 }; };
800  template<> struct SixthDensePt< 0, 3, 1> { enum { idx = 108 }; };
801 
802  template<> struct SixthDensePt< 0,-3,-1> { enum { idx = 109 }; };
803  template<> struct SixthDensePt< 0,-2,-1> { enum { idx = 110 }; };
804  template<> struct SixthDensePt< 0,-1,-1> { enum { idx = 111 }; };
805  template<> struct SixthDensePt< 0, 1,-1> { enum { idx = 112 }; };
806  template<> struct SixthDensePt< 0, 2,-1> { enum { idx = 113 }; };
807  template<> struct SixthDensePt< 0, 3,-1> { enum { idx = 114 }; };
808 
809  template<> struct SixthDensePt< 0,-3,-2> { enum { idx = 115 }; };
810  template<> struct SixthDensePt< 0,-2,-2> { enum { idx = 116 }; };
811  template<> struct SixthDensePt< 0,-1,-2> { enum { idx = 117 }; };
812  template<> struct SixthDensePt< 0, 1,-2> { enum { idx = 118 }; };
813  template<> struct SixthDensePt< 0, 2,-2> { enum { idx = 119 }; };
814  template<> struct SixthDensePt< 0, 3,-2> { enum { idx = 120 }; };
815 
816  template<> struct SixthDensePt< 0,-3,-3> { enum { idx = 121 }; };
817  template<> struct SixthDensePt< 0,-2,-3> { enum { idx = 122 }; };
818  template<> struct SixthDensePt< 0,-1,-3> { enum { idx = 123 }; };
819  template<> struct SixthDensePt< 0, 1,-3> { enum { idx = 124 }; };
820  template<> struct SixthDensePt< 0, 2,-3> { enum { idx = 125 }; };
821  template<> struct SixthDensePt< 0, 3,-3> { enum { idx = 126 }; };
822 
823 }
824 
825 
826 template<typename GridType>
827 class SixthOrderDenseStencil: public BaseStencil<GridType, SixthOrderDenseStencil<GridType> >
828 {
829 public:
832  typedef typename GridType::ValueType ValueType;
834 
835  static const int SIZE = 127;
836 
837  SixthOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
838 
840  template<int i, int j, int k>
841  unsigned int pos() const { return SixthDensePt<i,j,k>::idx; }
842 
843 private:
844  inline void init(const Coord& ijk)
845  {
846  mStencil[SixthDensePt<-3, 3, 0>::idx] = mCache.getValue(ijk.offsetBy(-3, 3, 0));
847  mStencil[SixthDensePt<-2, 3, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 3, 0));
848  mStencil[SixthDensePt<-1, 3, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 3, 0));
849  mStencil[SixthDensePt< 0, 3, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 3, 0));
850  mStencil[SixthDensePt< 1, 3, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 3, 0));
851  mStencil[SixthDensePt< 2, 3, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 3, 0));
852  mStencil[SixthDensePt< 3, 3, 0>::idx] = mCache.getValue(ijk.offsetBy( 3, 3, 0));
853 
854  mStencil[SixthDensePt<-3, 2, 0>::idx] = mCache.getValue(ijk.offsetBy(-3, 2, 0));
855  mStencil[SixthDensePt<-2, 2, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 2, 0));
856  mStencil[SixthDensePt<-1, 2, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 2, 0));
857  mStencil[SixthDensePt< 0, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 0));
858  mStencil[SixthDensePt< 1, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 2, 0));
859  mStencil[SixthDensePt< 2, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 2, 0));
860  mStencil[SixthDensePt< 3, 2, 0>::idx] = mCache.getValue(ijk.offsetBy( 3, 2, 0));
861 
862  mStencil[SixthDensePt<-3, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-3, 1, 0));
863  mStencil[SixthDensePt<-2, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 1, 0));
864  mStencil[SixthDensePt<-1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 1, 0));
865  mStencil[SixthDensePt< 0, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
866  mStencil[SixthDensePt< 1, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 1, 0));
867  mStencil[SixthDensePt< 2, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 1, 0));
868  mStencil[SixthDensePt< 3, 1, 0>::idx] = mCache.getValue(ijk.offsetBy( 3, 1, 0));
869 
870  mStencil[SixthDensePt<-3, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-3, 0, 0));
871  mStencil[SixthDensePt<-2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 0));
872  mStencil[SixthDensePt<-1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
873  mStencil[SixthDensePt< 1, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
874  mStencil[SixthDensePt< 2, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 0));
875  mStencil[SixthDensePt< 3, 0, 0>::idx] = mCache.getValue(ijk.offsetBy( 3, 0, 0));
876 
877  mStencil[SixthDensePt<-3,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-3,-1, 0));
878  mStencil[SixthDensePt<-2,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-2,-1, 0));
879  mStencil[SixthDensePt<-1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy(-1,-1, 0));
880  mStencil[SixthDensePt< 0,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 0));
881  mStencil[SixthDensePt< 1,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 1,-1, 0));
882  mStencil[SixthDensePt< 2,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 2,-1, 0));
883  mStencil[SixthDensePt< 3,-1, 0>::idx] = mCache.getValue(ijk.offsetBy( 3,-1, 0));
884 
885  mStencil[SixthDensePt<-3,-2, 0>::idx] = mCache.getValue(ijk.offsetBy(-3,-2, 0));
886  mStencil[SixthDensePt<-2,-2, 0>::idx] = mCache.getValue(ijk.offsetBy(-2,-2, 0));
887  mStencil[SixthDensePt<-1,-2, 0>::idx] = mCache.getValue(ijk.offsetBy(-1,-2, 0));
888  mStencil[SixthDensePt< 0,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 0));
889  mStencil[SixthDensePt< 1,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 1,-2, 0));
890  mStencil[SixthDensePt< 2,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 2,-2, 0));
891  mStencil[SixthDensePt< 3,-2, 0>::idx] = mCache.getValue(ijk.offsetBy( 3,-2, 0));
892 
893  mStencil[SixthDensePt<-3,-3, 0>::idx] = mCache.getValue(ijk.offsetBy(-3,-3, 0));
894  mStencil[SixthDensePt<-2,-3, 0>::idx] = mCache.getValue(ijk.offsetBy(-2,-3, 0));
895  mStencil[SixthDensePt<-1,-3, 0>::idx] = mCache.getValue(ijk.offsetBy(-1,-3, 0));
896  mStencil[SixthDensePt< 0,-3, 0>::idx] = mCache.getValue(ijk.offsetBy( 0,-3, 0));
897  mStencil[SixthDensePt< 1,-3, 0>::idx] = mCache.getValue(ijk.offsetBy( 1,-3, 0));
898  mStencil[SixthDensePt< 2,-3, 0>::idx] = mCache.getValue(ijk.offsetBy( 2,-3, 0));
899  mStencil[SixthDensePt< 3,-3, 0>::idx] = mCache.getValue(ijk.offsetBy( 3,-3, 0));
900 
901  mStencil[SixthDensePt<-3, 0, 3>::idx] = mCache.getValue(ijk.offsetBy(-3, 0, 3));
902  mStencil[SixthDensePt<-2, 0, 3>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 3));
903  mStencil[SixthDensePt<-1, 0, 3>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 3));
904  mStencil[SixthDensePt< 0, 0, 3>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 3));
905  mStencil[SixthDensePt< 1, 0, 3>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 3));
906  mStencil[SixthDensePt< 2, 0, 3>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 3));
907  mStencil[SixthDensePt< 3, 0, 3>::idx] = mCache.getValue(ijk.offsetBy( 3, 0, 3));
908 
909  mStencil[SixthDensePt<-3, 0, 2>::idx] = mCache.getValue(ijk.offsetBy(-3, 0, 2));
910  mStencil[SixthDensePt<-2, 0, 2>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 2));
911  mStencil[SixthDensePt<-1, 0, 2>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 2));
912  mStencil[SixthDensePt< 0, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 2));
913  mStencil[SixthDensePt< 1, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 2));
914  mStencil[SixthDensePt< 2, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 2));
915  mStencil[SixthDensePt< 3, 0, 2>::idx] = mCache.getValue(ijk.offsetBy( 3, 0, 2));
916 
917  mStencil[SixthDensePt<-3, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-3, 0, 1));
918  mStencil[SixthDensePt<-2, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-2, 0, 1));
919  mStencil[SixthDensePt<-1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0, 1));
920  mStencil[SixthDensePt< 0, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
921  mStencil[SixthDensePt< 1, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0, 1));
922  mStencil[SixthDensePt< 2, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 2, 0, 1));
923  mStencil[SixthDensePt< 3, 0, 1>::idx] = mCache.getValue(ijk.offsetBy( 3, 0, 1));
924 
925  mStencil[SixthDensePt<-3, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-3, 0,-1));
926  mStencil[SixthDensePt<-2, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-2, 0,-1));
927  mStencil[SixthDensePt<-1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy(-1, 0,-1));
928  mStencil[SixthDensePt< 0, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 0,-1));
929  mStencil[SixthDensePt< 1, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 1, 0,-1));
930  mStencil[SixthDensePt< 2, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 2, 0,-1));
931  mStencil[SixthDensePt< 3, 0,-1>::idx] = mCache.getValue(ijk.offsetBy( 3, 0,-1));
932 
933  mStencil[SixthDensePt<-3, 0,-2>::idx] = mCache.getValue(ijk.offsetBy(-3, 0,-2));
934  mStencil[SixthDensePt<-2, 0,-2>::idx] = mCache.getValue(ijk.offsetBy(-2, 0,-2));
935  mStencil[SixthDensePt<-1, 0,-2>::idx] = mCache.getValue(ijk.offsetBy(-1, 0,-2));
936  mStencil[SixthDensePt< 0, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 0,-2));
937  mStencil[SixthDensePt< 1, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 1, 0,-2));
938  mStencil[SixthDensePt< 2, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 2, 0,-2));
939  mStencil[SixthDensePt< 3, 0,-2>::idx] = mCache.getValue(ijk.offsetBy( 3, 0,-2));
940 
941  mStencil[SixthDensePt<-3, 0,-3>::idx] = mCache.getValue(ijk.offsetBy(-3, 0,-3));
942  mStencil[SixthDensePt<-2, 0,-3>::idx] = mCache.getValue(ijk.offsetBy(-2, 0,-3));
943  mStencil[SixthDensePt<-1, 0,-3>::idx] = mCache.getValue(ijk.offsetBy(-1, 0,-3));
944  mStencil[SixthDensePt< 0, 0,-3>::idx] = mCache.getValue(ijk.offsetBy( 0, 0,-3));
945  mStencil[SixthDensePt< 1, 0,-3>::idx] = mCache.getValue(ijk.offsetBy( 1, 0,-3));
946  mStencil[SixthDensePt< 2, 0,-3>::idx] = mCache.getValue(ijk.offsetBy( 2, 0,-3));
947  mStencil[SixthDensePt< 3, 0,-3>::idx] = mCache.getValue(ijk.offsetBy( 3, 0,-3));
948 
949  mStencil[SixthDensePt< 0,-3, 3>::idx] = mCache.getValue(ijk.offsetBy( 0,-3, 3));
950  mStencil[SixthDensePt< 0,-2, 3>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 3));
951  mStencil[SixthDensePt< 0,-1, 3>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 3));
952  mStencil[SixthDensePt< 0, 1, 3>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 3));
953  mStencil[SixthDensePt< 0, 2, 3>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 3));
954  mStencil[SixthDensePt< 0, 3, 3>::idx] = mCache.getValue(ijk.offsetBy( 0, 3, 3));
955 
956  mStencil[SixthDensePt< 0,-3, 2>::idx] = mCache.getValue(ijk.offsetBy( 0,-3, 2));
957  mStencil[SixthDensePt< 0,-2, 2>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 2));
958  mStencil[SixthDensePt< 0,-1, 2>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 2));
959  mStencil[SixthDensePt< 0, 1, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 2));
960  mStencil[SixthDensePt< 0, 2, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 2));
961  mStencil[SixthDensePt< 0, 3, 2>::idx] = mCache.getValue(ijk.offsetBy( 0, 3, 2));
962 
963  mStencil[SixthDensePt< 0,-3, 1>::idx] = mCache.getValue(ijk.offsetBy( 0,-3, 1));
964  mStencil[SixthDensePt< 0,-2, 1>::idx] = mCache.getValue(ijk.offsetBy( 0,-2, 1));
965  mStencil[SixthDensePt< 0,-1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0,-1, 1));
966  mStencil[SixthDensePt< 0, 1, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1, 1));
967  mStencil[SixthDensePt< 0, 2, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 2, 1));
968  mStencil[SixthDensePt< 0, 3, 1>::idx] = mCache.getValue(ijk.offsetBy( 0, 3, 1));
969 
970  mStencil[SixthDensePt< 0,-3,-1>::idx] = mCache.getValue(ijk.offsetBy( 0,-3,-1));
971  mStencil[SixthDensePt< 0,-2,-1>::idx] = mCache.getValue(ijk.offsetBy( 0,-2,-1));
972  mStencil[SixthDensePt< 0,-1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0,-1,-1));
973  mStencil[SixthDensePt< 0, 1,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 1,-1));
974  mStencil[SixthDensePt< 0, 2,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 2,-1));
975  mStencil[SixthDensePt< 0, 3,-1>::idx] = mCache.getValue(ijk.offsetBy( 0, 3,-1));
976 
977  mStencil[SixthDensePt< 0,-3,-2>::idx] = mCache.getValue(ijk.offsetBy( 0,-3,-2));
978  mStencil[SixthDensePt< 0,-2,-2>::idx] = mCache.getValue(ijk.offsetBy( 0,-2,-2));
979  mStencil[SixthDensePt< 0,-1,-2>::idx] = mCache.getValue(ijk.offsetBy( 0,-1,-2));
980  mStencil[SixthDensePt< 0, 1,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 1,-2));
981  mStencil[SixthDensePt< 0, 2,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 2,-2));
982  mStencil[SixthDensePt< 0, 3,-2>::idx] = mCache.getValue(ijk.offsetBy( 0, 3,-2));
983 
984  mStencil[SixthDensePt< 0,-3,-3>::idx] = mCache.getValue(ijk.offsetBy( 0,-3,-3));
985  mStencil[SixthDensePt< 0,-2,-3>::idx] = mCache.getValue(ijk.offsetBy( 0,-2,-3));
986  mStencil[SixthDensePt< 0,-1,-3>::idx] = mCache.getValue(ijk.offsetBy( 0,-1,-3));
987  mStencil[SixthDensePt< 0, 1,-3>::idx] = mCache.getValue(ijk.offsetBy( 0, 1,-3));
988  mStencil[SixthDensePt< 0, 2,-3>::idx] = mCache.getValue(ijk.offsetBy( 0, 2,-3));
989  mStencil[SixthDensePt< 0, 3,-3>::idx] = mCache.getValue(ijk.offsetBy( 0, 3,-3));
990  }
991 
992  template<typename, typename> friend class BaseStencil; // allow base class to call init()
993  using BaseType::mCache;
994  using BaseType::mStencil;
995 };
996 
997 
999 
1000 
1007 template<typename GridType>
1008 class GradStencil: public BaseStencil<GridType, GradStencil<GridType> >
1009 {
1010 public:
1013  typedef typename GridType::ValueType ValueType;
1015 
1016  static const int SIZE = 7;
1017 
1018  GradStencil(const GridType& grid):
1019  BaseType(grid, SIZE),
1020  mInv2Dx(ValueType(0.5 / grid.voxelSize()[0])),
1021  mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1022  {
1023  }
1024 
1025  GradStencil(const GridType& grid, Real dx):
1026  BaseType(grid, SIZE),
1027  mInv2Dx(ValueType(0.5 / dx)),
1028  mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1029  {
1030  }
1031 
1037  inline ValueType normSqGrad() const
1038  {
1039  return mInvDx2 * math::GudonovsNormSqrd(mStencil[0] > 0,
1040  mStencil[0] - mStencil[1],
1041  mStencil[2] - mStencil[0],
1042  mStencil[0] - mStencil[3],
1043  mStencil[4] - mStencil[0],
1044  mStencil[0] - mStencil[5],
1045  mStencil[6] - mStencil[0]);
1046  }
1047 
1053  inline Vec3Type gradient() const
1054  {
1055  return Vec3Type(mStencil[2] - mStencil[1],
1056  mStencil[4] - mStencil[3],
1057  mStencil[6] - mStencil[5])*mInv2Dx;
1058  }
1063  inline Vec3Type gradient(const Vec3Type& V) const
1064  {
1065  return Vec3Type(V[0]>0 ? mStencil[0] - mStencil[1] : mStencil[2] - mStencil[0],
1066  V[1]>0 ? mStencil[0] - mStencil[3] : mStencil[4] - mStencil[0],
1067  V[2]>0 ? mStencil[0] - mStencil[5] : mStencil[6] - mStencil[0])*2*mInv2Dx;
1068  }
1069 
1072  inline ValueType laplacian() const
1073  {
1074  return mInvDx2 * (mStencil[1] + mStencil[2] +
1075  mStencil[3] + mStencil[4] +
1076  mStencil[5] + mStencil[6] - 6*mStencil[0]);
1077  }
1078 
1081  inline bool zeroCrossing() const
1082  {
1083  const BufferType& v = mStencil;
1084  return (v[0]>0 ? (v[1]<0 || v[2]<0 || v[3]<0 || v[4]<0 || v[5]<0 || v[6]<0)
1085  : (v[1]>0 || v[2]>0 || v[3]>0 || v[4]>0 || v[5]>0 || v[6]>0));
1086  }
1087 
1095  inline Vec3Type cpt()
1096  {
1097  const Coord& ijk = BaseType::getCenterCoord();
1098  const ValueType d = ValueType(mStencil[0] * 0.5 * mInvDx2); // distance in voxels / (2dx^2)
1099  return Vec3Type(ijk[0] - d*(mStencil[2] - mStencil[1]),
1100  ijk[1] - d*(mStencil[4] - mStencil[3]),
1101  ijk[2] - d*(mStencil[6] - mStencil[5]));
1102  }
1103 
1104 private:
1105 
1106  inline void init(const Coord& ijk)
1107  {
1108  mStencil[1] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
1109  mStencil[2] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
1110 
1111  mStencil[3] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
1112  mStencil[4] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
1113 
1114  mStencil[5] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
1115  mStencil[6] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
1116  }
1117 
1118  template<typename, typename> friend class BaseStencil; // allow base class to call init()
1119  using BaseType::mCache;
1120  using BaseType::mStencil;
1121  const ValueType mInv2Dx, mInvDx2;
1122 }; // class GradStencil
1123 
1124 
1126 
1127 
1133 template<typename GridType>
1134 class WenoStencil: public BaseStencil<GridType, WenoStencil<GridType> >
1135 {
1136 public:
1139  typedef typename GridType::ValueType ValueType;
1141 
1142  static const int SIZE = 19;
1143 
1144  WenoStencil(const GridType& grid):
1145  BaseType(grid, SIZE),
1146  mDx2(ValueType(math::Pow2(grid.voxelSize()[0]))),
1147  mInv2Dx(ValueType(0.5 / grid.voxelSize()[0])),
1148  mInvDx2(ValueType(1.0 / mDx2))
1149  {
1150  }
1151 
1152  WenoStencil(const GridType& grid, Real dx):
1153  BaseType(grid, SIZE),
1154  mDx2(ValueType(dx * dx)),
1155  mInv2Dx(ValueType(0.5 / dx)),
1156  mInvDx2(ValueType(1.0 / mDx2))
1157  {
1158  }
1159 
1165  inline ValueType normSqGrad() const
1166  {
1167  const BufferType& v = mStencil;
1168 #ifdef DWA_OPENVDB
1169  // SSE optimized
1170  const simd::Float4
1171  v1(v[2]-v[1], v[ 8]-v[ 7], v[14]-v[13], 0),
1172  v2(v[3]-v[2], v[ 9]-v[ 8], v[15]-v[14], 0),
1173  v3(v[0]-v[3], v[ 0]-v[ 9], v[ 0]-v[15], 0),
1174  v4(v[4]-v[0], v[10]-v[ 0], v[16]-v[ 0], 0),
1175  v5(v[5]-v[4], v[11]-v[10], v[17]-v[16], 0),
1176  v6(v[6]-v[5], v[12]-v[11], v[18]-v[17], 0),
1177  dP_m = math::WENO5(v1, v2, v3, v4, v5, mDx2),
1178  dP_p = math::WENO5(v6, v5, v4, v3, v2, mDx2);
1179 
1180  return mInvDx2 * math::GudonovsNormSqrd(mStencil[0] > 0, dP_m, dP_p);
1181 #else
1182  const Real
1183  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),
1184  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),
1185  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),
1186  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),
1187  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),
1188  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);
1189  return mInvDx2*math::GudonovsNormSqrd(v[0]>0,dP_xm,dP_xp,dP_ym,dP_yp,dP_zm,dP_zp);
1190 #endif
1191  }
1192 
1198  inline Vec3Type gradient(const Vec3Type& V) const
1199  {
1200  const BufferType& v = mStencil;
1201  return 2*mInv2Dx * Vec3Type(
1202  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)
1203  : math::WENO5(v[ 6]-v[ 5],v[ 5]-v[ 4],v[ 4]-v[ 0], v[ 0]-v[ 3],v[ 3]-v[ 2],mDx2),
1204  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)
1205  : math::WENO5(v[12]-v[11],v[11]-v[10],v[10]-v[ 0], v[ 0]-v[ 9],v[ 9]-v[ 8],mDx2),
1206  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)
1207  : math::WENO5(v[18]-v[17],v[17]-v[16],v[16]-v[ 0], v[ 0]-v[15],v[15]-v[14],mDx2));
1208  }
1214  inline Vec3Type gradient() const
1215  {
1216  return mInv2Dx * Vec3Type(
1217  mStencil[ 4] - mStencil[ 3],
1218  mStencil[10] - mStencil[ 9],
1219  mStencil[16] - mStencil[15]);
1220  }
1221 
1227  inline ValueType laplacian() const
1228  {
1229  return mInvDx2 * (
1230  mStencil[ 3] + mStencil[ 4] +
1231  mStencil[ 9] + mStencil[10] +
1232  mStencil[15] + mStencil[16] - 6*mStencil[0]);
1233  }
1234 
1237  inline bool zeroCrossing() const
1238  {
1239  const BufferType& v = mStencil;
1240  return (v[ 0]>0 ? (v[ 3]<0 || v[ 4]<0 || v[ 9]<0 || v[10]<0 || v[15]<0 || v[16]<0)
1241  : (v[ 3]>0 || v[ 4]>0 || v[ 9]>0 || v[10]>0 || v[15]>0 || v[16]>0));
1242  }
1243 
1244 private:
1245  inline void init(const Coord& ijk)
1246  {
1247  mStencil[ 1] = mCache.getValue(ijk.offsetBy(-3, 0, 0));
1248  mStencil[ 2] = mCache.getValue(ijk.offsetBy(-2, 0, 0));
1249  mStencil[ 3] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
1250  mStencil[ 4] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
1251  mStencil[ 5] = mCache.getValue(ijk.offsetBy( 2, 0, 0));
1252  mStencil[ 6] = mCache.getValue(ijk.offsetBy( 3, 0, 0));
1253 
1254  mStencil[ 7] = mCache.getValue(ijk.offsetBy( 0, -3, 0));
1255  mStencil[ 8] = mCache.getValue(ijk.offsetBy( 0, -2, 0));
1256  mStencil[ 9] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
1257  mStencil[10] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
1258  mStencil[11] = mCache.getValue(ijk.offsetBy( 0, 2, 0));
1259  mStencil[12] = mCache.getValue(ijk.offsetBy( 0, 3, 0));
1260 
1261  mStencil[13] = mCache.getValue(ijk.offsetBy( 0, 0, -3));
1262  mStencil[14] = mCache.getValue(ijk.offsetBy( 0, 0, -2));
1263  mStencil[15] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
1264  mStencil[16] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
1265  mStencil[17] = mCache.getValue(ijk.offsetBy( 0, 0, 2));
1266  mStencil[18] = mCache.getValue(ijk.offsetBy( 0, 0, 3));
1267  }
1268 
1269  template<typename, typename> friend class BaseStencil; // allow base class to call init()
1270  using BaseType::mCache;
1271  using BaseType::mStencil;
1272  const ValueType mDx2, mInv2Dx, mInvDx2;
1273 }; // class WenoStencil
1274 
1275 
1277 
1278 
1279 template<typename GridType>
1280 class CurvatureStencil: public BaseStencil<GridType, CurvatureStencil<GridType> >
1281 {
1282 public:
1284  typedef typename GridType::ValueType ValueType;
1285 
1286  static const int SIZE = 19;
1287 
1289  BaseType(grid, SIZE),
1290  mInv2Dx(ValueType(0.5 / grid.voxelSize()[0])),
1291  mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1292  {
1293  }
1294 
1295  CurvatureStencil(const GridType& grid, Real dx):
1296  BaseType(grid, SIZE),
1297  mInv2Dx(ValueType(0.5 / dx)),
1298  mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1299  {
1300  }
1301 
1307  {
1308  Real alpha, beta;
1309  this->meanCurvature(alpha, beta);
1310  return ValueType(alpha*mInv2Dx/math::Pow3(beta));
1311  }
1312 
1320  {
1321  Real alpha, beta;
1322  this->meanCurvature(alpha, beta);
1323  return ValueType(alpha*mInvDx2/(2*math::Pow2(beta)));
1324  }
1325 
1331  inline ValueType laplacian() const
1332  {
1333  return mInvDx2 * (
1334  mStencil[1] + mStencil[2] +
1335  mStencil[3] + mStencil[4] +
1336  mStencil[5] + mStencil[6] - 6*mStencil[0]);
1337  }
1338 
1345  {
1346  return math::Vec3<ValueType>(
1347  mStencil[2] - mStencil[1],
1348  mStencil[4] - mStencil[3],
1349  mStencil[6] - mStencil[5])*mInv2Dx;
1350  }
1351 
1352 private:
1353  inline void init(const Coord &ijk)
1354  {
1355  mStencil[ 1] = mCache.getValue(ijk.offsetBy(-1, 0, 0));
1356  mStencil[ 2] = mCache.getValue(ijk.offsetBy( 1, 0, 0));
1357 
1358  mStencil[ 3] = mCache.getValue(ijk.offsetBy( 0, -1, 0));
1359  mStencil[ 4] = mCache.getValue(ijk.offsetBy( 0, 1, 0));
1360 
1361  mStencil[ 5] = mCache.getValue(ijk.offsetBy( 0, 0, -1));
1362  mStencil[ 6] = mCache.getValue(ijk.offsetBy( 0, 0, 1));
1363 
1364  mStencil[ 7] = mCache.getValue(ijk.offsetBy(-1, -1, 0));
1365  mStencil[ 8] = mCache.getValue(ijk.offsetBy( 1, -1, 0));
1366  mStencil[ 9] = mCache.getValue(ijk.offsetBy(-1, 1, 0));
1367  mStencil[10] = mCache.getValue(ijk.offsetBy( 1, 1, 0));
1368 
1369  mStencil[11] = mCache.getValue(ijk.offsetBy(-1, 0, -1));
1370  mStencil[12] = mCache.getValue(ijk.offsetBy( 1, 0, -1));
1371  mStencil[13] = mCache.getValue(ijk.offsetBy(-1, 0, 1));
1372  mStencil[14] = mCache.getValue(ijk.offsetBy( 1, 0, 1));
1373 
1374  mStencil[15] = mCache.getValue(ijk.offsetBy( 0, -1, -1));
1375  mStencil[16] = mCache.getValue(ijk.offsetBy( 0, 1, -1));
1376  mStencil[17] = mCache.getValue(ijk.offsetBy( 0, -1, 1));
1377  mStencil[18] = mCache.getValue(ijk.offsetBy( 0, 1, 1));
1378  }
1379 
1380  inline void meanCurvature(Real& alpha, Real& beta) const
1381  {
1382  // For performance all finite differences are unscaled wrt dx
1383  const Real
1384  Half(0.5), Quarter(0.25),
1385  Dx = Half * (mStencil[2] - mStencil[1]), Dx2 = Dx * Dx, // * 1/dx
1386  Dy = Half * (mStencil[4] - mStencil[3]), Dy2 = Dy * Dy, // * 1/dx
1387  Dz = Half * (mStencil[6] - mStencil[5]), Dz2 = Dz * Dz, // * 1/dx
1388  Dxx = mStencil[2] - 2 * mStencil[0] + mStencil[1], // * 1/dx2
1389  Dyy = mStencil[4] - 2 * mStencil[0] + mStencil[3], // * 1/dx2
1390  Dzz = mStencil[6] - 2 * mStencil[0] + mStencil[5], // * 1/dx2
1391  Dxy = Quarter * (mStencil[10] - mStencil[ 8] + mStencil[7] - mStencil[ 9]), // * 1/dx2
1392  Dxz = Quarter * (mStencil[14] - mStencil[12] + mStencil[11] - mStencil[13]), // * 1/dx2
1393  Dyz = Quarter * (mStencil[18] - mStencil[16] + mStencil[15] - mStencil[17]); // * 1/dx2
1394  alpha = (Dx2*(Dyy+Dzz)+Dy2*(Dxx+Dzz)+Dz2*(Dxx+Dyy)-2*(Dx*(Dy*Dxy+Dz*Dxz)+Dy*Dz*Dyz));
1395  beta = std::sqrt(Dx2 + Dy2 + Dz2); // * 1/dx
1396  }
1397 
1398  template<typename, typename> friend class BaseStencil; // allow base class to call init()
1399  using BaseType::mCache;
1400  using BaseType::mStencil;
1401  const ValueType mInv2Dx, mInvDx2;
1402 }; // class CurvatureStencil
1403 
1404 
1406 
1407 
1409 template<typename GridType>
1410 class DenseStencil: public BaseStencil<GridType, DenseStencil<GridType> >
1411 {
1412 public:
1414  typedef typename GridType::ValueType ValueType;
1415 
1416  DenseStencil(const GridType& grid, int halfWidth) :
1417  BaseType(grid, /*size=*/math::Pow3(2 * halfWidth + 1)),
1418  mHalfWidth(halfWidth)
1419  {
1420  assert(halfWidth>0);
1421  }
1422 
1423  inline const ValueType& getCenterValue() const { return mStencil[(mStencil.size()-1)>>1]; }
1424 
1427  inline void moveTo(const Coord& ijk)
1428  {
1429  BaseType::mCenter = ijk;
1430  this->init(ijk);
1431  }
1434  template<typename IterType>
1435  inline void moveTo(const IterType& iter)
1436  {
1437  BaseType::mCenter = iter.getCoord();
1438  this->init(BaseType::mCenter);
1439  }
1440 
1441 private:
1444  inline void init(const Coord& ijk)
1445  {
1446  int n = 0;
1447  for (Coord p=ijk.offsetBy(-mHalfWidth), q=ijk.offsetBy(mHalfWidth); p[0] <= q[0]; ++p[0]) {
1448  for (p[1] = ijk[1]-mHalfWidth; p[1] <= q[1]; ++p[1]) {
1449  for (p[2] = ijk[2]-mHalfWidth; p[2] <= q[2]; ++p[2]) {
1450  mStencil[n++] = mCache.getValue(p);
1451  }
1452  }
1453  }
1454  }
1455 
1456  template<typename, typename> friend class BaseStencil; // allow base class to call init()
1457  using BaseType::mCache;
1458  using BaseType::mStencil;
1459  const int mHalfWidth;
1460 };
1461 
1462 
1463 } // end math namespace
1464 } // namespace OPENVDB_VERSION_NAME
1465 } // end openvdb namespace
1466 
1467 #endif // OPENVDB_MATH_STENCILS_HAS_BEEN_INCLUDED
1468 
1469 // Copyright (c) 2012-2013 DreamWorks Animation LLC
1470 // All rights reserved. This software is distributed under the
1471 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )