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