39 #ifndef OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
40 #define OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
44 #include <openvdb/math/FiniteDifference.h>
65 template <
typename VelGr
idT,
typename Interpolator = BoxSampler>
72 DiscreteField(
const VelGridT &vel): mAccessor(vel.tree()), mTransform(&vel.transform()) {}
83 Interpolator::sample(mAccessor, mTransform->worldToIndex(xyz), result);
90 return mAccessor.getValue(ijk);
94 const typename VelGridT::ConstAccessor mAccessor;
105 template <
typename ScalarT =
float>
121 inline VectorType operator() (
const Vec3d& xyz, ScalarType time)
const;
126 return (*
this)(ijk.asVec3d(), time);
170 template<
typename GridT,
171 typename FieldT = EnrightField<typename GridT::ValueType>,
186 mTracker(grid, interrupt), mField(field),
188 mTemporalScheme(math::
TVD_RK2) {}
229 size_t advect(ScalarType time0, ScalarType time1);
242 LevelSetAdvect(
const LevelSetAdvect& other);
244 LevelSetAdvect(LevelSetAdvect& other, tbb::split);
246 virtual ~LevelSetAdvect() {
if (mIsMaster) this->clearField();};
249 size_t advect(ScalarType time0, ScalarType time1);
251 void operator()(
const RangeType& r)
const
253 if (mTask) mTask(const_cast<LevelSetAdvect*>(
this), r);
257 void operator()(
const RangeType& r)
259 if (mTask) mTask(
this, r);
263 void join(
const LevelSetAdvect& other) { mMaxAbsV =
math::Max(mMaxAbsV, other.mMaxAbsV); }
265 typedef typename boost::function<void (LevelSetAdvect*, const RangeType&)> FuncType;
266 LevelSetAdvection& mParent;
268 const ScalarType mMinAbsV;
272 const bool mIsMaster;
274 enum ThreadingMode { PARALLEL_FOR, PARALLEL_REDUCE };
276 void cook(ThreadingMode mode,
size_t swapBuffer = 0);
278 typename GridT::ValueType sampleField(ScalarType time0, ScalarType time1);
280 void sampleXformedField(
const RangeType& r, ScalarType time0, ScalarType time1);
281 void sampleAlignedField(
const RangeType& r, ScalarType time0, ScalarType time1);
283 void euler1(
const RangeType& r, ScalarType dt,
Index resultBuffer);
286 void euler2(
const RangeType& r, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer);
289 template<math::BiasedGradientScheme SpatialScheme>
290 size_t advect1(ScalarType time0, ScalarType time1);
294 size_t advect2(ScalarType time0, ScalarType time1);
299 size_t advect3(ScalarType time0, ScalarType time1);
308 void operator=(
const LevelSetAdvection& other) {}
312 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
316 switch (mSpatialScheme) {
318 return this->advect1<math::FIRST_BIAS >(time0, time1);
320 return this->advect1<math::SECOND_BIAS >(time0, time1);
322 return this->advect1<math::THIRD_BIAS >(time0, time1);
324 return this->advect1<math::WENO5_BIAS >(time0, time1);
326 return this->advect1<math::HJWENO5_BIAS>(time0, time1);
333 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
334 template<math::BiasedGradientScheme SpatialScheme>
338 switch (mTemporalScheme) {
340 return this->advect2<SpatialScheme, math::TVD_RK1>(time0, time1);
342 return this->advect2<SpatialScheme, math::TVD_RK2>(time0, time1);
344 return this->advect2<SpatialScheme, math::TVD_RK3>(time0, time1);
351 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
355 LevelSetAdvection<GridT, FieldT, InterruptT>::advect2(ScalarType time0, ScalarType time1)
358 if (trans.
mapType() == math::UniformScaleMap::mapType()) {
359 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleMap>(time0, time1);
360 }
else if (trans.
mapType() == math::UniformScaleTranslateMap::mapType()) {
361 return this->advect3<SpatialScheme, TemporalScheme, math::UniformScaleTranslateMap>(time0, time1);
362 }
else if (trans.
mapType() == math::UnitaryMap::mapType()) {
363 return this->advect3<SpatialScheme, TemporalScheme, math::UnitaryMap >(time0, time1);
364 }
else if (trans.
mapType() == math::TranslationMap::mapType()) {
365 return this->advect3<SpatialScheme, TemporalScheme, math::TranslationMap>(time0, time1);
372 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
377 LevelSetAdvection<GridT, FieldT, InterruptT>::advect3(ScalarType time0, ScalarType time1)
379 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme> tmp(*
this);
380 return tmp.advect(time0, time1);
385 template <
typename ScalarT>
386 inline math::Vec3<ScalarT>
389 static const ScalarT pi = ScalarT(3.1415926535897931), phase = pi / ScalarT(3.0);
390 const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
391 const ScalarT tr = cos(ScalarT(time) * phase);
392 const ScalarT a = sin(ScalarT(2.0)*Py);
393 const ScalarT b = -sin(ScalarT(2.0)*Px);
394 const ScalarT c = sin(ScalarT(2.0)*Pz);
396 tr * ( ScalarT(2) *
math::Pow2(sin(Px)) * a * c ),
405 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
415 mMap(parent.mTracker.grid().transform().template constMap<MapT>().get()),
421 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
425 LevelSetAdvection<GridT, FieldT, InterruptT>::
426 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
427 LevelSetAdvect(
const LevelSetAdvect& other):
428 mParent(other.mParent),
430 mMinAbsV(other.mMinAbsV),
431 mMaxAbsV(other.mMaxAbsV),
438 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
442 LevelSetAdvection<GridT, FieldT, InterruptT>::
443 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
444 LevelSetAdvect(LevelSetAdvect& other, tbb::split):
445 mParent(other.mParent),
447 mMinAbsV(other.mMinAbsV),
448 mMaxAbsV(other.mMaxAbsV),
455 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
461 advect(ScalarType time0, ScalarType time1)
465 const bool isForward = time0 < time1;
466 while ((isForward ? time0<time1 : time0>time1) && mParent.mTracker.checkInterrupter()) {
468 mParent.mTracker.leafs().rebuildAuxBuffers(TemporalScheme ==
math::TVD_RK3 ? 2 : 1);
470 const ScalarType dt = this->sampleField(time0, time1);
473 switch(TemporalScheme) {
477 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
479 this->cook(PARALLEL_FOR, 1);
484 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
486 this->cook(PARALLEL_FOR, 1);
490 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.5), 1, 1);
492 this->cook(PARALLEL_FOR, 1);
497 mTask = boost::bind(&LevelSetAdvect::euler1, _1, _2, dt, 1);
499 this->cook(PARALLEL_FOR, 1);
503 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(0.75), 1, 2);
505 this->cook(PARALLEL_FOR, 2);
509 mTask = boost::bind(&LevelSetAdvect::euler2, _1, _2, dt, ScalarType(1.0/3.0), 1, 2);
511 this->cook(PARALLEL_FOR, 2);
514 OPENVDB_THROW(ValueError,
"Temporal integration scheme not supported!");
516 time0 += isForward ? dt : -dt;
518 mParent.mTracker.leafs().removeAuxBuffers();
521 mParent.mTracker.track();
526 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
529 inline typename GridT::ValueType
530 LevelSetAdvection<GridT, FieldT, InterruptT>::
531 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
532 sampleField(ScalarType time0, ScalarType time1)
535 const size_t leafCount = mParent.mTracker.leafs().leafCount();
536 if (leafCount==0)
return ScalarType(0.0);
537 mVec =
new VectorType*[leafCount];
538 if (mParent.mField.transform() == mParent.mTracker.grid().transform()) {
539 mTask = boost::bind(&LevelSetAdvect::sampleAlignedField, _1, _2, time0, time1);
541 mTask = boost::bind(&LevelSetAdvect::sampleXformedField, _1, _2, time0, time1);
543 this->cook(PARALLEL_REDUCE);
545 static const ScalarType CFL = (TemporalScheme ==
math::TVD_RK1 ? ScalarType(0.3) :
546 TemporalScheme == math::
TVD_RK2 ? ScalarType(0.9) :
547 ScalarType(1.0))/math::
Sqrt(ScalarType(3.0));
548 const ScalarType dt =
math::Abs(time1 - time0), dx = mParent.mTracker.voxelSize();
552 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
556 LevelSetAdvection<GridT, FieldT, InterruptT>::
557 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
558 sampleXformedField(
const RangeType& range, ScalarType time0, ScalarType time1)
560 const bool isForward = time0 < time1;
561 typedef typename LeafType::ValueOnCIter VoxelIterT;
562 const MapT& map = *
mMap;
563 mParent.mTracker.checkInterrupter();
564 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
565 const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
566 VectorType* vec =
new VectorType[leaf.onVoxelCount()];
568 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
569 const VectorType V = mParent.mField(map.applyMap(iter.getCoord().asVec3d()), time0);
571 vec[m] = isForward ? V : -V;
577 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
581 LevelSetAdvection<GridT, FieldT, InterruptT>::
582 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
583 sampleAlignedField(
const RangeType& range, ScalarType time0, ScalarType time1)
585 const bool isForward = time0 < time1;
586 typedef typename LeafType::ValueOnCIter VoxelIterT;
587 mParent.mTracker.checkInterrupter();
588 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
589 const LeafType& leaf = mParent.mTracker.leafs().leaf(n);
590 VectorType* vec =
new VectorType[leaf.onVoxelCount()];
592 for (VoxelIterT iter = leaf.cbeginValueOn(); iter; ++iter, ++m) {
593 const VectorType V = mParent.mField(iter.getCoord(), time0);
595 vec[m] = isForward ? V : -V;
601 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
605 LevelSetAdvection<GridT, FieldT, InterruptT>::
606 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
609 if (mVec == NULL)
return;
610 for (
size_t n=0, e=mParent.mTracker.leafs().leafCount(); n<e; ++n)
delete [] mVec[n];
615 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
619 LevelSetAdvection<GridT, FieldT, InterruptT>::
620 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
621 cook(ThreadingMode mode,
size_t swapBuffer)
623 mParent.mTracker.startInterrupter(
"Advecting level set");
625 if (mParent.mTracker.getGrainSize()==0) {
626 (*this)(mParent.mTracker.leafs().getRange());
627 }
else if (mode == PARALLEL_FOR) {
628 tbb::parallel_for(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *
this);
629 }
else if (mode == PARALLEL_REDUCE) {
630 tbb::parallel_reduce(mParent.mTracker.leafs().getRange(mParent.mTracker.getGrainSize()), *
this);
632 throw std::runtime_error(
"Undefined threading mode");
635 mParent.mTracker.leafs().swapLeafBuffer(swapBuffer, mParent.mTracker.getGrainSize()==0);
637 mParent.mTracker.endInterrupter();
642 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
646 LevelSetAdvection<GridT, FieldT, InterruptT>::
647 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
648 euler1(
const RangeType& range, ScalarType dt,
Index resultBuffer)
650 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
651 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
652 typedef typename LeafType::ValueOnCIter VoxelIterT;
653 mParent.mTracker.checkInterrupter();
654 const MapT& map = *
mMap;
655 typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
656 Stencil stencil(mParent.mTracker.grid());
657 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
658 BufferType& result = leafs.getBuffer(n, resultBuffer);
659 const VectorType* vec = mVec[n];
661 for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
662 stencil.moveTo(iter);
664 result.setValue(iter.pos(), *iter - dt * V.dot(G));
671 template<
typename Gr
idT,
typename FieldT,
typename InterruptT>
675 LevelSetAdvection<GridT, FieldT, InterruptT>::
676 LevelSetAdvect<MapT, SpatialScheme, TemporalScheme>::
677 euler2(
const RangeType& range, ScalarType dt, ScalarType alpha,
Index phiBuffer,
Index resultBuffer)
679 typedef math::BIAS_SCHEME<SpatialScheme> Scheme;
680 typedef typename Scheme::template ISStencil<GridType>::StencilType Stencil;
681 typedef typename LeafType::ValueOnCIter VoxelIterT;
682 mParent.mTracker.checkInterrupter();
683 const MapT& map = *
mMap;
684 typename TrackerT::LeafManagerType& leafs = mParent.mTracker.leafs();
685 const ScalarType beta = ScalarType(1.0) - alpha;
686 Stencil stencil(mParent.mTracker.grid());
687 for (
size_t n=range.begin(), e=range.end(); n != e; ++n) {
688 const BufferType& phi = leafs.getBuffer(n, phiBuffer);
689 BufferType& result = leafs.getBuffer(n, resultBuffer);
690 const VectorType* vec = mVec[n];
692 for (VoxelIterT iter = leafs.leaf(n).cbeginValueOn(); iter; ++iter, ++m) {
693 stencil.moveTo(iter);
695 result.setValue(iter.pos(), alpha*phi[iter.pos()] + beta*(*iter - dt * V.dot(G)));
704 #endif // OPENVDB_TOOLS_LEVEL_SET_ADVECT_HAS_BEEN_INCLUDED
Definition: FiniteDifference.h:264
Definition: FiniteDifference.h:263
Definition: FiniteDifference.h:195
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
TemporalIntegrationScheme
Temporal integrations schemes.
Definition: FiniteDifference.h:261
const MapT & mMap
Definition: GridOperators.h:260
Definition: FiniteDifference.h:198
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:648
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Definition: FiniteDifference.h:194
Type Pow2(Type x)
Return .
Definition: Math.h:460
const Type & Max(const Type &a, const Type &b)
Return the maximum of two values.
Definition: Math.h:507
Definition: FiniteDifference.h:265
#define OPENVDB_VERSION_NAME
Definition: version.h:45
bool isZero(const Type &x)
Return true if x is exactly equal to zero.
Definition: Math.h:276
static math::Vec3< typename Accessor::ValueType > result(const MapType &map, const Accessor &grid, const Coord &ijk, const Vec3< typename Accessor::ValueType > &V)
Definition: Operators.h:798
Vec3< double > Vec3d
Definition: Vec3.h:605
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:353
Index32 Index
Definition: Types.h:56
Definition: FiniteDifference.h:197
Definition: Exceptions.h:88
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:246
BiasedGradientScheme
Biased Gradients are limited to non-centered differences.
Definition: FiniteDifference.h:192
Definition: FiniteDifference.h:196
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:56
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:568
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52