40 #ifndef OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 41 #define OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 43 #include <tbb/parallel_for.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);
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;
392 using LeafNodeT =
typename LeafManagerT::LeafNodeType;
393 using LeafRangeT =
typename LeafManagerT::LeafRange;
395 using RealT =
typename VelocityIntegratorT::ElementType;
396 using VoxelIterT =
typename TreeT::LeafNodeType::ValueOnIter;
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;
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()];
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;
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);
536 mVelocityInt.template rungeKutta<1, Vec3d>(dt, 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;
568 #endif // OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
Delta for small floating-point offsets.
Definition: Math.h:124
double max() const
Return the maximum value.
Definition: Stats.h:154
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:545
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...
This class computes the minimum and maximum values of a population of floating-point values...
Definition: Stats.h:118
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:108
static Coord floor(const Vec3< T > &xyz)
Return the largest integer coordinates that are not greater than xyz (node centered conversion)...
Definition: Coord.h:83
Functions to efficiently compute histograms, extremas (min/max) and statistics (mean, variance, etc.) of grid values.
float RoundUp(float x)
Return x rounded up to the nearest integer.
Definition: Math.h:737
Vec3< double > Vec3d
Definition: Vec3.h:679
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:354
math::Vec3< Real > Vec3R
Definition: Types.h:78
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:51
Coord Abs(const Coord &xyz)
Definition: Coord.h:513
Defined various multi-threaded utility functions for trees.
Definition: Exceptions.h:89
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:230
Index32 Index
Definition: Types.h:60
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:136
Definition: Exceptions.h:91
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:39
void add(double val)
Add a single sample.
Definition: Stats.h:132
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:110
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
Vec3SGrid Vec3fGrid
Definition: openvdb.h:83
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:188