10#ifndef OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED
11#define OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED
28template<
typename VecGr
idT>
31 typename VecGridT::template ValueConverter<typename VecGridT::ValueType::value_type>::Type;
32 using Ptr =
typename Type::Ptr;
42template<typename GridT, typename MaskT = typename GridT::template ValueConverter<ValueMask>::Type>
56template<
typename Vec3T,
typename Gr
idT,
typename MaskT>
57typename GridT::template ValueConverter<Vec3T>::Type::Ptr
59 const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
60 const Vec3T& backgroundVelocity);
74template<
typename Vec3Gr
idT,
typename MaskT,
typename InterrupterT = util::NullInterrupter>
77 InterrupterT* interrupter =
nullptr);
86template<
typename Vec3Gr
idT>
87typename Vec3GridT::Ptr
89 const Vec3GridT& neumann,
90 const typename Vec3GridT::ValueType backgroundVelocity =
91 zeroVal<typename Vec3GridT::TreeType::ValueType>());
98namespace potential_flow_internal {
103template<
typename Gr
idT>
104typename GridT::TreeType::template ValueConverter<ValueMask>::Type::Ptr
105extractOuterVoxelMask(GridT& inGrid)
107 using MaskTreeT =
typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
109 typename MaskTreeT::Ptr boundaryMask(
new MaskTreeT(inGrid.tree(),
false,
TopologyCopy()));
111 tools::erodeActiveValues(*interiorMask, 1, tools::NN_FACE, tools::IGNORE_TILES);
112 tools::pruneInactive(*interiorMask);
113 boundaryMask->topologyDifference(*interiorMask);
119template<
typename Vec3Gr
idT,
typename GradientT>
120struct ComputeNeumannVelocityOp
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;
128 ComputeNeumannVelocityOp(
const GradientT& gradient,
129 const Vec3GridT& velocity,
130 const ValueT& backgroundVelocity)
131 : mGradient(gradient)
132 , mVelocity(&velocity)
133 , mBackgroundVelocity(backgroundVelocity) { }
135 ComputeNeumannVelocityOp(
const GradientT& gradient,
136 const ValueT& backgroundVelocity)
137 : mGradient(gradient)
138 , mBackgroundVelocity(backgroundVelocity) { }
140 void operator()(
typename Vec3GridT::TreeType::LeafNodeType& leaf,
size_t)
const {
141 auto gradientAccessor = mGradient.getConstAccessor();
143 std::unique_ptr<VelocityAccessor> velocityAccessor;
144 std::unique_ptr<VelocitySamplerT> velocitySampler;
146 velocityAccessor.reset(
new VelocityAccessor(mVelocity->getConstAccessor()));
147 velocitySampler.reset(
new VelocitySamplerT(*velocityAccessor, mVelocity->transform()));
150 for (
auto it = leaf.beginValueOn(); it; ++it) {
151 Coord ijk = it.getCoord();
152 auto gradient = gradientAccessor.getValue(ijk);
154 const Vec3d xyz = mGradient.transform().indexToWorld(ijk);
155 const ValueT sampledVelocity = velocitySampler ?
156 velocitySampler->wsSample(xyz) : zeroVal<ValueT>();
157 auto velocity = sampledVelocity + mBackgroundVelocity;
168 const GradientT& mGradient;
169 const Vec3GridT* mVelocity =
nullptr;
170 const ValueT& mBackgroundVelocity;
175template<
typename Vec3Gr
idT,
typename MaskT>
176struct SolveBoundaryOp
178 SolveBoundaryOp(
const Vec3GridT& velGrid,
const MaskT& domainGrid)
179 : mVoxelSize(domainGrid.voxelSize()[0])
181 , mDomainGrid(domainGrid)
184 void operator()(
const Coord& ijk,
const Coord& neighbor,
185 double& source,
double& diagonal)
const {
187 typename Vec3GridT::ConstAccessor velGridAccessor = mVelGrid.getAccessor();
188 const Coord diff = (ijk - neighbor);
190 if (velGridAccessor.isValueOn(ijk)) {
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];
201 const double mVoxelSize;
202 const Vec3GridT& mVelGrid;
203 const MaskT& mDomainGrid;
213template<
typename Gr
idT,
typename MaskT>
217 using MaskTreeT =
typename MaskT::TreeType;
219 if (!grid.hasUniformVoxels()) {
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());
234 mask->tree().topologyDifference(interior->tree());
240template<
typename Vec3T,
typename Gr
idT,
typename MaskT>
242 const GridT& collider,
244 const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
245 const Vec3T& backgroundVelocity)
247 using Vec3GridT =
typename GridT::template ValueConverter<Vec3T>::Type;
248 using TreeT =
typename Vec3GridT::TreeType;
249 using ValueT =
typename TreeT::ValueType;
252 using potential_flow_internal::ComputeNeumannVelocityOp;
257 !std::is_floating_point<typename GridT::TreeType::ValueType>::value) {
262 if (backgroundVelocity == zeroVal<Vec3T>() &&
263 (!boundaryVelocity || boundaryVelocity->empty())) {
264 auto neumann = Vec3GridT::create();
265 neumann->setTransform(collider.transform().copy());
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());
274 typename TreeT::Ptr neumannTree(
new TreeT(*boundary, zeroVal<ValueT>(),
TopologyCopy()));
275 neumannTree->voxelizeActiveTiles();
278 const typename GradientT::Ptr
gradient = tools::gradient(collider);
282 if (boundaryVelocity && !boundaryVelocity->empty()) {
283 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
284 neumannOp(*
gradient, *boundaryVelocity, backgroundVelocity);
285 leafManager.
foreach(neumannOp,
false);
288 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
289 neumannOp(*
gradient, backgroundVelocity);
290 leafManager.
foreach(neumannOp,
false);
294 tools::pruneInactive(*neumannTree);
296 typename Vec3GridT::Ptr neumann(Vec3GridT::create(neumannTree));
297 neumann->setTransform(collider.transform().copy());
303template<
typename Vec3Gr
idT,
typename MaskT,
typename InterrupterT>
304typename VectorToScalarGrid<Vec3GridT>::Ptr
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;
312 using potential_flow_internal::SolveBoundaryOp;
315 ScalarTreeT solveTree(domain.tree(), zeroVal<ScalarT>(),
TopologyCopy());
316 solveTree.voxelizeActiveTiles();
319 if (!interrupter) interrupter = &nullInterrupt;
322 SolveBoundaryOp<Vec3GridT, MaskT> solve(neumann, domain);
323 typename ScalarTreeT::Ptr potentialTree =
324 poisson::solveWithBoundaryConditions(solveTree, solve, state, *interrupter,
true);
326 auto potential = ScalarGridT::create(potentialTree);
327 potential->setTransform(domain.transform().copy());
333template<
typename Vec3Gr
idT>
334typename Vec3GridT::Ptr
336 const Vec3GridT& neumann,
337 const typename Vec3GridT::ValueType backgroundVelocity)
339 using Vec3T =
const typename Vec3GridT::ValueType;
340 using potential_flow_internal::extractOuterVoxelMask;
351 auto gradient = tools::gradient(potential);
355 auto applyNeumann = [&
gradient, &neumann] (
356 const MaskGrid::TreeType::LeafNodeType& leaf, size_t)
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);
371 leafManager.
foreach(applyNeumann);
375 if (backgroundVelocity != zeroVal<Vec3T>()) {
376 auto applyBackgroundVelocity = [&backgroundVelocity] (
377 typename Vec3GridT::TreeType::LeafNodeType& leaf, size_t)
379 for (
auto it = leaf.beginValueOn(); it; ++it) {
380 it.setValue(it.getValue() - backgroundVelocity);
385 leafManager2.
foreach(applyBackgroundVelocity);
397#ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
399#ifdef OPENVDB_INSTANTIATE_POTENTIALFLOW
403#define _FUNCTION(TreeT) \
404 Grid<TreeT>::Ptr createPotentialFlowNeumannVelocities(const FloatGrid&, const MaskGrid&, \
405 const Grid<TreeT>::ConstPtr, const TreeT::ValueType&)
409#define _FUNCTION(TreeT) \
410 VectorToScalarGrid<Grid<TreeT>>::Ptr computeScalarPotential(const MaskGrid&, const Grid<TreeT>&, \
411 math::pcg::State&, util::NullInterrupter*)
415#define _FUNCTION(TreeT) \
416 Grid<TreeT>::Ptr computePotentialFlow( \
417 const VectorToScalarGrid<Grid<TreeT>>::Type&, const Grid<TreeT>&, const TreeT::ValueType)
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:573
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition Types.h:683
Definition Exceptions.h:64
Definition Exceptions.h:65
Signed (x, y, z) 32-bit integer coordinates.
Definition Coord.h:25
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
Vec3< double > Vec3d
Definition Vec3.h:664
@ GRID_LEVEL_SET
Definition Types.h:455
Definition Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition Exceptions.h:74
Information about the state of a conjugate gradient solution.
Definition ConjGradient.h:46
Base class for interrupters.
Definition NullInterrupter.h:26
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212
#define OPENVDB_VEC3_TREE_INSTANTIATE(Function)
Definition version.h.in:159