OpenVDB  9.0.1
VelocityFields.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 ///////////////////////////////////////////////////////////////////////////
5 //
6 /// @author Ken Museth
7 ///
8 /// @file VelocityFields.h
9 ///
10 /// @brief Defines two simple wrapper classes for advection velocity
11 /// fields as well as VelocitySampler and VelocityIntegrator
12 ///
13 ///
14 /// @details DiscreteField wraps a velocity grid and EnrightField is mostly
15 /// intended for debugging (it's an analytical divergence free and
16 /// periodic field). They both share the same API required by the
17 /// LevelSetAdvection class defined in LevelSetAdvect.h. Thus, any
18 /// class with this API should work with LevelSetAdvection.
19 ///
20 /// @warning Note the Field wrapper classes below always assume the velocity
21 /// is represented in the world-frame of reference. For DiscreteField
22 /// this implies the input grid must contain velocities in world
23 /// coordinates.
24 
25 #ifndef OPENVDB_TOOLS_VELOCITY_FIELDS_HAS_BEEN_INCLUDED
26 #define OPENVDB_TOOLS_VELOCITY_FIELDS_HAS_BEEN_INCLUDED
27 
28 #include <tbb/parallel_reduce.h>
29 #include <openvdb/Platform.h>
30 #include <openvdb/openvdb.h>
31 #include "Interpolation.h" // for Sampler, etc.
33 
34 namespace openvdb {
36 namespace OPENVDB_VERSION_NAME {
37 namespace tools {
38 
39 /// @brief Thin wrapper class for a velocity grid
40 /// @note Consider replacing BoxSampler with StaggeredBoxSampler
41 template <typename VelGridT, typename Interpolator = BoxSampler>
43 {
44 public:
45  typedef typename VelGridT::ValueType VectorType;
46  typedef typename VectorType::ValueType ValueType;
48  "DiscreteField requires a floating point grid.");
49 
50  DiscreteField(const VelGridT &vel)
51  : mAccessor(vel.tree())
52  , mTransform(&vel.transform())
53  {
54  }
55 
56  /// @brief Copy constructor
58  : mAccessor(other.mAccessor.tree())
59  , mTransform(other.mTransform)
60  {
61  }
62 
63  /// @return const reference to the transform between world and index space
64  /// @note Use this method to determine if a client grid is
65  /// aligned with the coordinate space of the velocity grid.
66  const math::Transform& transform() const { return *mTransform; }
67 
68  /// @return the interpolated velocity at the world space position xyz
69  ///
70  /// @warning Not threadsafe since it uses a ValueAccessor! So use
71  /// one instance per thread (which is fine since its lightweight).
72  inline VectorType operator() (const Vec3d& xyz, ValueType/*dummy time*/) const
73  {
74  return Interpolator::sample(mAccessor, mTransform->worldToIndex(xyz));
75  }
76 
77  /// @return the velocity at the coordinate space position ijk
78  ///
79  /// @warning Not threadsafe since it uses a ValueAccessor! So use
80  /// one instance per thread (which is fine since its lightweight).
81  inline VectorType operator() (const Coord& ijk, ValueType/*dummy time*/) const
82  {
83  return mAccessor.getValue(ijk);
84  }
85 
86 private:
87  const typename VelGridT::ConstAccessor mAccessor;//Not thread-safe
88  const math::Transform* mTransform;
89 
90 }; // end of DiscreteField
91 
92 ///////////////////////////////////////////////////////////////////////
93 
94 /// @brief Analytical, divergence-free and periodic velocity field
95 /// @note Primarily intended for debugging!
96 /// @warning This analytical velocity only produce meaningful values
97 /// in the unit box in world space. In other words make sure any level
98 /// set surface is fully enclosed in the axis aligned bounding box
99 /// spanning 0->1 in world units.
100 template <typename ScalarT = float>
102 {
103 public:
104  typedef ScalarT ValueType;
107  "EnrightField requires a floating point grid.");
108 
110 
111  /// @return const reference to the identity transform between world and index space
112  /// @note Use this method to determine if a client grid is
113  /// aligned with the coordinate space of this velocity field
115 
116  /// @return the velocity in world units, evaluated at the world
117  /// position xyz and at the specified time
118  inline VectorType operator() (const Vec3d& xyz, ValueType time) const;
119 
120  /// @return the velocity at the coordinate space position ijk
121  inline VectorType operator() (const Coord& ijk, ValueType time) const
122  {
123  return (*this)(ijk.asVec3d(), time);
124  }
125 }; // end of EnrightField
126 
127 template <typename ScalarT>
128 inline math::Vec3<ScalarT>
129 EnrightField<ScalarT>::operator() (const Vec3d& xyz, ValueType time) const
130 {
131  const ScalarT pi = math::pi<ScalarT>();
132  const ScalarT phase = pi / ScalarT(3);
133  const ScalarT Px = pi * ScalarT(xyz[0]), Py = pi * ScalarT(xyz[1]), Pz = pi * ScalarT(xyz[2]);
134  const ScalarT tr = math::Cos(ScalarT(time) * phase);
135  const ScalarT a = math::Sin(ScalarT(2)*Py);
136  const ScalarT b = -math::Sin(ScalarT(2)*Px);
137  const ScalarT c = math::Sin(ScalarT(2)*Pz);
138  return math::Vec3<ScalarT>(
139  tr * ( ScalarT(2) * math::Pow2(math::Sin(Px)) * a * c ),
140  tr * ( b * math::Pow2(math::Sin(Py)) * c ),
141  tr * ( b * a * math::Pow2(math::Sin(Pz)) ));
142 }
143 
144 
145 ///////////////////////////////////////////////////////////////////////
146 
147 /// Class to hold a Vec3 field interpreted as a velocity field.
148 /// Primarily exists to provide a method(s) that integrate a passive
149 /// point forward in the velocity field for a single time-step (dt)
150 template<typename GridT = Vec3fGrid,
151  bool Staggered = false,
152  size_t Order = 1>
154 {
155 public:
156  typedef typename GridT::ConstAccessor AccessorType;
157  typedef typename GridT::ValueType ValueType;
158 
159  /// @brief Constructor from a grid
160  VelocitySampler(const GridT& grid):
161  mGrid(&grid),
162  mAcc(grid.getAccessor())
163  {
164  }
165  /// @brief Copy-constructor
167  mGrid(other.mGrid),
168  mAcc(mGrid->getAccessor())
169  {
170  }
171  /// @brief Samples the velocity at world position onto result. Supports both
172  /// staggered (i.e. MAC) and collocated velocity grids.
173  ///
174  /// @return @c true if any one of the sampled values is active.
175  ///
176  /// @warning Not threadsafe since it uses a ValueAccessor! So use
177  /// one instance per thread (which is fine since its lightweight).
178  template <typename LocationType>
179  inline bool sample(const LocationType& world, ValueType& result) const
180  {
181  const Vec3R xyz = mGrid->worldToIndex(Vec3R(world[0], world[1], world[2]));
182  bool active = Sampler<Order, Staggered>::sample(mAcc, xyz, result);
183  return active;
184  }
185 
186  /// @brief Samples the velocity at world position onto result. Supports both
187  /// staggered (i.e. MAC) and co-located velocity grids.
188  ///
189  /// @warning Not threadsafe since it uses a ValueAccessor! So use
190  /// one instance per thread (which is fine since its lightweight).
191  template <typename LocationType>
192  inline ValueType sample(const LocationType& world) const
193  {
194  const Vec3R xyz = mGrid->worldToIndex(Vec3R(world[0], world[1], world[2]));
195  return Sampler<Order, Staggered>::sample(mAcc, xyz);
196  }
197 
198 private:
199  // holding the Grids for the transforms
200  const GridT* mGrid; // Velocity vector field
201  AccessorType mAcc;
202 };// end of VelocitySampler class
203 
204 ///////////////////////////////////////////////////////////////////////
205 
206 /// @brief Performs Runge-Kutta time integration of variable order in
207 /// a static velocity field.
208 ///
209 /// @note Note that the order of the velocity sampling is controlled
210 /// with the SampleOrder template parameter, which defaults
211 /// to one, i.e. a tri-linear interpolation kernel.
212 template<typename GridT = Vec3fGrid,
213  bool Staggered = false,
214  size_t SampleOrder = 1>
216 {
217 public:
218  typedef typename GridT::ValueType VecType;
219  typedef typename VecType::ValueType ElementType;
220 
221  VelocityIntegrator(const GridT& velGrid):
222  mVelSampler(velGrid)
223  {
224  }
225  /// @brief Variable order Runge-Kutta time integration for a single time step
226  ///
227  /// @param dt Time sub-step for the Runge-Kutte integrator of order OrderRK
228  /// @param world Location in world space coordinates (both input and output)
229  template<size_t OrderRK, typename LocationType>
230  inline void rungeKutta(const ElementType dt, LocationType& world) const
231  {
232  BOOST_STATIC_ASSERT(OrderRK <= 4);
233  VecType P(static_cast<ElementType>(world[0]),
234  static_cast<ElementType>(world[1]),
235  static_cast<ElementType>(world[2]));
236  // Note the if-branching below is optimized away at compile time
237  if (OrderRK == 0) {
238  return;// do nothing
239  } else if (OrderRK == 1) {
240  VecType V0;
241  mVelSampler.sample(P, V0);
242  P = dt * V0;
243  } else if (OrderRK == 2) {
244  VecType V0, V1;
245  mVelSampler.sample(P, V0);
246  mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
247  P = dt * V1;
248  } else if (OrderRK == 3) {
249  VecType V0, V1, V2;
250  mVelSampler.sample(P, V0);
251  mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
252  mVelSampler.sample(P + dt * (ElementType(2.0) * V1 - V0), V2);
253  P = dt * (V0 + ElementType(4.0) * V1 + V2) * ElementType(1.0 / 6.0);
254  } else if (OrderRK == 4) {
255  VecType V0, V1, V2, V3;
256  mVelSampler.sample(P, V0);
257  mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
258  mVelSampler.sample(P + ElementType(0.5) * dt * V1, V2);
259  mVelSampler.sample(P + dt * V2, V3);
260  P = dt * (V0 + ElementType(2.0) * (V1 + V2) + V3) * ElementType(1.0 / 6.0);
261  }
262  typedef typename LocationType::ValueType OutType;
263  world += LocationType(static_cast<OutType>(P[0]),
264  static_cast<OutType>(P[1]),
265  static_cast<OutType>(P[2]));
266  }
267 private:
269 };// end of VelocityIntegrator class
270 
271 
272 } // namespace tools
273 } // namespace OPENVDB_VERSION_NAME
274 } // namespace openvdb
275 
276 #endif // OPENVDB_TOOLS_VELOCITY_FIELDS_HAS_BEEN_INCLUDED
GridT::ValueType ValueType
Definition: VelocityFields.h:157
math::Transform transform() const
Definition: VelocityFields.h:114
Definition: VelocityFields.h:153
const math::Transform & transform() const
Definition: VelocityFields.h:66
void rungeKutta(const ElementType dt, LocationType &world) const
Variable order Runge-Kutta time integration for a single time step.
Definition: VelocityFields.h:230
VectorType::ValueType ValueType
Definition: VelocityFields.h:46
Definition: Transform.h:39
Vec3SGrid Vec3fGrid
Definition: openvdb.h:54
VelocityIntegrator(const GridT &velGrid)
Definition: VelocityFields.h:221
DiscreteField(const VelGridT &vel)
Definition: VelocityFields.h:50
GridT::ValueType VecType
Definition: VelocityFields.h:218
float Cos(const float &x)
Return cos x.
Definition: Math.h:728
math::Vec3< Real > Vec3R
Definition: Types.h:72
Analytical, divergence-free and periodic velocity field.
Definition: VelocityFields.h:101
Definition: Exceptions.h:13
ScalarT ValueType
Definition: VelocityFields.h:104
ValueT value
Definition: GridBuilder.h:1287
math::Vec3< ScalarT > VectorType
Definition: VelocityFields.h:105
Performs Runge-Kutta time integration of variable order in a static velocity field.
Definition: VelocityFields.h:215
Definition: Mat.h:187
constexpr T pi()
Pi constant taken from Boost to match old behaviour.
Definition: Math.h:118
VecType::ValueType ElementType
Definition: VelocityFields.h:219
float Sin(const float &x)
Return sin x.
Definition: Math.h:719
ValueType sample(const LocationType &world) const
Samples the velocity at world position onto result. Supports both staggered (i.e. MAC) and co-located...
Definition: VelocityFields.h:192
Provises a unified interface for sampling, i.e. interpolation.
Definition: Interpolation.h:63
DiscreteField(const DiscreteField &other)
Copy constructor.
Definition: VelocityFields.h:57
VelocitySampler(const GridT &grid)
Constructor from a grid.
Definition: VelocityFields.h:160
GridT::ConstAccessor AccessorType
Definition: VelocityFields.h:156
VelGridT::ValueType VectorType
Definition: VelocityFields.h:45
VelocitySampler(const VelocitySampler &other)
Copy-constructor.
Definition: VelocityFields.h:166
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
bool sample(const LocationType &world, ValueType &result) const
Samples the velocity at world position onto result. Supports both staggered (i.e. MAC) and collocated...
Definition: VelocityFields.h:179
EnrightField()
Definition: VelocityFields.h:109
Thin wrapper class for a velocity grid.
Definition: VelocityFields.h:42
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:202