40 #ifndef OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED
41 #define OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED
43 #include <tbb/parallel_for.h>
44 #include <openvdb/Types.h>
45 #include <openvdb/math/Math.h>
46 #include <openvdb/util/NullInterrupter.h>
96 template<
typename VelocityGridT =
Vec3fGrid,
97 bool StaggeredVelocity =
false,
112 , mInterrupter(interrupter)
113 , mIntegrator( Scheme::
SEMI )
114 , mLimiter( Scheme::
CLAMP )
119 e.
add(velGrid.background().length());
120 mMaxVelocity = e.
max();
141 switch (mIntegrator) {
204 template<
typename VolumeGr
idT>
207 if (!inGrid.hasUniformVoxels()) {
210 const double d = mMaxVelocity*
math::Abs(dt)/inGrid.voxelSize()[0];
232 template<
typename VolumeGridT,
233 typename VolumeSamplerT>
234 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
double timeStep)
236 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
237 const double dt = timeStep/mSubSteps;
238 const int n = this->getMaxDistance(inGrid, dt);
240 this->
template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
241 for (
int step = 1; step < mSubSteps; ++step) {
242 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
244 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
245 outGrid.swap( tmpGrid );
278 template<
typename VolumeGridT,
280 typename VolumeSamplerT>
281 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
const MaskGridT& mask,
double timeStep)
283 if (inGrid.transform() != mask.transform()) {
285 "resampling either of the two grids into the index space of the other.");
287 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
288 const double dt = timeStep/mSubSteps;
289 const int n = this->getMaxDistance(inGrid, dt);
291 outGrid->topologyIntersection( mask );
293 this->
template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
294 outGrid->topologyUnion( inGrid );
296 for (
int step = 1; step < mSubSteps; ++step) {
297 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
299 tmpGrid->topologyIntersection( mask );
301 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
302 tmpGrid->topologyUnion( inGrid );
303 outGrid.swap( tmpGrid );
313 void start(
const char* str)
const
315 if (mInterrupter) mInterrupter->start(str);
319 if (mInterrupter) mInterrupter->end();
321 bool interrupt()
const
324 tbb::task::self().cancel_group_execution();
330 template<
typename VolumeGr
idT,
typename VolumeSamplerT>
331 void cook(VolumeGridT& outGrid,
const VolumeGridT& inGrid,
double dt)
333 switch (mIntegrator) {
335 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
336 adv.cook(outGrid, dt);
340 Advect<VolumeGridT, 2, VolumeSamplerT> adv(inGrid, *
this);
341 adv.cook(outGrid, dt);
345 Advect<VolumeGridT, 3, VolumeSamplerT> adv(inGrid, *
this);
346 adv.cook(outGrid, dt);
350 Advect<VolumeGridT, 4, VolumeSamplerT> adv(inGrid, *
this);
351 adv.cook(outGrid, dt);
355 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
356 adv.cook(outGrid, dt);
360 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
361 adv.cook(outGrid, dt);
365 OPENVDB_THROW(ValueError,
"Spatial difference scheme not supported!");
371 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
struct Advect;
374 const VelocityGridT& mVelGrid;
376 InterrupterType* mInterrupter;
384 template<
typename VelocityGr
idT,
bool StaggeredVelocity,
typename InterrupterType>
385 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
386 struct VolumeAdvection<VelocityGridT, StaggeredVelocity, InterrupterType>::Advect
388 using TreeT =
typename VolumeGridT::TreeType;
389 using AccT =
typename VolumeGridT::ConstAccessor;
390 using ValueT =
typename TreeT::ValueType;
391 using LeafManagerT =
typename tree::LeafManager<TreeT>;
392 using LeafNodeT =
typename LeafManagerT::LeafNodeType;
393 using LeafRangeT =
typename LeafManagerT::LeafRange;
394 using VelocityIntegratorT = VelocityIntegrator<VelocityGridT, StaggeredVelocity>;
395 using RealT =
typename VelocityIntegratorT::ElementType;
396 using VoxelIterT =
typename TreeT::LeafNodeType::ValueOnIter;
398 Advect(
const VolumeGridT& inGrid,
const VolumeAdvection& parent)
401 , mVelocityInt(parent.mVelGrid)
405 inline void cook(
const LeafRangeT& range)
407 if (mParent->mGrainSize > 0) {
408 tbb::parallel_for(range, *
this);
413 void operator()(
const LeafRangeT& range)
const
416 mTask(const_cast<Advect*>(
this), range);
418 void cook(VolumeGridT& outGrid,
double time_step)
420 namespace ph = std::placeholders;
422 mParent->start(
"Advecting volume");
423 LeafManagerT manager(outGrid.tree(), mParent->spatialOrder()==2 ? 1 : 0);
424 const LeafRangeT range = manager.leafRange(mParent->mGrainSize);
425 const RealT dt = static_cast<RealT>(-time_step);
427 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
429 mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);
431 mTask = std::bind(&Advect::mac, ph::_1, ph::_2);
434 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
436 mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);
438 mTask = std::bind(&Advect::bfecc, ph::_1, ph::_2);
440 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 1, &outGrid);
442 manager.swapLeafBuffer(1);
444 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
448 if (mParent->spatialOrder()==2) manager.removeAuxBuffers();
450 mTask = std::bind(&Advect::limiter, ph::_1, ph::_2, dt);
456 void mac(
const LeafRangeT& range)
const
458 if (mParent->interrupt())
return;
460 AccT acc = mInGrid->getAccessor();
461 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
462 ValueT* out0 = leafIter.buffer( 0 ).data();
463 const ValueT* out1 = leafIter.buffer( 1 ).data();
464 const LeafNodeT* leaf = acc.probeConstLeaf( leafIter->origin() );
465 if (leaf !=
nullptr) {
466 const ValueT* in0 = leaf->buffer().data();
467 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
468 const Index i = voxelIter.pos();
469 out0[i] += RealT(0.5) * ( in0[i] - out1[i] );
472 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
473 const Index i = voxelIter.pos();
474 out0[i] += RealT(0.5) * ( acc.getValue(voxelIter.getCoord()) - out1[i] );
480 void bfecc(
const LeafRangeT& range)
const
482 if (mParent->interrupt())
return;
484 AccT acc = mInGrid->getAccessor();
485 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
486 ValueT* out0 = leafIter.buffer( 0 ).data();
487 const ValueT* out1 = leafIter.buffer( 1 ).data();
488 const LeafNodeT* leaf = acc.probeConstLeaf(leafIter->origin());
489 if (leaf !=
nullptr) {
490 const ValueT* in0 = leaf->buffer().data();
491 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
492 const Index i = voxelIter.pos();
493 out0[i] = RealT(0.5)*( RealT(3)*in0[i] - out1[i] );
496 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
497 const Index i = voxelIter.pos();
498 out0[i] = RealT(0.5)*( RealT(3)*acc.getValue(voxelIter.getCoord()) - out1[i] );
504 void rk(
const LeafRangeT& range, RealT dt,
size_t n,
const VolumeGridT* grid)
const
506 if (mParent->interrupt())
return;
507 const math::Transform& xform = mInGrid->transform();
508 AccT acc = grid->getAccessor();
509 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
510 ValueT* phi = leafIter.buffer( n ).data();
511 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
512 ValueT& value = phi[voxelIter.pos()];
513 Vec3d wPos = xform.indexToWorld(voxelIter.getCoord());
514 mVelocityInt.template rungeKutta<OrderRK, Vec3d>(dt, wPos);
515 value = SamplerT::sample(acc, xform.worldToIndex(wPos));
519 void limiter(
const LeafRangeT& range, RealT dt)
const
521 if (mParent->interrupt())
return;
522 const bool doLimiter = mParent->isLimiterOn();
524 ValueT data[2][2][2], vMin, vMax;
525 const math::Transform& xform = mInGrid->transform();
526 AccT acc = mInGrid->getAccessor();
527 const ValueT backg = mInGrid->background();
528 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
529 ValueT* phi = leafIter.buffer( 0 ).data();
530 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
531 ValueT& value = phi[voxelIter.pos()];
534 assert(OrderRK == 1);
535 Vec3d wPos = xform.indexToWorld(voxelIter.getCoord());
536 mVelocityInt.template rungeKutta<1, Vec3d>(dt, wPos);
537 Vec3d iPos = xform.worldToIndex(wPos);
539 BoxSampler::getValues(data, acc, ijk);
543 }
else if (value < vMin || value > vMax ) {
544 iPos -=
Vec3R(ijk[0], ijk[1], ijk[2]);
545 value = BoxSampler::trilinearInterpolation( data, iPos );
551 leafIter->setValueOff( voxelIter.pos() );
558 typename std::function<void (Advect*,
const LeafRangeT&)> mTask;
559 const VolumeGridT* mInGrid;
560 const VelocityIntegratorT mVelocityInt;
561 const VolumeAdvection* mParent;
568 #endif // OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED