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