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>
58 #include <openvdb/math/FiniteDifference.h>
59 #include <boost/math/constants/constants.hpp>
68 template <
typename VelGr
idT,
typename Interpolator = BoxSampler>
74 BOOST_STATIC_ASSERT(boost::is_floating_point<ValueType>::value);
76 DiscreteField(
const VelGridT &vel): mAccessor(vel.tree()), mTransform(&vel.transform())
86 inline VectorType operator() (
const Vec3d& xyz, ValueType)
const
88 return Interpolator::sample(mAccessor, mTransform->worldToIndex(xyz));
92 inline VectorType operator() (
const Coord& ijk, ValueType)
const
94 return mAccessor.getValue(ijk);
98 const typename VelGridT::ConstAccessor mAccessor;
111 template <
typename ScalarT =
float>
117 BOOST_STATIC_ASSERT(boost::is_floating_point<ScalarT>::value);
128 inline VectorType operator() (
const Vec3d& xyz, ValueType time)
const;
131 inline VectorType operator() (
const Coord& ijk, ValueType time)
const
133 return (*
this)(ijk.asVec3d(), time);
137 template <
typename ScalarT>
141 const ScalarT pi = boost::math::constants::pi<ScalarT>();
142 const ScalarT phase = pi / ScalarT(3.0);
143 const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
144 const ScalarT tr = cos(ScalarT(time) * phase);
145 const ScalarT a = sin(ScalarT(2.0)*Py);
146 const ScalarT b = -sin(ScalarT(2.0)*Px);
147 const ScalarT c = sin(ScalarT(2.0)*Pz);
149 tr * ( ScalarT(2) *
math::Pow2(sin(Px)) * a * c ),
161 bool Staggered =
false,
172 mAcc(grid.getAccessor())
178 mAcc(mGrid->getAccessor())
184 template <
typename LocationType>
185 inline bool sample(
const LocationType& world, ValueType& result)
const
187 const Vec3R xyz = mGrid->worldToIndex(
Vec3R(world[0], world[1], world[2]));
194 template <
typename LocationType>
195 inline ValueType
sample(
const LocationType& world)
const
197 const Vec3R xyz = mGrid->worldToIndex(
Vec3R(world[0], world[1], world[2]));
216 bool Staggered =
false,
217 size_t SampleOrder = 1>
232 template<
size_t OrderRK,
typename LocationType>
233 inline void rungeKutta(
const ElementType dt, LocationType& world)
const
235 BOOST_STATIC_ASSERT(OrderRK <= 4);
236 VecType P(static_cast<ElementType>(world[0]),
237 static_cast<ElementType>(world[1]),
238 static_cast<ElementType>(world[2]));
242 }
else if (OrderRK == 1) {
244 mVelSampler.sample(P, V0);
246 }
else if (OrderRK == 2) {
248 mVelSampler.sample(P, V0);
249 mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
251 }
else if (OrderRK == 3) {
253 mVelSampler.sample(P, V0);
254 mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
255 mVelSampler.sample(P + dt * (ElementType(2.0) * V1 - V0), V2);
256 P = dt * (V0 + ElementType(4.0) * V1 + V2) * ElementType(1.0 / 6.0);
257 }
else if (OrderRK == 4) {
258 VecType V0, V1, V2, V3;
259 mVelSampler.sample(P, V0);
260 mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
261 mVelSampler.sample(P + ElementType(0.5) * dt * V1, V2);
262 mVelSampler.sample(P + dt * V2, V3);
263 P = dt * (V0 + ElementType(2.0) * (V1 + V2) + V3) * ElementType(1.0 / 6.0);
265 typedef typename LocationType::ValueType OutType;
266 world += LocationType(static_cast<OutType>(P[0]),
267 static_cast<OutType>(P[1]),
268 static_cast<OutType>(P[2]));
279 #endif // OPENVDB_TOOLS_VELOCITY_FIELDS_HAS_BEEN_INCLUDED
math::Vec3< Real > Vec3R
Definition: Types.h:76
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Vec3SGrid Vec3fGrid
Definition: openvdb.h:77
Definition: Exceptions.h:39
Vec3< double > Vec3d
Definition: Vec3.h:643
Type Pow2(Type x)
Return .
Definition: Math.h:514
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71