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