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 <boost/bind.hpp> 45 #include <boost/function.hpp> 97 template<
typename VelocityGridT =
Vec3fGrid,
98 bool StaggeredVelocity =
false,
113 , mInterrupter(interrupter)
114 , mIntegrator( Scheme::
SEMI )
115 , mLimiter( Scheme::
CLAMP )
120 e.
add(velGrid.background().length());
121 mMaxVelocity = e.
max();
142 switch (mIntegrator) {
205 template<
typename VolumeGr
idT>
208 if (!inGrid.hasUniformVoxels()) {
211 const double d = mMaxVelocity*
math::Abs(dt)/inGrid.voxelSize()[0];
233 template<
typename VolumeGridT,
234 typename VolumeSamplerT>
235 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
double timeStep)
237 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
238 const double dt = timeStep/mSubSteps;
239 const int n = this->getMaxDistance(inGrid, dt);
241 this->
template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
242 for (
int step = 1; step < mSubSteps; ++step) {
243 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
245 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
246 outGrid.swap( tmpGrid );
279 template<
typename VolumeGridT,
281 typename VolumeSamplerT>
282 typename VolumeGridT::Ptr
advect(
const VolumeGridT& inGrid,
const MaskGridT& mask,
double timeStep)
284 if (inGrid.transform() != mask.transform()) {
286 "resampling either of the two grids into the index space of the other.");
288 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy();
289 const double dt = timeStep/mSubSteps;
290 const int n = this->getMaxDistance(inGrid, dt);
292 outGrid->topologyIntersection( mask );
294 this->
template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt);
295 outGrid->topologyUnion( inGrid );
297 for (
int step = 1; step < mSubSteps; ++step) {
298 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy();
300 tmpGrid->topologyIntersection( mask );
302 this->
template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt);
303 tmpGrid->topologyUnion( inGrid );
304 outGrid.swap( tmpGrid );
314 void start(
const char* str)
const 316 if (mInterrupter) mInterrupter->start(str);
320 if (mInterrupter) mInterrupter->end();
322 bool interrupt()
const 325 tbb::task::self().cancel_group_execution();
331 template<
typename VolumeGr
idT,
typename VolumeSamplerT>
332 void cook(VolumeGridT& outGrid,
const VolumeGridT& inGrid,
double dt)
334 switch (mIntegrator) {
336 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
337 adv.cook(outGrid, dt);
341 Advect<VolumeGridT, 2, VolumeSamplerT> adv(inGrid, *
this);
342 adv.cook(outGrid, dt);
346 Advect<VolumeGridT, 3, VolumeSamplerT> adv(inGrid, *
this);
347 adv.cook(outGrid, dt);
351 Advect<VolumeGridT, 4, VolumeSamplerT> adv(inGrid, *
this);
352 adv.cook(outGrid, dt);
356 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
357 adv.cook(outGrid, dt);
361 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *
this);
362 adv.cook(outGrid, dt);
372 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
struct Advect;
375 const VelocityGridT& mVelGrid;
377 InterrupterType* mInterrupter;
385 template<
typename VelocityGr
idT,
bool StaggeredVelocity,
typename InterrupterType>
386 template<
typename VolumeGr
idT,
size_t OrderRK,
typename SamplerT>
387 struct VolumeAdvection<VelocityGridT, StaggeredVelocity, InterrupterType>::Advect
389 typedef typename VolumeGridT::TreeType TreeT;
390 typedef typename VolumeGridT::ConstAccessor AccT;
391 typedef typename TreeT::ValueType ValueT;
393 typedef typename LeafManagerT::LeafNodeType LeafNodeT;
394 typedef typename LeafManagerT::LeafRange LeafRangeT;
396 typedef typename VelocityIntegratorT::ElementType RealT;
397 typedef typename TreeT::LeafNodeType::ValueOnIter VoxelIterT;
402 , mVelocityInt(parent.mVelGrid)
406 inline void cook(
const LeafRangeT& range)
408 if (mParent->mGrainSize > 0) {
409 tbb::parallel_for(range, *
this);
414 void operator()(
const LeafRangeT& range)
const 417 mTask(const_cast<Advect*>(
this), range);
419 void cook(VolumeGridT& outGrid,
double time_step)
421 mParent->start(
"Advecting volume");
422 LeafManagerT manager(outGrid.tree(), mParent->spatialOrder()==2 ? 1 : 0);
423 const LeafRangeT range = manager.leafRange(mParent->mGrainSize);
424 const RealT dt =
static_cast<RealT
>(-time_step);
426 mTask = boost::bind(&Advect::rk, _1, _2, dt, 0, mInGrid);
428 mTask = boost::bind(&Advect::rk, _1, _2,-dt, 1, &outGrid);
430 mTask = boost::bind(&Advect::mac, _1, _2);
433 mTask = boost::bind(&Advect::rk, _1, _2, dt, 0, mInGrid);
435 mTask = boost::bind(&Advect::rk, _1, _2,-dt, 1, &outGrid);
437 mTask = boost::bind(&Advect::bfecc, _1, _2);
439 mTask = boost::bind(&Advect::rk, _1, _2, dt, 1, &outGrid);
441 manager.swapLeafBuffer(1);
443 mTask = boost::bind(&Advect::rk, _1, _2, dt, 0, mInGrid);
447 if (mParent->spatialOrder()==2) manager.removeAuxBuffers();
449 mTask = boost::bind(&Advect::limiter, _1, _2, dt);
455 void mac(
const LeafRangeT& range)
const 457 if (mParent->interrupt())
return;
459 AccT acc = mInGrid->getAccessor();
460 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
461 ValueT* out0 = leafIter.buffer( 0 ).data();
462 const ValueT* out1 = leafIter.buffer( 1 ).data();
463 const LeafNodeT* leaf = acc.probeConstLeaf( leafIter->origin() );
465 const ValueT* in0 = leaf->buffer().data();
466 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
467 const Index i = voxelIter.pos();
468 out0[i] += RealT(0.5) * ( in0[i] - out1[i] );
471 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
472 const Index i = voxelIter.pos();
473 out0[i] += RealT(0.5) * ( acc.getValue(voxelIter.getCoord()) - out1[i] );
479 void bfecc(
const LeafRangeT& range)
const 481 if (mParent->interrupt())
return;
483 AccT acc = mInGrid->getAccessor();
484 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
485 ValueT* out0 = leafIter.buffer( 0 ).data();
486 const ValueT* out1 = leafIter.buffer( 1 ).data();
487 const LeafNodeT* leaf = acc.probeConstLeaf(leafIter->origin());
489 const ValueT* in0 = leaf->buffer().data();
490 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
491 const Index i = voxelIter.pos();
492 out0[i] = RealT(0.5)*( RealT(3)*in0[i] - out1[i] );
495 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
496 const Index i = voxelIter.pos();
497 out0[i] = RealT(0.5)*( RealT(3)*acc.getValue(voxelIter.getCoord()) - out1[i] );
503 void rk(
const LeafRangeT& range, RealT dt,
size_t n,
const VolumeGridT* grid)
const 505 if (mParent->interrupt())
return;
507 AccT acc = grid->getAccessor();
508 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
509 ValueT* phi = leafIter.buffer( n ).data();
510 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
511 ValueT& value = phi[voxelIter.pos()];
513 mVelocityInt.template rungeKutta<OrderRK, Vec3d>(dt, wPos);
514 value = SamplerT::sample(acc, xform.
worldToIndex(wPos));
518 void limiter(
const LeafRangeT& range, RealT dt)
const 520 if (mParent->interrupt())
return;
521 const bool doLimiter = mParent->isLimiterOn();
523 ValueT data[2][2][2], vMin, vMax;
525 AccT acc = mInGrid->getAccessor();
526 const ValueT backg = mInGrid->background();
527 for (
typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) {
528 ValueT* phi = leafIter.buffer( 0 ).data();
529 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) {
530 ValueT& value = phi[voxelIter.pos()];
533 assert(OrderRK == 1);
535 mVelocityInt.template rungeKutta<1, Vec3d>(dt, wPos);
538 BoxSampler::getValues(data, acc, ijk);
542 }
else if (value < vMin || value > vMax ) {
543 iPos -=
Vec3R(ijk[0], ijk[1], ijk[2]);
544 value = BoxSampler::trilinearInterpolation( data, iPos );
550 leafIter->setValueOff( voxelIter.pos() );
557 typename boost::function<void (Advect*, const LeafRangeT&)> mTask;
558 const VolumeGridT* mInGrid;
559 const VelocityIntegratorT mVelocityInt;
567 #endif // OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 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:370
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...
Coord Abs(const Coord &xyz)
Definition: Coord.h:254
static Coord floor(const Vec3< T > &xyz)
Return the largest integer coordinates that are not greater than xyz (node centered conversion)...
Definition: Coord.h:80
Delta for small floating-point offsets.
Definition: Math.h:132
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:101
Functions to efficiently compute histograms, extremas (min/max) and statistics (mean, variance, etc.) of grid values.
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:115
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:48
Defined various multi-threaded utility functions for trees.
float RoundUp(float x)
Return x rounded up to the nearest integer.
Definition: Math.h:753
Definition: Exceptions.h:90
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Vec3< double > Vec3d
Definition: Vec3.h:708
This class computes the minimum and maximum values of a population of floating-point values...
Definition: Stats.h:105
Vec3SGrid Vec3fGrid
Definition: openvdb.h:83
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:39
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:246
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
void add(double val)
Add a single sample.
Definition: Stats.h:119
Definition: Exceptions.h:92
math::Vec3< Real > Vec3R
Definition: Types.h:75
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:561
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Index32 Index
Definition: Types.h:57
double max() const
Return the maximum value.
Definition: Stats.h:141