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