13 #ifndef OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 14 #define OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 16 #include <tbb/parallel_for.h> 69 template<
typename VelocityGridT =
Vec3fGrid,
70 bool StaggeredVelocity =
false,
83 VolumeAdvection(
const VelocityGridT& velGrid, InterrupterType* interrupter =
nullptr)
85 , mInterrupter(interrupter)
86 , mIntegrator( Scheme::
SEMI )
87 , mLimiter( Scheme::
CLAMP )
92 e.
add(velGrid.background().length());
93 mMaxVelocity = e.
max();
114 switch (mIntegrator) {
177 template<
typename VolumeGr
idT>
180 if (!inGrid.hasUniformVoxels()) {
183 const double d = mMaxVelocity*
math::Abs(dt)/inGrid.voxelSize()[0];
205 template<
typename VolumeGridT,
206 typename VolumeSamplerT>
207 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
double timeStep)
209 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
210 const double dt = timeStep/mSubSteps;
211 const int n = this->getMaxDistance(inGrid, dt);
213 this->
template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
214 for (
int step = 1; step < mSubSteps; ++step) {
215 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
217 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
218 outGrid.swap( tmpGrid );
251 template<
typename VolumeGridT,
253 typename VolumeSamplerT>
254 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
const MaskGridT& mask,
double timeStep)
256 if (inGrid.transform() != mask.transform()) {
258 "resampling either of the two grids into the index space of the other.");
260 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
261 const double dt = timeStep/mSubSteps;
262 const int n = this->getMaxDistance(inGrid, dt);
264 outGrid->topologyIntersection( mask );
266 this->
template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
267 outGrid->topologyUnion( inGrid );
269 for (
int step = 1; step < mSubSteps; ++step) {
270 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
272 tmpGrid->topologyIntersection( mask );
274 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
275 tmpGrid->topologyUnion( inGrid );
276 outGrid.swap( tmpGrid );
286 void start(
const char* str)
const 288 if (mInterrupter) mInterrupter->start(str);
292 if (mInterrupter) mInterrupter->end();
294 bool interrupt()
const 297 tbb::task::self().cancel_group_execution();
303 template<
typename VolumeGr
idT,
typename VolumeSamplerT>
304 void cook(VolumeGridT& outGrid,
const VolumeGridT& inGrid,
double dt)
306 switch (mIntegrator) {
308 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
309 adv.cook(outGrid, dt);
313 Advect<VolumeGridT, 2, VolumeSamplerT> adv(inGrid, *
this);
314 adv.cook(outGrid, dt);
318 Advect<VolumeGridT, 3, VolumeSamplerT> adv(inGrid, *
this);
319 adv.cook(outGrid, dt);
323 Advect<VolumeGridT, 4, VolumeSamplerT> adv(inGrid, *
this);
324 adv.cook(outGrid, dt);
328 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
329 adv.cook(outGrid, dt);
333 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
334 adv.cook(outGrid, dt);
338 OPENVDB_THROW(ValueError,
"Spatial difference scheme not supported!");
344 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
struct Advect;
347 const VelocityGridT& mVelGrid;
349 InterrupterType* mInterrupter;
357 template<
typename VelocityGr
idT,
bool StaggeredVelocity,
typename InterrupterType>
358 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
359 struct VolumeAdvection<VelocityGridT, StaggeredVelocity, InterrupterType>::Advect
361 using TreeT =
typename VolumeGridT::TreeType;
362 using AccT =
typename VolumeGridT::ConstAccessor;
363 using ValueT =
typename TreeT::ValueType;
364 using LeafManagerT =
typename tree::LeafManager<TreeT>;
365 using LeafNodeT =
typename LeafManagerT::LeafNodeType;
366 using LeafRangeT =
typename LeafManagerT::LeafRange;
367 using VelocityIntegratorT = VelocityIntegrator<VelocityGridT, StaggeredVelocity>;
368 using RealT =
typename VelocityIntegratorT::ElementType;
369 using VoxelIterT =
typename TreeT::LeafNodeType::ValueOnIter;
371 Advect(
const VolumeGridT& inGrid,
const VolumeAdvection& parent)
374 , mVelocityInt(parent.mVelGrid)
378 inline void cook(
const LeafRangeT& range)
380 if (mParent->mGrainSize > 0) {
381 tbb::parallel_for(range, *
this);
386 void operator()(
const LeafRangeT& range)
const 389 mTask(const_cast<Advect*>(
this), range);
391 void cook(VolumeGridT& outGrid,
double time_step)
393 namespace ph = std::placeholders;
395 mParent->start(
"Advecting volume");
396 LeafManagerT manager(outGrid.tree(), mParent->spatialOrder()==2 ? 1 : 0);
397 const LeafRangeT range = manager.leafRange(mParent->mGrainSize);
398 const RealT dt = static_cast<RealT>(-time_step);
400 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
402 mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);
404 mTask = std::bind(&Advect::mac, ph::_1, ph::_2);
407 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
409 mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);
411 mTask = std::bind(&Advect::bfecc, ph::_1, ph::_2);
413 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 1, &outGrid);
415 manager.swapLeafBuffer(1);
417 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
421 if (mParent->spatialOrder()==2) manager.removeAuxBuffers();
423 mTask = std::bind(&Advect::limiter, ph::_1, ph::_2, dt);
429 void mac(
const LeafRangeT& range)
const 431 if (mParent->interrupt())
return;
433 AccT acc = mInGrid->getAccessor();
434 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
435 ValueT* out0 = leafIter.buffer( 0 ).data();
436 const ValueT* out1 = leafIter.buffer( 1 ).data();
437 const LeafNodeT* leaf = acc.probeConstLeaf( leafIter->origin() );
438 if (leaf !=
nullptr) {
439 const ValueT* in0 = leaf->buffer().data();
440 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
441 const Index i = voxelIter.pos();
442 out0[i] += RealT(0.5) * ( in0[i] - out1[i] );
445 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
446 const Index i = voxelIter.pos();
447 out0[i] += RealT(0.5) * ( acc.getValue(voxelIter.getCoord()) - out1[i] );
453 void bfecc(
const LeafRangeT& range)
const 455 if (mParent->interrupt())
return;
457 AccT acc = mInGrid->getAccessor();
458 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
459 ValueT* out0 = leafIter.buffer( 0 ).data();
460 const ValueT* out1 = leafIter.buffer( 1 ).data();
461 const LeafNodeT* leaf = acc.probeConstLeaf(leafIter->origin());
462 if (leaf !=
nullptr) {
463 const ValueT* in0 = leaf->buffer().data();
464 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
465 const Index i = voxelIter.pos();
466 out0[i] = RealT(0.5)*( RealT(3)*in0[i] - out1[i] );
469 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
470 const Index i = voxelIter.pos();
471 out0[i] = RealT(0.5)*( RealT(3)*acc.getValue(voxelIter.getCoord()) - out1[i] );
477 void rk(
const LeafRangeT& range, RealT dt,
size_t n,
const VolumeGridT* grid)
const 479 if (mParent->interrupt())
return;
480 const math::Transform& xform = mInGrid->transform();
481 AccT acc = grid->getAccessor();
482 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
483 ValueT* phi = leafIter.buffer( n ).data();
484 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
485 ValueT& value = phi[voxelIter.pos()];
486 Vec3d wPos = xform.indexToWorld(voxelIter.getCoord());
487 mVelocityInt.template rungeKutta<OrderRK, Vec3d>(dt, wPos);
488 value = SamplerT::sample(acc, xform.worldToIndex(wPos));
492 void limiter(
const LeafRangeT& range, RealT dt)
const 494 if (mParent->interrupt())
return;
495 const bool doLimiter = mParent->isLimiterOn();
497 ValueT data[2][2][2], vMin, vMax;
498 const math::Transform& xform = mInGrid->transform();
499 AccT acc = mInGrid->getAccessor();
500 const ValueT backg = mInGrid->background();
501 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
502 ValueT* phi = leafIter.buffer( 0 ).data();
503 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
504 ValueT& value = phi[voxelIter.pos()];
507 assert(OrderRK == 1);
508 Vec3d wPos = xform.indexToWorld(voxelIter.getCoord());
509 mVelocityInt.template rungeKutta<1, Vec3d>(dt, wPos);
510 Vec3d iPos = xform.worldToIndex(wPos);
512 BoxSampler::getValues(data, acc, ijk);
516 }
else if (value < vMin || value > vMax ) {
517 iPos -=
Vec3R(ijk[0], ijk[1], ijk[2]);
518 value = BoxSampler::trilinearInterpolation( data, iPos );
524 leafIter->setValueOff( voxelIter.pos() );
531 typename std::function<void (Advect*,
const LeafRangeT&)> mTask;
532 const VolumeGridT* mInGrid;
533 const VelocityIntegratorT mVelocityInt;
534 const VolumeAdvection* mParent;
541 #endif // OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED
double max() const
Return the maximum value.
Definition: Stats.h:128
bool isApproxEqual(const Type &a, const Type &b)
Return true if a is equal to b to within the default floating-point comparison tolerance.
Definition: Math.h:351
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Defines two simple wrapper classes for advection velocity fields as well as VelocitySampler and Veloc...
Vec3< double > Vec3d
Definition: Vec3.h:662
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:82
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
Functions to efficiently compute histograms, extremas (min/max) and statistics (mean,...
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:542
Defined various multi-threaded utility functions for trees.
Vec3SGrid Vec3fGrid
Definition: openvdb.h:56
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:13
void add(double val)
Add a single sample.
Definition: Stats.h:106
math::Vec3< Real > Vec3R
Definition: Types.h:49
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:25
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
Definition: Exceptions.h:63
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:203
This class computes the minimum and maximum values of a population of floating-point values.
Definition: Stats.h:92
float RoundUp(float x)
Return x rounded up to the nearest integer.
Definition: Math.h:734
Index32 Index
Definition: Types.h:31
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
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154