1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3
4 /// @author Dan Bailey
5 ///
6 /// @file points/PointAdvect.h
7 ///
8 /// @brief Ability to advect VDB Points through a velocity field.
9
10 #ifndef OPENVDB_POINTS_POINT_ADVECT_HAS_BEEN_INCLUDED
11 #define OPENVDB_POINTS_POINT_ADVECT_HAS_BEEN_INCLUDED
12
13 #include <openvdb/openvdb.h>
14 #include <openvdb/tools/Prune.h>
15 #include <openvdb/tools/VelocityFields.h>
16
17 #include <openvdb/points/AttributeGroup.h>
18 #include <openvdb/points/PointDataGrid.h>
19 #include <openvdb/points/PointGroup.h>
20 #include <openvdb/points/PointMove.h>
21
22 #include <memory>
23
24
25 namespace openvdb {
26 OPENVDB_USE_VERSION_NAMESPACE
27 namespace OPENVDB_VERSION_NAME {
28 namespace points {
29
30
31 /// @brief Advect points in a PointDataGrid through a velocity grid
32 /// @param points the PointDataGrid containing the points to be advected.
33 /// @param velocity a velocity grid to be sampled.
34 /// @param integrationOrder the integration scheme to use (1 is forward euler, 4 is runge-kutta 4th)
35 /// @param dt delta time.
36 /// @param timeSteps number of advection steps to perform.
37 /// @param advectFilter an optional advection index filter (moves a subset of the points)
38 /// @param filter an optional index filter (deletes a subset of the points)
39 /// @param cached caches velocity interpolation for faster performance, disable to use
40 /// less memory (default is on).
41 template <typename PointDataGridT, typename VelGridT,
42 typename AdvectFilterT = NullFilter, typename FilterT = NullFilter>
43 inline void advectPoints(PointDataGridT& points, const VelGridT& velocity,
44 const Index integrationOrder, const double dt, const Index timeSteps,
45 const AdvectFilterT& advectFilter = NullFilter(),
46 const FilterT& filter = NullFilter(),
47 const bool cached = true);
48
49
50 ////////////////////////////////////////
51
52 /// @cond OPENVDB_DOCS_INTERNAL
53
54 namespace point_advect_internal {
55
56 enum IntegrationOrder {
57 INTEGRATION_ORDER_FWD_EULER = 1,
58 INTEGRATION_ORDER_RK_2ND,
59 INTEGRATION_ORDER_RK_3RD,
60 INTEGRATION_ORDER_RK_4TH
61 };
62
63 template <typename VelGridT, Index IntegrationOrder, bool Staggered, typename FilterT>
64 class AdvectionDeformer
65 {
66 public:
67 using IntegratorT = openvdb::tools::VelocityIntegrator<VelGridT, Staggered>;
68
AdvectionDeformer(const VelGridT & velocityGrid,const double timeStep,const int steps,const FilterT & filter)69 AdvectionDeformer(const VelGridT& velocityGrid, const double timeStep, const int steps,
70 const FilterT& filter)
71 : mIntegrator(velocityGrid)
72 , mTimeStep(timeStep)
73 , mSteps(steps)
74 , mFilter(filter) { }
75
76 template <typename LeafT>
reset(const LeafT & leaf,size_t)77 void reset(const LeafT& leaf, size_t /*idx*/)
78 {
79 mFilter.reset(leaf);
80 }
81
82 template <typename IndexIterT>
apply(Vec3d & position,const IndexIterT & iter)83 void apply(Vec3d& position, const IndexIterT& iter) const
84 {
85 if (mFilter.valid(iter)) {
86 for (int n = 0; n < mSteps; ++n) {
87 mIntegrator.template rungeKutta<IntegrationOrder, openvdb::Vec3d>(
88 static_cast<typename IntegratorT::ElementType>(mTimeStep), position);
89 }
90 }
91 }
92
93 private:
94 IntegratorT mIntegrator;
95 double mTimeStep;
96 const int mSteps;
97 FilterT mFilter;
98 }; // class AdvectionDeformer
99
100
101 template <typename PointDataGridT, typename VelGridT, typename AdvectFilterT, typename FilterT>
102 struct AdvectionOp
103 {
104 using CachedDeformerT = CachedDeformer<double>;
105
AdvectionOpAdvectionOp106 AdvectionOp(PointDataGridT& points, const VelGridT& velocity,
107 const Index integrationOrder, const double timeStep, const Index steps,
108 const AdvectFilterT& advectFilter,
109 const FilterT& filter)
110 : mPoints(points)
111 , mVelocity(velocity)
112 , mIntegrationOrder(integrationOrder)
113 , mTimeStep(timeStep)
114 , mSteps(steps)
115 , mAdvectFilter(advectFilter)
116 , mFilter(filter) { }
117
cacheAdvectionOp118 void cache()
119 {
120 mCachedDeformer.reset(new CachedDeformerT(mCache));
121 (*this)(true);
122 }
123
advectAdvectionOp124 void advect()
125 {
126 (*this)(false);
127 }
128
129 private:
130 template <int IntegrationOrder, bool Staggered>
resolveIntegrationOrderAdvectionOp131 void resolveIntegrationOrder(bool buildCache)
132 {
133 const auto leaf = mPoints.constTree().cbeginLeaf();
134 if (!leaf) return;
135
136 // move points according to the pre-computed cache
137 if (!buildCache && mCachedDeformer) {
138 movePoints(mPoints, *mCachedDeformer, mFilter);
139 return;
140 }
141
142 NullFilter nullFilter;
143
144 if (buildCache) {
145 // disable group filtering from the advection deformer and perform group filtering
146 // in the cache deformer instead, this restricts the cache to just containing
147 // positions from points which are both deforming *and* are not being deleted
148 AdvectionDeformer<VelGridT, IntegrationOrder, Staggered, NullFilter> deformer(
149 mVelocity, mTimeStep, mSteps, nullFilter);
150 if (mFilter.state() == index::ALL && mAdvectFilter.state() == index::ALL) {
151 mCachedDeformer->evaluate(mPoints, deformer, nullFilter);
152 } else {
153 BinaryFilter<AdvectFilterT, FilterT, /*And=*/true> binaryFilter(
154 mAdvectFilter, mFilter);
155 mCachedDeformer->evaluate(mPoints, deformer, binaryFilter);
156 }
157 }
158 else {
159 // revert to NullFilter if all points are being evaluated
160 if (mAdvectFilter.state() == index::ALL) {
161 AdvectionDeformer<VelGridT, IntegrationOrder, Staggered, NullFilter> deformer(
162 mVelocity, mTimeStep, mSteps, nullFilter);
163 movePoints(mPoints, deformer, mFilter);
164 }
165 else {
166 AdvectionDeformer<VelGridT, IntegrationOrder, Staggered, AdvectFilterT> deformer(
167 mVelocity, mTimeStep, mSteps, mAdvectFilter);
168 movePoints(mPoints, deformer, mFilter);
169 }
170 }
171 }
172
173 template <bool Staggered>
resolveStaggeredAdvectionOp174 void resolveStaggered(bool buildCache)
175 {
176 if (mIntegrationOrder == INTEGRATION_ORDER_FWD_EULER) {
177 resolveIntegrationOrder<1, Staggered>(buildCache);
178 } else if (mIntegrationOrder == INTEGRATION_ORDER_RK_2ND) {
179 resolveIntegrationOrder<2, Staggered>(buildCache);
180 } else if (mIntegrationOrder == INTEGRATION_ORDER_RK_3RD) {
181 resolveIntegrationOrder<3, Staggered>(buildCache);
182 } else if (mIntegrationOrder == INTEGRATION_ORDER_RK_4TH) {
183 resolveIntegrationOrder<4, Staggered>(buildCache);
184 }
185 }
186
operatorAdvectionOp187 void operator()(bool buildCache)
188 {
189 // early-exit if no leafs
190 if (mPoints.constTree().leafCount() == 0) return;
191
192 if (mVelocity.getGridClass() == openvdb::GRID_STAGGERED) {
193 resolveStaggered<true>(buildCache);
194 } else {
195 resolveStaggered<false>(buildCache);
196 }
197 }
198
199 PointDataGridT& mPoints;
200 const VelGridT& mVelocity;
201 const Index mIntegrationOrder;
202 const double mTimeStep;
203 const Index mSteps;
204 const AdvectFilterT& mAdvectFilter;
205 const FilterT& mFilter;
206 CachedDeformerT::Cache mCache;
207 std::unique_ptr<CachedDeformerT> mCachedDeformer;
208 }; // struct AdvectionOp
209
210 } // namespace point_advect_internal
211
212 /// @endcond
213
214 ////////////////////////////////////////
215
216
217 template <typename PointDataGridT, typename VelGridT, typename AdvectFilterT, typename FilterT>
advectPoints(PointDataGridT & points,const VelGridT & velocity,const Index integrationOrder,const double timeStep,const Index steps,const AdvectFilterT & advectFilter,const FilterT & filter,const bool cached)218 inline void advectPoints(PointDataGridT& points, const VelGridT& velocity,
219 const Index integrationOrder, const double timeStep, const Index steps,
220 const AdvectFilterT& advectFilter,
221 const FilterT& filter,
222 const bool cached)
223 {
224 using namespace point_advect_internal;
225
226 if (steps == 0) return;
227
228 if (integrationOrder > 4) {
229 throw ValueError{"Unknown integration order for advecting points."};
230 }
231
232 AdvectionOp<PointDataGridT, VelGridT, AdvectFilterT, FilterT> op(
233 points, velocity, integrationOrder, timeStep, steps,
234 advectFilter, filter);
235
236 // if caching is enabled, sample the velocity field using a CachedDeformer to store the
237 // intermediate positions before moving the points, this uses more memory but typically
238 // results in faster overall performance
239 if (cached) op.cache();
240
241 // advect the points
242 op.advect();
243 }
244
245 } // namespace points
246 } // namespace OPENVDB_VERSION_NAME
247 } // namespace openvdb
248
249 #endif // OPENVDB_POINTS_POINT_ADVECT_HAS_BEEN_INCLUDED
250