10 #ifndef OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED
11 #define OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED
28 template<
typename VecGr
idT>
32 using Ptr =
typename Type::Ptr;
56 template<
typename Vec3T,
typename Gr
idT,
typename MaskT>
57 typename GridT::template ValueConverter<Vec3T>::Type::Ptr
59 const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
60 const Vec3T& backgroundVelocity);
74 template<
typename Vec3Gr
idT,
typename MaskT,
typename InterrupterT = util::NullInterrupter>
77 InterrupterT* interrupter =
nullptr);
86 template<
typename Vec3Gr
idT>
87 typename Vec3GridT::Ptr
89 const Vec3GridT& neumann,
91 zeroVal<typename Vec3GridT::TreeType::ValueType>());
98 namespace potential_flow_internal {
103 template<
typename Gr
idT>
104 typename GridT::TreeType::template ValueConverter<ValueMask>::Type::Ptr
105 extractOuterVoxelMask(GridT& inGrid)
109 typename MaskTreeT::Ptr boundaryMask(
new MaskTreeT(inGrid.tree(),
false,
TopologyCopy()));
119 template<
typename Vec3Gr
idT,
typename GradientT>
120 struct ComputeNeumannVelocityOp
123 using VelocityAccessor =
typename Vec3GridT::ConstAccessor;
124 using VelocitySamplerT = GridSampler<
125 typename Vec3GridT::ConstAccessor, BoxSampler>;
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;
175 template<
typename Vec3Gr
idT,
typename MaskT>
176 struct 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)) {
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;
213 template<
typename Gr
idT,
typename MaskT>
217 using MaskTreeT =
typename MaskT::TreeType;
219 if (!grid.hasUniformVoxels()) {
220 OPENVDB_THROW(ValueError,
"Transform must have uniform voxels for Potential Flow mask.");
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());
240 template<
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)
248 using TreeT =
typename Vec3GridT::TreeType;
252 using potential_flow_internal::ComputeNeumannVelocityOp;
258 OPENVDB_THROW(TypeError,
"Potential Flow expecting the collider to be a level set.");
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();
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);
296 typename Vec3GridT::Ptr neumann(Vec3GridT::create(neumannTree));
297 neumann->setTransform(collider.transform().copy());
303 template<
typename Vec3Gr
idT,
typename MaskT,
typename InterrupterT>
304 typename VectorToScalarGrid<Vec3GridT>::Ptr
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 =
326 auto potential = ScalarGridT::create(potentialTree);
327 potential->setTransform(domain.transform().copy());
333 template<
typename Vec3Gr
idT>
334 typename Vec3GridT::Ptr
336 const Vec3GridT& neumann,
340 using potential_flow_internal::extractOuterVoxelMask;
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();
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)
421 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
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...
GLsizei const GLfloat * value
#define OPENVDB_USE_VERSION_NAMESPACE
Information about the state of a conjugate gradient solution.
Base class for interrupters.
GLsizei GLsizei GLchar * source
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
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.
#define OPENVDB_VEC3_TREE_INSTANTIATE(Function)
Construct boolean mask grids from grids of arbitrary type.
Implementation of morphological dilation and erosion.
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Apply an operator to an input grid to produce an output grid with the same active voxel topology but ...
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
#define OPENVDB_THROW(exception, message)