OpenVDB  12.0.0
PotentialFlow.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: Apache-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 
22 namespace openvdb {
24 namespace OPENVDB_VERSION_NAME {
25 namespace tools {
26 
27 /// @brief Metafunction to convert a vector-valued grid type to a scalar grid type
28 template<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
42 template<typename GridT, typename MaskT = typename GridT::template ValueConverter<ValueMask>::Type>
43 typename MaskT::Ptr
44 createPotentialFlowMask(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.
56 template<typename Vec3T, typename GridT, typename MaskT>
57 typename GridT::template ValueConverter<Vec3T>::Type::Ptr
58 createPotentialFlowNeumannVelocities(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.
74 template<typename Vec3GridT, typename MaskT, typename InterrupterT = util::NullInterrupter>
76 computeScalarPotential(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
86 template<typename Vec3GridT>
87 typename 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 
98 namespace potential_flow_internal {
99 
100 
101 /// @private
102 // helper function for retrieving a mask that comprises the outer-most layer of voxels
103 template<typename GridT>
104 typename GridT::TreeType::template ValueConverter<ValueMask>::Type::Ptr
105 extractOuterVoxelMask(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 
113  boundaryMask->topologyDifference(*interiorMask);
114  return boundaryMask;
115 }
116 
117 
118 // computes Neumann velocities through sampling the gradient and velocities
119 template<typename Vec3GridT, typename GradientT>
120 struct 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 
167 private:
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
175 template<typename Vec3GridT, typename MaskT>
176 struct 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 
213 template<typename GridT, typename MaskT>
214 typename MaskT::Ptr
215 createPotentialFlowMask(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 
240 template<typename Vec3T, typename GridT, typename MaskT>
241 typename 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 
303 template<typename Vec3GridT, typename MaskT, typename InterrupterT>
305 computeScalarPotential(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 
333 template<typename Vec3GridT>
334 typename 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
Solve Poisson&#39;s equation ∇2x = b for x, where b is a vector comprising the values of all of the acti...
Definition: Morphology.h:60
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:1134
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Apply an operator to an input grid to produce an output grid with the same active voxel topology but ...
Class that provides the interface for continuous sampling of values in a tree.
Definition: Interpolation.h:284
typename Type::Ptr Ptr
Definition: PotentialFlow.h:32
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
typename Type::ConstPtr ConstPtr
Definition: PotentialFlow.h:33
Construct boolean mask grids from grids of arbitrary type.
ScalarToVectorConverter< GridType >::Type::Ptr gradient(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the gradient of the given scalar grid.
Definition: GridOperators.h:1001
Base class for interrupters.
Definition: NullInterrupter.h:25
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:1608
Definition: Exceptions.h:65
Definition: Interpolation.h:120
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
Definition: Types.h:455
Definition: Morphology.h:82
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:769
Implementation of morphological dilation and erosion.
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:1057
SharedPtr< Grid > Ptr
Definition: Grid.h:573
Definition: Exceptions.h:13
ScalarGridType::template ValueConverter< VectorValueT >::Type Type
Definition: GridOperators.h:44
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
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
Definition: Morphology.h:60
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:484
Definition: Exceptions.h:64
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
Metafunction to convert a vector-valued grid type to a scalar grid type.
Definition: PotentialFlow.h:29
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
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:683
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
typename VecGridT::template ValueConverter< typename VecGridT::ValueType::value_type >::Type Type
Definition: PotentialFlow.h:31
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
#define OPENVDB_VEC3_TREE_INSTANTIATE(Function)
Definition: version.h.in:164
MaskT::Ptr createPotentialFlowMask(const GridT &grid, int dilation=5)
Construct a mask for the Potential Flow domain.
Definition: PotentialFlow.h:215
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218