1 // Copyright Contributors to the OpenVDB Project 2 // SPDX-License-Identifier: MPL-2.0 3 // 4 /////////////////////////////////////////////////////////////////////////// 5 // 6 /// @author Ken Museth 7 /// 8 /// @file tools/VolumeAdvect.h 9 /// 10 /// @brief Sparse hyperbolic advection of volumes, e.g. a density or 11 /// velocity (vs a level set interface). 12 13 #ifndef OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 14 #define OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 15 16 #include "openvdb/Types.h" 17 #include "openvdb/math/Math.h" 18 #include "openvdb/util/NullInterrupter.h" 19 #include "openvdb/thread/Threading.h" 20 #include "Interpolation.h"// for Sampler 21 #include "VelocityFields.h" // for VelocityIntegrator 22 #include "Morphology.h"//for dilateActiveValues 23 #include "Prune.h"// for prune 24 #include "Statistics.h" // for extrema 25 26 #include <tbb/parallel_for.h> 27 28 #include <functional> 29 30 31 namespace openvdb { 32 OPENVDB_USE_VERSION_NAMESPACE 33 namespace OPENVDB_VERSION_NAME { 34 namespace tools { 35 36 namespace Scheme { 37 /// @brief Numerical advections schemes. 38 enum SemiLagrangian { SEMI, MID, RK3, RK4, MAC, BFECC }; 39 /// @brief Flux-limiters employed to stabilize the second-order 40 /// advection schemes MacCormack and BFECC. 41 enum Limiter { NO_LIMITER, CLAMP, REVERT }; 42 } 43 44 /// @brief Performs advections of an arbitrary type of volume in a 45 /// static velocity field. The advections are performed by means 46 /// of various derivatives of Semi-Lagrangian integration, i.e. 47 /// backwards tracking along the hyperbolic characteristics 48 /// followed by interpolation. 49 /// 50 /// @note Optionally a limiter can be combined with the higher-order 51 /// integration schemes MacCormack and BFECC. There are two 52 /// types of limiters (CLAMP and REVERT) that suppress 53 /// non-physical oscillations by means of either claminging or 54 /// reverting to a first-order schemes when the function is not 55 /// bounded by the cell values used for tri-linear interpolation. 56 /// 57 /// @verbatim The supported integrations schemes: 58 /// 59 /// ================================================================ 60 /// | Lable | Accuracy | Integration Scheme | Interpolations | 61 /// | |Time/Space| | velocity/volume | 62 /// ================================================================ 63 /// | SEMI | 1/1 | Semi-Lagrangian | 1/1 | 64 /// | MID | 2/1 | Mid-Point | 2/1 | 65 /// | RK3 | 3/1 | 3rd Order Runge-Kutta | 3/1 | 66 /// | RK4 | 4/1 | 4th Order Runge-Kutta | 4/1 | 67 /// | MAC | 2/2 | MacCormack | 2/2 | 68 /// | BFECC | 2/2 | BFECC | 3/2 | 69 /// ================================================================ 70 /// @endverbatim 71 72 template<typename VelocityGridT = Vec3fGrid, 73 bool StaggeredVelocity = false, 74 typename InterrupterType = util::NullInterrupter> 75 class VolumeAdvection 76 { 77 public: 78 79 /// @brief Constructor 80 /// 81 /// @param velGrid Velocity grid responsible for the (passive) advection. 82 /// @param interrupter Optional interrupter used to prematurely end computations. 83 /// 84 /// @note The velocity field is assumed to be constant for the duration of the 85 /// advection. 86 VolumeAdvection(const VelocityGridT& velGrid, InterrupterType* interrupter = nullptr) mVelGrid(velGrid)87 : mVelGrid(velGrid) 88 , mInterrupter(interrupter) 89 , mIntegrator( Scheme::SEMI ) 90 , mLimiter( Scheme::CLAMP ) 91 , mGrainSize( 128 ) 92 , mSubSteps( 1 ) 93 { 94 math::Extrema e = extrema(velGrid.cbeginValueAll(), /*threading*/true); 95 e.add(velGrid.background().length()); 96 mMaxVelocity = e.max(); 97 } 98 ~VolumeAdvection()99 virtual ~VolumeAdvection() 100 { 101 } 102 103 /// @brief Return the spatial order of accuracy of the advection scheme 104 /// 105 /// @note This is the optimal order in smooth regions. In 106 /// non-smooth regions the flux-limiter will drop the order of 107 /// accuracy to add numerical dissipation. spatialOrder()108 int spatialOrder() const { return (mIntegrator == Scheme::MAC || 109 mIntegrator == Scheme::BFECC) ? 2 : 1; } 110 111 /// @brief Return the temporal order of accuracy of the advection scheme 112 /// 113 /// @note This is the optimal order in smooth regions. In 114 /// non-smooth regions the flux-limiter will drop the order of 115 /// accuracy to add numerical dissipation. temporalOrder()116 int temporalOrder() const { 117 switch (mIntegrator) { 118 case Scheme::SEMI: return 1; 119 case Scheme::MID: return 2; 120 case Scheme::RK3: return 3; 121 case Scheme::RK4: return 4; 122 case Scheme::BFECC:return 2; 123 case Scheme::MAC: return 2; 124 } 125 return 0;//should never reach this point 126 } 127 128 /// @brief Set the integrator (see details in the table above) setIntegrator(Scheme::SemiLagrangian integrator)129 void setIntegrator(Scheme::SemiLagrangian integrator) { mIntegrator = integrator; } 130 131 /// @brief Return the integrator (see details in the table above) getIntegrator()132 Scheme::SemiLagrangian getIntegrator() const { return mIntegrator; } 133 134 /// @brief Set the limiter (see details above) setLimiter(Scheme::Limiter limiter)135 void setLimiter(Scheme::Limiter limiter) { mLimiter = limiter; } 136 137 /// @brief Retrun the limiter (see details above) getLimiter()138 Scheme::Limiter getLimiter() const { return mLimiter; } 139 140 /// @brief Return @c true if a limiter will be applied based on 141 /// the current settings. isLimiterOn()142 bool isLimiterOn() const { return this->spatialOrder()>1 && 143 mLimiter != Scheme::NO_LIMITER; } 144 145 /// @return the grain-size used for multi-threading 146 /// @note A grainsize of 0 implies serial execution getGrainSize()147 size_t getGrainSize() const { return mGrainSize; } 148 149 /// @brief Set the grain-size used for multi-threading 150 /// @note A grainsize of 0 disables multi-threading 151 /// @warning A small grainsize can degrade performance, 152 /// both in terms of time and memory footprint! setGrainSize(size_t grainsize)153 void setGrainSize(size_t grainsize) { mGrainSize = grainsize; } 154 155 /// @return the number of sub-steps per integration (always larger 156 /// than or equal to 1). getSubSteps()157 int getSubSteps() const { return mSubSteps; } 158 159 /// @brief Set the number of sub-steps per integration. 160 /// @note The only reason to increase the sub-step above its 161 /// default value of one is to reduce the memory footprint 162 /// due to significant dilation. Values smaller than 1 will 163 /// be clamped to 1! setSubSteps(int substeps)164 void setSubSteps(int substeps) { mSubSteps = math::Max(1, substeps); } 165 166 /// @brief Return the maximum magnitude of the velocity in the 167 /// advection velocity field defined during construction. getMaxVelocity()168 double getMaxVelocity() const { return mMaxVelocity; } 169 170 /// @return Returns the maximum distance in voxel units of @a inGrid 171 /// that a particle can travel in the time-step @a dt when advected 172 /// in the velocity field defined during construction. 173 /// 174 /// @details This method is useful when dilating sparse volume 175 /// grids to pad boundary regions. Excessive dilation can be 176 /// computationally expensive so use this method to prevent 177 /// or warn against run-away computation. 178 /// 179 /// @throw RuntimeError if @a inGrid does not have uniform voxels. 180 template<typename VolumeGridT> getMaxDistance(const VolumeGridT & inGrid,double dt)181 int getMaxDistance(const VolumeGridT& inGrid, double dt) const 182 { 183 if (!inGrid.hasUniformVoxels()) { 184 OPENVDB_THROW(RuntimeError, "Volume grid does not have uniform voxels!"); 185 } 186 const double d = mMaxVelocity*math::Abs(dt)/inGrid.voxelSize()[0]; 187 return static_cast<int>( math::RoundUp(d) ); 188 } 189 190 /// @return Returns a new grid that is the result of passive advection 191 /// of all the active values the input grid by @a timeStep. 192 /// 193 /// @param inGrid The input grid to be advected (unmodified) 194 /// @param timeStep Time-step of the Runge-Kutta integrator. 195 /// 196 /// @details This method will advect all of the active values in 197 /// the input @a inGrid. To achieve this a 198 /// deep-copy is dilated to account for the material 199 /// transport. This dilation step can be slow for large 200 /// time steps @a dt or a velocity field with large magnitudes. 201 /// 202 /// @warning If the VolumeSamplerT is of higher order than one 203 /// (i.e. tri-linear interpolation) instabilities are 204 /// known to occur. To suppress those monotonicity 205 /// constrains or flux-limiters need to be applies. 206 /// 207 /// @throw RuntimeError if @a inGrid does not have uniform voxels. 208 template<typename VolumeGridT, 209 typename VolumeSamplerT>//only C++11 allows for a default argument advect(const VolumeGridT & inGrid,double timeStep)210 typename VolumeGridT::Ptr advect(const VolumeGridT& inGrid, double timeStep) 211 { 212 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy(); 213 const double dt = timeStep/mSubSteps; 214 const int n = this->getMaxDistance(inGrid, dt); 215 dilateActiveValues( outGrid->tree(), n, NN_FACE, EXPAND_TILES); 216 this->template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt); 217 for (int step = 1; step < mSubSteps; ++step) { 218 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy(); 219 dilateActiveValues( tmpGrid->tree(), n, NN_FACE, EXPAND_TILES); 220 this->template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt); 221 outGrid.swap( tmpGrid ); 222 } 223 224 return outGrid; 225 } 226 227 /// @return Returns a new grid that is the result of 228 /// passive advection of the active values in @a inGrid 229 /// that intersect the active values in @c mask. The time 230 /// of the output grid is incremented by @a timeStep. 231 /// 232 /// @param inGrid The input grid to be advected (unmodified). 233 /// @param mask The mask of active values defining the active voxels 234 /// in @c inGrid on which to perform advection. Only 235 /// if a value is active in both grids will it be modified. 236 /// @param timeStep Time-step for a single Runge-Kutta integration step. 237 /// 238 /// 239 /// @details This method will advect all of the active values in 240 /// the input @a inGrid that intersects with the 241 /// active values in @a mask. To achieve this a 242 /// deep-copy is dilated to account for the material 243 /// transport and finally cropped to the intersection 244 /// with @a mask. The dilation step can be slow for large 245 /// time steps @a dt or fast moving velocity fields. 246 /// 247 /// @warning If the VolumeSamplerT is of higher order the one 248 /// (i.e. tri-linear interpolation) instabilities are 249 /// known to occur. To suppress those monotonicity 250 /// constrains or flux-limiters need to be applies. 251 /// 252 /// @throw RuntimeError if @a inGrid is not aligned with @a mask 253 /// or if its voxels are not uniform. 254 template<typename VolumeGridT, 255 typename MaskGridT, 256 typename VolumeSamplerT>//only C++11 allows for a default argument advect(const VolumeGridT & inGrid,const MaskGridT & mask,double timeStep)257 typename VolumeGridT::Ptr advect(const VolumeGridT& inGrid, const MaskGridT& mask, double timeStep) 258 { 259 if (inGrid.transform() != mask.transform()) { 260 OPENVDB_THROW(RuntimeError, "Volume grid and mask grid are misaligned! Consider " 261 "resampling either of the two grids into the index space of the other."); 262 } 263 typename VolumeGridT::Ptr outGrid = inGrid.deepCopy(); 264 const double dt = timeStep/mSubSteps; 265 const int n = this->getMaxDistance(inGrid, dt); 266 dilateActiveValues( outGrid->tree(), n, NN_FACE, EXPAND_TILES); 267 outGrid->topologyIntersection( mask ); 268 pruneInactive( outGrid->tree(), mGrainSize>0, mGrainSize ); 269 this->template cook<VolumeGridT, VolumeSamplerT>(*outGrid, inGrid, dt); 270 outGrid->topologyUnion( inGrid ); 271 272 for (int step = 1; step < mSubSteps; ++step) { 273 typename VolumeGridT::Ptr tmpGrid = outGrid->deepCopy(); 274 dilateActiveValues( tmpGrid->tree(), n, NN_FACE, EXPAND_TILES); 275 tmpGrid->topologyIntersection( mask ); 276 pruneInactive( tmpGrid->tree(), mGrainSize>0, mGrainSize ); 277 this->template cook<VolumeGridT, VolumeSamplerT>(*tmpGrid, *outGrid, dt); 278 tmpGrid->topologyUnion( inGrid ); 279 outGrid.swap( tmpGrid ); 280 } 281 return outGrid; 282 } 283 284 private: 285 // disallow copy construction and copy by assignment! 286 VolumeAdvection(const VolumeAdvection&);// not implemented 287 VolumeAdvection& operator=(const VolumeAdvection&);// not implemented 288 start(const char * str)289 void start(const char* str) const 290 { 291 if (mInterrupter) mInterrupter->start(str); 292 } stop()293 void stop() const 294 { 295 if (mInterrupter) mInterrupter->end(); 296 } interrupt()297 bool interrupt() const 298 { 299 if (mInterrupter && util::wasInterrupted(mInterrupter)) { 300 thread::cancelGroupExecution(); 301 return true; 302 } 303 return false; 304 } 305 306 template<typename VolumeGridT, typename VolumeSamplerT> cook(VolumeGridT & outGrid,const VolumeGridT & inGrid,double dt)307 void cook(VolumeGridT& outGrid, const VolumeGridT& inGrid, double dt) 308 { 309 switch (mIntegrator) { 310 case Scheme::SEMI: { 311 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *this); 312 adv.cook(outGrid, dt); 313 break; 314 } 315 case Scheme::MID: { 316 Advect<VolumeGridT, 2, VolumeSamplerT> adv(inGrid, *this); 317 adv.cook(outGrid, dt); 318 break; 319 } 320 case Scheme::RK3: { 321 Advect<VolumeGridT, 3, VolumeSamplerT> adv(inGrid, *this); 322 adv.cook(outGrid, dt); 323 break; 324 } 325 case Scheme::RK4: { 326 Advect<VolumeGridT, 4, VolumeSamplerT> adv(inGrid, *this); 327 adv.cook(outGrid, dt); 328 break; 329 } 330 case Scheme::BFECC: { 331 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *this); 332 adv.cook(outGrid, dt); 333 break; 334 } 335 case Scheme::MAC: { 336 Advect<VolumeGridT, 1, VolumeSamplerT> adv(inGrid, *this); 337 adv.cook(outGrid, dt); 338 break; 339 } 340 default: 341 OPENVDB_THROW(ValueError, "Spatial difference scheme not supported!"); 342 } 343 pruneInactive(outGrid.tree(), mGrainSize>0, mGrainSize); 344 } 345 346 // Private class that implements the multi-threaded advection 347 template<typename VolumeGridT, size_t OrderRK, typename SamplerT> struct Advect; 348 349 // Private member data of VolumeAdvection 350 const VelocityGridT& mVelGrid; 351 double mMaxVelocity; 352 InterrupterType* mInterrupter; 353 Scheme::SemiLagrangian mIntegrator; 354 Scheme::Limiter mLimiter; 355 size_t mGrainSize; 356 int mSubSteps; 357 };//end of VolumeAdvection class 358 359 // Private class that implements the multi-threaded advection 360 template<typename VelocityGridT, bool StaggeredVelocity, typename InterrupterType> 361 template<typename VolumeGridT, size_t OrderRK, typename SamplerT> 362 struct VolumeAdvection<VelocityGridT, StaggeredVelocity, InterrupterType>::Advect 363 { 364 using TreeT = typename VolumeGridT::TreeType; 365 using AccT = typename VolumeGridT::ConstAccessor; 366 using ValueT = typename TreeT::ValueType; 367 using LeafManagerT = typename tree::LeafManager<TreeT>; 368 using LeafNodeT = typename LeafManagerT::LeafNodeType; 369 using LeafRangeT = typename LeafManagerT::LeafRange; 370 using VelocityIntegratorT = VelocityIntegrator<VelocityGridT, StaggeredVelocity>; 371 using RealT = typename VelocityIntegratorT::ElementType; 372 using VoxelIterT = typename TreeT::LeafNodeType::ValueOnIter; 373 374 Advect(const VolumeGridT& inGrid, const VolumeAdvection& parent) 375 : mTask(nullptr) 376 , mInGrid(&inGrid) 377 , mVelocityInt(parent.mVelGrid) 378 , mParent(&parent) 379 { 380 } 381 inline void cook(const LeafRangeT& range) 382 { 383 if (mParent->mGrainSize > 0) { 384 tbb::parallel_for(range, *this); 385 } else { 386 (*this)(range); 387 } 388 } 389 void operator()(const LeafRangeT& range) const 390 { 391 assert(mTask); 392 mTask(const_cast<Advect*>(this), range); 393 } 394 void cook(VolumeGridT& outGrid, double time_step) 395 { 396 namespace ph = std::placeholders; 397 398 mParent->start("Advecting volume"); 399 LeafManagerT manager(outGrid.tree(), mParent->spatialOrder()==2 ? 1 : 0); 400 const LeafRangeT range = manager.leafRange(mParent->mGrainSize); 401 const RealT dt = static_cast<RealT>(-time_step);//method of characteristics backtracks 402 if (mParent->mIntegrator == Scheme::MAC) { 403 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);//out[0]=forward 404 this->cook(range); 405 mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);//out[1]=backward 406 this->cook(range); 407 mTask = std::bind(&Advect::mac, ph::_1, ph::_2);//out[0] = out[0] + (in[0] - out[1])/2 408 this->cook(range); 409 } else if (mParent->mIntegrator == Scheme::BFECC) { 410 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);//out[0]=forward 411 this->cook(range); 412 mTask = std::bind(&Advect::rk, ph::_1, ph::_2,-dt, 1, &outGrid);//out[1]=backward 413 this->cook(range); 414 mTask = std::bind(&Advect::bfecc, ph::_1, ph::_2);//out[0] = (3*in[0] - out[1])/2 415 this->cook(range); 416 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 1, &outGrid);//out[1]=forward 417 this->cook(range); 418 manager.swapLeafBuffer(1);// out[0] = out[1] 419 } else {// SEMI, MID, RK3 and RK4 420 mTask = std::bind(&Advect::rk, ph::_1, ph::_2, dt, 0, mInGrid);//forward 421 this->cook(range); 422 } 423 424 if (mParent->spatialOrder()==2) manager.removeAuxBuffers(); 425 426 mTask = std::bind(&Advect::limiter, ph::_1, ph::_2, dt);// out[0] = limiter( out[0] ) 427 this->cook(range); 428 429 mParent->stop(); 430 } 431 // Last step of the MacCormack scheme: out[0] = out[0] + (in[0] - out[1])/2 432 void mac(const LeafRangeT& range) const 433 { 434 if (mParent->interrupt()) return; 435 assert( mParent->mIntegrator == Scheme::MAC ); 436 AccT acc = mInGrid->getAccessor(); 437 for (typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) { 438 ValueT* out0 = leafIter.buffer( 0 ).data();// forward 439 const ValueT* out1 = leafIter.buffer( 1 ).data();// backward 440 const LeafNodeT* leaf = acc.probeConstLeaf( leafIter->origin() ); 441 if (leaf != nullptr) { 442 const ValueT* in0 = leaf->buffer().data(); 443 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) { 444 const Index i = voxelIter.pos(); 445 out0[i] += RealT(0.5) * ( in0[i] - out1[i] ); 446 } 447 } else { 448 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) { 449 const Index i = voxelIter.pos(); 450 out0[i] += RealT(0.5) * ( acc.getValue(voxelIter.getCoord()) - out1[i] ); 451 }//loop over active voxels 452 } 453 }//loop over leaf nodes 454 } 455 // Intermediate step in the BFECC scheme: out[0] = (3*in[0] - out[1])/2 456 void bfecc(const LeafRangeT& range) const 457 { 458 if (mParent->interrupt()) return; 459 assert( mParent->mIntegrator == Scheme::BFECC ); 460 AccT acc = mInGrid->getAccessor(); 461 for (typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) { 462 ValueT* out0 = leafIter.buffer( 0 ).data();// forward 463 const ValueT* out1 = leafIter.buffer( 1 ).data();// backward 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)*( RealT(3)*in0[i] - out1[i] ); 470 }//loop over active voxels 471 } else { 472 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) { 473 const Index i = voxelIter.pos(); 474 out0[i] = RealT(0.5)*( RealT(3)*acc.getValue(voxelIter.getCoord()) - out1[i] ); 475 }//loop over active voxels 476 } 477 }//loop over leaf nodes 478 } 479 // Semi-Lagrangian integration with Runge-Kutta of various orders (1->4) 480 void rk(const LeafRangeT& range, RealT dt, size_t n, const VolumeGridT* grid) const 481 { 482 if (mParent->interrupt()) return; 483 const math::Transform& xform = mInGrid->transform(); 484 AccT acc = grid->getAccessor(); 485 for (typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) { 486 ValueT* phi = leafIter.buffer( n ).data(); 487 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) { 488 ValueT& value = phi[voxelIter.pos()]; 489 Vec3d wPos = xform.indexToWorld(voxelIter.getCoord()); 490 mVelocityInt.template rungeKutta<OrderRK, Vec3d>(dt, wPos); 491 value = SamplerT::sample(acc, xform.worldToIndex(wPos)); 492 }//loop over active voxels 493 }//loop over leaf nodes 494 } 495 void limiter(const LeafRangeT& range, RealT dt) const 496 { 497 if (mParent->interrupt()) return; 498 const bool doLimiter = mParent->isLimiterOn(); 499 const bool doClamp = mParent->mLimiter == Scheme::CLAMP; 500 ValueT data[2][2][2], vMin, vMax; 501 const math::Transform& xform = mInGrid->transform(); 502 AccT acc = mInGrid->getAccessor(); 503 const ValueT backg = mInGrid->background(); 504 for (typename LeafRangeT::Iterator leafIter = range.begin(); leafIter; ++leafIter) { 505 ValueT* phi = leafIter.buffer( 0 ).data(); 506 for (VoxelIterT voxelIter = leafIter->beginValueOn(); voxelIter; ++voxelIter) { 507 ValueT& value = phi[voxelIter.pos()]; 508 509 if ( doLimiter ) { 510 assert(OrderRK == 1); 511 Vec3d wPos = xform.indexToWorld(voxelIter.getCoord()); 512 mVelocityInt.template rungeKutta<1, Vec3d>(dt, wPos);// Explicit Euler 513 Vec3d iPos = xform.worldToIndex(wPos); 514 Coord ijk = Coord::floor( iPos ); 515 BoxSampler::getValues(data, acc, ijk); 516 BoxSampler::extrema(data, vMin, vMax); 517 if ( doClamp ) { 518 value = math::Clamp( value, vMin, vMax); 519 } else if (value < vMin || value > vMax ) { 520 iPos -= Vec3R(ijk[0], ijk[1], ijk[2]);//unit coordinates 521 value = BoxSampler::trilinearInterpolation( data, iPos ); 522 } 523 } 524 525 if (math::isApproxEqual(value, backg, math::Delta<ValueT>::value())) { 526 value = backg; 527 leafIter->setValueOff( voxelIter.pos() ); 528 } 529 }//loop over active voxels 530 }//loop over leaf nodes 531 } 532 // Public member data of the private Advect class 533 534 typename std::function<void (Advect*, const LeafRangeT&)> mTask; 535 const VolumeGridT* mInGrid; 536 const VelocityIntegratorT mVelocityInt;// lightweight! 537 const VolumeAdvection* mParent; 538 };// end of private member class Advect 539 540 541 //////////////////////////////////////// 542 543 544 // Explicit Template Instantiation 545 546 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION 547 548 #ifdef OPENVDB_INSTANTIATE_VOLUMEADVECT 549 #include <openvdb/util/ExplicitInstantiation.h> 550 #endif 551 552 OPENVDB_INSTANTIATE_CLASS VolumeAdvection<Vec3SGrid, true, util::NullInterrupter>; 553 OPENVDB_INSTANTIATE_CLASS VolumeAdvection<Vec3SGrid, false, util::NullInterrupter>; 554 555 OPENVDB_INSTANTIATE FloatGrid::Ptr VolumeAdvection<Vec3SGrid, true, util::NullInterrupter>::advect<FloatGrid, Sampler<1, false>>(const FloatGrid&, double); 556 OPENVDB_INSTANTIATE DoubleGrid::Ptr VolumeAdvection<Vec3SGrid, true, util::NullInterrupter>::advect<DoubleGrid, Sampler<1, false>>(const DoubleGrid&, double); 557 OPENVDB_INSTANTIATE Vec3SGrid::Ptr VolumeAdvection<Vec3SGrid, true, util::NullInterrupter>::advect<Vec3SGrid, Sampler<1, false>>(const Vec3SGrid&, double); 558 559 OPENVDB_INSTANTIATE FloatGrid::Ptr VolumeAdvection<Vec3SGrid, false, util::NullInterrupter>::advect<FloatGrid, Sampler<1, false>>(const FloatGrid&, double); 560 OPENVDB_INSTANTIATE DoubleGrid::Ptr VolumeAdvection<Vec3SGrid, false, util::NullInterrupter>::advect<DoubleGrid, Sampler<1, false>>(const DoubleGrid&, double); 561 OPENVDB_INSTANTIATE Vec3SGrid::Ptr VolumeAdvection<Vec3SGrid, false, util::NullInterrupter>::advect<Vec3SGrid, Sampler<1, false>>(const Vec3SGrid&, double); 562 563 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION 564 565 566 } // namespace tools 567 } // namespace OPENVDB_VERSION_NAME 568 } // namespace openvdb 569 570 #endif // OPENVDB_TOOLS_VOLUME_ADVECT_HAS_BEEN_INCLUDED 571