52 #ifndef OPENVDB_TOOLS_VELOCITY_FIELDS_HAS_BEEN_INCLUDED 53 #define OPENVDB_TOOLS_VELOCITY_FIELDS_HAS_BEEN_INCLUDED 55 #include <tbb/parallel_reduce.h> 56 #include <openvdb/Platform.h> 57 #include <openvdb/openvdb.h> 59 #include <openvdb/math/FiniteDifference.h> 60 #include <boost/math/constants/constants.hpp> 69 template <
typename VelGr
idT,
typename Interpolator = BoxSampler>
75 BOOST_STATIC_ASSERT(boost::is_floating_point<ValueType>::value);
78 : mAccessor(vel.tree())
79 , mTransform(&vel.transform())
85 : mAccessor(other.mAccessor.tree())
86 , mTransform(other.mTransform)
99 inline VectorType operator() (
const Vec3d& xyz, ValueType)
const 101 return Interpolator::sample(mAccessor, mTransform->worldToIndex(xyz));
108 inline VectorType operator() (
const Coord& ijk, ValueType)
const 110 return mAccessor.getValue(ijk);
114 const typename VelGridT::ConstAccessor mAccessor;
127 template <
typename ScalarT =
float>
133 BOOST_STATIC_ASSERT(boost::is_floating_point<ScalarT>::value);
144 inline VectorType operator() (
const Vec3d& xyz, ValueType time)
const;
147 inline VectorType operator() (
const Coord& ijk, ValueType time)
const 149 return (*
this)(ijk.
asVec3d(), time);
153 template <
typename ScalarT>
157 const ScalarT pi = boost::math::constants::pi<ScalarT>();
158 const ScalarT phase = pi / ScalarT(3);
159 const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
160 const ScalarT tr =
math::Cos(ScalarT(time) * phase);
161 const ScalarT a =
math::Sin(ScalarT(2)*Py);
162 const ScalarT b = -
math::Sin(ScalarT(2)*Px);
163 const ScalarT c =
math::Sin(ScalarT(2)*Pz);
177 bool Staggered =
false,
188 mAcc(grid.getAccessor())
194 mAcc(mGrid->getAccessor())
204 template <
typename LocationType>
205 inline bool sample(
const LocationType& world, ValueType& result)
const 207 const Vec3R xyz = mGrid->worldToIndex(
Vec3R(world[0], world[1], world[2]));
217 template <
typename LocationType>
218 inline ValueType
sample(
const LocationType& world)
const 220 const Vec3R xyz = mGrid->worldToIndex(
Vec3R(world[0], world[1], world[2]));
239 bool Staggered =
false,
240 size_t SampleOrder = 1>
255 template<
size_t OrderRK,
typename LocationType>
256 inline void rungeKutta(
const ElementType dt, LocationType& world)
const 258 BOOST_STATIC_ASSERT(OrderRK <= 4);
259 VecType P(static_cast<ElementType>(world[0]),
260 static_cast<ElementType>(world[1]),
261 static_cast<ElementType>(world[2]));
265 }
else if (OrderRK == 1) {
267 mVelSampler.sample(P, V0);
269 }
else if (OrderRK == 2) {
271 mVelSampler.sample(P, V0);
272 mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
274 }
else if (OrderRK == 3) {
276 mVelSampler.sample(P, V0);
277 mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
278 mVelSampler.sample(P + dt * (ElementType(2.0) * V1 - V0), V2);
279 P = dt * (V0 + ElementType(4.0) * V1 + V2) * ElementType(1.0 / 6.0);
280 }
else if (OrderRK == 4) {
281 VecType V0, V1, V2, V3;
282 mVelSampler.sample(P, V0);
283 mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
284 mVelSampler.sample(P + ElementType(0.5) * dt * V1, V2);
285 mVelSampler.sample(P + dt * V2, V3);
286 P = dt * (V0 + ElementType(2.0) * (V1 + V2) + V3) * ElementType(1.0 / 6.0);
288 typedef typename LocationType::ValueType OutType;
289 world += LocationType(static_cast<OutType>(P[0]),
290 static_cast<OutType>(P[1]),
291 static_cast<OutType>(P[2]));
302 #endif // OPENVDB_TOOLS_VELOCITY_FIELDS_HAS_BEEN_INCLUDED
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:51
Vec3SGrid Vec3fGrid
Definition: openvdb.h:83
Vec3< double > Vec3d
Definition: Vec3.h:679
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:136
math::Vec3< Real > Vec3R
Definition: Types.h:79
Definition: Exceptions.h:40
Type Pow2(Type x)
Return x2.
Definition: Math.h:502
float Cos(const float &x)
Return cos x.
Definition: Math.h:679
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:188
Vec3d asVec3d() const
Definition: Coord.h:170
float Sin(const float &x)
Return sin x.
Definition: Math.h:670