36 #ifndef OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
37 #define OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
39 #include <openvdb/openvdb.h>
40 #include <openvdb/math/Math.h>
41 #include <openvdb/Types.h>
42 #include <openvdb/Grid.h>
43 #include <openvdb/util/NullInterrupter.h>
46 #include <boost/static_assert.hpp>
47 #include <tbb/blocked_range.h>
48 #include <tbb/parallel_for.h>
60 template<
typename CptGr
idT = Vec3fGr
id>
74 mCptAccessor(cptGrid.getAccessor()),
79 mCptGrid(other.mCptGrid),
80 mCptAccessor(mCptGrid->getAccessor()),
81 mCptIterations(other.mCptIterations)
88 template <
typename LocationType>
96 for (
unsigned int i = 0; i < mCptIterations; ++i) {
97 const Vec3R location = mCptGrid->worldToIndex(
Vec3R(result[0], result[1], result[2]));
98 BoxSampler::sample<CptAccessor>(mCptAccessor, location, result);
106 const CptGridType* mCptGrid;
107 CptAccessor mCptAccessor;
108 unsigned int mCptIterations;
115 template<
typename Gr
idT = Vec3fGr
id,
bool StaggeredVelocity = false>
124 mVelAccessor(mVelGrid->getAccessor())
128 mVelGrid(other.mVelGrid),
129 mVelAccessor(mVelGrid->getAccessor())
137 template <
typename LocationType>
140 const Vec3R location = mVelGrid->worldToIndex(
Vec3R(W[0], W[1], W[2]));
142 if (StaggeredVelocity) {
144 StaggeredBoxSampler::sample<VelAccessor>(mVelAccessor, location, result);
147 BoxSampler::sample<VelAccessor>(mVelAccessor, location, result);
153 const GridT* mVelGrid;
154 VelAccessor mVelAccessor;
160 template<
typename Gr
idT = Vec3fGr
id,
bool StaggeredVelocity = false>
172 template<
int Order,
typename LocationType>
174 VecType P(loc[0],loc[1],loc[2]), V0, V1, V2, V3;
176 BOOST_STATIC_ASSERT((Order < 5) && (Order > -1));
181 }
else if (Order == 1) {
182 mVelField.sample(P, V0);
185 }
else if (Order == 2) {
186 mVelField.sample(P, V0);
190 }
else if (Order == 3) {
191 mVelField.sample(P, V0);
193 mVelField.sample(P+dt*(
ElementType(2.0)*V1-V0), V2);
196 }
else if (Order == 4) {
197 mVelField.sample(P, V0);
200 mVelField.sample(P+ dt*V2, V3);
204 loc += LocationType(P[0], P[1], P[2]);
237 typename PointListT = std::vector<typename GridT::ValueType>,
238 bool StaggeredVelocity =
false,
248 PointAdvect(
const GridT& velGrid, InterrupterType* interrupter=NULL) :
251 mIntegrationOrder(1),
253 mInterrupter(interrupter)
257 mVelGrid(other.mVelGrid),
258 mPoints(other.mPoints),
260 mAdvIterations(other.mAdvIterations),
261 mIntegrationOrder(other.mIntegrationOrder),
262 mThreaded(other.mThreaded),
263 mInterrupter(other.mInterrupter)
270 bool earlyOut()
const {
return (mIntegrationOrder==0);}
277 void advect(PointListT& points,
float dt,
unsigned int advIterations = 1)
279 if (this->earlyOut())
return;
282 mAdvIterations = advIterations;
284 if (mInterrupter) mInterrupter->start(
"Advecting points by OpenVDB velocity field: ");
286 tbb::parallel_for(tbb::blocked_range<size_t>(0, mPoints->size()), *
this);
288 (*this)(tbb::blocked_range<size_t>(0, mPoints->size()));
290 if (mInterrupter) mInterrupter->end();
294 void operator() (
const tbb::blocked_range<size_t> &range)
const
296 if (mInterrupter && mInterrupter->wasInterrupted()) {
297 tbb::task::self().cancel_group_execution();
301 switch (mIntegrationOrder) {
304 for (
size_t n = range.begin(); n != range.end(); ++n) {
307 for (
unsigned int i = 0; i < mAdvIterations; ++i) {
308 velField.template rungeKutta<1>(mDt, X0);
315 for (
size_t n = range.begin(); n != range.end(); ++n) {
318 for (
unsigned int i = 0; i < mAdvIterations; ++i) {
319 velField.template rungeKutta<2>(mDt, X0);
326 for (
size_t n = range.begin(); n != range.end(); ++n) {
329 for (
unsigned int i = 0; i < mAdvIterations; ++i) {
330 velField.template rungeKutta<3>(mDt, X0);
337 for (
size_t n = range.begin(); n != range.end(); ++n) {
340 for (
unsigned int i = 0; i < mAdvIterations; ++i) {
341 velField.template rungeKutta<4>(mDt, X0);
351 const GridType* mVelGrid;
358 unsigned int mAdvIterations;
359 unsigned int mIntegrationOrder;
363 InterrupterType* mInterrupter;
369 typename PointListT = std::vector<typename GridT::ValueType>,
370 bool StaggeredVelocity =
false,
371 typename CptGridType = GridT,
383 const GridType& cptGrid,
int cptn, InterrupterType* interrupter = NULL):
387 mInterrupter(interrupter)
391 mVelGrid(other.mVelGrid),
392 mCptGrid(other.mCptGrid),
393 mCptIter(other.mCptIter),
394 mPoints(other.mPoints),
396 mAdvIterations(other.mAdvIterations),
397 mIntegrationOrder(other.mIntegrationOrder),
398 mThreaded(other.mThreaded),
399 mInterrupter(other.mInterrupter)
411 void advect(PointListT& points,
float dt,
unsigned int advIterations = 1)
416 if (mIntegrationOrder==0 && mCptIter == 0) {
419 (mIntegrationOrder>0) ? mAdvIterations = advIterations : mAdvIterations = 1;
421 if (mInterrupter) mInterrupter->start(
"Advecting points by OpenVDB velocity field: ");
422 const int N = mPoints->size();
425 tbb::parallel_for(tbb::blocked_range<size_t>(0, N), *
this);
427 (*this)(tbb::blocked_range<size_t>(0, N));
429 if (mInterrupter) mInterrupter->end();
434 void operator() (
const tbb::blocked_range<size_t> &range)
const
436 if (mInterrupter && mInterrupter->wasInterrupted()) {
437 tbb::task::self().cancel_group_execution();
442 switch (mIntegrationOrder) {
445 for (
size_t n = range.begin(); n != range.end(); ++n) {
447 for (
unsigned int i = 0; i < mAdvIterations; ++i) {
455 for (
size_t n = range.begin(); n != range.end(); ++n) {
457 for (
unsigned int i = 0; i < mAdvIterations; ++i) {
458 velField.template rungeKutta<1>(mDt, X0);
466 for (
size_t n = range.begin(); n != range.end(); ++n) {
468 for (
unsigned int i = 0; i < mAdvIterations; ++i) {
469 velField.template rungeKutta<2>(mDt, X0);
478 for (
size_t n = range.begin(); n != range.end(); ++n) {
480 for (
unsigned int i = 0; i < mAdvIterations; ++i) {
481 velField.template rungeKutta<3>(mDt, X0);
489 for (
size_t n = range.begin(); n != range.end(); ++n) {
491 for (
unsigned int i = 0; i < mAdvIterations; ++i) {
492 velField.template rungeKutta<4>(mDt, X0);
502 const GridType* mVelGrid;
503 const GridType* mCptGrid;
509 unsigned int mAdvIterations;
510 unsigned int mIntegrationOrder;
513 InterrupterType* mInterrupter;
520 #endif // OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
math::Vec3< Real > Vec3R
Definition: Types.h:74
Vec3SGrid Vec3fGrid
Definition: openvdb.h:79
#define OPENVDB_VERSION_NAME
Definition: version.h:45
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67