OpenVDB  8.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 
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 
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 
42 template<typename GridT, typename MaskT = typename GridT::template ValueConverter<ValueMask>::Type>
43 inline typename MaskT::Ptr
44 createPotentialFlowMask(const GridT& grid, int dilation = 5);
45 
46 
56 template<typename Vec3T, typename GridT, typename MaskT>
57 inline 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 
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 
86 template<typename Vec3GridT>
87 inline typename Vec3GridT::Ptr
89  const Vec3GridT& neumann,
90  const typename Vec3GridT::ValueType backgroundVelocity =
91  zeroVal<typename Vec3GridT::TreeType::ValueType>());
92 
93 
95 
96 
97 namespace potential_flow_internal {
98 
99 
101 // helper function for retrieving a mask that comprises the outer-most layer of voxels
102 template<typename GridT>
103 inline typename GridT::TreeType::template ValueConverter<ValueMask>::Type::Ptr
104 extractOuterVoxelMask(GridT& inGrid)
105 {
106  using MaskTreeT = typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
107  typename MaskTreeT::Ptr interiorMask(new MaskTreeT(inGrid.tree(), false, TopologyCopy()));
108  typename MaskTreeT::Ptr boundaryMask(new MaskTreeT(inGrid.tree(), false, TopologyCopy()));
109 
112  boundaryMask->topologyDifference(*interiorMask);
113  return boundaryMask;
114 }
115 
116 
117 // computes Neumann velocities through sampling the gradient and velocities
118 template<typename Vec3GridT, typename GradientT>
120 {
121  using ValueT = typename Vec3GridT::ValueType;
122  using VelocityAccessor = typename Vec3GridT::ConstAccessor;
124  typename Vec3GridT::ConstAccessor, BoxSampler>;
125  using GradientValueT = typename GradientT::TreeType::ValueType;
126 
128  const Vec3GridT& velocity,
129  const ValueT& backgroundVelocity)
130  : mGradient(gradient)
131  , mVelocity(&velocity)
132  , mBackgroundVelocity(backgroundVelocity) { }
133 
135  const ValueT& backgroundVelocity)
136  : mGradient(gradient)
137  , mBackgroundVelocity(backgroundVelocity) { }
138 
139  void operator()(typename Vec3GridT::TreeType::LeafNodeType& leaf, size_t) const {
140  auto gradientAccessor = mGradient.getConstAccessor();
141 
142  std::unique_ptr<VelocityAccessor> velocityAccessor;
143  std::unique_ptr<VelocitySamplerT> velocitySampler;
144  if (mVelocity) {
145  velocityAccessor.reset(new VelocityAccessor(mVelocity->getConstAccessor()));
146  velocitySampler.reset(new VelocitySamplerT(*velocityAccessor, mVelocity->transform()));
147  }
148 
149  for (auto it = leaf.beginValueOn(); it; ++it) {
150  Coord ijk = it.getCoord();
151  auto gradient = gradientAccessor.getValue(ijk);
152  if (gradient.normalize()) {
153  const Vec3d xyz = mGradient.transform().indexToWorld(ijk);
154  const ValueT sampledVelocity = velocitySampler ?
155  velocitySampler->wsSample(xyz) : zeroVal<ValueT>();
156  auto velocity = sampledVelocity + mBackgroundVelocity;
157  auto value = gradient.dot(velocity) * gradient;
158  it.setValue(value);
159  }
160  else {
161  it.setValueOff();
162  }
163  }
164  }
165 
166 private:
167  const GradientT& mGradient;
168  const Vec3GridT* mVelocity = nullptr;
169  const ValueT& mBackgroundVelocity;
170 }; // struct ComputeNeumannVelocityOp
171 
172 
173 // initalizes the boundary conditions for use in the Poisson Solver
174 template<typename Vec3GridT, typename MaskT>
176 {
177  SolveBoundaryOp(const Vec3GridT& velGrid, const MaskT& domainGrid)
178  : mVoxelSize(domainGrid.voxelSize()[0])
179  , mVelGrid(velGrid)
180  , mDomainGrid(domainGrid)
181  { }
182 
183  void operator()(const Coord& ijk, const Coord& neighbor,
184  double& source, double& diagonal) const {
185 
186  typename Vec3GridT::ConstAccessor velGridAccessor = mVelGrid.getAccessor();
187  const Coord diff = (ijk - neighbor);
188 
189  if (velGridAccessor.isValueOn(ijk)) { // Neumann
190  const typename Vec3GridT::ValueType& sampleVel = velGridAccessor.getValue(ijk);
191  source += mVoxelSize*diff[0]*sampleVel[0];
192  source += mVoxelSize*diff[1]*sampleVel[1];
193  source += mVoxelSize*diff[2]*sampleVel[2];
194  } else {
195  diagonal -= 1; // Zero Dirichlet
196  }
197 
198  }
199 
200  const double mVoxelSize;
201  const Vec3GridT& mVelGrid;
202  const MaskT& mDomainGrid;
203 }; // struct SolveBoundaryOp
204 
205 
206 } // namespace potential_flow_internal
207 
208 
210 
211 template<typename GridT, typename MaskT>
212 inline typename MaskT::Ptr
213 createPotentialFlowMask(const GridT& grid, int dilation)
214 {
215  using MaskTreeT = typename MaskT::TreeType;
216 
217  if (!grid.hasUniformVoxels()) {
218  OPENVDB_THROW(ValueError, "Transform must have uniform voxels for Potential Flow mask.");
219  }
220 
221  // construct a new mask grid representing the interior region
222  auto interior = interiorMask(grid);
223 
224  // create a new mask grid from the interior topology
225  typename MaskTreeT::Ptr maskTree(new MaskTreeT(interior->tree(), false, TopologyCopy()));
226  typename MaskT::Ptr mask = MaskT::create(maskTree);
227  mask->setTransform(grid.transform().copy());
228 
229  dilateActiveValues(*maskTree, dilation, NN_FACE_EDGE);
230 
231  // subtract the interior region from the mask to leave just the exterior narrow band
232  mask->tree().topologyDifference(interior->tree());
233 
234  return mask;
235 }
236 
237 
238 template<typename Vec3T, typename GridT, typename MaskT>
239 typename GridT::template ValueConverter<Vec3T>::Type::Ptr createPotentialFlowNeumannVelocities(
240  const GridT& collider,
241  const MaskT& domain,
242  const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
243  const Vec3T& backgroundVelocity)
244 {
245  using Vec3GridT = typename GridT::template ValueConverter<Vec3T>::Type;
246  using TreeT = typename Vec3GridT::TreeType;
247  using ValueT = typename TreeT::ValueType;
248  using GradientT = typename ScalarToVectorConverter<GridT>::Type;
249 
251 
252  // this method requires the collider to be a level set to generate the gradient
253  // use the tools::topologyToLevelset() method if you need to convert a mask into a level set
254  if (collider.getGridClass() != GRID_LEVEL_SET ||
255  !std::is_floating_point<typename GridT::TreeType::ValueType>::value) {
256  OPENVDB_THROW(TypeError, "Potential Flow expecting the collider to be a level set.");
257  }
258 
259  // empty grid if there are no velocities
260  if (backgroundVelocity == zeroVal<Vec3T>() &&
261  (!boundaryVelocity || boundaryVelocity->empty())) {
262  auto neumann = Vec3GridT::create();
263  neumann->setTransform(collider.transform().copy());
264  return neumann;
265  }
266 
267  // extract the intersection between the collider and the domain
268  using MaskTreeT = typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
269  typename MaskTreeT::Ptr boundary(new MaskTreeT(domain.tree(), false, TopologyCopy()));
270  boundary->topologyIntersection(collider.tree());
271 
272  typename TreeT::Ptr neumannTree(new TreeT(*boundary, zeroVal<ValueT>(), TopologyCopy()));
273  neumannTree->voxelizeActiveTiles();
274 
275  // compute the gradient from the collider
276  const typename GradientT::Ptr gradient = tools::gradient(collider);
277 
278  typename tree::LeafManager<TreeT> leafManager(*neumannTree);
279 
280  if (boundaryVelocity && !boundaryVelocity->empty()) {
281  ComputeNeumannVelocityOp<Vec3GridT, GradientT>
282  neumannOp(*gradient, *boundaryVelocity, backgroundVelocity);
283  leafManager.foreach(neumannOp, false);
284  }
285  else {
286  ComputeNeumannVelocityOp<Vec3GridT, GradientT>
287  neumannOp(*gradient, backgroundVelocity);
288  leafManager.foreach(neumannOp, false);
289  }
290 
291  // prune any inactive values
292  tools::pruneInactive(*neumannTree);
293 
294  typename Vec3GridT::Ptr neumann(Vec3GridT::create(neumannTree));
295  neumann->setTransform(collider.transform().copy());
296 
297  return neumann;
298 }
299 
300 
301 template<typename Vec3GridT, typename MaskT, typename InterrupterT>
302 inline typename VectorToScalarGrid<Vec3GridT>::Ptr
303 computeScalarPotential(const MaskT& domain, const Vec3GridT& neumann,
304  math::pcg::State& state, InterrupterT* interrupter)
305 {
306  using ScalarT = typename Vec3GridT::ValueType::value_type;
307  using ScalarTreeT = typename Vec3GridT::TreeType::template ValueConverter<ScalarT>::Type;
308  using ScalarGridT = typename Vec3GridT::template ValueConverter<ScalarT>::Type;
309 
311 
312  // create the solution tree and activate using domain topology
313  ScalarTreeT solveTree(domain.tree(), zeroVal<ScalarT>(), TopologyCopy());
314  solveTree.voxelizeActiveTiles();
315 
316  util::NullInterrupter nullInterrupt;
317  if (!interrupter) interrupter = &nullInterrupt;
318 
319  // solve for scalar potential
320  SolveBoundaryOp<Vec3GridT, MaskT> solve(neumann, domain);
321  typename ScalarTreeT::Ptr potentialTree =
322  poisson::solveWithBoundaryConditions(solveTree, solve, state, *interrupter, true);
323 
324  auto potential = ScalarGridT::create(potentialTree);
325  potential->setTransform(domain.transform().copy());
326 
327  return potential;
328 }
329 
330 
331 template<typename Vec3GridT>
332 inline typename Vec3GridT::Ptr
334  const Vec3GridT& neumann,
335  const typename Vec3GridT::ValueType backgroundVelocity)
336 {
337  using Vec3T = const typename Vec3GridT::ValueType;
338  using potential_flow_internal::extractOuterVoxelMask;
339 
340  // The VDB gradient op uses the background grid value, which is zero by default, when
341  // computing the gradient at the boundary. This works at the zero-dirichlet boundaries, but
342  // give spurious values at Neumann ones as the potential should be non-zero there. To avoid
343  // the extra error, we just substitute the Neumann condition on the boundaries.
344  // Technically, we should allow for some tangential velocity, coming from the gradient of
345  // potential. However, considering the voxelized nature of our solve, a decent approximation
346  // to a tangential derivative isn't probably worth our time. Any tangential component will be
347  // found in the next interior ring of voxels.
348 
349  auto gradient = tools::gradient(potential);
350 
351  // apply Neumann values to the gradient
352 
353  auto applyNeumann = [&gradient, &neumann] (
354  const MaskGrid::TreeType::LeafNodeType& leaf, size_t)
355  {
356  typename Vec3GridT::Accessor gradientAccessor = gradient->getAccessor();
357  typename Vec3GridT::ConstAccessor neumannAccessor = neumann.getAccessor();
358  for (auto it = leaf.beginValueOn(); it; ++it) {
359  const Coord ijk = it.getCoord();
360  typename Vec3GridT::ValueType value;
361  if (neumannAccessor.probeValue(ijk, value)) {
362  gradientAccessor.setValue(ijk, value);
363  }
364  }
365  };
366 
367  const MaskGrid::TreeType::Ptr boundary = extractOuterVoxelMask(*gradient);
368  typename tree::LeafManager<const typename MaskGrid::TreeType> leafManager(*boundary);
369  leafManager.foreach(applyNeumann);
370 
371  // apply the background value to the gradient if supplied
372 
373  if (backgroundVelocity != zeroVal<Vec3T>()) {
374  auto applyBackgroundVelocity = [&backgroundVelocity] (
375  typename Vec3GridT::TreeType::LeafNodeType& leaf, size_t)
376  {
377  for (auto it = leaf.beginValueOn(); it; ++it) {
378  it.setValue(it.getValue() - backgroundVelocity);
379  }
380  };
381 
382  typename tree::LeafManager<typename Vec3GridT::TreeType> leafManager2(gradient->tree());
383  leafManager2.foreach(applyBackgroundVelocity);
384  }
385 
386  return gradient;
387 }
388 
389 
391 
392 } // namespace tools
393 } // namespace OPENVDB_VERSION_NAME
394 } // namespace openvdb
395 
396 #endif // OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED
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:579
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: openvdb/Types.h:560
Definition: openvdb/Exceptions.h:64
Definition: openvdb/Exceptions.h:65
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:26
Class that provides the interface for continuous sampling of values in a tree.
Definition: Interpolation.h:284
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:757
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:239
MaskT::Ptr createPotentialFlowMask(const GridT &grid, int dilation=5)
Construct a mask for the Potential Flow domain.
Definition: PotentialFlow.h:213
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:303
@ NN_FACE_EDGE
Definition: Morphology.h:56
@ NN_FACE
Definition: Morphology.h:56
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:111
ScalarToVectorConverter< GridType >::Type::Ptr gradient(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the gradient of the given scalar grid.
Definition: GridOperators.h:996
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:1126
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:1049
@ IGNORE_TILES
Definition: Morphology.h:78
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:333
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:354
@ GRID_LEVEL_SET
Definition: openvdb/Types.h:333
Definition: openvdb/Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition: openvdb/Exceptions.h:74
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
Definition: Interpolation.h:120
ScalarGridType::template ValueConverter< VectorValueT >::Type Type
Definition: GridOperators.h:41
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
typename Vec3GridT::ValueType ValueT
Definition: PotentialFlow.h:121
ComputeNeumannVelocityOp(const GradientT &gradient, const ValueT &backgroundVelocity)
Definition: PotentialFlow.h:134
void operator()(typename Vec3GridT::TreeType::LeafNodeType &leaf, size_t) const
Definition: PotentialFlow.h:139
typename Vec3GridT::ConstAccessor VelocityAccessor
Definition: PotentialFlow.h:122
typename GradientT::TreeType::ValueType GradientValueT
Definition: PotentialFlow.h:125
ComputeNeumannVelocityOp(const GradientT &gradient, const Vec3GridT &velocity, const ValueT &backgroundVelocity)
Definition: PotentialFlow.h:127
const Vec3GridT & mVelGrid
Definition: PotentialFlow.h:201
const double mVoxelSize
Definition: PotentialFlow.h:200
const MaskT & mDomainGrid
Definition: PotentialFlow.h:202
void operator()(const Coord &ijk, const Coord &neighbor, double &source, double &diagonal) const
Definition: PotentialFlow.h:183
SolveBoundaryOp(const Vec3GridT &velGrid, const MaskT &domainGrid)
Definition: PotentialFlow.h:177
Dummy NOOP interrupter class defining interface.
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:178