HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VelocityFields.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 VelocityFields.h
36 ///
37 /// @brief Defines two simple wrapper classes for advection velocity
38 /// fields as well as VelocitySampler and VelocityIntegrator
39 ///
40 ///
41 /// @details DiscreteField wraps a velocity grid and EnrightField is mostly
42 /// intended for debugging (it's an analytical divergence free and
43 /// periodic field). They both share the same API required by the
44 /// LevelSetAdvection class defined in LevelSetAdvect.h. Thus, any
45 /// class with this API should work with LevelSetAdvection.
46 ///
47 /// @warning Note the Field wrapper classes below always assume the velocity
48 /// is represented in the world-frame of reference. For DiscreteField
49 /// this implies the input grid must contain velocities in world
50 /// coordinates.
51 
52 #ifndef OPENVDB_TOOLS_VELOCITY_FIELDS_HAS_BEEN_INCLUDED
53 #define OPENVDB_TOOLS_VELOCITY_FIELDS_HAS_BEEN_INCLUDED
54 
55 #include <tbb/parallel_reduce.h>
56 #include <openvdb/Platform.h>
57 #include <openvdb/openvdb.h>
58 #include "Interpolation.h" // for Sampler, etc.
60 #include <hboost/math/constants/constants.hpp>
61 
62 namespace openvdb {
64 namespace OPENVDB_VERSION_NAME {
65 namespace tools {
66 
67 /// @brief Thin wrapper class for a velocity grid
68 /// @note Consider replacing BoxSampler with StaggeredBoxSampler
69 template <typename VelGridT, typename Interpolator = BoxSampler>
71 {
72 public:
73  typedef typename VelGridT::ValueType VectorType;
74  typedef typename VectorType::ValueType ValueType;
76 
77  DiscreteField(const VelGridT &vel)
78  : mAccessor(vel.tree())
79  , mTransform(&vel.transform())
80  {
81  }
82 
83  /// @brief Copy constructor
85  : mAccessor(other.mAccessor.tree())
86  , mTransform(other.mTransform)
87  {
88  }
89 
90  /// @return const reference to the transform between world and index space
91  /// @note Use this method to determine if a client grid is
92  /// aligned with the coordinate space of the velocity grid.
93  const math::Transform& transform() const { return *mTransform; }
94 
95  /// @return the interpolated velocity at the world space position xyz
96  ///
97  /// @warning Not threadsafe since it uses a ValueAccessor! So use
98  /// one instance per thread (which is fine since its lightweight).
99  inline VectorType operator() (const Vec3d& xyz, ValueType/*dummy time*/) const
100  {
101  return Interpolator::sample(mAccessor, mTransform->worldToIndex(xyz));
102  }
103 
104  /// @return the velocity at the coordinate space position ijk
105  ///
106  /// @warning Not threadsafe since it uses a ValueAccessor! So use
107  /// one instance per thread (which is fine since its lightweight).
108  inline VectorType operator() (const Coord& ijk, ValueType/*dummy time*/) const
109  {
110  return mAccessor.getValue(ijk);
111  }
112 
113 private:
114  const typename VelGridT::ConstAccessor mAccessor;//Not thread-safe
115  const math::Transform* mTransform;
116 
117 }; // end of DiscreteField
118 
119 ///////////////////////////////////////////////////////////////////////
120 
121 /// @brief Analytical, divergence-free and periodic velocity field
122 /// @note Primarily intended for debugging!
123 /// @warning This analytical velocity only produce meaningful values
124 /// in the unit box in world space. In other words make sure any level
125 /// set surface is fully enclosed in the axis aligned bounding box
126 /// spanning 0->1 in world units.
127 template <typename ScalarT = float>
129 {
130 public:
131  typedef ScalarT ValueType;
134 
136 
137  /// @return const reference to the identity transform between world and index space
138  /// @note Use this method to determine if a client grid is
139  /// aligned with the coordinate space of this velocity field
141 
142  /// @return the velocity in world units, evaluated at the world
143  /// position xyz and at the specified time
144  inline VectorType operator() (const Vec3d& xyz, ValueType time) const;
145 
146  /// @return the velocity at the coordinate space position ijk
147  inline VectorType operator() (const Coord& ijk, ValueType time) const
148  {
149  return (*this)(ijk.asVec3d(), time);
150  }
151 }; // end of EnrightField
152 
153 template <typename ScalarT>
154 inline math::Vec3<ScalarT>
156 {
157  const ScalarT pi = hboost::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);
164  return math::Vec3<ScalarT>(
165  tr * ( ScalarT(2) * math::Pow2(math::Sin(Px)) * a * c ),
166  tr * ( b * math::Pow2(math::Sin(Py)) * c ),
167  tr * ( b * a * math::Pow2(math::Sin(Pz)) ));
168 }
169 
170 
171 ///////////////////////////////////////////////////////////////////////
172 
173 /// Class to hold a Vec3 field interpreted as a velocity field.
174 /// Primarily exists to provide a method(s) that integrate a passive
175 /// point forward in the velocity field for a single time-step (dt)
176 template<typename GridT = Vec3fGrid,
177  bool Staggered = false,
178  size_t Order = 1>
180 {
181 public:
182  typedef typename GridT::ConstAccessor AccessorType;
183  typedef typename GridT::ValueType ValueType;
184 
185  /// @brief Constructor from a grid
186  VelocitySampler(const GridT& grid):
187  mGrid(&grid),
188  mAcc(grid.getAccessor())
189  {
190  }
191  /// @brief Copy-constructor
193  mGrid(other.mGrid),
194  mAcc(mGrid->getAccessor())
195  {
196  }
197  /// @brief Samples the velocity at world position onto result. Supports both
198  /// staggered (i.e. MAC) and collocated velocity grids.
199  ///
200  /// @return @c true if any one of the sampled values is active.
201  ///
202  /// @warning Not threadsafe since it uses a ValueAccessor! So use
203  /// one instance per thread (which is fine since its lightweight).
204  template <typename LocationType>
205  inline bool sample(const LocationType& world, ValueType& result) const
206  {
207  const Vec3R xyz = mGrid->worldToIndex(Vec3R(world[0], world[1], world[2]));
208  bool active = Sampler<Order, Staggered>::sample(mAcc, xyz, result);
209  return active;
210  }
211 
212  /// @brief Samples the velocity at world position onto result. Supports both
213  /// staggered (i.e. MAC) and co-located velocity grids.
214  ///
215  /// @warning Not threadsafe since it uses a ValueAccessor! So use
216  /// one instance per thread (which is fine since its lightweight).
217  template <typename LocationType>
218  inline ValueType sample(const LocationType& world) const
219  {
220  const Vec3R xyz = mGrid->worldToIndex(Vec3R(world[0], world[1], world[2]));
221  return Sampler<Order, Staggered>::sample(mAcc, xyz);
222  }
223 
224 private:
225  // holding the Grids for the transforms
226  const GridT* mGrid; // Velocity vector field
227  AccessorType mAcc;
228 };// end of VelocitySampler class
229 
230 ///////////////////////////////////////////////////////////////////////
231 
232 /// @brief Performs Runge-Kutta time integration of variable order in
233 /// a static velocity field.
234 ///
235 /// @note Note that the order of the velocity sampling is controlled
236 /// with the SampleOrder template parameter, which defaults
237 /// to one, i.e. a tri-linear interpolation kernel.
238 template<typename GridT = Vec3fGrid,
239  bool Staggered = false,
240  size_t SampleOrder = 1>
242 {
243 public:
244  typedef typename GridT::ValueType VecType;
245  typedef typename VecType::ValueType ElementType;
246 
247  VelocityIntegrator(const GridT& velGrid):
248  mVelSampler(velGrid)
249  {
250  }
251  /// @brief Variable order Runge-Kutta time integration for a single time step
252  ///
253  /// @param dt Time sub-step for the Runge-Kutte integrator of order OrderRK
254  /// @param world Location in world space coordinates (both input and output)
255  template<size_t OrderRK, typename LocationType>
256  inline void rungeKutta(const ElementType dt, LocationType& world) const
257  {
258  HBOOST_STATIC_ASSERT(OrderRK <= 4);
259  VecType P(static_cast<ElementType>(world[0]),
260  static_cast<ElementType>(world[1]),
261  static_cast<ElementType>(world[2]));
262  // Note the if-branching below is optimized away at compile time
263  if (OrderRK == 0) {
264  return;// do nothing
265  } else if (OrderRK == 1) {
266  VecType V0;
267  mVelSampler.sample(P, V0);
268  P = dt * V0;
269  } else if (OrderRK == 2) {
270  VecType V0, V1;
271  mVelSampler.sample(P, V0);
272  mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
273  P = dt * V1;
274  } else if (OrderRK == 3) {
275  VecType V0, V1, V2;
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);
287  }
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]));
292  }
293 private:
295 };// end of VelocityIntegrator class
296 
297 
298 } // namespace tools
299 } // namespace OPENVDB_VERSION_NAME
300 } // namespace openvdb
301 
302 #endif // OPENVDB_TOOLS_VELOCITY_FIELDS_HAS_BEEN_INCLUDED
303 
304 // Copyright (c) 2012-2018 DreamWorks Animation LLC
305 // All rights reserved. This software is distributed under the
306 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
void rungeKutta(const ElementType dt, LocationType &world) const
Variable order Runge-Kutta time integration for a single time step.
VectorType operator()(const Vec3d &xyz, ValueType) const
HBOOST_STATIC_ASSERT(hboost::is_floating_point< ScalarT >::value)
DiscreteField(const DiscreteField &other)
Copy constructor.
math::Vec3< Real > Vec3R
Definition: Types.h:79
Type Pow2(Type x)
Return x2.
Definition: Math.h:502
GT_API const UT_StringHolder time
static bool sample(const TreeT &inTree, const Vec3R &inCoord, typename TreeT::ValueType &result)
Sample inTree at the floating-point index coordinate inCoord and store the result in result...
Thin wrapper class for a velocity grid.
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1221
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:189
float Cos(const float &x)
Return cos x.
Definition: Math.h:679
VelocitySampler(const GridT &grid)
Constructor from a grid.
VectorType operator()(const Vec3d &xyz, ValueType time) const
Performs Runge-Kutta time integration of variable order in a static velocity field.
HBOOST_STATIC_ASSERT(hboost::is_floating_point< ValueType >::value)
VelocitySampler(const VelocitySampler &other)
Copy-constructor.
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1221
ValueType sample(const LocationType &world) const
Samples the velocity at world position onto result. Supports both staggered (i.e. MAC) and co-located...
GLsizei const GLfloat * value
Definition: glcorearb.h:823
float Sin(const float &x)
Return sin x.
Definition: Math.h:670
Vec3d worldToIndex(const Vec3d &xyz) const
Apply this transformation to the given coordinates.
Definition: Transform.h:137
bool sample(const LocationType &world, ValueType &result) const
Samples the velocity at world position onto result. Supports both staggered (i.e. MAC) and collocated...
Analytical, divergence-free and periodic velocity field.
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:135
const math::Transform & transform() const