GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/tools/PotentialFlow.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 88 90 97.8%
Functions: 13 18 72.2%
Branches: 97 208 46.6%
Line Branch Exec Source
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3
4 /// @file tools/PotentialFlow.h
5 ///
6 /// @brief Tools for creating potential flow fields through solving Laplace's equation
7 ///
8 /// @authors Todd Keeler, Dan Bailey
9
10 #ifndef OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED
11 #define OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED
12
13 #include <openvdb/openvdb.h>
14
15 #include "GridOperators.h"
16 #include "GridTransformer.h"
17 #include "Mask.h" // interiorMask
18 #include "Morphology.h" // erodeActiveValues
19 #include "PoissonSolver.h"
20 #include <openvdb/openvdb.h>
21
22
23 namespace openvdb {
24 OPENVDB_USE_VERSION_NAMESPACE
25 namespace OPENVDB_VERSION_NAME {
26 namespace tools {
27
28 /// @brief Metafunction to convert a vector-valued grid type to a scalar grid type
29 template<typename VecGridT>
30 struct VectorToScalarGrid {
31 using Type =
32 typename VecGridT::template ValueConverter<typename VecGridT::ValueType::value_type>::Type;
33 using Ptr = typename Type::Ptr;
34 using ConstPtr = typename Type::ConstPtr;
35 };
36
37
38 /// @brief Construct a mask for the Potential Flow domain.
39 /// @details For a level set, this represents a rebuilt exterior narrow band.
40 /// For any other grid it is a new region that surrounds the active voxels.
41 /// @param grid source grid to use for computing the mask
42 /// @param dilation dilation in voxels of the source grid to form the new potential flow mask
43 template<typename GridT, typename MaskT = typename GridT::template ValueConverter<ValueMask>::Type>
44 typename MaskT::Ptr
45 createPotentialFlowMask(const GridT& grid, int dilation = 5);
46
47
48 /// @brief Create a Potential Flow velocities grid for the Neumann boundary.
49 /// @param collider a level set that represents the boundary
50 /// @param domain a mask to represent the potential flow domain
51 /// @param boundaryVelocity an optional grid pointer to stores the velocities of the boundary
52 /// @param backgroundVelocity a background velocity value
53 /// @details Typically this method involves supplying a velocity grid for the
54 /// collider boundary, however it can also be used for a global wind field
55 /// around the collider by supplying an empty boundary Velocity and a
56 /// non-zero background velocity.
57 template<typename Vec3T, typename GridT, typename MaskT>
58 typename GridT::template ValueConverter<Vec3T>::Type::Ptr
59 createPotentialFlowNeumannVelocities(const GridT& collider, const MaskT& domain,
60 const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
61 const Vec3T& backgroundVelocity);
62
63
64 /// @brief Compute the Potential on the domain using the Neumann boundary conditions on
65 /// solid boundaries
66 /// @param domain a mask to represent the domain in which to perform the solve
67 /// @param neumann the topology of this grid defines where the solid boundaries are and grid
68 /// values give the Neumann boundaries that should be applied there
69 /// @param state the solver parameters for computing the solution
70 /// @param interrupter pointer to an optional interrupter adhering to the
71 /// util::NullInterrupter interface
72 /// @details On input, the State object should specify convergence criteria
73 /// (minimum error and maximum number of iterations); on output, it gives
74 /// the actual termination conditions.
75 template<typename Vec3GridT, typename MaskT, typename InterrupterT = util::NullInterrupter>
76 typename VectorToScalarGrid<Vec3GridT>::Ptr
77 computeScalarPotential(const MaskT& domain, const Vec3GridT& neumann, math::pcg::State& state,
78 InterrupterT* interrupter = nullptr);
79
80
81 /// @brief Compute a vector Flow Field comprising the gradient of the potential with Neumann
82 /// boundary conditions applied
83 /// @param potential scalar potential, typically computed from computeScalarPotential()
84 /// @param neumann the topology of this grid defines where the solid boundaries are and grid
85 /// values give the Neumann boundaries that should be applied there
86 /// @param backgroundVelocity a background velocity value
87 template<typename Vec3GridT>
88 typename Vec3GridT::Ptr
89 computePotentialFlow(const typename VectorToScalarGrid<Vec3GridT>::Type& potential,
90 const Vec3GridT& neumann,
91 const typename Vec3GridT::ValueType backgroundVelocity =
92 zeroVal<typename Vec3GridT::TreeType::ValueType>());
93
94
95 //////////////////////////////////////////////////////////
96
97 /// @cond OPENVDB_DOCS_INTERNAL
98
99 namespace potential_flow_internal {
100
101
102 /// @private
103 // helper function for retrieving a mask that comprises the outer-most layer of voxels
104 template<typename GridT>
105 typename GridT::TreeType::template ValueConverter<ValueMask>::Type::Ptr
106
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 extractOuterVoxelMask(GridT& inGrid)
107 {
108 using MaskTreeT = typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
109
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 16 typename MaskTreeT::Ptr interiorMask(new MaskTreeT(inGrid.tree(), false, TopologyCopy()));
110
3/8
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 16 typename MaskTreeT::Ptr boundaryMask(new MaskTreeT(inGrid.tree(), false, TopologyCopy()));
111
112
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 tools::erodeActiveValues(*interiorMask, /*iterations=*/1, tools::NN_FACE, tools::IGNORE_TILES);
113
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 tools::pruneInactive(*interiorMask);
114 boundaryMask->topologyDifference(*interiorMask);
115 8 return boundaryMask;
116 }
117
118
119 // computes Neumann velocities through sampling the gradient and velocities
120 template<typename Vec3GridT, typename GradientT>
121 struct ComputeNeumannVelocityOp
122 {
123 using ValueT = typename Vec3GridT::ValueType;
124 using VelocityAccessor = typename Vec3GridT::ConstAccessor;
125 using VelocitySamplerT = GridSampler<
126 typename Vec3GridT::ConstAccessor, BoxSampler>;
127 using GradientValueT = typename GradientT::TreeType::ValueType;
128
129 3 ComputeNeumannVelocityOp( const GradientT& gradient,
130 const Vec3GridT& velocity,
131 const ValueT& backgroundVelocity)
132 : mGradient(gradient)
133 , mVelocity(&velocity)
134 3 , mBackgroundVelocity(backgroundVelocity) { }
135
136 4 ComputeNeumannVelocityOp( const GradientT& gradient,
137 const ValueT& backgroundVelocity)
138 : mGradient(gradient)
139 4 , mBackgroundVelocity(backgroundVelocity) { }
140
141 240 void operator()(typename Vec3GridT::TreeType::LeafNodeType& leaf, size_t) const {
142 240 auto gradientAccessor = mGradient.getConstAccessor();
143
144 240 std::unique_ptr<VelocityAccessor> velocityAccessor;
145 240 std::unique_ptr<VelocitySamplerT> velocitySampler;
146
2/2
✓ Branch 0 taken 60 times.
✓ Branch 1 taken 80 times.
 240 if (mVelocity) {
147
2/4
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
 120 velocityAccessor.reset(new VelocityAccessor(mVelocity->getConstAccessor()));
148
1/2
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
 120 velocitySampler.reset(new VelocitySamplerT(*velocityAccessor, mVelocity->transform()));
149 }
150
151
2/2
✓ Branch 0 taken 14518 times.
✓ Branch 1 taken 140 times.
 25128 for (auto it = leaf.beginValueOn(); it; ++it) {
152
1/2
✓ Branch 1 taken 14518 times.
✗ Branch 2 not taken.
 24888 Coord ijk = it.getCoord();
153
1/2
✓ Branch 1 taken 14518 times.
✗ Branch 2 not taken.
 24888 auto gradient = gradientAccessor.getValue(ijk);
154
1/2
✓ Branch 1 taken 14518 times.
✗ Branch 2 not taken.
 24888 if (gradient.normalize()) {
155
1/2
✓ Branch 1 taken 14518 times.
✗ Branch 2 not taken.
 24888 const Vec3d xyz = mGradient.transform().indexToWorld(ijk);
156
3/4
✓ Branch 0 taken 6222 times.
✓ Branch 1 taken 8296 times.
✓ Branch 3 taken 6222 times.
✗ Branch 4 not taken.
 24888 const ValueT sampledVelocity = velocitySampler ?
157 4148 velocitySampler->wsSample(xyz) : zeroVal<ValueT>();
158
1/2
✓ Branch 1 taken 14518 times.
✗ Branch 2 not taken.
 24888 auto velocity = sampledVelocity + mBackgroundVelocity;
159 4148 auto value = gradient.dot(velocity) * gradient;
160
1/4
✓ Branch 1 taken 14518 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 24888 it.setValue(value);
161 }
162 else {
163 ✗ it.setValueOff();
164 }
165 }
166 240 }
167
168 private:
169 const GradientT& mGradient;
170 const Vec3GridT* mVelocity = nullptr;
171 const ValueT& mBackgroundVelocity;
172 }; // struct ComputeNeumannVelocityOp
173
174
175 // initializes the boundary conditions for use in the Poisson Solver
176 template<typename Vec3GridT, typename MaskT>
177 struct SolveBoundaryOp
178 {
179
2/6
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 3 SolveBoundaryOp(const Vec3GridT& velGrid, const MaskT& domainGrid)
180 3 : mVoxelSize(domainGrid.voxelSize()[0])
181 , mVelGrid(velGrid)
182 3 , mDomainGrid(domainGrid)
183 { }
184
185 418032 void operator()(const Coord& ijk, const Coord& neighbor,
186 double& source, double& diagonal) const {
187
188 418032 typename Vec3GridT::ConstAccessor velGridAccessor = mVelGrid.getAccessor();
189 const Coord diff = (ijk - neighbor);
190
191
2/2
✓ Branch 1 taken 61308 times.
✓ Branch 2 taken 147708 times.
 418032 if (velGridAccessor.isValueOn(ijk)) { // Neumann
192
1/2
✓ Branch 1 taken 61308 times.
✗ Branch 2 not taken.
 122616 const typename Vec3GridT::ValueType& sampleVel = velGridAccessor.getValue(ijk);
193 122616 source += mVoxelSize*diff[0]*sampleVel[0];
194 122616 source += mVoxelSize*diff[1]*sampleVel[1];
195 122616 source += mVoxelSize*diff[2]*sampleVel[2];
196 } else {
197 295416 diagonal -= 1; // Zero Dirichlet
198 }
199
200 }
201
202 const double mVoxelSize;
203 const Vec3GridT& mVelGrid;
204 const MaskT& mDomainGrid;
205 }; // struct SolveBoundaryOp
206
207
208 } // namespace potential_flow_internal
209
210 /// @endcond
211
212 ////////////////////////////////////////////////////////////////////////////
213
214 template<typename GridT, typename MaskT>
215 typename MaskT::Ptr
216 9 createPotentialFlowMask(const GridT& grid, int dilation)
217 {
218 using MaskTreeT = typename MaskT::TreeType;
219
220
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 8 times.
 9 if (!grid.hasUniformVoxels()) {
221
2/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 4 OPENVDB_THROW(ValueError, "Transform must have uniform voxels for Potential Flow mask.");
222 }
223
224 // construct a new mask grid representing the interior region
225 8 auto interior = interiorMask(grid);
226
227 // create a new mask grid from the interior topology
228
2/6
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 16 typename MaskTreeT::Ptr maskTree(new MaskTreeT(interior->tree(), false, TopologyCopy()));
229
2/6
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 16 typename MaskT::Ptr mask = MaskT::create(maskTree);
230
3/8
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 16 mask->setTransform(grid.transform().copy());
231
232
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 dilateActiveValues(*maskTree, dilation, NN_FACE_EDGE);
233
234 // subtract the interior region from the mask to leave just the exterior narrow band
235 mask->tree().topologyDifference(interior->tree());
236
237 8 return mask;
238 }
239
240
241 template<typename Vec3T, typename GridT, typename MaskT>
242 16 typename GridT::template ValueConverter<Vec3T>::Type::Ptr createPotentialFlowNeumannVelocities(
243 const GridT& collider,
244 const MaskT& domain,
245 const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
246 const Vec3T& backgroundVelocity)
247 {
248 using Vec3GridT = typename GridT::template ValueConverter<Vec3T>::Type;
249 using TreeT = typename Vec3GridT::TreeType;
250 using ValueT = typename TreeT::ValueType;
251 using GradientT = typename ScalarToVectorConverter<GridT>::Type;
252
253 using potential_flow_internal::ComputeNeumannVelocityOp;
254
255 // this method requires the collider to be a level set to generate the gradient
256 // use the tools::topologyToLevelset() method if you need to convert a mask into a level set
257
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 8 times.
 16 if (collider.getGridClass() != GRID_LEVEL_SET ||
258 !std::is_floating_point<typename GridT::TreeType::ValueType>::value) {
259
2/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 8 OPENVDB_THROW(TypeError, "Potential Flow expecting the collider to be a level set.");
260 }
261
262 // empty grid if there are no velocities
263
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 6 if (backgroundVelocity == zeroVal<Vec3T>() &&
264
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 2 times.
 6 (!boundaryVelocity || boundaryVelocity->empty())) {
265 ✗ auto neumann = Vec3GridT::create();
266
3/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 4 neumann->setTransform(collider.transform().copy());
267 return neumann;
268 }
269
270 // extract the intersection between the collider and the domain
271 using MaskTreeT = typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
272
3/6
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
 24 typename MaskTreeT::Ptr boundary(new MaskTreeT(domain.tree(), false, TopologyCopy()));
273 boundary->topologyIntersection(collider.tree());
274
275
3/8
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 24 typename TreeT::Ptr neumannTree(new TreeT(*boundary, zeroVal<ValueT>(), TopologyCopy()));
276 neumannTree->voxelizeActiveTiles();
277
278 // compute the gradient from the collider
279
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 12 const typename GradientT::Ptr gradient = tools::gradient(collider);
280
281
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 24 typename tree::LeafManager<TreeT> leafManager(*neumannTree);
282
283
5/6
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 3 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 3 times.
 12 if (boundaryVelocity && !boundaryVelocity->empty()) {
284 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
285 neumannOp(*gradient, *boundaryVelocity, backgroundVelocity);
286
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 leafManager.foreach(neumannOp, false);
287 }
288 else {
289 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
290 neumannOp(*gradient, backgroundVelocity);
291
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 6 leafManager.foreach(neumannOp, false);
292 }
293
294 // prune any inactive values
295
1/2
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
 12 tools::pruneInactive(*neumannTree);
296
297
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
 24 typename Vec3GridT::Ptr neumann(Vec3GridT::create(neumannTree));
298
3/8
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 7 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 24 neumann->setTransform(collider.transform().copy());
299
300 return neumann;
301 }
302
303
304 template<typename Vec3GridT, typename MaskT, typename InterrupterT>
305 typename VectorToScalarGrid<Vec3GridT>::Ptr
306 6 computeScalarPotential(const MaskT& domain, const Vec3GridT& neumann,
307 math::pcg::State& state, InterrupterT* interrupter)
308 {
309 using ScalarT = typename Vec3GridT::ValueType::value_type;
310 using ScalarTreeT = typename Vec3GridT::TreeType::template ValueConverter<ScalarT>::Type;
311 using ScalarGridT = typename Vec3GridT::template ValueConverter<ScalarT>::Type;
312
313 using potential_flow_internal::SolveBoundaryOp;
314
315 // create the solution tree and activate using domain topology
316
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 12 ScalarTreeT solveTree(domain.tree(), zeroVal<ScalarT>(), TopologyCopy());
317 solveTree.voxelizeActiveTiles();
318
319 6 util::NullInterrupter nullInterrupt;
320
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 6 if (!interrupter) interrupter = &nullInterrupt;
321
322 // solve for scalar potential
323 SolveBoundaryOp<Vec3GridT, MaskT> solve(neumann, domain);
324
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 typename ScalarTreeT::Ptr potentialTree =
325 poisson::solveWithBoundaryConditions(solveTree, solve, state, *interrupter, true);
326
327
2/6
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 12 auto potential = ScalarGridT::create(potentialTree);
328
3/8
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 12 potential->setTransform(domain.transform().copy());
329
330 6 return potential;
331 }
332
333
334 template<typename Vec3GridT>
335 typename Vec3GridT::Ptr
336 4 computePotentialFlow(const typename VectorToScalarGrid<Vec3GridT>::Type& potential,
337 const Vec3GridT& neumann,
338 const typename Vec3GridT::ValueType backgroundVelocity)
339 {
340 using Vec3T = const typename Vec3GridT::ValueType;
341 using potential_flow_internal::extractOuterVoxelMask;
342
343 // The VDB gradient op uses the background grid value, which is zero by default, when
344 // computing the gradient at the boundary. This works at the zero-dirichlet boundaries, but
345 // give spurious values at Neumann ones as the potential should be non-zero there. To avoid
346 // the extra error, we just substitute the Neumann condition on the boundaries.
347 // Technically, we should allow for some tangential velocity, coming from the gradient of
348 // potential. However, considering the voxelized nature of our solve, a decent approximation
349 // to a tangential derivative isn't probably worth our time. Any tangential component will be
350 // found in the next interior ring of voxels.
351
352 4 auto gradient = tools::gradient(potential);
353
354 // apply Neumann values to the gradient
355
356 20276 auto applyNeumann = [&gradient, &neumann] (
357 const MaskGrid::TreeType::LeafNodeType& leaf, size_t)
358 {
359 typename Vec3GridT::Accessor gradientAccessor = gradient->getAccessor();
360 typename Vec3GridT::ConstAccessor neumannAccessor = neumann.getAccessor();
361
4/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 107552 times.
✓ Branch 3 taken 5208 times.
✓ Branch 4 taken 97488 times.
✓ Branch 5 taken 4928 times.
 215176 for (auto it = leaf.beginValueOn(); it; ++it) {
362
2/6
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 4 taken 107552 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 97488 times.
✗ Branch 8 not taken.
 205040 const Coord ijk = it.getCoord();
363 typename Vec3GridT::ValueType value;
364
6/12
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 6 taken 107552 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 59222 times.
✓ Branch 9 taken 48330 times.
✓ Branch 11 taken 97488 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 828 times.
✓ Branch 14 taken 96660 times.
 205040 if (neumannAccessor.probeValue(ijk, value)) {
365
2/6
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 4 taken 59222 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 828 times.
✗ Branch 8 not taken.
 60050 gradientAccessor.setValue(ijk, value);
366 }
367 }
368 };
369
370 4 const MaskGrid::TreeType::Ptr boundary = extractOuterVoxelMask(*gradient);
371 8 typename tree::LeafManager<const typename MaskGrid::TreeType> leafManager(*boundary);
372 4 leafManager.foreach(applyNeumann);
373
374 // apply the background value to the gradient if supplied
375
376 if (backgroundVelocity != zeroVal<Vec3T>()) {
377 2465 auto applyBackgroundVelocity = [&backgroundVelocity] (
378 typename Vec3GridT::TreeType::LeafNodeType& leaf, size_t)
379 {
380
2/6
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1131300 times.
✓ Branch 5 taken 2464 times.
 1133764 for (auto it = leaf.beginValueOn(); it; ++it) {
381 1131300 it.setValue(it.getValue() - backgroundVelocity);
382 }
383 };
384
385 2 typename tree::LeafManager<typename Vec3GridT::TreeType> leafManager2(gradient->tree());
386 1 leafManager2.foreach(applyBackgroundVelocity);
387 }
388
389 4 return gradient;
390 }
391
392
393 ////////////////////////////////////////
394
395
396 // Explicit Template Instantiation
397
398 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
399
400 #ifdef OPENVDB_INSTANTIATE_POTENTIALFLOW
401 #include <openvdb/util/ExplicitInstantiation.h>
402 #endif
403
404 #define _FUNCTION(TreeT) \
405 Grid<TreeT>::Ptr createPotentialFlowNeumannVelocities(const FloatGrid&, const MaskGrid&, \
406 const Grid<TreeT>::ConstPtr, const TreeT::ValueType&)
407 OPENVDB_VEC3_TREE_INSTANTIATE(_FUNCTION)
408 #undef _FUNCTION
409
410 #define _FUNCTION(TreeT) \
411 VectorToScalarGrid<Grid<TreeT>>::Ptr computeScalarPotential(const MaskGrid&, const Grid<TreeT>&, \
412 math::pcg::State&, util::NullInterrupter*)
413 OPENVDB_VEC3_TREE_INSTANTIATE(_FUNCTION)
414 #undef _FUNCTION
415
416 #define _FUNCTION(TreeT) \
417 Grid<TreeT>::Ptr computePotentialFlow( \
418 const VectorToScalarGrid<Grid<TreeT>>::Type&, const Grid<TreeT>&, const TreeT::ValueType)
419 OPENVDB_VEC3_TREE_INSTANTIATE(_FUNCTION)
420 #undef _FUNCTION
421
422 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
423
424
425 } // namespace tools
426 } // namespace OPENVDB_VERSION_NAME
427 } // namespace openvdb
428
429 #endif // OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED
430