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