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