HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VolumeAdvect.h
Go to the documentation of this file.
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>
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 {
33 namespace OPENVDB_VERSION_NAME {
34 namespace tools {
35 
36 namespace Scheme {
37  /// @brief Numerical advections schemes.
39  /// @brief Flux-limiters employed to stabilize the second-order
40  /// advection schemes MacCormack and BFECC.
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>
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)
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 
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.
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.
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)
129  void setIntegrator(Scheme::SemiLagrangian integrator) { mIntegrator = integrator; }
130 
131  /// @brief Return the integrator (see details in the table above)
132  Scheme::SemiLagrangian getIntegrator() const { return mIntegrator; }
133 
134  /// @brief Set the limiter (see details above)
135  void setLimiter(Scheme::Limiter limiter) { mLimiter = limiter; }
136 
137  /// @brief Retrun the limiter (see details above)
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.
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
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!
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).
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!
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.
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>
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
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
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 
289  void start(const char* str) const
290  {
291  if (mInterrupter) mInterrupter->start(str);
292  }
293  void stop() const
294  {
295  if (mInterrupter) mInterrupter->end();
296  }
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>
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
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
#define CLAMP(EXPR, TYPE)
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
Definition: parallel.h:127
GLenum GLint * range
Definition: glcorearb.h:1925
math::Vec3< Real > Vec3R
Definition: Types.h:72
SemiLagrangian
Numerical advections schemes.
Definition: VolumeAdvect.h:38
int getMaxDistance(const VolumeGridT &inGrid, double dt) const
Definition: VolumeAdvect.h:181
IMATH_HOSTDEVICE constexpr int floor(T x) IMATH_NOEXCEPT
Definition: ImathFun.h:112
void setGrainSize(size_t grainsize)
Set the grain-size used for multi-threading.
Definition: VolumeAdvect.h:153
GLuint start
Definition: glcorearb.h:475
GLsizei const GLfloat * value
Definition: glcorearb.h:824
static ValueT trilinearInterpolation(ValueT(&data)[N][N][N], const Vec3R &uvw)
int temporalOrder() const
Return the temporal order of accuracy of the advection scheme.
Definition: VolumeAdvect.h:116
void setSubSteps(int substeps)
Set the number of sub-steps per integration.
Definition: VolumeAdvect.h:164
static void extrema(ValueT(&data)[N][N][N], ValueT &vMin, ValueT &vMax)
Find the minimum and maximum values of the eight cell values in @ data.
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:239
float RoundUp(float x)
Return x rounded up to the nearest integer.
Definition: Math.h:787
Performs advections of an arbitrary type of volume in a static velocity field. The advections are per...
Definition: VolumeAdvect.h:75
VolumeGridT::Ptr advect(const VolumeGridT &inGrid, const MaskGridT &mask, double timeStep)
Definition: VolumeAdvect.h:257
Limiter
Flux-limiters employed to stabilize the second-order advection schemes MacCormack and BFECC...
Definition: VolumeAdvect.h:41
void dilateActiveValues(TreeOrLeafManagerT &tree, const int iterations=1, const NearestNeighbors nn=NN_FACE, const TilePolicy mode=PRESERVE_TILES, const bool threaded=true)
Topologically dilate all active values (i.e. both voxels and tiles) in a tree using one of three near...
Definition: Morphology.h:1056
Vec3< double > Vec3d
Definition: NanoVDB.h:1685
GLdouble n
Definition: glcorearb.h:2008
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
#define OPENVDB_INSTANTIATE_CLASS
Coord Abs(const Coord &xyz)
Definition: Coord.h:517
Scheme::SemiLagrangian getIntegrator() const
Return the integrator (see details in the table above)
Definition: VolumeAdvect.h:132
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
GLint GLuint mask
Definition: glcorearb.h:124
Defined various multi-threaded utility functions for trees.
VolumeAdvection(const VelocityGridT &velGrid, InterrupterType *interrupter=nullptr)
Constructor.
Definition: VolumeAdvect.h:86
static void getValues(ValueT(&data)[N][N][N], const TreeT &inTree, Coord ijk)
Import all eight values from inTree to support tri-linear interpolation.
Type Clamp(Type x, Type min, Type max)
Return x clamped to [min, max].
Definition: Math.h:261
#define OPENVDB_INSTANTIATE
bool isLimiterOn() const
Return true if a limiter will be applied based on the current settings.
Definition: VolumeAdvect.h:142
double getMaxVelocity() const
Return the maximum magnitude of the velocity in the advection velocity field defined during construct...
Definition: VolumeAdvect.h:168
void setIntegrator(Scheme::SemiLagrangian integrator)
Set the integrator (see details in the table above)
Definition: VolumeAdvect.h:129
Defines two simple wrapper classes for advection velocity fields as well as VelocitySampler and Veloc...
int spatialOrder() const
Return the spatial order of accuracy of the advection scheme.
Definition: VolumeAdvect.h:108
VolumeGridT::Ptr advect(const VolumeGridT &inGrid, double timeStep)
Definition: VolumeAdvect.h:210
Implementation of morphological dilation and erosion.
Functions to efficiently compute histograms, extrema (min/max) and statistics (mean, variance, etc.) of grid values.
Scheme::Limiter getLimiter() const
Retrun the limiter (see details above)
Definition: VolumeAdvect.h:138
Definition: core.h:1131
Grid< DoubleTree > DoubleGrid
Definition: openvdb.h:74
Grid< FloatTree > FloatGrid
Definition: openvdb.h:75
bool wasInterrupted(T *i, int percent=-1)
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:119
void setLimiter(Scheme::Limiter limiter)
Set the limiter (see details above)
Definition: VolumeAdvect.h:135
bool ValueType
Definition: NanoVDB.h:5729
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:595
Grid< Vec3STree > Vec3SGrid
Definition: openvdb.h:81
double max() const
Return the maximum value.
Definition: Stats.h:124
math::Extrema extrema(const IterT &iter, bool threaded=true)
Iterate over a scalar grid and compute extrema (min/max) of the values of the voxels that are visited...
Definition: Statistics.h:354
Definition: format.h:895
void add(double val)
Add a single sample.
Definition: Stats.h:102
void pruneInactive(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with background tiles any nodes whose values are a...
Definition: Prune.h:355
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
This class computes the minimum and maximum values of a population of floating-point values...
Definition: Stats.h:88