OpenVDB  9.1.0
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 #include <openvdb/openvdb.h>
21 
22 
23 namespace openvdb {
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>
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>
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
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 extractOuterVoxelMask(GridT& inGrid)
107 {
108  using MaskTreeT = typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
109  typename MaskTreeT::Ptr interiorMask(new MaskTreeT(inGrid.tree(), false, TopologyCopy()));
110  typename MaskTreeT::Ptr boundaryMask(new MaskTreeT(inGrid.tree(), false, TopologyCopy()));
111 
114  boundaryMask->topologyDifference(*interiorMask);
115  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  ComputeNeumannVelocityOp( const GradientT& gradient,
130  const Vec3GridT& velocity,
131  const ValueT& backgroundVelocity)
132  : mGradient(gradient)
133  , mVelocity(&velocity)
134  , mBackgroundVelocity(backgroundVelocity) { }
135 
136  ComputeNeumannVelocityOp( const GradientT& gradient,
137  const ValueT& backgroundVelocity)
138  : mGradient(gradient)
139  , mBackgroundVelocity(backgroundVelocity) { }
140 
141  void operator()(typename Vec3GridT::TreeType::LeafNodeType& leaf, size_t) const {
142  auto gradientAccessor = mGradient.getConstAccessor();
143 
144  std::unique_ptr<VelocityAccessor> velocityAccessor;
145  std::unique_ptr<VelocitySamplerT> velocitySampler;
146  if (mVelocity) {
147  velocityAccessor.reset(new VelocityAccessor(mVelocity->getConstAccessor()));
148  velocitySampler.reset(new VelocitySamplerT(*velocityAccessor, mVelocity->transform()));
149  }
150 
151  for (auto it = leaf.beginValueOn(); it; ++it) {
152  Coord ijk = it.getCoord();
153  auto gradient = gradientAccessor.getValue(ijk);
154  if (gradient.normalize()) {
155  const Vec3d xyz = mGradient.transform().indexToWorld(ijk);
156  const ValueT sampledVelocity = velocitySampler ?
157  velocitySampler->wsSample(xyz) : zeroVal<ValueT>();
158  auto velocity = sampledVelocity + mBackgroundVelocity;
159  auto value = gradient.dot(velocity) * gradient;
160  it.setValue(value);
161  }
162  else {
163  it.setValueOff();
164  }
165  }
166  }
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  SolveBoundaryOp(const Vec3GridT& velGrid, const MaskT& domainGrid)
180  : mVoxelSize(domainGrid.voxelSize()[0])
181  , mVelGrid(velGrid)
182  , mDomainGrid(domainGrid)
183  { }
184 
185  void operator()(const Coord& ijk, const Coord& neighbor,
186  double& source, double& diagonal) const {
187 
188  typename Vec3GridT::ConstAccessor velGridAccessor = mVelGrid.getAccessor();
189  const Coord diff = (ijk - neighbor);
190 
191  if (velGridAccessor.isValueOn(ijk)) { // Neumann
192  const typename Vec3GridT::ValueType& sampleVel = velGridAccessor.getValue(ijk);
193  source += mVoxelSize*diff[0]*sampleVel[0];
194  source += mVoxelSize*diff[1]*sampleVel[1];
195  source += mVoxelSize*diff[2]*sampleVel[2];
196  } else {
197  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 createPotentialFlowMask(const GridT& grid, int dilation)
217 {
218  using MaskTreeT = typename MaskT::TreeType;
219 
220  if (!grid.hasUniformVoxels()) {
221  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  auto interior = interiorMask(grid);
226 
227  // create a new mask grid from the interior topology
228  typename MaskTreeT::Ptr maskTree(new MaskTreeT(interior->tree(), false, TopologyCopy()));
229  typename MaskT::Ptr mask = MaskT::create(maskTree);
230  mask->setTransform(grid.transform().copy());
231 
232  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  return mask;
238 }
239 
240 
241 template<typename Vec3T, typename GridT, typename MaskT>
242 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  if (collider.getGridClass() != GRID_LEVEL_SET ||
259  OPENVDB_THROW(TypeError, "Potential Flow expecting the collider to be a level set.");
260  }
261 
262  // empty grid if there are no velocities
263  if (backgroundVelocity == zeroVal<Vec3T>() &&
264  (!boundaryVelocity || boundaryVelocity->empty())) {
265  auto neumann = Vec3GridT::create();
266  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  typename MaskTreeT::Ptr boundary(new MaskTreeT(domain.tree(), false, TopologyCopy()));
273  boundary->topologyIntersection(collider.tree());
274 
275  typename TreeT::Ptr neumannTree(new TreeT(*boundary, zeroVal<ValueT>(), TopologyCopy()));
276  neumannTree->voxelizeActiveTiles();
277 
278  // compute the gradient from the collider
279  const typename GradientT::Ptr gradient = tools::gradient(collider);
280 
281  typename tree::LeafManager<TreeT> leafManager(*neumannTree);
282 
283  if (boundaryVelocity && !boundaryVelocity->empty()) {
284  ComputeNeumannVelocityOp<Vec3GridT, GradientT>
285  neumannOp(*gradient, *boundaryVelocity, backgroundVelocity);
286  leafManager.foreach(neumannOp, false);
287  }
288  else {
289  ComputeNeumannVelocityOp<Vec3GridT, GradientT>
290  neumannOp(*gradient, backgroundVelocity);
291  leafManager.foreach(neumannOp, false);
292  }
293 
294  // prune any inactive values
295  tools::pruneInactive(*neumannTree);
296 
297  typename Vec3GridT::Ptr neumann(Vec3GridT::create(neumannTree));
298  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 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  ScalarTreeT solveTree(domain.tree(), zeroVal<ScalarT>(), TopologyCopy());
317  solveTree.voxelizeActiveTiles();
318 
319  util::NullInterrupter nullInterrupt;
320  if (!interrupter) interrupter = &nullInterrupt;
321 
322  // solve for scalar potential
323  SolveBoundaryOp<Vec3GridT, MaskT> solve(neumann, domain);
324  typename ScalarTreeT::Ptr potentialTree =
325  poisson::solveWithBoundaryConditions(solveTree, solve, state, *interrupter, true);
326 
327  auto potential = ScalarGridT::create(potentialTree);
328  potential->setTransform(domain.transform().copy());
329 
330  return potential;
331 }
332 
333 
334 template<typename Vec3GridT>
335 typename Vec3GridT::Ptr
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  auto gradient = tools::gradient(potential);
353 
354  // apply Neumann values to the gradient
355 
356  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  for (auto it = leaf.beginValueOn(); it; ++it) {
362  const Coord ijk = it.getCoord();
363  typename Vec3GridT::ValueType value;
364  if (neumannAccessor.probeValue(ijk, value)) {
365  gradientAccessor.setValue(ijk, value);
366  }
367  }
368  };
369 
370  const MaskGrid::TreeType::Ptr boundary = extractOuterVoxelMask(*gradient);
371  typename tree::LeafManager<const typename MaskGrid::TreeType> leafManager(*boundary);
372  leafManager.foreach(applyNeumann);
373 
374  // apply the background value to the gradient if supplied
375 
376  if (backgroundVelocity != zeroVal<Vec3T>()) {
377  auto applyBackgroundVelocity = [&backgroundVelocity] (
378  typename Vec3GridT::TreeType::LeafNodeType& leaf, size_t)
379  {
380  for (auto it = leaf.beginValueOn(); it; ++it) {
381  it.setValue(it.getValue() - backgroundVelocity);
382  }
383  };
384 
385  typename tree::LeafManager<typename Vec3GridT::TreeType> leafManager2(gradient->tree());
386  leafManager2.foreach(applyBackgroundVelocity);
387  }
388 
389  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
402 #endif
403 
404 #define _FUNCTION(TreeT) \
405  Grid<TreeT>::Ptr createPotentialFlowNeumannVelocities(const FloatGrid&, const MaskGrid&, \
406  const Grid<TreeT>::ConstPtr, const TreeT::ValueType&)
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*)
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)
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
ValueT value
Definition: GridBuilder.h:1287
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:575
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
State solve(const PositiveDefMatrix &A, const Vector< typename PositiveDefMatrix::ValueType > &b, Vector< typename PositiveDefMatrix::ValueType > &x, Preconditioner< typename PositiveDefMatrix::ValueType > &preconditioner, const State &termination=terminationDefaults< typename PositiveDefMatrix::ValueType >())
Solve Ax = b via the preconditioned conjugate gradient method.
Definition: ConjGradient.h:1612
Vec3< double > Vec3d
Definition: Vec3.h:668
TreeType::Ptr solveWithBoundaryConditions(const TreeType &, const BoundaryOp &, math::pcg::State &, Interrupter &, bool staggered=false)
Solve ∇2x = b for x with user-specified boundary conditions, where b is a vector comprising the value...
Definition: PoissonSolver.h:767
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:242
MaskT::Ptr createPotentialFlowMask(const GridT &grid, int dilation=5)
Construct a mask for the Potential Flow domain.
Definition: PotentialFlow.h:216
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:306
@ NN_FACE_EDGE
Definition: Morphology.h:58
@ NN_FACE
Definition: Morphology.h:58
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:1000
void erodeActiveValues(TreeOrLeafManagerT &tree, const int iterations=1, const NearestNeighbors nn=NN_FACE, const TilePolicy mode=PRESERVE_TILES, const bool threaded=true)
Topologically erode all active values (i.e. both voxels and tiles) in a tree using one of three neare...
Definition: Morphology.h:1132
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:1055
@ IGNORE_TILES
Definition: Morphology.h:80
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:336
void pruneInactive(TreeT &tree, bool threaded=true, size_t grainSize=1)
Reduce the memory footprint of a tree by replacing with background tiles any nodes whose values are a...
Definition: Prune.h:355
@ 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:43
Metafunction to convert a vector-valued grid type to a scalar grid type.
Definition: PotentialFlow.h:30
typename Type::Ptr Ptr
Definition: PotentialFlow.h:33
typename Type::ConstPtr ConstPtr
Definition: PotentialFlow.h:34
typename VecGridT::template ValueConverter< typename VecGridT::ValueType::value_type >::Type Type
Definition: PotentialFlow.h:32
Base class for interrupters.
Definition: NullInterrupter.h:26
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:202
#define OPENVDB_VEC3_TREE_INSTANTIATE(Function)
Definition: version.h.in:149