OpenVDB 10.0.1
Loading...
Searching...
No Matches
Stencils.h
Go to the documentation of this file.
1// Copyright Contributors to the OpenVDB Project
2// SPDX-License-Identifier: MPL-2.0
3//
4/// @author Ken Museth
5///
6/// @file Stencils.h
7///
8/// @brief Defines various finite difference stencils by means of the
9/// "curiously recurring template pattern" on a BaseStencil
10/// that caches stencil values and stores a ValueAccessor for
11/// fast lookup.
12
13#ifndef OPENVDB_MATH_STENCILS_HAS_BEEN_INCLUDED
14#define OPENVDB_MATH_STENCILS_HAS_BEEN_INCLUDED
15
16#include <algorithm>
17#include <vector> // for std::vector
18#include <bitset> // for std::bitset
19#include <openvdb/Types.h> // for Real
21
22#include "Math.h" // for Pow2, needed by WENO and Godunov
23#include "Coord.h" // for Coord
24#include "FiniteDifference.h" // for WENO5 and GodunovsNormSqrd
25
26namespace openvdb {
28namespace OPENVDB_VERSION_NAME {
29namespace math {
30
31
32////////////////////////////////////////
33
34template<typename DerivedType, typename GridT, bool IsSafe>
36{
37public:
38 typedef GridT GridType;
39 typedef typename GridT::TreeType TreeType;
40 typedef typename GridT::ValueType ValueType;
42 typedef std::vector<ValueType> BufferType;
43 typedef typename BufferType::iterator IterType;
44
45 /// @brief Initialize the stencil buffer with the values of voxel (i, j, k)
46 /// and its neighbors.
47 /// @param ijk Index coordinates of stencil center
48 inline void moveTo(const Coord& ijk)
49 {
50 mCenter = ijk;
51 mValues[0] = mAcc.getValue(ijk);
52 static_cast<DerivedType&>(*this).init(mCenter);
53 }
54
55 /// @brief Initialize the stencil buffer with the values of voxel (i, j, k)
56 /// and its neighbors. The method also takes a value of the center
57 /// element of the stencil, assuming it is already known.
58 /// @param ijk Index coordinates of stnecil center
59 /// @param centerValue Value of the center element of the stencil
60 inline void moveTo(const Coord& ijk, const ValueType& centerValue)
61 {
62 mCenter = ijk;
63 mValues[0] = centerValue;
64 static_cast<DerivedType&>(*this).init(mCenter);
65 }
66
67 /// @brief Initialize the stencil buffer with the values of voxel
68 /// (x, y, z) and its neighbors.
69 ///
70 /// @note This version is slightly faster than the one above, since
71 /// the center voxel's value is read directly from the iterator.
72 template<typename IterType>
73 inline void moveTo(const IterType& iter)
74 {
75 mCenter = iter.getCoord();
76 mValues[0] = *iter;
77 static_cast<DerivedType&>(*this).init(mCenter);
78 }
79
80 /// @brief Initialize the stencil buffer with the values of voxel (x, y, z)
81 /// and its neighbors.
82 /// @param xyz Floating point voxel coordinates of stencil center
83 /// @details This method will check to see if it is necessary to
84 /// update the stencil based on the cached index coordinates of
85 /// the center point.
86 template<typename RealType>
87 inline void moveTo(const Vec3<RealType>& xyz)
88 {
89 Coord ijk = Coord::floor(xyz);
90 if (ijk != mCenter) this->moveTo(ijk);
91 }
92
93 /// @brief Return the value from the stencil buffer with linear
94 /// offset pos.
95 ///
96 /// @note The default (@a pos = 0) corresponds to the first element
97 /// which is typically the center point of the stencil.
98 inline const ValueType& getValue(unsigned int pos = 0) const
99 {
100 assert(pos < mValues.size());
101 return mValues[pos];
102 }
103
104 /// @brief Return the value at the specified location relative to the center of the stencil
105 template<int i, int j, int k>
106 inline const ValueType& getValue() const
107 {
108 return mValues[static_cast<const DerivedType&>(*this).template pos<i,j,k>()];
109 }
110
111 /// @brief Set the value at the specified location relative to the center of the stencil
112 template<int i, int j, int k>
113 inline void setValue(const ValueType& value)
114 {
115 mValues[static_cast<const DerivedType&>(*this).template pos<i,j,k>()] = value;
116 }
117
118 /// @brief Return the size of the stencil buffer.
119 inline int size() { return mValues.size(); }
120
121 /// @brief Return the median value of the current stencil.
122 inline ValueType median() const
123 {
124 BufferType tmp(mValues);//local copy
125 assert(!tmp.empty());
126 size_t midpoint = (tmp.size() - 1) >> 1;
127 // Partially sort the vector until the median value is at the midpoint.
128#if !defined(_MSC_VER) || _MSC_VER < 1924
129 std::nth_element(tmp.begin(), tmp.begin() + midpoint, tmp.end());
130#else
131 // Workaround MSVC bool warning C4804 unsafe use of type 'bool'
132 std::nth_element(tmp.begin(), tmp.begin() + midpoint, tmp.end(),
133 std::less<ValueType>());
134#endif
135 return tmp[midpoint];
136 }
137
138 /// @brief Return the mean value of the current stencil.
139 inline ValueType mean() const
140 {
141 ValueType sum = 0.0;
142 for (int n = 0, s = int(mValues.size()); n < s; ++n) sum += mValues[n];
143 return sum / ValueType(mValues.size());
144 }
145
146 /// @brief Return the smallest value in the stencil buffer.
147 inline ValueType min() const
148 {
149 IterType iter = std::min_element(mValues.begin(), mValues.end());
150 return *iter;
151 }
152
153 /// @brief Return the largest value in the stencil buffer.
154 inline ValueType max() const
155 {
156 IterType iter = std::max_element(mValues.begin(), mValues.end());
157 return *iter;
158 }
159
160 /// @brief Return the coordinates of the center point of the stencil.
161 inline const Coord& getCenterCoord() const { return mCenter; }
162
163 /// @brief Return the value at the center of the stencil
164 inline const ValueType& getCenterValue() const { return mValues[0]; }
165
166 /// @brief Return true if the center of the stencil intersects the
167 /// iso-contour specified by the isoValue
168 inline bool intersects(const ValueType &isoValue = zeroVal<ValueType>()) const
169 {
170 const bool less = this->getValue< 0, 0, 0>() < isoValue;
171 return (less ^ (this->getValue<-1, 0, 0>() < isoValue)) ||
172 (less ^ (this->getValue< 1, 0, 0>() < isoValue)) ||
173 (less ^ (this->getValue< 0,-1, 0>() < isoValue)) ||
174 (less ^ (this->getValue< 0, 1, 0>() < isoValue)) ||
175 (less ^ (this->getValue< 0, 0,-1>() < isoValue)) ||
176 (less ^ (this->getValue< 0, 0, 1>() < isoValue)) ;
177 }
178
179 /// @brief Return true a bit-mask where the 6 bits indicates if the
180 /// center of the stencil intersects the iso-contour specified by the isoValue.
181 ///
182 /// @note There are 2^6 = 64 different possible cases, including no intersections!
183 ///
184 /// @details The ordering of bit mask is ( -x, +x, -y, +y, -z, +z ), so to
185 /// check if there is an intersection in -y use mask.test(2) where mask is
186 /// ther return value from this function. To check if there are any
187 /// intersections use mask.any(), and for no intersections use mask.none().
188 /// To count the number of intersections use mask.count().
189 inline std::bitset<6> intersectionMask(const ValueType &isoValue = zeroVal<ValueType>()) const
190 {
191 std::bitset<6> mask;
192 const bool less = this->getValue< 0, 0, 0>() < isoValue;
193 mask[0] = less ^ (this->getValue<-1, 0, 0>() < isoValue);
194 mask[1] = less ^ (this->getValue< 1, 0, 0>() < isoValue);
195 mask[2] = less ^ (this->getValue< 0,-1, 0>() < isoValue);
196 mask[3] = less ^ (this->getValue< 0, 1, 0>() < isoValue);
197 mask[4] = less ^ (this->getValue< 0, 0,-1>() < isoValue);
198 mask[5] = less ^ (this->getValue< 0, 0, 1>() < isoValue);
199 return mask;
200 }
201
202 /// @brief Return a const reference to the grid from which this
203 /// stencil was constructed.
204 inline const GridType& grid() const { return *mGrid; }
205
206 /// @brief Return a const reference to the ValueAccessor
207 /// associated with this Stencil.
208 inline const AccessorType& accessor() const { return mAcc; }
209
210protected:
211 // Constructor is protected to prevent direct instantiation.
212 BaseStencil(const GridType& grid, int size)
213 : mGrid(&grid)
214 , mAcc(grid.tree())
215 , mValues(size)
216 , mCenter(Coord::max())
217 {
218 }
219
224
225}; // BaseStencil class
226
227
228////////////////////////////////////////
229
230
231namespace { // anonymous namespace for stencil-layout map
232
233 // the seven point stencil
234 template<int i, int j, int k> struct SevenPt {};
235 template<> struct SevenPt< 0, 0, 0> { enum { idx = 0 }; };
236 template<> struct SevenPt< 1, 0, 0> { enum { idx = 1 }; };
237 template<> struct SevenPt< 0, 1, 0> { enum { idx = 2 }; };
238 template<> struct SevenPt< 0, 0, 1> { enum { idx = 3 }; };
239 template<> struct SevenPt<-1, 0, 0> { enum { idx = 4 }; };
240 template<> struct SevenPt< 0,-1, 0> { enum { idx = 5 }; };
241 template<> struct SevenPt< 0, 0,-1> { enum { idx = 6 }; };
242
243}
244
245
246template<typename GridT, bool IsSafe = true>
247class SevenPointStencil: public BaseStencil<SevenPointStencil<GridT, IsSafe>, GridT, IsSafe>
248{
251public:
252 typedef GridT GridType;
253 typedef typename GridT::TreeType TreeType;
254 typedef typename GridT::ValueType ValueType;
255
256 static const int SIZE = 7;
257
258 SevenPointStencil(const GridT& grid): BaseType(grid, SIZE) {}
259
260 /// Return linear offset for the specified stencil point relative to its center
261 template<int i, int j, int k>
262 unsigned int pos() const { return SevenPt<i,j,k>::idx; }
263
264private:
265 inline void init(const Coord& ijk)
266 {
267 BaseType::template setValue<-1, 0, 0>(mAcc.getValue(ijk.offsetBy(-1, 0, 0)));
268 BaseType::template setValue< 1, 0, 0>(mAcc.getValue(ijk.offsetBy( 1, 0, 0)));
269
270 BaseType::template setValue< 0,-1, 0>(mAcc.getValue(ijk.offsetBy( 0,-1, 0)));
271 BaseType::template setValue< 0, 1, 0>(mAcc.getValue(ijk.offsetBy( 0, 1, 0)));
272
273 BaseType::template setValue< 0, 0,-1>(mAcc.getValue(ijk.offsetBy( 0, 0,-1)));
274 BaseType::template setValue< 0, 0, 1>(mAcc.getValue(ijk.offsetBy( 0, 0, 1)));
275 }
276
277 template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
278 using BaseType::mAcc;
279 using BaseType::mValues;
280};// SevenPointStencil class
281
282
283////////////////////////////////////////
284
285
286namespace { // anonymous namespace for stencil-layout map
287
288 // the eight point box stencil
289 template<int i, int j, int k> struct BoxPt {};
290 template<> struct BoxPt< 0, 0, 0> { enum { idx = 0 }; };
291 template<> struct BoxPt< 0, 0, 1> { enum { idx = 1 }; };
292 template<> struct BoxPt< 0, 1, 1> { enum { idx = 2 }; };
293 template<> struct BoxPt< 0, 1, 0> { enum { idx = 3 }; };
294 template<> struct BoxPt< 1, 0, 0> { enum { idx = 4 }; };
295 template<> struct BoxPt< 1, 0, 1> { enum { idx = 5 }; };
296 template<> struct BoxPt< 1, 1, 1> { enum { idx = 6 }; };
297 template<> struct BoxPt< 1, 1, 0> { enum { idx = 7 }; };
298}
299
300template<typename GridT, bool IsSafe = true>
301class BoxStencil: public BaseStencil<BoxStencil<GridT, IsSafe>, GridT, IsSafe>
302{
305public:
306 typedef GridT GridType;
307 typedef typename GridT::TreeType TreeType;
308 typedef typename GridT::ValueType ValueType;
309
310 static const int SIZE = 8;
311
312 BoxStencil(const GridType& grid): BaseType(grid, SIZE) {}
313
314 /// Return linear offset for the specified stencil point relative to its center
315 template<int i, int j, int k>
316 unsigned int pos() const { return BoxPt<i,j,k>::idx; }
317
318 /// @brief Return true if the center of the stencil intersects the
319 /// iso-contour specified by the isoValue
320 inline bool intersects(const ValueType &isoValue = zeroVal<ValueType>()) const
321 {
322 const bool less = mValues[0] < isoValue;
323 return (less ^ (mValues[1] < isoValue)) ||
324 (less ^ (mValues[2] < isoValue)) ||
325 (less ^ (mValues[3] < isoValue)) ||
326 (less ^ (mValues[4] < isoValue)) ||
327 (less ^ (mValues[5] < isoValue)) ||
328 (less ^ (mValues[6] < isoValue)) ||
329 (less ^ (mValues[7] < isoValue)) ;
330 }
331
332 /// @brief Return the trilinear interpolation at the normalized position.
333 /// @param xyz Floating point coordinate position.
334 /// @warning It is assumed that the stencil has already been moved
335 /// to the relevant voxel position, e.g. using moveTo(xyz).
336 /// @note Trilinear interpolation kernal reads as:
337 /// v000 (1-u)(1-v)(1-w) + v001 (1-u)(1-v)w + v010 (1-u)v(1-w) + v011 (1-u)vw
338 /// + v100 u(1-v)(1-w) + v101 u(1-v)w + v110 uv(1-w) + v111 uvw
340 {
342 const ValueType u = xyz[0] - BaseType::mCenter[0];
343 const ValueType v = xyz[1] - BaseType::mCenter[1];
344 const ValueType w = xyz[2] - BaseType::mCenter[2];
346
347 assert(u>=0 && u<=1);
348 assert(v>=0 && v<=1);
349 assert(w>=0 && w<=1);
350
351 ValueType V = BaseType::template getValue<0,0,0>();
352 ValueType A = static_cast<ValueType>(V + (BaseType::template getValue<0,0,1>() - V) * w);
353 V = BaseType::template getValue< 0, 1, 0>();
354 ValueType B = static_cast<ValueType>(V + (BaseType::template getValue<0,1,1>() - V) * w);
355 ValueType C = static_cast<ValueType>(A + (B - A) * v);
356
357 V = BaseType::template getValue<1,0,0>();
358 A = static_cast<ValueType>(V + (BaseType::template getValue<1,0,1>() - V) * w);
359 V = BaseType::template getValue<1,1,0>();
360 B = static_cast<ValueType>(V + (BaseType::template getValue<1,1,1>() - V) * w);
361 ValueType D = static_cast<ValueType>(A + (B - A) * v);
362
363 return static_cast<ValueType>(C + (D - C) * u);
364 }
365
366 /// @brief Return the gradient in world space of the trilinear interpolation kernel.
367 /// @param xyz Floating point coordinate position.
368 /// @warning It is assumed that the stencil has already been moved
369 /// to the relevant voxel position, e.g. using moveTo(xyz).
370 /// @note Computed as partial derivatives of the trilinear interpolation kernel:
371 /// v000 (1-u)(1-v)(1-w) + v001 (1-u)(1-v)w + v010 (1-u)v(1-w) + v011 (1-u)vw
372 /// + v100 u(1-v)(1-w) + v101 u(1-v)w + v110 uv(1-w) + v111 uvw
374 {
376 const ValueType u = xyz[0] - BaseType::mCenter[0];
377 const ValueType v = xyz[1] - BaseType::mCenter[1];
378 const ValueType w = xyz[2] - BaseType::mCenter[2];
380
381 assert(u>=0 && u<=1);
382 assert(v>=0 && v<=1);
383 assert(w>=0 && w<=1);
384
385 ValueType D[4]={BaseType::template getValue<0,0,1>()-BaseType::template getValue<0,0,0>(),
386 BaseType::template getValue<0,1,1>()-BaseType::template getValue<0,1,0>(),
387 BaseType::template getValue<1,0,1>()-BaseType::template getValue<1,0,0>(),
388 BaseType::template getValue<1,1,1>()-BaseType::template getValue<1,1,0>()};
389
390 // Z component
391 ValueType A = static_cast<ValueType>(D[0] + (D[1]- D[0]) * v);
392 ValueType B = static_cast<ValueType>(D[2] + (D[3]- D[2]) * v);
393 math::Vec3<ValueType> grad(zeroVal<ValueType>(),
394 zeroVal<ValueType>(),
395 static_cast<ValueType>(A + (B - A) * u));
396
397 D[0] = static_cast<ValueType>(BaseType::template getValue<0,0,0>() + D[0] * w);
398 D[1] = static_cast<ValueType>(BaseType::template getValue<0,1,0>() + D[1] * w);
399 D[2] = static_cast<ValueType>(BaseType::template getValue<1,0,0>() + D[2] * w);
400 D[3] = static_cast<ValueType>(BaseType::template getValue<1,1,0>() + D[3] * w);
401
402 // X component
403 A = static_cast<ValueType>(D[0] + (D[1] - D[0]) * v);
404 B = static_cast<ValueType>(D[2] + (D[3] - D[2]) * v);
405
406 grad[0] = B - A;
407
408 // Y component
409 A = D[1] - D[0];
410 B = D[3] - D[2];
411
412 grad[1] = static_cast<ValueType>(A + (B - A) * u);
413
414 return BaseType::mGrid->transform().baseMap()->applyIJT(grad, xyz);
415 }
416
417private:
418 inline void init(const Coord& ijk)
419 {
420 BaseType::template setValue< 0, 0, 1>(mAcc.getValue(ijk.offsetBy( 0, 0, 1)));
421 BaseType::template setValue< 0, 1, 1>(mAcc.getValue(ijk.offsetBy( 0, 1, 1)));
422 BaseType::template setValue< 0, 1, 0>(mAcc.getValue(ijk.offsetBy( 0, 1, 0)));
423 BaseType::template setValue< 1, 0, 0>(mAcc.getValue(ijk.offsetBy( 1, 0, 0)));
424 BaseType::template setValue< 1, 0, 1>(mAcc.getValue(ijk.offsetBy( 1, 0, 1)));
425 BaseType::template setValue< 1, 1, 1>(mAcc.getValue(ijk.offsetBy( 1, 1, 1)));
426 BaseType::template setValue< 1, 1, 0>(mAcc.getValue(ijk.offsetBy( 1, 1, 0)));
427 }
428
429 template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
430 using BaseType::mAcc;
431 using BaseType::mValues;
432};// BoxStencil class
433
434
435////////////////////////////////////////
436
437
438namespace { // anonymous namespace for stencil-layout map
439
440 // the dense point stencil
441 template<int i, int j, int k> struct DensePt {};
442 template<> struct DensePt< 0, 0, 0> { enum { idx = 0 }; };
443
444 template<> struct DensePt< 1, 0, 0> { enum { idx = 1 }; };
445 template<> struct DensePt< 0, 1, 0> { enum { idx = 2 }; };
446 template<> struct DensePt< 0, 0, 1> { enum { idx = 3 }; };
447
448 template<> struct DensePt<-1, 0, 0> { enum { idx = 4 }; };
449 template<> struct DensePt< 0,-1, 0> { enum { idx = 5 }; };
450 template<> struct DensePt< 0, 0,-1> { enum { idx = 6 }; };
451
452 template<> struct DensePt<-1,-1, 0> { enum { idx = 7 }; };
453 template<> struct DensePt< 0,-1,-1> { enum { idx = 8 }; };
454 template<> struct DensePt<-1, 0,-1> { enum { idx = 9 }; };
455
456 template<> struct DensePt< 1,-1, 0> { enum { idx = 10 }; };
457 template<> struct DensePt< 0, 1,-1> { enum { idx = 11 }; };
458 template<> struct DensePt<-1, 0, 1> { enum { idx = 12 }; };
459
460 template<> struct DensePt<-1, 1, 0> { enum { idx = 13 }; };
461 template<> struct DensePt< 0,-1, 1> { enum { idx = 14 }; };
462 template<> struct DensePt< 1, 0,-1> { enum { idx = 15 }; };
463
464 template<> struct DensePt< 1, 1, 0> { enum { idx = 16 }; };
465 template<> struct DensePt< 0, 1, 1> { enum { idx = 17 }; };
466 template<> struct DensePt< 1, 0, 1> { enum { idx = 18 }; };
467}
468
469
470template<typename GridT, bool IsSafe = true>
472 : public BaseStencil<SecondOrderDenseStencil<GridT, IsSafe>, GridT, IsSafe >
473{
476public:
477 typedef GridT GridType;
478 typedef typename GridT::TreeType TreeType;
479 typedef typename GridType::ValueType ValueType;
480
481 static const int SIZE = 19;
482
483 SecondOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
484
485 /// Return linear offset for the specified stencil point relative to its center
486 template<int i, int j, int k>
487 unsigned int pos() const { return DensePt<i,j,k>::idx; }
488
489private:
490 inline void init(const Coord& ijk)
491 {
492 mValues[DensePt< 1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
493 mValues[DensePt< 0, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
494 mValues[DensePt< 0, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
495
496 mValues[DensePt<-1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
497 mValues[DensePt< 0,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
498 mValues[DensePt< 0, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
499
500 mValues[DensePt<-1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, -1, 0));
501 mValues[DensePt< 1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, -1, 0));
502 mValues[DensePt<-1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 1, 0));
503 mValues[DensePt< 1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 1, 0));
504
505 mValues[DensePt<-1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, -1));
506 mValues[DensePt< 1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, -1));
507 mValues[DensePt<-1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 1));
508 mValues[DensePt< 1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 1));
509
510 mValues[DensePt< 0,-1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, -1, -1));
511 mValues[DensePt< 0, 1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, -1));
512 mValues[DensePt< 0,-1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, -1, 1));
513 mValues[DensePt< 0, 1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 1));
514 }
515
516 template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
517 using BaseType::mAcc;
518 using BaseType::mValues;
519};// SecondOrderDenseStencil class
520
521
522////////////////////////////////////////
523
524
525namespace { // anonymous namespace for stencil-layout map
526
527 // the dense point stencil
528 template<int i, int j, int k> struct ThirteenPt {};
529 template<> struct ThirteenPt< 0, 0, 0> { enum { idx = 0 }; };
530
531 template<> struct ThirteenPt< 1, 0, 0> { enum { idx = 1 }; };
532 template<> struct ThirteenPt< 0, 1, 0> { enum { idx = 2 }; };
533 template<> struct ThirteenPt< 0, 0, 1> { enum { idx = 3 }; };
534
535 template<> struct ThirteenPt<-1, 0, 0> { enum { idx = 4 }; };
536 template<> struct ThirteenPt< 0,-1, 0> { enum { idx = 5 }; };
537 template<> struct ThirteenPt< 0, 0,-1> { enum { idx = 6 }; };
538
539 template<> struct ThirteenPt< 2, 0, 0> { enum { idx = 7 }; };
540 template<> struct ThirteenPt< 0, 2, 0> { enum { idx = 8 }; };
541 template<> struct ThirteenPt< 0, 0, 2> { enum { idx = 9 }; };
542
543 template<> struct ThirteenPt<-2, 0, 0> { enum { idx = 10 }; };
544 template<> struct ThirteenPt< 0,-2, 0> { enum { idx = 11 }; };
545 template<> struct ThirteenPt< 0, 0,-2> { enum { idx = 12 }; };
546
547}
548
549
550template<typename GridT, bool IsSafe = true>
552 : public BaseStencil<ThirteenPointStencil<GridT, IsSafe>, GridT, IsSafe>
553{
556public:
557 typedef GridT GridType;
558 typedef typename GridT::TreeType TreeType;
559 typedef typename GridType::ValueType ValueType;
560
561 static const int SIZE = 13;
562
563 ThirteenPointStencil(const GridType& grid): BaseType(grid, SIZE) {}
564
565 /// Return linear offset for the specified stencil point relative to its center
566 template<int i, int j, int k>
567 unsigned int pos() const { return ThirteenPt<i,j,k>::idx; }
568
569private:
570 inline void init(const Coord& ijk)
571 {
572 mValues[ThirteenPt< 2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 0));
573 mValues[ThirteenPt< 1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
574 mValues[ThirteenPt<-1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
575 mValues[ThirteenPt<-2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 0));
576
577 mValues[ThirteenPt< 0, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 0));
578 mValues[ThirteenPt< 0, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
579 mValues[ThirteenPt< 0,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
580 mValues[ThirteenPt< 0,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -2, 0));
581
582 mValues[ThirteenPt< 0, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 2));
583 mValues[ThirteenPt< 0, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
584 mValues[ThirteenPt< 0, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
585 mValues[ThirteenPt< 0, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -2));
586 }
587
588 template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
589 using BaseType::mAcc;
590 using BaseType::mValues;
591};// ThirteenPointStencil class
592
593
594////////////////////////////////////////
595
596
597namespace { // anonymous namespace for stencil-layout map
598
599 // the 4th-order dense point stencil
600 template<int i, int j, int k> struct FourthDensePt {};
601 template<> struct FourthDensePt< 0, 0, 0> { enum { idx = 0 }; };
602
603 template<> struct FourthDensePt<-2, 2, 0> { enum { idx = 1 }; };
604 template<> struct FourthDensePt<-1, 2, 0> { enum { idx = 2 }; };
605 template<> struct FourthDensePt< 0, 2, 0> { enum { idx = 3 }; };
606 template<> struct FourthDensePt< 1, 2, 0> { enum { idx = 4 }; };
607 template<> struct FourthDensePt< 2, 2, 0> { enum { idx = 5 }; };
608
609 template<> struct FourthDensePt<-2, 1, 0> { enum { idx = 6 }; };
610 template<> struct FourthDensePt<-1, 1, 0> { enum { idx = 7 }; };
611 template<> struct FourthDensePt< 0, 1, 0> { enum { idx = 8 }; };
612 template<> struct FourthDensePt< 1, 1, 0> { enum { idx = 9 }; };
613 template<> struct FourthDensePt< 2, 1, 0> { enum { idx = 10 }; };
614
615 template<> struct FourthDensePt<-2, 0, 0> { enum { idx = 11 }; };
616 template<> struct FourthDensePt<-1, 0, 0> { enum { idx = 12 }; };
617 template<> struct FourthDensePt< 1, 0, 0> { enum { idx = 13 }; };
618 template<> struct FourthDensePt< 2, 0, 0> { enum { idx = 14 }; };
619
620 template<> struct FourthDensePt<-2,-1, 0> { enum { idx = 15 }; };
621 template<> struct FourthDensePt<-1,-1, 0> { enum { idx = 16 }; };
622 template<> struct FourthDensePt< 0,-1, 0> { enum { idx = 17 }; };
623 template<> struct FourthDensePt< 1,-1, 0> { enum { idx = 18 }; };
624 template<> struct FourthDensePt< 2,-1, 0> { enum { idx = 19 }; };
625
626 template<> struct FourthDensePt<-2,-2, 0> { enum { idx = 20 }; };
627 template<> struct FourthDensePt<-1,-2, 0> { enum { idx = 21 }; };
628 template<> struct FourthDensePt< 0,-2, 0> { enum { idx = 22 }; };
629 template<> struct FourthDensePt< 1,-2, 0> { enum { idx = 23 }; };
630 template<> struct FourthDensePt< 2,-2, 0> { enum { idx = 24 }; };
631
632
633 template<> struct FourthDensePt<-2, 0, 2> { enum { idx = 25 }; };
634 template<> struct FourthDensePt<-1, 0, 2> { enum { idx = 26 }; };
635 template<> struct FourthDensePt< 0, 0, 2> { enum { idx = 27 }; };
636 template<> struct FourthDensePt< 1, 0, 2> { enum { idx = 28 }; };
637 template<> struct FourthDensePt< 2, 0, 2> { enum { idx = 29 }; };
638
639 template<> struct FourthDensePt<-2, 0, 1> { enum { idx = 30 }; };
640 template<> struct FourthDensePt<-1, 0, 1> { enum { idx = 31 }; };
641 template<> struct FourthDensePt< 0, 0, 1> { enum { idx = 32 }; };
642 template<> struct FourthDensePt< 1, 0, 1> { enum { idx = 33 }; };
643 template<> struct FourthDensePt< 2, 0, 1> { enum { idx = 34 }; };
644
645 template<> struct FourthDensePt<-2, 0,-1> { enum { idx = 35 }; };
646 template<> struct FourthDensePt<-1, 0,-1> { enum { idx = 36 }; };
647 template<> struct FourthDensePt< 0, 0,-1> { enum { idx = 37 }; };
648 template<> struct FourthDensePt< 1, 0,-1> { enum { idx = 38 }; };
649 template<> struct FourthDensePt< 2, 0,-1> { enum { idx = 39 }; };
650
651 template<> struct FourthDensePt<-2, 0,-2> { enum { idx = 40 }; };
652 template<> struct FourthDensePt<-1, 0,-2> { enum { idx = 41 }; };
653 template<> struct FourthDensePt< 0, 0,-2> { enum { idx = 42 }; };
654 template<> struct FourthDensePt< 1, 0,-2> { enum { idx = 43 }; };
655 template<> struct FourthDensePt< 2, 0,-2> { enum { idx = 44 }; };
656
657
658 template<> struct FourthDensePt< 0,-2, 2> { enum { idx = 45 }; };
659 template<> struct FourthDensePt< 0,-1, 2> { enum { idx = 46 }; };
660 template<> struct FourthDensePt< 0, 1, 2> { enum { idx = 47 }; };
661 template<> struct FourthDensePt< 0, 2, 2> { enum { idx = 48 }; };
662
663 template<> struct FourthDensePt< 0,-2, 1> { enum { idx = 49 }; };
664 template<> struct FourthDensePt< 0,-1, 1> { enum { idx = 50 }; };
665 template<> struct FourthDensePt< 0, 1, 1> { enum { idx = 51 }; };
666 template<> struct FourthDensePt< 0, 2, 1> { enum { idx = 52 }; };
667
668 template<> struct FourthDensePt< 0,-2,-1> { enum { idx = 53 }; };
669 template<> struct FourthDensePt< 0,-1,-1> { enum { idx = 54 }; };
670 template<> struct FourthDensePt< 0, 1,-1> { enum { idx = 55 }; };
671 template<> struct FourthDensePt< 0, 2,-1> { enum { idx = 56 }; };
672
673 template<> struct FourthDensePt< 0,-2,-2> { enum { idx = 57 }; };
674 template<> struct FourthDensePt< 0,-1,-2> { enum { idx = 58 }; };
675 template<> struct FourthDensePt< 0, 1,-2> { enum { idx = 59 }; };
676 template<> struct FourthDensePt< 0, 2,-2> { enum { idx = 60 }; };
677
678}
679
680
681template<typename GridT, bool IsSafe = true>
683 : public BaseStencil<FourthOrderDenseStencil<GridT, IsSafe>, GridT, IsSafe>
684{
687public:
688 typedef GridT GridType;
689 typedef typename GridT::TreeType TreeType;
690 typedef typename GridType::ValueType ValueType;
691
692 static const int SIZE = 61;
693
694 FourthOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
695
696 /// Return linear offset for the specified stencil point relative to its center
697 template<int i, int j, int k>
698 unsigned int pos() const { return FourthDensePt<i,j,k>::idx; }
699
700private:
701 inline void init(const Coord& ijk)
702 {
703 mValues[FourthDensePt<-2, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 2, 0));
704 mValues[FourthDensePt<-1, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 2, 0));
705 mValues[FourthDensePt< 0, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 0));
706 mValues[FourthDensePt< 1, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 2, 0));
707 mValues[FourthDensePt< 2, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 2, 0));
708
709 mValues[FourthDensePt<-2, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 1, 0));
710 mValues[FourthDensePt<-1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 1, 0));
711 mValues[FourthDensePt< 0, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
712 mValues[FourthDensePt< 1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 1, 0));
713 mValues[FourthDensePt< 2, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 1, 0));
714
715 mValues[FourthDensePt<-2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 0));
716 mValues[FourthDensePt<-1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
717 mValues[FourthDensePt< 1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
718 mValues[FourthDensePt< 2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 0));
719
720 mValues[FourthDensePt<-2,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2,-1, 0));
721 mValues[FourthDensePt<-1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1,-1, 0));
722 mValues[FourthDensePt< 0,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 0));
723 mValues[FourthDensePt< 1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1,-1, 0));
724 mValues[FourthDensePt< 2,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2,-1, 0));
725
726 mValues[FourthDensePt<-2,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2,-2, 0));
727 mValues[FourthDensePt<-1,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1,-2, 0));
728 mValues[FourthDensePt< 0,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 0));
729 mValues[FourthDensePt< 1,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1,-2, 0));
730 mValues[FourthDensePt< 2,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2,-2, 0));
731
732 mValues[FourthDensePt<-2, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 2));
733 mValues[FourthDensePt<-1, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 2));
734 mValues[FourthDensePt< 0, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 2));
735 mValues[FourthDensePt< 1, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 2));
736 mValues[FourthDensePt< 2, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 2));
737
738 mValues[FourthDensePt<-2, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 1));
739 mValues[FourthDensePt<-1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 1));
740 mValues[FourthDensePt< 0, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
741 mValues[FourthDensePt< 1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 1));
742 mValues[FourthDensePt< 2, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 1));
743
744 mValues[FourthDensePt<-2, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0,-1));
745 mValues[FourthDensePt<-1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0,-1));
746 mValues[FourthDensePt< 0, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0,-1));
747 mValues[FourthDensePt< 1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0,-1));
748 mValues[FourthDensePt< 2, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0,-1));
749
750 mValues[FourthDensePt<-2, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0,-2));
751 mValues[FourthDensePt<-1, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0,-2));
752 mValues[FourthDensePt< 0, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0,-2));
753 mValues[FourthDensePt< 1, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0,-2));
754 mValues[FourthDensePt< 2, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0,-2));
755
756
757 mValues[FourthDensePt< 0,-2, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 2));
758 mValues[FourthDensePt< 0,-1, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 2));
759 mValues[FourthDensePt< 0, 1, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 2));
760 mValues[FourthDensePt< 0, 2, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 2));
761
762 mValues[FourthDensePt< 0,-2, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 1));
763 mValues[FourthDensePt< 0,-1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 1));
764 mValues[FourthDensePt< 0, 1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 1));
765 mValues[FourthDensePt< 0, 2, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 1));
766
767 mValues[FourthDensePt< 0,-2,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2,-1));
768 mValues[FourthDensePt< 0,-1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1,-1));
769 mValues[FourthDensePt< 0, 1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1,-1));
770 mValues[FourthDensePt< 0, 2,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2,-1));
771
772 mValues[FourthDensePt< 0,-2,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2,-2));
773 mValues[FourthDensePt< 0,-1,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1,-2));
774 mValues[FourthDensePt< 0, 1,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1,-2));
775 mValues[FourthDensePt< 0, 2,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2,-2));
776 }
777
778 template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
779 using BaseType::mAcc;
780 using BaseType::mValues;
781};// FourthOrderDenseStencil class
782
783
784////////////////////////////////////////
785
786
787namespace { // anonymous namespace for stencil-layout map
788
789 // the dense point stencil
790 template<int i, int j, int k> struct NineteenPt {};
791 template<> struct NineteenPt< 0, 0, 0> { enum { idx = 0 }; };
792
793 template<> struct NineteenPt< 1, 0, 0> { enum { idx = 1 }; };
794 template<> struct NineteenPt< 0, 1, 0> { enum { idx = 2 }; };
795 template<> struct NineteenPt< 0, 0, 1> { enum { idx = 3 }; };
796
797 template<> struct NineteenPt<-1, 0, 0> { enum { idx = 4 }; };
798 template<> struct NineteenPt< 0,-1, 0> { enum { idx = 5 }; };
799 template<> struct NineteenPt< 0, 0,-1> { enum { idx = 6 }; };
800
801 template<> struct NineteenPt< 2, 0, 0> { enum { idx = 7 }; };
802 template<> struct NineteenPt< 0, 2, 0> { enum { idx = 8 }; };
803 template<> struct NineteenPt< 0, 0, 2> { enum { idx = 9 }; };
804
805 template<> struct NineteenPt<-2, 0, 0> { enum { idx = 10 }; };
806 template<> struct NineteenPt< 0,-2, 0> { enum { idx = 11 }; };
807 template<> struct NineteenPt< 0, 0,-2> { enum { idx = 12 }; };
808
809 template<> struct NineteenPt< 3, 0, 0> { enum { idx = 13 }; };
810 template<> struct NineteenPt< 0, 3, 0> { enum { idx = 14 }; };
811 template<> struct NineteenPt< 0, 0, 3> { enum { idx = 15 }; };
812
813 template<> struct NineteenPt<-3, 0, 0> { enum { idx = 16 }; };
814 template<> struct NineteenPt< 0,-3, 0> { enum { idx = 17 }; };
815 template<> struct NineteenPt< 0, 0,-3> { enum { idx = 18 }; };
816
817}
818
819
820template<typename GridT, bool IsSafe = true>
822 : public BaseStencil<NineteenPointStencil<GridT, IsSafe>, GridT, IsSafe>
823{
826public:
827 typedef GridT GridType;
828 typedef typename GridT::TreeType TreeType;
829 typedef typename GridType::ValueType ValueType;
830
831 static const int SIZE = 19;
832
833 NineteenPointStencil(const GridType& grid): BaseType(grid, SIZE) {}
834
835 /// Return linear offset for the specified stencil point relative to its center
836 template<int i, int j, int k>
837 unsigned int pos() const { return NineteenPt<i,j,k>::idx; }
838
839private:
840 inline void init(const Coord& ijk)
841 {
842 mValues[NineteenPt< 3, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0, 0));
843 mValues[NineteenPt< 2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 0));
844 mValues[NineteenPt< 1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
845 mValues[NineteenPt<-1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
846 mValues[NineteenPt<-2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 0));
847 mValues[NineteenPt<-3, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0, 0));
848
849 mValues[NineteenPt< 0, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3, 0));
850 mValues[NineteenPt< 0, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 0));
851 mValues[NineteenPt< 0, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
852 mValues[NineteenPt< 0,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
853 mValues[NineteenPt< 0,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -2, 0));
854 mValues[NineteenPt< 0,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, -3, 0));
855
856 mValues[NineteenPt< 0, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 3));
857 mValues[NineteenPt< 0, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 2));
858 mValues[NineteenPt< 0, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
859 mValues[NineteenPt< 0, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
860 mValues[NineteenPt< 0, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -2));
861 mValues[NineteenPt< 0, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, -3));
862 }
863
864 template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
865 using BaseType::mAcc;
866 using BaseType::mValues;
867};// NineteenPointStencil class
868
869
870////////////////////////////////////////
871
872
873namespace { // anonymous namespace for stencil-layout map
874
875 // the 4th-order dense point stencil
876 template<int i, int j, int k> struct SixthDensePt { };
877 template<> struct SixthDensePt< 0, 0, 0> { enum { idx = 0 }; };
878
879 template<> struct SixthDensePt<-3, 3, 0> { enum { idx = 1 }; };
880 template<> struct SixthDensePt<-2, 3, 0> { enum { idx = 2 }; };
881 template<> struct SixthDensePt<-1, 3, 0> { enum { idx = 3 }; };
882 template<> struct SixthDensePt< 0, 3, 0> { enum { idx = 4 }; };
883 template<> struct SixthDensePt< 1, 3, 0> { enum { idx = 5 }; };
884 template<> struct SixthDensePt< 2, 3, 0> { enum { idx = 6 }; };
885 template<> struct SixthDensePt< 3, 3, 0> { enum { idx = 7 }; };
886
887 template<> struct SixthDensePt<-3, 2, 0> { enum { idx = 8 }; };
888 template<> struct SixthDensePt<-2, 2, 0> { enum { idx = 9 }; };
889 template<> struct SixthDensePt<-1, 2, 0> { enum { idx = 10 }; };
890 template<> struct SixthDensePt< 0, 2, 0> { enum { idx = 11 }; };
891 template<> struct SixthDensePt< 1, 2, 0> { enum { idx = 12 }; };
892 template<> struct SixthDensePt< 2, 2, 0> { enum { idx = 13 }; };
893 template<> struct SixthDensePt< 3, 2, 0> { enum { idx = 14 }; };
894
895 template<> struct SixthDensePt<-3, 1, 0> { enum { idx = 15 }; };
896 template<> struct SixthDensePt<-2, 1, 0> { enum { idx = 16 }; };
897 template<> struct SixthDensePt<-1, 1, 0> { enum { idx = 17 }; };
898 template<> struct SixthDensePt< 0, 1, 0> { enum { idx = 18 }; };
899 template<> struct SixthDensePt< 1, 1, 0> { enum { idx = 19 }; };
900 template<> struct SixthDensePt< 2, 1, 0> { enum { idx = 20 }; };
901 template<> struct SixthDensePt< 3, 1, 0> { enum { idx = 21 }; };
902
903 template<> struct SixthDensePt<-3, 0, 0> { enum { idx = 22 }; };
904 template<> struct SixthDensePt<-2, 0, 0> { enum { idx = 23 }; };
905 template<> struct SixthDensePt<-1, 0, 0> { enum { idx = 24 }; };
906 template<> struct SixthDensePt< 1, 0, 0> { enum { idx = 25 }; };
907 template<> struct SixthDensePt< 2, 0, 0> { enum { idx = 26 }; };
908 template<> struct SixthDensePt< 3, 0, 0> { enum { idx = 27 }; };
909
910
911 template<> struct SixthDensePt<-3,-1, 0> { enum { idx = 28 }; };
912 template<> struct SixthDensePt<-2,-1, 0> { enum { idx = 29 }; };
913 template<> struct SixthDensePt<-1,-1, 0> { enum { idx = 30 }; };
914 template<> struct SixthDensePt< 0,-1, 0> { enum { idx = 31 }; };
915 template<> struct SixthDensePt< 1,-1, 0> { enum { idx = 32 }; };
916 template<> struct SixthDensePt< 2,-1, 0> { enum { idx = 33 }; };
917 template<> struct SixthDensePt< 3,-1, 0> { enum { idx = 34 }; };
918
919
920 template<> struct SixthDensePt<-3,-2, 0> { enum { idx = 35 }; };
921 template<> struct SixthDensePt<-2,-2, 0> { enum { idx = 36 }; };
922 template<> struct SixthDensePt<-1,-2, 0> { enum { idx = 37 }; };
923 template<> struct SixthDensePt< 0,-2, 0> { enum { idx = 38 }; };
924 template<> struct SixthDensePt< 1,-2, 0> { enum { idx = 39 }; };
925 template<> struct SixthDensePt< 2,-2, 0> { enum { idx = 40 }; };
926 template<> struct SixthDensePt< 3,-2, 0> { enum { idx = 41 }; };
927
928
929 template<> struct SixthDensePt<-3,-3, 0> { enum { idx = 42 }; };
930 template<> struct SixthDensePt<-2,-3, 0> { enum { idx = 43 }; };
931 template<> struct SixthDensePt<-1,-3, 0> { enum { idx = 44 }; };
932 template<> struct SixthDensePt< 0,-3, 0> { enum { idx = 45 }; };
933 template<> struct SixthDensePt< 1,-3, 0> { enum { idx = 46 }; };
934 template<> struct SixthDensePt< 2,-3, 0> { enum { idx = 47 }; };
935 template<> struct SixthDensePt< 3,-3, 0> { enum { idx = 48 }; };
936
937
938 template<> struct SixthDensePt<-3, 0, 3> { enum { idx = 49 }; };
939 template<> struct SixthDensePt<-2, 0, 3> { enum { idx = 50 }; };
940 template<> struct SixthDensePt<-1, 0, 3> { enum { idx = 51 }; };
941 template<> struct SixthDensePt< 0, 0, 3> { enum { idx = 52 }; };
942 template<> struct SixthDensePt< 1, 0, 3> { enum { idx = 53 }; };
943 template<> struct SixthDensePt< 2, 0, 3> { enum { idx = 54 }; };
944 template<> struct SixthDensePt< 3, 0, 3> { enum { idx = 55 }; };
945
946
947 template<> struct SixthDensePt<-3, 0, 2> { enum { idx = 56 }; };
948 template<> struct SixthDensePt<-2, 0, 2> { enum { idx = 57 }; };
949 template<> struct SixthDensePt<-1, 0, 2> { enum { idx = 58 }; };
950 template<> struct SixthDensePt< 0, 0, 2> { enum { idx = 59 }; };
951 template<> struct SixthDensePt< 1, 0, 2> { enum { idx = 60 }; };
952 template<> struct SixthDensePt< 2, 0, 2> { enum { idx = 61 }; };
953 template<> struct SixthDensePt< 3, 0, 2> { enum { idx = 62 }; };
954
955 template<> struct SixthDensePt<-3, 0, 1> { enum { idx = 63 }; };
956 template<> struct SixthDensePt<-2, 0, 1> { enum { idx = 64 }; };
957 template<> struct SixthDensePt<-1, 0, 1> { enum { idx = 65 }; };
958 template<> struct SixthDensePt< 0, 0, 1> { enum { idx = 66 }; };
959 template<> struct SixthDensePt< 1, 0, 1> { enum { idx = 67 }; };
960 template<> struct SixthDensePt< 2, 0, 1> { enum { idx = 68 }; };
961 template<> struct SixthDensePt< 3, 0, 1> { enum { idx = 69 }; };
962
963
964 template<> struct SixthDensePt<-3, 0,-1> { enum { idx = 70 }; };
965 template<> struct SixthDensePt<-2, 0,-1> { enum { idx = 71 }; };
966 template<> struct SixthDensePt<-1, 0,-1> { enum { idx = 72 }; };
967 template<> struct SixthDensePt< 0, 0,-1> { enum { idx = 73 }; };
968 template<> struct SixthDensePt< 1, 0,-1> { enum { idx = 74 }; };
969 template<> struct SixthDensePt< 2, 0,-1> { enum { idx = 75 }; };
970 template<> struct SixthDensePt< 3, 0,-1> { enum { idx = 76 }; };
971
972
973 template<> struct SixthDensePt<-3, 0,-2> { enum { idx = 77 }; };
974 template<> struct SixthDensePt<-2, 0,-2> { enum { idx = 78 }; };
975 template<> struct SixthDensePt<-1, 0,-2> { enum { idx = 79 }; };
976 template<> struct SixthDensePt< 0, 0,-2> { enum { idx = 80 }; };
977 template<> struct SixthDensePt< 1, 0,-2> { enum { idx = 81 }; };
978 template<> struct SixthDensePt< 2, 0,-2> { enum { idx = 82 }; };
979 template<> struct SixthDensePt< 3, 0,-2> { enum { idx = 83 }; };
980
981
982 template<> struct SixthDensePt<-3, 0,-3> { enum { idx = 84 }; };
983 template<> struct SixthDensePt<-2, 0,-3> { enum { idx = 85 }; };
984 template<> struct SixthDensePt<-1, 0,-3> { enum { idx = 86 }; };
985 template<> struct SixthDensePt< 0, 0,-3> { enum { idx = 87 }; };
986 template<> struct SixthDensePt< 1, 0,-3> { enum { idx = 88 }; };
987 template<> struct SixthDensePt< 2, 0,-3> { enum { idx = 89 }; };
988 template<> struct SixthDensePt< 3, 0,-3> { enum { idx = 90 }; };
989
990
991 template<> struct SixthDensePt< 0,-3, 3> { enum { idx = 91 }; };
992 template<> struct SixthDensePt< 0,-2, 3> { enum { idx = 92 }; };
993 template<> struct SixthDensePt< 0,-1, 3> { enum { idx = 93 }; };
994 template<> struct SixthDensePt< 0, 1, 3> { enum { idx = 94 }; };
995 template<> struct SixthDensePt< 0, 2, 3> { enum { idx = 95 }; };
996 template<> struct SixthDensePt< 0, 3, 3> { enum { idx = 96 }; };
997
998 template<> struct SixthDensePt< 0,-3, 2> { enum { idx = 97 }; };
999 template<> struct SixthDensePt< 0,-2, 2> { enum { idx = 98 }; };
1000 template<> struct SixthDensePt< 0,-1, 2> { enum { idx = 99 }; };
1001 template<> struct SixthDensePt< 0, 1, 2> { enum { idx = 100 }; };
1002 template<> struct SixthDensePt< 0, 2, 2> { enum { idx = 101 }; };
1003 template<> struct SixthDensePt< 0, 3, 2> { enum { idx = 102 }; };
1004
1005 template<> struct SixthDensePt< 0,-3, 1> { enum { idx = 103 }; };
1006 template<> struct SixthDensePt< 0,-2, 1> { enum { idx = 104 }; };
1007 template<> struct SixthDensePt< 0,-1, 1> { enum { idx = 105 }; };
1008 template<> struct SixthDensePt< 0, 1, 1> { enum { idx = 106 }; };
1009 template<> struct SixthDensePt< 0, 2, 1> { enum { idx = 107 }; };
1010 template<> struct SixthDensePt< 0, 3, 1> { enum { idx = 108 }; };
1011
1012 template<> struct SixthDensePt< 0,-3,-1> { enum { idx = 109 }; };
1013 template<> struct SixthDensePt< 0,-2,-1> { enum { idx = 110 }; };
1014 template<> struct SixthDensePt< 0,-1,-1> { enum { idx = 111 }; };
1015 template<> struct SixthDensePt< 0, 1,-1> { enum { idx = 112 }; };
1016 template<> struct SixthDensePt< 0, 2,-1> { enum { idx = 113 }; };
1017 template<> struct SixthDensePt< 0, 3,-1> { enum { idx = 114 }; };
1018
1019 template<> struct SixthDensePt< 0,-3,-2> { enum { idx = 115 }; };
1020 template<> struct SixthDensePt< 0,-2,-2> { enum { idx = 116 }; };
1021 template<> struct SixthDensePt< 0,-1,-2> { enum { idx = 117 }; };
1022 template<> struct SixthDensePt< 0, 1,-2> { enum { idx = 118 }; };
1023 template<> struct SixthDensePt< 0, 2,-2> { enum { idx = 119 }; };
1024 template<> struct SixthDensePt< 0, 3,-2> { enum { idx = 120 }; };
1025
1026 template<> struct SixthDensePt< 0,-3,-3> { enum { idx = 121 }; };
1027 template<> struct SixthDensePt< 0,-2,-3> { enum { idx = 122 }; };
1028 template<> struct SixthDensePt< 0,-1,-3> { enum { idx = 123 }; };
1029 template<> struct SixthDensePt< 0, 1,-3> { enum { idx = 124 }; };
1030 template<> struct SixthDensePt< 0, 2,-3> { enum { idx = 125 }; };
1031 template<> struct SixthDensePt< 0, 3,-3> { enum { idx = 126 }; };
1032
1033}
1034
1035
1036template<typename GridT, bool IsSafe = true>
1038 : public BaseStencil<SixthOrderDenseStencil<GridT, IsSafe>, GridT, IsSafe>
1039{
1042public:
1043 typedef GridT GridType;
1044 typedef typename GridT::TreeType TreeType;
1045 typedef typename GridType::ValueType ValueType;
1046
1047 static const int SIZE = 127;
1048
1049 SixthOrderDenseStencil(const GridType& grid): BaseType(grid, SIZE) {}
1050
1051 /// Return linear offset for the specified stencil point relative to its center
1052 template<int i, int j, int k>
1053 unsigned int pos() const { return SixthDensePt<i,j,k>::idx; }
1054
1055private:
1056 inline void init(const Coord& ijk)
1057 {
1058 mValues[SixthDensePt<-3, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3, 3, 0));
1059 mValues[SixthDensePt<-2, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 3, 0));
1060 mValues[SixthDensePt<-1, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 3, 0));
1061 mValues[SixthDensePt< 0, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3, 0));
1062 mValues[SixthDensePt< 1, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 3, 0));
1063 mValues[SixthDensePt< 2, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 3, 0));
1064 mValues[SixthDensePt< 3, 3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3, 3, 0));
1065
1066 mValues[SixthDensePt<-3, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3, 2, 0));
1067 mValues[SixthDensePt<-2, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 2, 0));
1068 mValues[SixthDensePt<-1, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 2, 0));
1069 mValues[SixthDensePt< 0, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 0));
1070 mValues[SixthDensePt< 1, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 2, 0));
1071 mValues[SixthDensePt< 2, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 2, 0));
1072 mValues[SixthDensePt< 3, 2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3, 2, 0));
1073
1074 mValues[SixthDensePt<-3, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3, 1, 0));
1075 mValues[SixthDensePt<-2, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 1, 0));
1076 mValues[SixthDensePt<-1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 1, 0));
1077 mValues[SixthDensePt< 0, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
1078 mValues[SixthDensePt< 1, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 1, 0));
1079 mValues[SixthDensePt< 2, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 1, 0));
1080 mValues[SixthDensePt< 3, 1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3, 1, 0));
1081
1082 mValues[SixthDensePt<-3, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0, 0));
1083 mValues[SixthDensePt<-2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 0));
1084 mValues[SixthDensePt<-1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
1085 mValues[SixthDensePt< 1, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
1086 mValues[SixthDensePt< 2, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 0));
1087 mValues[SixthDensePt< 3, 0, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0, 0));
1088
1089 mValues[SixthDensePt<-3,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3,-1, 0));
1090 mValues[SixthDensePt<-2,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2,-1, 0));
1091 mValues[SixthDensePt<-1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1,-1, 0));
1092 mValues[SixthDensePt< 0,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 0));
1093 mValues[SixthDensePt< 1,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1,-1, 0));
1094 mValues[SixthDensePt< 2,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2,-1, 0));
1095 mValues[SixthDensePt< 3,-1, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3,-1, 0));
1096
1097 mValues[SixthDensePt<-3,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3,-2, 0));
1098 mValues[SixthDensePt<-2,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2,-2, 0));
1099 mValues[SixthDensePt<-1,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1,-2, 0));
1100 mValues[SixthDensePt< 0,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 0));
1101 mValues[SixthDensePt< 1,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1,-2, 0));
1102 mValues[SixthDensePt< 2,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2,-2, 0));
1103 mValues[SixthDensePt< 3,-2, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3,-2, 0));
1104
1105 mValues[SixthDensePt<-3,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-3,-3, 0));
1106 mValues[SixthDensePt<-2,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-2,-3, 0));
1107 mValues[SixthDensePt<-1,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy(-1,-3, 0));
1108 mValues[SixthDensePt< 0,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3, 0));
1109 mValues[SixthDensePt< 1,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 1,-3, 0));
1110 mValues[SixthDensePt< 2,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 2,-3, 0));
1111 mValues[SixthDensePt< 3,-3, 0>::idx] = mAcc.getValue(ijk.offsetBy( 3,-3, 0));
1112
1113 mValues[SixthDensePt<-3, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0, 3));
1114 mValues[SixthDensePt<-2, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 3));
1115 mValues[SixthDensePt<-1, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 3));
1116 mValues[SixthDensePt< 0, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 3));
1117 mValues[SixthDensePt< 1, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 3));
1118 mValues[SixthDensePt< 2, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 3));
1119 mValues[SixthDensePt< 3, 0, 3>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0, 3));
1120
1121 mValues[SixthDensePt<-3, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0, 2));
1122 mValues[SixthDensePt<-2, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 2));
1123 mValues[SixthDensePt<-1, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 2));
1124 mValues[SixthDensePt< 0, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 2));
1125 mValues[SixthDensePt< 1, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 2));
1126 mValues[SixthDensePt< 2, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 2));
1127 mValues[SixthDensePt< 3, 0, 2>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0, 2));
1128
1129 mValues[SixthDensePt<-3, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0, 1));
1130 mValues[SixthDensePt<-2, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0, 1));
1131 mValues[SixthDensePt<-1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0, 1));
1132 mValues[SixthDensePt< 0, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
1133 mValues[SixthDensePt< 1, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0, 1));
1134 mValues[SixthDensePt< 2, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0, 1));
1135 mValues[SixthDensePt< 3, 0, 1>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0, 1));
1136
1137 mValues[SixthDensePt<-3, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0,-1));
1138 mValues[SixthDensePt<-2, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0,-1));
1139 mValues[SixthDensePt<-1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0,-1));
1140 mValues[SixthDensePt< 0, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0,-1));
1141 mValues[SixthDensePt< 1, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0,-1));
1142 mValues[SixthDensePt< 2, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0,-1));
1143 mValues[SixthDensePt< 3, 0,-1>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0,-1));
1144
1145 mValues[SixthDensePt<-3, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0,-2));
1146 mValues[SixthDensePt<-2, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0,-2));
1147 mValues[SixthDensePt<-1, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0,-2));
1148 mValues[SixthDensePt< 0, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0,-2));
1149 mValues[SixthDensePt< 1, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0,-2));
1150 mValues[SixthDensePt< 2, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0,-2));
1151 mValues[SixthDensePt< 3, 0,-2>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0,-2));
1152
1153 mValues[SixthDensePt<-3, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy(-3, 0,-3));
1154 mValues[SixthDensePt<-2, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy(-2, 0,-3));
1155 mValues[SixthDensePt<-1, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy(-1, 0,-3));
1156 mValues[SixthDensePt< 0, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 0,-3));
1157 mValues[SixthDensePt< 1, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy( 1, 0,-3));
1158 mValues[SixthDensePt< 2, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy( 2, 0,-3));
1159 mValues[SixthDensePt< 3, 0,-3>::idx] = mAcc.getValue(ijk.offsetBy( 3, 0,-3));
1160
1161 mValues[SixthDensePt< 0,-3, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3, 3));
1162 mValues[SixthDensePt< 0,-2, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 3));
1163 mValues[SixthDensePt< 0,-1, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 3));
1164 mValues[SixthDensePt< 0, 1, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 3));
1165 mValues[SixthDensePt< 0, 2, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 3));
1166 mValues[SixthDensePt< 0, 3, 3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3, 3));
1167
1168 mValues[SixthDensePt< 0,-3, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3, 2));
1169 mValues[SixthDensePt< 0,-2, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 2));
1170 mValues[SixthDensePt< 0,-1, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 2));
1171 mValues[SixthDensePt< 0, 1, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 2));
1172 mValues[SixthDensePt< 0, 2, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 2));
1173 mValues[SixthDensePt< 0, 3, 2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3, 2));
1174
1175 mValues[SixthDensePt< 0,-3, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3, 1));
1176 mValues[SixthDensePt< 0,-2, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2, 1));
1177 mValues[SixthDensePt< 0,-1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1, 1));
1178 mValues[SixthDensePt< 0, 1, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1, 1));
1179 mValues[SixthDensePt< 0, 2, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2, 1));
1180 mValues[SixthDensePt< 0, 3, 1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3, 1));
1181
1182 mValues[SixthDensePt< 0,-3,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3,-1));
1183 mValues[SixthDensePt< 0,-2,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2,-1));
1184 mValues[SixthDensePt< 0,-1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1,-1));
1185 mValues[SixthDensePt< 0, 1,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1,-1));
1186 mValues[SixthDensePt< 0, 2,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2,-1));
1187 mValues[SixthDensePt< 0, 3,-1>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3,-1));
1188
1189 mValues[SixthDensePt< 0,-3,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3,-2));
1190 mValues[SixthDensePt< 0,-2,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2,-2));
1191 mValues[SixthDensePt< 0,-1,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1,-2));
1192 mValues[SixthDensePt< 0, 1,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1,-2));
1193 mValues[SixthDensePt< 0, 2,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2,-2));
1194 mValues[SixthDensePt< 0, 3,-2>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3,-2));
1195
1196 mValues[SixthDensePt< 0,-3,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-3,-3));
1197 mValues[SixthDensePt< 0,-2,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-2,-3));
1198 mValues[SixthDensePt< 0,-1,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0,-1,-3));
1199 mValues[SixthDensePt< 0, 1,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 1,-3));
1200 mValues[SixthDensePt< 0, 2,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 2,-3));
1201 mValues[SixthDensePt< 0, 3,-3>::idx] = mAcc.getValue(ijk.offsetBy( 0, 3,-3));
1202 }
1203
1204 template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
1205 using BaseType::mAcc;
1206 using BaseType::mValues;
1207};// SixthOrderDenseStencil class
1208
1209
1210//////////////////////////////////////////////////////////////////////
1211
1212namespace { // anonymous namespace for stencil-layout map
1213
1214 // the seven point stencil with a different layout from SevenPt
1215 template<int i, int j, int k> struct GradPt {};
1216 template<> struct GradPt< 0, 0, 0> { enum { idx = 0 }; };
1217 template<> struct GradPt< 1, 0, 0> { enum { idx = 2 }; };
1218 template<> struct GradPt< 0, 1, 0> { enum { idx = 4 }; };
1219 template<> struct GradPt< 0, 0, 1> { enum { idx = 6 }; };
1220 template<> struct GradPt<-1, 0, 0> { enum { idx = 1 }; };
1221 template<> struct GradPt< 0,-1, 0> { enum { idx = 3 }; };
1222 template<> struct GradPt< 0, 0,-1> { enum { idx = 5 }; };
1223}
1224
1225/// This is a simple 7-point nearest neighbor stencil that supports
1226/// gradient by second-order central differencing, first-order upwinding,
1227/// Laplacian, closest-point transform and zero-crossing test.
1228///
1229/// @note For optimal random access performance this class
1230/// includes its own grid accessor.
1231template<typename GridT, bool IsSafe = true>
1232class GradStencil : public BaseStencil<GradStencil<GridT, IsSafe>, GridT, IsSafe>
1233{
1236public:
1237 typedef GridT GridType;
1238 typedef typename GridT::TreeType TreeType;
1239 typedef typename GridType::ValueType ValueType;
1240
1241 static const int SIZE = 7;
1242
1244 : BaseType(grid, SIZE)
1245 , mInv2Dx(ValueType(0.5 / grid.voxelSize()[0]))
1246 , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1247 {
1248 }
1249
1250 GradStencil(const GridType& grid, Real dx)
1251 : BaseType(grid, SIZE)
1252 , mInv2Dx(ValueType(0.5 / dx))
1253 , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1254 {
1255 }
1256
1257 /// @brief Return the norm square of the single-sided upwind gradient
1258 /// (computed via Godunov's scheme) at the previously buffered location.
1259 ///
1260 /// @note This method should not be called until the stencil
1261 /// buffer has been populated via a call to moveTo(ijk).
1262 inline ValueType normSqGrad() const
1263 {
1264 return mInvDx2 * math::GodunovsNormSqrd(mValues[0] > zeroVal<ValueType>(),
1265 mValues[0] - mValues[1],
1266 mValues[2] - mValues[0],
1267 mValues[0] - mValues[3],
1268 mValues[4] - mValues[0],
1269 mValues[0] - mValues[5],
1270 mValues[6] - mValues[0]);
1271 }
1272
1273 /// @brief Return the gradient computed at the previously buffered
1274 /// location by second order central differencing.
1275 ///
1276 /// @note This method should not be called until the stencil
1277 /// buffer has been populated via a call to moveTo(ijk).
1279 {
1280 return math::Vec3<ValueType>(mValues[2] - mValues[1],
1281 mValues[4] - mValues[3],
1282 mValues[6] - mValues[5])*mInv2Dx;
1283 }
1284 /// @brief Return the first-order upwind gradient corresponding to the direction V.
1285 ///
1286 /// @note This method should not be called until the stencil
1287 /// buffer has been populated via a call to moveTo(ijk).
1289 {
1290 return math::Vec3<ValueType>(
1291 V[0]>0 ? mValues[0] - mValues[1] : mValues[2] - mValues[0],
1292 V[1]>0 ? mValues[0] - mValues[3] : mValues[4] - mValues[0],
1293 V[2]>0 ? mValues[0] - mValues[5] : mValues[6] - mValues[0])*2*mInv2Dx;
1294 }
1295
1296 /// Return the Laplacian computed at the previously buffered
1297 /// location by second-order central differencing.
1298 inline ValueType laplacian() const
1299 {
1300 return mInvDx2 * (mValues[1] + mValues[2] +
1301 mValues[3] + mValues[4] +
1302 mValues[5] + mValues[6] - 6*mValues[0]);
1303 }
1304
1305 /// Return @c true if the sign of the value at the center point of the stencil
1306 /// is different from the signs of any of its six nearest neighbors.
1307 inline bool zeroCrossing() const
1308 {
1309 const typename BaseType::BufferType& v = mValues;
1310 return (v[0]>0 ? (v[1]<0 || v[2]<0 || v[3]<0 || v[4]<0 || v[5]<0 || v[6]<0)
1311 : (v[1]>0 || v[2]>0 || v[3]>0 || v[4]>0 || v[5]>0 || v[6]>0));
1312 }
1313
1314 /// @brief Compute the closest-point transform to a level set.
1315 /// @return the closest point in index space to the surface
1316 /// from which the level set was derived.
1317 ///
1318 /// @note This method assumes that the grid represents a level set
1319 /// with distances in world units and a simple affine transfrom
1320 /// with uniform scaling.
1322 {
1323 const Coord& ijk = BaseType::getCenterCoord();
1324 const ValueType d = ValueType(mValues[0] * 0.5 * mInvDx2); // distance in voxels / (2dx^2)
1326 const auto value = math::Vec3<ValueType>( ijk[0] - d*(mValues[2] - mValues[1]),
1327 ijk[1] - d*(mValues[4] - mValues[3]),
1328 ijk[2] - d*(mValues[6] - mValues[5]));
1330 return value;
1331 }
1332
1333 /// Return linear offset for the specified stencil point relative to its center
1334 template<int i, int j, int k>
1335 unsigned int pos() const { return GradPt<i,j,k>::idx; }
1336
1337private:
1338
1339 inline void init(const Coord& ijk)
1340 {
1341 BaseType::template setValue<-1, 0, 0>(mAcc.getValue(ijk.offsetBy(-1, 0, 0)));
1342 BaseType::template setValue< 1, 0, 0>(mAcc.getValue(ijk.offsetBy( 1, 0, 0)));
1343
1344 BaseType::template setValue< 0,-1, 0>(mAcc.getValue(ijk.offsetBy( 0,-1, 0)));
1345 BaseType::template setValue< 0, 1, 0>(mAcc.getValue(ijk.offsetBy( 0, 1, 0)));
1346
1347 BaseType::template setValue< 0, 0,-1>(mAcc.getValue(ijk.offsetBy( 0, 0,-1)));
1348 BaseType::template setValue< 0, 0, 1>(mAcc.getValue(ijk.offsetBy( 0, 0, 1)));
1349 }
1350
1351 template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
1352 using BaseType::mAcc;
1353 using BaseType::mValues;
1354 const ValueType mInv2Dx, mInvDx2;
1355}; // GradStencil class
1356
1357////////////////////////////////////////
1358
1359
1360/// @brief This is a special 19-point stencil that supports optimal fifth-order WENO
1361/// upwinding, second-order central differencing, Laplacian, and zero-crossing test.
1362///
1363/// @note For optimal random access performance this class
1364/// includes its own grid accessor.
1365template<typename GridT, bool IsSafe = true>
1366class WenoStencil: public BaseStencil<WenoStencil<GridT, IsSafe>, GridT, IsSafe>
1367{
1370public:
1371 typedef GridT GridType;
1372 typedef typename GridT::TreeType TreeType;
1373 typedef typename GridType::ValueType ValueType;
1374
1375 static const int SIZE = 19;
1376
1378 : BaseType(grid, SIZE)
1379 , _mDx2(ValueType(math::Pow2(grid.voxelSize()[0])))
1380 , mInv2Dx(ValueType(0.5 / grid.voxelSize()[0]))
1381 , mInvDx2(ValueType(1.0 / _mDx2))
1382 , mDx2(static_cast<float>(_mDx2))
1383 {
1384 }
1385
1386 WenoStencil(const GridType& grid, Real dx)
1387 : BaseType(grid, SIZE)
1388 , _mDx2(ValueType(dx * dx))
1389 , mInv2Dx(ValueType(0.5 / dx))
1390 , mInvDx2(ValueType(1.0 / _mDx2))
1391 , mDx2(static_cast<float>(_mDx2))
1392 {
1393 }
1394
1395 /// @brief Return the norm-square of the WENO upwind gradient (computed via
1396 /// WENO upwinding and Godunov's scheme) at the previously buffered location.
1397 ///
1398 /// @note This method should not be called until the stencil
1399 /// buffer has been populated via a call to moveTo(ijk).
1400 inline ValueType normSqGrad(const ValueType &isoValue = zeroVal<ValueType>()) const
1401 {
1402 const typename BaseType::BufferType& v = mValues;
1403#ifdef DWA_OPENVDB
1404 // SSE optimized
1405 const simd::Float4
1406 v1(v[2]-v[1], v[ 8]-v[ 7], v[14]-v[13], 0),
1407 v2(v[3]-v[2], v[ 9]-v[ 8], v[15]-v[14], 0),
1408 v3(v[0]-v[3], v[ 0]-v[ 9], v[ 0]-v[15], 0),
1409 v4(v[4]-v[0], v[10]-v[ 0], v[16]-v[ 0], 0),
1410 v5(v[5]-v[4], v[11]-v[10], v[17]-v[16], 0),
1411 v6(v[6]-v[5], v[12]-v[11], v[18]-v[17], 0),
1412 dP_m = math::WENO5(v1, v2, v3, v4, v5, mDx2),
1413 dP_p = math::WENO5(v6, v5, v4, v3, v2, mDx2);
1414
1415 return mInvDx2 * math::GodunovsNormSqrd(mValues[0] > isoValue, dP_m, dP_p);
1416#else
1417 const Real
1418 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),
1419 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),
1420 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),
1421 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),
1422 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),
1423 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);
1424 return static_cast<ValueType>(
1425 mInvDx2*math::GodunovsNormSqrd(v[0]>isoValue, dP_xm, dP_xp, dP_ym, dP_yp, dP_zm, dP_zp));
1426#endif
1427 }
1428
1429 /// Return the optimal fifth-order upwind gradient corresponding to the
1430 /// direction V.
1431 ///
1432 /// @note This method should not be called until the stencil
1433 /// buffer has been populated via a call to moveTo(ijk).
1435 {
1436 const typename BaseType::BufferType& v = mValues;
1437 return 2*mInv2Dx * math::Vec3<ValueType>(
1438 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)
1439 : math::WENO5(v[ 6]-v[ 5],v[ 5]-v[ 4],v[ 4]-v[ 0], v[ 0]-v[ 3],v[ 3]-v[ 2],mDx2),
1440 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)
1441 : math::WENO5(v[12]-v[11],v[11]-v[10],v[10]-v[ 0], v[ 0]-v[ 9],v[ 9]-v[ 8],mDx2),
1442 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)
1443 : math::WENO5(v[18]-v[17],v[17]-v[16],v[16]-v[ 0], v[ 0]-v[15],v[15]-v[14],mDx2));
1444 }
1445 /// Return the gradient computed at the previously buffered
1446 /// location by second-order central differencing.
1447 ///
1448 /// @note This method should not be called until the stencil
1449 /// buffer has been populated via a call to moveTo(ijk).
1451 {
1452 return mInv2Dx * math::Vec3<ValueType>(mValues[ 4] - mValues[ 3],
1453 mValues[10] - mValues[ 9],
1454 mValues[16] - mValues[15]);
1455 }
1456
1457 /// Return the Laplacian computed at the previously buffered
1458 /// location by second-order central differencing.
1459 ///
1460 /// @note This method should not be called until the stencil
1461 /// buffer has been populated via a call to moveTo(ijk).
1462 inline ValueType laplacian() const
1463 {
1464 return mInvDx2 * (
1465 mValues[ 3] + mValues[ 4] +
1466 mValues[ 9] + mValues[10] +
1467 mValues[15] + mValues[16] - 6*mValues[0]);
1468 }
1469
1470 /// Return @c true if the sign of the value at the center point of the stencil
1471 /// differs from the sign of any of its six nearest neighbors
1472 inline bool zeroCrossing() const
1473 {
1474 const typename BaseType::BufferType& v = mValues;
1475 return (v[ 0]>0 ? (v[ 3]<0 || v[ 4]<0 || v[ 9]<0 || v[10]<0 || v[15]<0 || v[16]<0)
1476 : (v[ 3]>0 || v[ 4]>0 || v[ 9]>0 || v[10]>0 || v[15]>0 || v[16]>0));
1477 }
1478
1479private:
1480 inline void init(const Coord& ijk)
1481 {
1482 mValues[ 1] = mAcc.getValue(ijk.offsetBy(-3, 0, 0));
1483 mValues[ 2] = mAcc.getValue(ijk.offsetBy(-2, 0, 0));
1484 mValues[ 3] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
1485 mValues[ 4] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
1486 mValues[ 5] = mAcc.getValue(ijk.offsetBy( 2, 0, 0));
1487 mValues[ 6] = mAcc.getValue(ijk.offsetBy( 3, 0, 0));
1488
1489 mValues[ 7] = mAcc.getValue(ijk.offsetBy( 0, -3, 0));
1490 mValues[ 8] = mAcc.getValue(ijk.offsetBy( 0, -2, 0));
1491 mValues[ 9] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
1492 mValues[10] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
1493 mValues[11] = mAcc.getValue(ijk.offsetBy( 0, 2, 0));
1494 mValues[12] = mAcc.getValue(ijk.offsetBy( 0, 3, 0));
1495
1496 mValues[13] = mAcc.getValue(ijk.offsetBy( 0, 0, -3));
1497 mValues[14] = mAcc.getValue(ijk.offsetBy( 0, 0, -2));
1498 mValues[15] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
1499 mValues[16] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
1500 mValues[17] = mAcc.getValue(ijk.offsetBy( 0, 0, 2));
1501 mValues[18] = mAcc.getValue(ijk.offsetBy( 0, 0, 3));
1502 }
1503
1504 template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
1505 using BaseType::mAcc;
1506 using BaseType::mValues;
1507 const ValueType _mDx2, mInv2Dx, mInvDx2;
1508 const float mDx2;
1509}; // WenoStencil class
1510
1511
1512//////////////////////////////////////////////////////////////////////
1513
1514
1515template<typename GridT, bool IsSafe = true>
1516class CurvatureStencil: public BaseStencil<CurvatureStencil<GridT, IsSafe>, GridT, IsSafe>
1517{
1520public:
1521 typedef GridT GridType;
1522 typedef typename GridT::TreeType TreeType;
1523 typedef typename GridT::ValueType ValueType;
1524
1525 static const int SIZE = 19;
1526
1528 : BaseType(grid, SIZE)
1529 , mInv2Dx(ValueType(0.5 / grid.voxelSize()[0]))
1530 , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1531 {
1532 }
1533
1535 : BaseType(grid, SIZE)
1536 , mInv2Dx(ValueType(0.5 / dx))
1537 , mInvDx2(ValueType(4.0 * mInv2Dx * mInv2Dx))
1538 {
1539 }
1540
1541 /// @brief Return the mean curvature at the previously buffered location.
1542 ///
1543 /// @note This method should not be called until the stencil
1544 /// buffer has been populated via a call to moveTo(ijk).
1546 {
1547 Real alpha, normGrad;
1548 return this->meanCurvature(alpha, normGrad) ?
1549 ValueType(alpha*mInv2Dx/math::Pow3(normGrad)) : 0;
1550 }
1551
1552 /// @brief Return the Gaussian curvature at the previously buffered location.
1553 ///
1554 /// @note This method should not be called until the stencil
1555 /// buffer has been populated via a call to moveTo(ijk).
1557 {
1558 Real alpha, normGrad;
1559 return this->gaussianCurvature(alpha, normGrad) ?
1560 ValueType(alpha*mInvDx2/math::Pow4(normGrad)) : 0;
1561 }
1562
1563 /// @brief Return both the mean and the Gaussian curvature at the
1564 /// previously buffered location.
1565 ///
1566 /// @note This method should not be called until the stencil
1567 /// buffer has been populated via a call to moveTo(ijk).
1568 inline void curvatures(ValueType &mean, ValueType& gauss) const
1569 {
1570 Real alphaM, alphaG, normGrad;
1571 if (this->curvatures(alphaM, alphaG, normGrad)) {
1572 mean = ValueType(alphaM*mInv2Dx/math::Pow3(normGrad));
1573 gauss = ValueType(alphaG*mInvDx2/math::Pow4(normGrad));
1574 } else {
1575 mean = gauss = 0;
1576 }
1577 }
1578
1579 /// Return the mean curvature multiplied by the norm of the
1580 /// central-difference gradient. This method is very useful for
1581 /// mean-curvature flow of level sets!
1582 ///
1583 /// @note This method should not be called until the stencil
1584 /// buffer has been populated via a call to moveTo(ijk).
1586 {
1587 Real alpha, normGrad;
1588 return this->meanCurvature(alpha, normGrad) ?
1589 ValueType(alpha*mInvDx2/(2*math::Pow2(normGrad))) : 0;
1590 }
1591
1592 /// Return the mean Gaussian multiplied by the norm of the
1593 /// central-difference gradient.
1594 ///
1595 /// @note This method should not be called until the stencil
1596 /// buffer has been populated via a call to moveTo(ijk).
1598 {
1599 Real alpha, normGrad;
1600 return this->gaussianCurvature(alpha, normGrad) ?
1601 ValueType(2*alpha*mInv2Dx*mInvDx2/math::Pow3(normGrad)) : 0;
1602 }
1603
1604 /// @brief Return both the mean and the Gaussian curvature at the
1605 /// previously buffered location.
1606 ///
1607 /// @note This method should not be called until the stencil
1608 /// buffer has been populated via a call to moveTo(ijk).
1609 inline void curvaturesNormGrad(ValueType &mean, ValueType& gauss) const
1610 {
1611 Real alphaM, alphaG, normGrad;
1612 if (this->curvatures(alphaM, alphaG, normGrad)) {
1613 mean = ValueType(alphaM*mInvDx2/(2*math::Pow2(normGrad)));
1614 gauss = ValueType(2*alphaG*mInv2Dx*mInvDx2/math::Pow3(normGrad));
1615 } else {
1616 mean = gauss = 0;
1617 }
1618 }
1619
1620 /// @brief Return the pair (minimum, maximum) principal curvature at the
1621 /// previously buffered location.
1622 ///
1623 /// @note This method should not be called until the stencil
1624 /// buffer has been populated via a call to moveTo(ijk).
1625 inline std::pair<ValueType, ValueType> principalCurvatures() const
1626 {
1627 std::pair<ValueType, ValueType> pair(0, 0);// min, max
1628 Real alphaM, alphaG, normGrad;
1629 if (this->curvatures(alphaM, alphaG, normGrad)) {
1630 const Real mean = alphaM*mInv2Dx/math::Pow3(normGrad);
1631 const Real tmp = std::sqrt(mean*mean - alphaG*mInvDx2/math::Pow4(normGrad));
1632 pair.first = ValueType(mean - tmp);
1633 pair.second = ValueType(mean + tmp);
1634 }
1635 return pair;// min, max
1636 }
1637
1638 /// Return the Laplacian computed at the previously buffered
1639 /// location by second-order central differencing.
1640 ///
1641 /// @note This method should not be called until the stencil
1642 /// buffer has been populated via a call to moveTo(ijk).
1643 inline ValueType laplacian() const
1644 {
1645 return mInvDx2 * (
1646 mValues[1] + mValues[2] +
1647 mValues[3] + mValues[4] +
1648 mValues[5] + mValues[6] - 6*mValues[0]);
1649 }
1650
1651 /// Return the gradient computed at the previously buffered
1652 /// location by second-order central differencing.
1653 ///
1654 /// @note This method should not be called until the stencil
1655 /// buffer has been populated via a call to moveTo(ijk).
1657 {
1658 return math::Vec3<ValueType>(
1659 mValues[2] - mValues[1],
1660 mValues[4] - mValues[3],
1661 mValues[6] - mValues[5])*mInv2Dx;
1662 }
1663
1664private:
1665 inline void init(const Coord &ijk)
1666 {
1667 mValues[ 1] = mAcc.getValue(ijk.offsetBy(-1, 0, 0));
1668 mValues[ 2] = mAcc.getValue(ijk.offsetBy( 1, 0, 0));
1669
1670 mValues[ 3] = mAcc.getValue(ijk.offsetBy( 0, -1, 0));
1671 mValues[ 4] = mAcc.getValue(ijk.offsetBy( 0, 1, 0));
1672
1673 mValues[ 5] = mAcc.getValue(ijk.offsetBy( 0, 0, -1));
1674 mValues[ 6] = mAcc.getValue(ijk.offsetBy( 0, 0, 1));
1675
1676 mValues[ 7] = mAcc.getValue(ijk.offsetBy(-1, -1, 0));
1677 mValues[ 8] = mAcc.getValue(ijk.offsetBy( 1, -1, 0));
1678 mValues[ 9] = mAcc.getValue(ijk.offsetBy(-1, 1, 0));
1679 mValues[10] = mAcc.getValue(ijk.offsetBy( 1, 1, 0));
1680
1681 mValues[11] = mAcc.getValue(ijk.offsetBy(-1, 0, -1));
1682 mValues[12] = mAcc.getValue(ijk.offsetBy( 1, 0, -1));
1683 mValues[13] = mAcc.getValue(ijk.offsetBy(-1, 0, 1));
1684 mValues[14] = mAcc.getValue(ijk.offsetBy( 1, 0, 1));
1685
1686 mValues[15] = mAcc.getValue(ijk.offsetBy( 0, -1, -1));
1687 mValues[16] = mAcc.getValue(ijk.offsetBy( 0, 1, -1));
1688 mValues[17] = mAcc.getValue(ijk.offsetBy( 0, -1, 1));
1689 mValues[18] = mAcc.getValue(ijk.offsetBy( 0, 1, 1));
1690 }
1691
1692 inline Real Dx() const { return 0.5*(mValues[2] - mValues[1]); }// * 1/dx
1693 inline Real Dy() const { return 0.5*(mValues[4] - mValues[3]); }// * 1/dx
1694 inline Real Dz() const { return 0.5*(mValues[6] - mValues[5]); }// * 1/dx
1695 inline Real Dxx() const { return mValues[2] - 2 * mValues[0] + mValues[1]; }// * 1/dx2
1696 inline Real Dyy() const { return mValues[4] - 2 * mValues[0] + mValues[3]; }// * 1/dx2}
1697 inline Real Dzz() const { return mValues[6] - 2 * mValues[0] + mValues[5]; }// * 1/dx2
1698 inline Real Dxy() const { return 0.25 * (mValues[10] - mValues[ 8] + mValues[ 7] - mValues[ 9]); }// * 1/dx2
1699 inline Real Dxz() const { return 0.25 * (mValues[14] - mValues[12] + mValues[11] - mValues[13]); }// * 1/dx2
1700 inline Real Dyz() const { return 0.25 * (mValues[18] - mValues[16] + mValues[15] - mValues[17]); }// * 1/dx2
1701
1702 inline bool meanCurvature(Real& alpha, Real& normGrad) const
1703 {
1704 // For performance all finite differences are unscaled wrt dx
1705 const Real Dx = this->Dx(), Dy = this->Dy(), Dz = this->Dz(),
1706 Dx2 = Dx*Dx, Dy2 = Dy*Dy, Dz2 = Dz*Dz, normGrad2 = Dx2 + Dy2 + Dz2;
1707 if (normGrad2 <= math::Tolerance<Real>::value()) {
1708 alpha = normGrad = 0;
1709 return false;
1710 }
1711 const Real Dxx = this->Dxx(), Dyy = this->Dyy(), Dzz = this->Dzz();
1712 alpha = Dx2*(Dyy + Dzz) + Dy2*(Dxx + Dzz) + Dz2*(Dxx + Dyy) -
1713 2*(Dx*(Dy*this->Dxy() + Dz*this->Dxz()) + Dy*Dz*this->Dyz());// * 1/dx^4
1714 normGrad = std::sqrt(normGrad2); // * 1/dx
1715 return true;
1716 }
1717
1718 inline bool gaussianCurvature(Real& alpha, Real& normGrad) const
1719 {
1720 // For performance all finite differences are unscaled wrt dx
1721 const Real Dx = this->Dx(), Dy = this->Dy(), Dz = this->Dz(),
1722 Dx2 = Dx*Dx, Dy2 = Dy*Dy, Dz2 = Dz*Dz, normGrad2 = Dx2 + Dy2 + Dz2;
1723 if (normGrad2 <= math::Tolerance<Real>::value()) {
1724 alpha = normGrad = 0;
1725 return false;
1726 }
1727 const Real Dxx = this->Dxx(), Dyy = this->Dyy(), Dzz = this->Dzz(),
1728 Dxy = this->Dxy(), Dxz = this->Dxz(), Dyz = this->Dyz();
1729 alpha = Dx2*(Dyy*Dzz - Dyz*Dyz) + Dy2*(Dxx*Dzz - Dxz*Dxz) + Dz2*(Dxx*Dyy - Dxy*Dxy) +
1730 2*( Dy*Dz*(Dxy*Dxz - Dyz*Dxx) + Dx*Dz*(Dxy*Dyz - Dxz*Dyy) + Dx*Dy*(Dxz*Dyz - Dxy*Dzz) );// * 1/dx^6
1731 normGrad = std::sqrt(normGrad2); // * 1/dx
1732 return true;
1733 }
1734 inline bool curvatures(Real& alphaM, Real& alphaG, Real& normGrad) const
1735 {
1736 // For performance all finite differences are unscaled wrt dx
1737 const Real Dx = this->Dx(), Dy = this->Dy(), Dz = this->Dz(),
1738 Dx2 = Dx*Dx, Dy2 = Dy*Dy, Dz2 = Dz*Dz, normGrad2 = Dx2 + Dy2 + Dz2;
1739 if (normGrad2 <= math::Tolerance<Real>::value()) {
1740 alphaM = alphaG =normGrad = 0;
1741 return false;
1742 }
1743 const Real Dxx = this->Dxx(), Dyy = this->Dyy(), Dzz = this->Dzz(),
1744 Dxy = this->Dxy(), Dxz = this->Dxz(), Dyz = this->Dyz();
1745 alphaM = Dx2*(Dyy + Dzz) + Dy2*(Dxx + Dzz) + Dz2*(Dxx + Dyy) -
1746 2*(Dx*(Dy*Dxy + Dz*Dxz) + Dy*Dz*Dyz);// *1/dx^4
1747 alphaG = Dx2*(Dyy*Dzz - Dyz*Dyz) + Dy2*(Dxx*Dzz - Dxz*Dxz) + Dz2*(Dxx*Dyy - Dxy*Dxy) +
1748 2*( Dy*Dz*(Dxy*Dxz - Dyz*Dxx) + Dx*Dz*(Dxy*Dyz - Dxz*Dyy) + Dx*Dy*(Dxz*Dyz - Dxy*Dzz) );// *1/dx^6
1749 normGrad = std::sqrt(normGrad2); // * 1/dx
1750 return true;
1751 }
1752
1753 template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
1754 using BaseType::mAcc;
1755 using BaseType::mValues;
1756 const ValueType mInv2Dx, mInvDx2;
1757}; // CurvatureStencil class
1758
1759
1760//////////////////////////////////////////////////////////////////////
1761
1762
1763/// @brief Dense stencil of a given width
1764template<typename GridT, bool IsSafe = true>
1765class DenseStencil: public BaseStencil<DenseStencil<GridT, IsSafe>, GridT, IsSafe>
1766{
1769public:
1770 typedef GridT GridType;
1771 typedef typename GridT::TreeType TreeType;
1772 typedef typename GridType::ValueType ValueType;
1773
1774 DenseStencil(const GridType& grid, int halfWidth)
1775 : BaseType(grid, /*size=*/math::Pow3(2 * halfWidth + 1))
1776 , mHalfWidth(halfWidth)
1777 {
1778 assert(halfWidth>0);
1779 }
1780
1781 inline const ValueType& getCenterValue() const { return mValues[(mValues.size()-1)>>1]; }
1782
1783 /// @brief Initialize the stencil buffer with the values of voxel (x, y, z)
1784 /// and its neighbors.
1785 inline void moveTo(const Coord& ijk)
1786 {
1787 BaseType::mCenter = ijk;
1788 this->init(ijk);
1789 }
1790 /// @brief Initialize the stencil buffer with the values of voxel
1791 /// (x, y, z) and its neighbors.
1792 template<typename IterType>
1793 inline void moveTo(const IterType& iter)
1794 {
1795 BaseType::mCenter = iter.getCoord();
1796 this->init(BaseType::mCenter);
1797 }
1798
1799private:
1800 /// Initialize the stencil buffer centered at (i, j, k).
1801 /// @warning The center point is NOT at mValues[0] for this DenseStencil!
1802 inline void init(const Coord& ijk)
1803 {
1804 int n = 0;
1805 for (Coord p=ijk.offsetBy(-mHalfWidth), q=ijk.offsetBy(mHalfWidth); p[0] <= q[0]; ++p[0]) {
1806 for (p[1] = ijk[1]-mHalfWidth; p[1] <= q[1]; ++p[1]) {
1807 for (p[2] = ijk[2]-mHalfWidth; p[2] <= q[2]; ++p[2]) {
1808 mValues[n++] = mAcc.getValue(p);
1809 }
1810 }
1811 }
1812 }
1813
1814 template<typename, typename, bool> friend class BaseStencil; // allow base class to call init()
1815 using BaseType::mAcc;
1816 using BaseType::mValues;
1817 const int mHalfWidth;
1818};// DenseStencil class
1819
1820
1821} // end math namespace
1822} // namespace OPENVDB_VERSION_NAME
1823} // end openvdb namespace
1824
1825#endif // OPENVDB_MATH_STENCILS_HAS_BEEN_INCLUDED
ValueT value
Definition: GridBuilder.h:1290
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
Bracket code with OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN/_END, to inhibit warnings about type conve...
Definition: Platform.h:204
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_END
Definition: Platform.h:205
Definition: Stencils.h:36
ValueType min() const
Return the smallest value in the stencil buffer.
Definition: Stencils.h:147
void moveTo(const Vec3< RealType > &xyz)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:87
BufferType mValues
Definition: Stencils.h:222
const AccessorType & accessor() const
Return a const reference to the ValueAccessor associated with this Stencil.
Definition: Stencils.h:208
void moveTo(const IterType &iter)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:73
const ValueType & getValue(unsigned int pos=0) const
Return the value from the stencil buffer with linear offset pos.
Definition: Stencils.h:98
ValueType mean() const
Return the mean value of the current stencil.
Definition: Stencils.h:139
void setValue(const ValueType &value)
Set the value at the specified location relative to the center of the stencil.
Definition: Stencils.h:113
void moveTo(const Coord &ijk, const ValueType &centerValue)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors....
Definition: Stencils.h:60
BaseStencil(const GridType &grid, int size)
Definition: Stencils.h:212
GridT::ValueType ValueType
Definition: Stencils.h:40
std::bitset< 6 > intersectionMask(const ValueType &isoValue=zeroVal< ValueType >()) const
Return true a bit-mask where the 6 bits indicates if the center of the stencil intersects the iso-con...
Definition: Stencils.h:189
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:48
Coord mCenter
Definition: Stencils.h:223
const GridType * mGrid
Definition: Stencils.h:220
const ValueType & getCenterValue() const
Return the value at the center of the stencil.
Definition: Stencils.h:164
const ValueType & getValue() const
Return the value at the specified location relative to the center of the stencil.
Definition: Stencils.h:106
const Coord & getCenterCoord() const
Return the coordinates of the center point of the stencil.
Definition: Stencils.h:161
bool intersects(const ValueType &isoValue=zeroVal< ValueType >()) const
Return true if the center of the stencil intersects the iso-contour specified by the isoValue.
Definition: Stencils.h:168
tree::ValueAccessor< const TreeType, IsSafe > AccessorType
Definition: Stencils.h:41
BufferType::iterator IterType
Definition: Stencils.h:43
ValueType median() const
Return the median value of the current stencil.
Definition: Stencils.h:122
const GridType & grid() const
Return a const reference to the grid from which this stencil was constructed.
Definition: Stencils.h:204
ValueType max() const
Return the largest value in the stencil buffer.
Definition: Stencils.h:154
AccessorType mAcc
Definition: Stencils.h:221
GridT::TreeType TreeType
Definition: Stencils.h:39
std::vector< ValueType > BufferType
Definition: Stencils.h:42
GridT GridType
Definition: Stencils.h:38
int size()
Return the size of the stencil buffer.
Definition: Stencils.h:119
Definition: Stencils.h:302
math::Vec3< ValueType > gradient(const math::Vec3< ValueType > &xyz) const
Return the gradient in world space of the trilinear interpolation kernel.
Definition: Stencils.h:373
BoxStencil(const GridType &grid)
Definition: Stencils.h:312
GridT::ValueType ValueType
Definition: Stencils.h:308
bool intersects(const ValueType &isoValue=zeroVal< ValueType >()) const
Return true if the center of the stencil intersects the.
Definition: Stencils.h:320
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:316
GridT::TreeType TreeType
Definition: Stencils.h:307
GridT GridType
Definition: Stencils.h:306
ValueType interpolation(const math::Vec3< ValueType > &xyz) const
Return the trilinear interpolation at the normalized position.
Definition: Stencils.h:339
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
Coord offsetBy(Int32 dx, Int32 dy, Int32 dz) const
Definition: Coord.h:91
Definition: Stencils.h:1517
void curvaturesNormGrad(ValueType &mean, ValueType &gauss) const
Return both the mean and the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:1609
CurvatureStencil(const GridType &grid, Real dx)
Definition: Stencils.h:1534
CurvatureStencil(const GridType &grid)
Definition: Stencils.h:1527
ValueType meanCurvature() const
Return the mean curvature at the previously buffered location.
Definition: Stencils.h:1545
ValueType gaussianCurvatureNormGrad() const
Definition: Stencils.h:1597
void curvatures(ValueType &mean, ValueType &gauss) const
Return both the mean and the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:1568
math::Vec3< ValueType > gradient() const
Definition: Stencils.h:1656
GridT::ValueType ValueType
Definition: Stencils.h:1523
ValueType laplacian() const
Definition: Stencils.h:1643
std::pair< ValueType, ValueType > principalCurvatures() const
Return the pair (minimum, maximum) principal curvature at the previously buffered location.
Definition: Stencils.h:1625
ValueType meanCurvatureNormGrad() const
Definition: Stencils.h:1585
ValueType gaussianCurvature() const
Return the Gaussian curvature at the previously buffered location.
Definition: Stencils.h:1556
GridT::TreeType TreeType
Definition: Stencils.h:1522
GridT GridType
Definition: Stencils.h:1521
Dense stencil of a given width.
Definition: Stencils.h:1766
void moveTo(const IterType &iter)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:1793
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (x, y, z) and its neighbors.
Definition: Stencils.h:1785
const ValueType & getCenterValue() const
Definition: Stencils.h:1781
DenseStencil(const GridType &grid, int halfWidth)
Definition: Stencils.h:1774
GridType::ValueType ValueType
Definition: Stencils.h:1772
GridT::TreeType TreeType
Definition: Stencils.h:1771
GridT GridType
Definition: Stencils.h:1770
FourthOrderDenseStencil(const GridType &grid)
Definition: Stencils.h:694
GridType::ValueType ValueType
Definition: Stencils.h:690
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:698
GridT::TreeType TreeType
Definition: Stencils.h:689
GridT GridType
Definition: Stencils.h:688
Definition: Stencils.h:1233
GradStencil(const GridType &grid)
Definition: Stencils.h:1243
math::Vec3< ValueType > gradient(const math::Vec3< ValueType > &V) const
Return the first-order upwind gradient corresponding to the direction V.
Definition: Stencils.h:1288
GradStencil(const GridType &grid, Real dx)
Definition: Stencils.h:1250
math::Vec3< ValueType > gradient() const
Return the gradient computed at the previously buffered location by second order central differencing...
Definition: Stencils.h:1278
ValueType laplacian() const
Definition: Stencils.h:1298
bool zeroCrossing() const
Definition: Stencils.h:1307
math::Vec3< ValueType > cpt()
Compute the closest-point transform to a level set.
Definition: Stencils.h:1321
GridType::ValueType ValueType
Definition: Stencils.h:1239
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:1335
ValueType normSqGrad() const
Return the norm square of the single-sided upwind gradient (computed via Godunov's scheme) at the pre...
Definition: Stencils.h:1262
GridT::TreeType TreeType
Definition: Stencils.h:1238
GridT GridType
Definition: Stencils.h:1237
GridType::ValueType ValueType
Definition: Stencils.h:829
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:837
GridT::TreeType TreeType
Definition: Stencils.h:828
GridT GridType
Definition: Stencils.h:827
NineteenPointStencil(const GridType &grid)
Definition: Stencils.h:833
SecondOrderDenseStencil(const GridType &grid)
Definition: Stencils.h:483
GridType::ValueType ValueType
Definition: Stencils.h:479
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:487
GridT::TreeType TreeType
Definition: Stencils.h:478
GridT GridType
Definition: Stencils.h:477
Definition: Stencils.h:248
SevenPointStencil(const GridT &grid)
Definition: Stencils.h:258
GridT::ValueType ValueType
Definition: Stencils.h:254
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:262
GridT::TreeType TreeType
Definition: Stencils.h:253
GridT GridType
Definition: Stencils.h:252
SixthOrderDenseStencil(const GridType &grid)
Definition: Stencils.h:1049
GridType::ValueType ValueType
Definition: Stencils.h:1045
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:1053
GridT::TreeType TreeType
Definition: Stencils.h:1044
GridT GridType
Definition: Stencils.h:1043
ThirteenPointStencil(const GridType &grid)
Definition: Stencils.h:563
GridType::ValueType ValueType
Definition: Stencils.h:559
unsigned int pos() const
Return linear offset for the specified stencil point relative to its center.
Definition: Stencils.h:567
GridT::TreeType TreeType
Definition: Stencils.h:558
GridT GridType
Definition: Stencils.h:557
Definition: Vec3.h:24
This is a special 19-point stencil that supports optimal fifth-order WENO upwinding,...
Definition: Stencils.h:1367
WenoStencil(const GridType &grid, Real dx)
Definition: Stencils.h:1386
math::Vec3< ValueType > gradient(const math::Vec3< ValueType > &V) const
Definition: Stencils.h:1434
WenoStencil(const GridType &grid)
Definition: Stencils.h:1377
ValueType normSqGrad(const ValueType &isoValue=zeroVal< ValueType >()) const
Return the norm-square of the WENO upwind gradient (computed via WENO upwinding and Godunov's scheme)...
Definition: Stencils.h:1400
math::Vec3< ValueType > gradient() const
Definition: Stencils.h:1450
ValueType laplacian() const
Definition: Stencils.h:1462
bool zeroCrossing() const
Definition: Stencils.h:1472
GridType::ValueType ValueType
Definition: Stencils.h:1373
GridT::TreeType TreeType
Definition: Stencils.h:1372
GridT GridType
Definition: Stencils.h:1371
Definition: ValueAccessor.h:191
Type Pow3(Type x)
Return x3.
Definition: Math.h:552
Type Pow2(Type x)
Return x2.
Definition: Math.h:548
GridType::Ptr meanCurvature(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the mean curvature of the given grid.
Definition: GridOperators.h:1035
double Real
Definition: Types.h:60
Definition: Exceptions.h:13
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:212