OpenVDB  3.1.0
VelocityFields.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2015 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 //
30 //
32 //
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 "Interpolation.h" // for Sampler, etc.
59 #include <boost/math/constants/constants.hpp>
60 
61 namespace openvdb {
63 namespace OPENVDB_VERSION_NAME {
64 namespace tools {
65 
68 template <typename VelGridT, typename Interpolator = BoxSampler>
70 {
71 public:
72  typedef typename VelGridT::ValueType VectorType;
73  typedef typename VectorType::ValueType ValueType;
74  BOOST_STATIC_ASSERT(boost::is_floating_point<ValueType>::value);
75 
76  DiscreteField(const VelGridT &vel): mAccessor(vel.tree()), mTransform(&vel.transform())
77  {
78  }
79 
83  const math::Transform& transform() const { return *mTransform; }
84 
86  inline VectorType operator() (const Vec3d& xyz, ValueType) const
87  {
88  return Interpolator::sample(mAccessor, mTransform->worldToIndex(xyz));
89  }
90 
92  inline VectorType operator() (const Coord& ijk, ValueType) const
93  {
94  return mAccessor.getValue(ijk);
95  }
96 
97 private:
98  const typename VelGridT::ConstAccessor mAccessor;//Not thread-safe
99  const math::Transform* mTransform;
100 
101 }; // end of DiscreteField
102 
104 
111 template <typename ScalarT = float>
113 {
114 public:
115  typedef ScalarT ValueType;
117  BOOST_STATIC_ASSERT(boost::is_floating_point<ScalarT>::value);
118 
120 
125 
128  inline VectorType operator() (const Vec3d& xyz, ValueType time) const;
129 
131  inline VectorType operator() (const Coord& ijk, ValueType time) const
132  {
133  return (*this)(ijk.asVec3d(), time);
134  }
135 }; // end of EnrightField
136 
137 template <typename ScalarT>
138 inline math::Vec3<ScalarT>
140 {
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);
148  return math::Vec3<ScalarT>(
149  tr * ( ScalarT(2) * math::Pow2(sin(Px)) * a * c ),
150  tr * ( b * math::Pow2(sin(Py)) * c ),
151  tr * ( b * a * math::Pow2(sin(Pz)) ));
152 }
153 
154 
156 
160 template<typename GridT = Vec3fGrid,
161  bool Staggered = false,
162  size_t Order = 1>
164 {
165 public:
166  typedef typename GridT::ConstAccessor AccessorType;
167  typedef typename GridT::ValueType ValueType;
168 
170  VelocitySampler(const GridT& grid):
171  mGrid(&grid),
172  mAcc(grid.getAccessor())
173  {
174  }
177  mGrid(other.mGrid),
178  mAcc(mGrid->getAccessor())
179  {
180  }
184  template <typename LocationType>
185  inline bool sample(const LocationType& world, ValueType& result) const
186  {
187  const Vec3R xyz = mGrid->worldToIndex(Vec3R(world[0], world[1], world[2]));
188  bool active = Sampler<Order, Staggered>::sample(mAcc, xyz, result);
189  return active;
190  }
191 
194  template <typename LocationType>
195  inline ValueType sample(const LocationType& world) const
196  {
197  const Vec3R xyz = mGrid->worldToIndex(Vec3R(world[0], world[1], world[2]));
198  return Sampler<Order, Staggered>::sample(mAcc, xyz);
199  }
200 
201 private:
202  // holding the Grids for the transforms
203  const GridT* mGrid; // Velocity vector field
204  AccessorType mAcc;
205 };// end of VelocitySampler class
206 
208 
215 template<typename GridT = Vec3fGrid,
216  bool Staggered = false,
217  size_t SampleOrder = 1>
219 {
220 public:
221  typedef typename GridT::ValueType VecType;
222  typedef typename VecType::ValueType ElementType;
223 
224  VelocityIntegrator(const GridT& velGrid):
225  mVelSampler(velGrid)
226  {
227  }
232  template<size_t OrderRK, typename LocationType>
233  inline void rungeKutta(const ElementType dt, LocationType& world) const
234  {
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]));
239  // Note the if-branching below is optimized away at compile time
240  if (OrderRK == 0) {
241  return;// do nothing
242  } else if (OrderRK == 1) {
243  VecType V0;
244  mVelSampler.sample(P, V0);
245  P = dt * V0;
246  } else if (OrderRK == 2) {
247  VecType V0, V1;
248  mVelSampler.sample(P, V0);
249  mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
250  P = dt * V1;
251  } else if (OrderRK == 3) {
252  VecType V0, V1, V2;
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);
264  }
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]));
269  }
270 private:
272 };// end of VelocityIntegrator class
273 
274 
275 } // namespace tools
276 } // namespace OPENVDB_VERSION_NAME
277 } // namespace openvdb
278 
279 #endif // OPENVDB_TOOLS_VELOCITY_FIELDS_HAS_BEEN_INCLUDED
280 
281 // Copyright (c) 2012-2015 DreamWorks Animation LLC
282 // All rights reserved. This software is distributed under the
283 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
math::Vec3< ScalarT > VectorType
Definition: VelocityFields.h:116
DiscreteField(const VelGridT &vel)
Definition: VelocityFields.h:76
void rungeKutta(const ElementType dt, LocationType &world) const
Variable order Runge-Kutta time integration for a single time step.
Definition: VelocityFields.h:233
VelGridT::ValueType VectorType
Definition: VelocityFields.h:72
VelocitySampler(const VelocitySampler &other)
Copy-constructor.
Definition: VelocityFields.h:176
math::Vec3< Real > Vec3R
Definition: Types.h:76
VecType::ValueType ElementType
Definition: VelocityFields.h:222
Analytical, divergence-free and periodic velocity field.
Definition: VelocityFields.h:112
const math::Transform & transform() const
Definition: VelocityFields.h:83
math::Transform transform() const
Definition: VelocityFields.h:124
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:195
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Thin wrapper class for a velocity grid.
Definition: VelocityFields.h:69
Vec3SGrid Vec3fGrid
Definition: openvdb.h:77
GridT::ValueType VecType
Definition: VelocityFields.h:221
VelocityIntegrator(const GridT &velGrid)
Definition: VelocityFields.h:224
Definition: Mat.h:146
Definition: Exceptions.h:39
Vec3< double > Vec3d
Definition: Vec3.h:643
Type Pow2(Type x)
Return .
Definition: Math.h:514
Performs Runge-Kutta time integration of variable order in a static velocity field.
Definition: VelocityFields.h:218
Calculate an axis-aligned bounding box in index space from a bounding sphere in world space...
Definition: Transform.h:66
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:185
VectorType::ValueType ValueType
Definition: VelocityFields.h:73
ScalarT ValueType
Definition: VelocityFields.h:115
GridT::ConstAccessor AccessorType
Definition: VelocityFields.h:166
Definition: VelocityFields.h:163
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
VelocitySampler(const GridT &grid)
Constructor from a grid.
Definition: VelocityFields.h:170
Provises a unified interface for sampling, i.e. interpolation.
Definition: Interpolation.h:90
EnrightField()
Definition: VelocityFields.h:119
GridT::ValueType ValueType
Definition: VelocityFields.h:167