OpenVDB  9.1.1
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>
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
Solve Poisson&#39;s equation ∇2x = b for x, where b is a vector comprising the values of all of the acti...
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 valu...
Definition: PoissonSolver.h:767
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
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
typename VecGridT::template ValueConverter< typename VecGridT::ValueType::value_type >::Type Type
Definition: PotentialFlow.h:32
Base class for interrupters.
Definition: NullInterrupter.h:25
Definition: Morphology.h:80
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
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 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
Apply an operator to an input grid to produce an output grid with the same active voxel topology but ...
Definition: Exceptions.h:64
Definition: Interpolation.h:119
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
Definition: Exceptions.h:65
Construct boolean mask grids from grids of arbitrary type.
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
Definition: Types.h:416
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:84
typename Type::ConstPtr ConstPtr
Definition: PotentialFlow.h:34
MaskT::Ptr createPotentialFlowMask(const GridT &grid, int dilation=5)
Construct a mask for the Potential Flow domain.
Definition: PotentialFlow.h:216
typename Type::Ptr Ptr
Definition: PotentialFlow.h:33
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:13
ValueT value
Definition: GridBuilder.h:1287
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
Definition: Morphology.h:58
Definition: Morphology.h:58
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
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
Metafunction to convert a vector-valued grid type to a scalar grid type.
Definition: PotentialFlow.h:30
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
Class that provides the interface for continuous sampling of values in a tree.
Definition: Interpolation.h:283
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
SharedPtr< Grid > Ptr
Definition: Grid.h:575
ScalarGridType::template ValueConverter< VectorValueT >::Type Type
Definition: GridOperators.h:43
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:644
#define OPENVDB_VEC3_TREE_INSTANTIATE(Function)
Definition: version.h.in:149
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:202