HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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 =
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
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 ||
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>
304 typename VectorToScalarGrid<Vec3GridT>::Ptr
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's equation ∇2x = b for x, where b is a vector comprising the values of all of the activ...
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...
GLsizei const GLfloat * value
Definition: glcorearb.h:824
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:239
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:1133
typename VecGridT::template ValueConverter< typename VecGridT::ValueType::value_type >::Type Type
Definition: PotentialFlow.h:31
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
uint64 value_type
Definition: GA_PrimCompat.h:29
Metafunction to convert a vector-valued grid type to a scalar grid type.
Definition: PotentialFlow.h:29
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...
GLsizei GLsizei GLchar * source
Definition: glcorearb.h:803
GLint GLuint mask
Definition: glcorearb.h:124
ScalarGridType::template ValueConverter< VectorValueT >::Type Type
Definition: GridOperators.h:44
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:84
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
#define OPENVDB_VEC3_TREE_INSTANTIATE(Function)
Definition: version.h:165
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
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.
MaskT::Ptr createPotentialFlowMask(const GridT &grid, int dilation=5)
Construct a mask for the Potential Flow domain.
ScalarToVectorConverter< GridType >::Type::Ptr gradient(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the gradient of the given scalar grid.
Implementation of morphological dilation and erosion.
GU_API void solve(const GU_Detail &gdp_a, const GA_Range &pts_a, const GU_Detail &gdp_b, const GA_Range &pts_b, Method method, bool compute_distortion, Result &result)
Definition: core.h:1131
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:683
Apply an operator to an input grid to produce an output grid with the same active voxel topology but ...
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...
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.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
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74