GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/tools/VelocityFields.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 26 34 76.5%
Functions: 2 5 40.0%
Branches: 2 530 0.4%
Line Branch Exec Source
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.
32 #include <openvdb/math/FiniteDifference.h>
33
34 namespace openvdb {
35 OPENVDB_USE_VERSION_NAMESPACE
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>
42 ✗ class DiscreteField
43 {
44 public:
45 typedef typename VelGridT::ValueType VectorType;
46 typedef typename VectorType::ValueType ValueType;
47 static_assert(std::is_floating_point<ValueType>::value,
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
57 ✗ DiscreteField(const DiscreteField& other)
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>
101 class EnrightField
102 {
103 public:
104 typedef ScalarT ValueType;
105 typedef math::Vec3<ScalarT> VectorType;
106 static_assert(std::is_floating_point<ScalarT>::value,
107 "EnrightField requires a floating point grid.");
108
109 EnrightField() {}
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
114 math::Transform transform() const { return math::Transform(); }
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>
153 class VelocitySampler
154 {
155 public:
156 typedef typename GridT::ConstAccessor AccessorType;
157 typedef typename GridT::ValueType ValueType;
158
159 /// @brief Constructor from a grid
160 1257 VelocitySampler(const GridT& grid):
161 mGrid(&grid),
162
2/48
✓ Branch 1 taken 240 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 10 taken 480 times.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
 1257 mAcc(grid.getAccessor())
163 {
164 }
165 /// @brief Copy-constructor
166 2327 VelocitySampler(const VelocitySampler& other):
167 2327 mGrid(other.mGrid),
168 2327 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 202321038 inline bool sample(const LocationType& world, ValueType& result) const
180 {
181 202321038 const Vec3R xyz = mGrid->worldToIndex(Vec3R(world[0], world[1], world[2]));
182 202321038 bool active = Sampler<Order, Staggered>::sample(mAcc, xyz, result);
183 202321038 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>
215 ✗ class VelocityIntegrator
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 30006456 inline void rungeKutta(const ElementType dt, LocationType& world) const
231 {
232 BOOST_STATIC_ASSERT(OrderRK <= 4);
233 30006456 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 4288427 mVelSampler.sample(P, V0);
242 4288427 P = dt * V0;
243 } else if (OrderRK == 2) {
244 VecType V0, V1;
245 2000008 mVelSampler.sample(P, V0);
246 2000008 mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
247 2000008 P = dt * V1;
248 } else if (OrderRK == 3) {
249 VecType V0, V1, V2;
250 2000008 mVelSampler.sample(P, V0);
251 2000008 mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
252 2000008 mVelSampler.sample(P + dt * (ElementType(2.0) * V1 - V0), V2);
253 2000008 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 21718013 mVelSampler.sample(P, V0);
257 21718013 mVelSampler.sample(P + ElementType(0.5) * dt * V0, V1);
258 21718013 mVelSampler.sample(P + ElementType(0.5) * dt * V1, V2);
259 21718013 mVelSampler.sample(P + dt * V2, V3);
260 21718013 P = dt * (V0 + ElementType(2.0) * (V1 + V2) + V3) * ElementType(1.0 / 6.0);
261 }
262 typedef typename LocationType::ValueType OutType;
263 22006440 world += LocationType(static_cast<OutType>(P[0]),
264 static_cast<OutType>(P[1]),
265 static_cast<OutType>(P[2]));
266 }
267 private:
268 VelocitySampler<GridT, Staggered, SampleOrder> mVelSampler;
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
277