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