13 #ifndef OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 14 #define OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 20 #include <openvdb/thread/Threading.h> 27 #include <tbb/parallel_for.h> 73 template<
typename VelocityGridT =
Vec3fGrid,
74 bool StaggeredVelocity =
false,
87 VolumeAdvection(
const VelocityGridT& velGrid, InterrupterType* interrupter =
nullptr)
89 , mInterrupter(interrupter)
90 , mIntegrator( Scheme::
SEMI )
91 , mLimiter( Scheme::
CLAMP )
96 e.
add(velGrid.background().length());
97 mMaxVelocity = e.
max();
118 switch (mIntegrator) {
181 template<
typename VolumeGr
idT>
184 if (!inGrid.hasUniformVoxels()) {
187 const double d = mMaxVelocity*
math::Abs(dt)/inGrid.voxelSize()[0];
209 template<
typename VolumeGridT,
210 typename VolumeSamplerT>
211 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
double timeStep)
213 typename VolumeGridT::Ptr
outGrid = inGrid.deepCopy();
214 const double dt = timeStep/mSubSteps;
215 const int n = this->getMaxDistance(inGrid, dt);
217 this->
template cook<VolumeGridT, VolumeSamplerT>(*
outGrid, inGrid, dt);
218 for (
int step = 1; step < mSubSteps; ++step) {
219 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
221 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *
outGrid, dt);
222 outGrid.swap( tmpGrid );
255 template<
typename VolumeGridT,
257 typename VolumeSamplerT>
258 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
const MaskGridT& mask,
double timeStep)
260 if (inGrid.transform() != mask.transform()) {
262 "resampling either of the two grids into the index space of the other.");
264 typename VolumeGridT::Ptr
outGrid = inGrid.deepCopy();
265 const double dt = timeStep/mSubSteps;
266 const int n = this->getMaxDistance(inGrid, dt);
268 outGrid->topologyIntersection( mask );
270 this->
template cook<VolumeGridT, VolumeSamplerT>(*
outGrid, inGrid, dt);
271 outGrid->topologyUnion( inGrid );
273 for (
int step = 1; step < mSubSteps; ++step) {
274 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
276 tmpGrid->topologyIntersection( mask );
278 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *
outGrid, dt);
279 tmpGrid->topologyUnion( inGrid );
280 outGrid.swap( tmpGrid );
290 void start(
const char* str)
const 292 if (mInterrupter) mInterrupter->start(str);
296 if (mInterrupter) mInterrupter->end();
298 bool interrupt()
const 301 thread::cancelGroupExecution();
307 template<
typename VolumeGr
idT,
typename VolumeSamplerT>
308 void cook(VolumeGridT&
outGrid,
const VolumeGridT& inGrid,
double dt)
310 switch (mIntegrator) {
312 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
313 adv.cook(outGrid, dt);
317 Advect<VolumeGridT, 2, VolumeSamplerT> adv(inGrid, *
this);
318 adv.cook(outGrid, dt);
322 Advect<VolumeGridT, 3, VolumeSamplerT> adv(inGrid, *
this);
323 adv.cook(outGrid, dt);
327 Advect<VolumeGridT, 4, VolumeSamplerT> adv(inGrid, *
this);
328 adv.cook(outGrid, dt);
332 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
333 adv.cook(outGrid, dt);
337 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
338 adv.cook(outGrid, dt);
348 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
struct Advect;
351 const VelocityGridT& mVelGrid;
353 InterrupterType* mInterrupter;
361 template<
typename VelocityGr
idT,
bool StaggeredVelocity,
typename InterrupterType>
362 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
363 struct VolumeAdvection<VelocityGridT, StaggeredVelocity, InterrupterType>::Advect
365 using TreeT =
typename VolumeGridT::TreeType;
366 using AccT =
typename VolumeGridT::ConstAccessor;
367 using ValueT =
typename TreeT::ValueType;
369 using LeafNodeT =
typename LeafManagerT::LeafNodeType;
370 using LeafRangeT =
typename LeafManagerT::LeafRange;
372 using RealT =
typename VelocityIntegratorT::ElementType;
373 using VoxelIterT =
typename TreeT::LeafNodeType::ValueOnIter;
378 , mVelocityInt(parent.mVelGrid)
382 inline void cook(
const LeafRangeT& range)
384 if (mParent->mGrainSize > 0) {
385 tbb::parallel_for(range, *
this);
390 void operator()(
const LeafRangeT& range)
const 393 mTask(const_cast<Advect*>(
this), range);
395 void cook(VolumeGridT&
outGrid,
double time_step)
397 namespace ph = std::placeholders;
399 mParent->start(
"Advecting volume");
400 LeafManagerT manager(outGrid.tree(), mParent->spatialOrder()==2 ? 1 : 0);
401 const LeafRangeT range = manager.leafRange(mParent->mGrainSize);
402 const RealT dt =
static_cast<RealT
>(-time_step);
404 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
406 mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);
408 mTask = std::bind(&Advect::mac, ph::_1, ph::_2);
411 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
413 mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);
415 mTask = std::bind(&Advect::bfecc, ph::_1, ph::_2);
417 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 1, &outGrid);
419 manager.swapLeafBuffer(1);
421 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);
425 if (mParent->spatialOrder()==2) manager.removeAuxBuffers();
427 mTask = std::bind(&Advect::limiter, ph::_1, ph::_2, dt);
433 void mac(
const LeafRangeT& range)
const 435 if (mParent->interrupt())
return;
437 AccT acc = mInGrid->getAccessor();
438 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
439 ValueT* out0 = leafIter.buffer( 0 ).data();
440 const ValueT* out1 = leafIter.buffer( 1 ).data();
441 const LeafNodeT* leaf = acc.probeConstLeaf( leafIter->origin() );
442 if (leaf !=
nullptr) {
443 const ValueT* in0 = leaf->buffer().data();
444 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
445 const Index i = voxelIter.pos();
446 out0[i] += RealT(0.5) * ( in0[i] - out1[i] );
449 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
450 const Index i = voxelIter.pos();
451 out0[i] += RealT(0.5) * ( acc.getValue(voxelIter.getCoord()) - out1[i] );
457 void bfecc(
const LeafRangeT& range)
const 459 if (mParent->interrupt())
return;
461 AccT acc = mInGrid->getAccessor();
462 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
463 ValueT* out0 = leafIter.buffer( 0 ).data();
464 const ValueT* out1 = leafIter.buffer( 1 ).data();
465 const LeafNodeT* leaf = acc.probeConstLeaf(leafIter->origin());
466 if (leaf !=
nullptr) {
467 const ValueT* in0 = leaf->buffer().data();
468 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
469 const Index i = voxelIter.pos();
470 out0[i] = RealT(0.5)*( RealT(3)*in0[i] - out1[i] );
473 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
474 const Index i = voxelIter.pos();
475 out0[i] = RealT(0.5)*( RealT(3)*acc.getValue(voxelIter.getCoord()) - out1[i] );
481 void rk(
const LeafRangeT& range, RealT dt,
size_t n,
const VolumeGridT* grid)
const 483 if (mParent->interrupt())
return;
485 AccT acc = grid->getAccessor();
486 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
487 ValueT* phi = leafIter.buffer( n ).data();
488 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
489 ValueT& value = phi[voxelIter.pos()];
491 mVelocityInt.template rungeKutta<OrderRK, Vec3d>(dt, wPos);
492 value = SamplerT::sample(acc, xform.
worldToIndex(wPos));
496 void limiter(
const LeafRangeT& range, RealT dt)
const 498 if (mParent->interrupt())
return;
499 const bool doLimiter = mParent->isLimiterOn();
501 ValueT data[2][2][2], vMin, vMax;
503 AccT acc = mInGrid->getAccessor();
504 const ValueT backg = mInGrid->background();
505 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
506 ValueT* phi = leafIter.buffer( 0 ).data();
507 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
508 ValueT& value = phi[voxelIter.pos()];
513 mVelocityInt.template rungeKutta<1, Vec3d>(dt, wPos);
515 Coord ijk = Coord::floor( iPos );
516 BoxSampler::getValues(data, acc, ijk);
520 }
else if (value < vMin || value > vMax ) {
521 iPos -=
Vec3R(ijk[0], ijk[1], ijk[2]);
522 value = BoxSampler::trilinearInterpolation( data, iPos );
528 leafIter->setValueOff( voxelIter.pos() );
535 typename std::function<void (Advect*, const LeafRangeT&)> mTask;
536 const VolumeGridT* mInGrid;
537 const VelocityIntegratorT mVelocityInt;
547 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION 549 #ifdef OPENVDB_INSTANTIATE_VOLUMEADVECT 564 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION 571 #endif // OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Coord Abs(const Coord &xyz)
Definition: Coord.h:518
void add(double val)
Add a single sample.
Definition: Stats.h:103
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...
#define OPENVDB_INSTANTIATE
Definition: version.h.in:157
Functions to efficiently compute histograms, extrema (min/max) and statistics (mean, variance, etc.) of grid values.
bool isApproxEqual(const Type &a, const Type &b, const Type &tolerance)
Return true if a is equal to b to within the given tolerance.
Definition: Math.h:406
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:595
Index32 Index
Definition: Types.h:54
Vec3< double > Vec3d
Definition: Vec3.h:665
Base class for interrupters.
Definition: NullInterrupter.h:25
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
Defined various multi-threaded utility functions for trees.
Definition: Exceptions.h:65
#define OPENVDB_ASSERT(X)
Definition: Assert.h:41
Grid< DoubleTree > DoubleGrid
Definition: openvdb.h:74
double max() const
Return the maximum value.
Definition: Stats.h:125
Implementation of morphological dilation and erosion.
SharedPtr< Grid > Ptr
Definition: Grid.h:573
Definition: Exceptions.h:13
#define OPENVDB_INSTANTIATE_CLASS
Definition: version.h.in:158
This class computes the minimum and maximum values of a population of floating-point values...
Definition: Stats.h:89
Vec3SGrid Vec3fGrid
Definition: openvdb.h:85
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
float RoundUp(float x)
Return x rounded up to the nearest integer.
Definition: Math.h:787
Grid< Vec3STree > Vec3SGrid
Definition: openvdb.h:81
Delta for small floating-point offsets.
Definition: Math.h:155
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:261
math::Vec3< Real > Vec3R
Definition: Types.h:72
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218
Grid< FloatTree > FloatGrid
Definition: openvdb.h:75
Definition: Exceptions.h:63