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