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