10 #ifndef OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED 11 #define OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED 28 template<
typename VecGr
idT>
31 typename VecGridT::template ValueConverter<typename VecGridT::ValueType::value_type>::Type;
32 using Ptr =
typename Type::Ptr;
42 template<typename GridT, typename MaskT = typename GridT::template ValueConverter<ValueMask>::Type>
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,
90 const typename Vec3GridT::ValueType backgroundVelocity =
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)
107 using MaskTreeT =
typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
109 typename MaskTreeT::Ptr boundaryMask(
new MaskTreeT(inGrid.tree(),
false,
TopologyCopy()));
119 template<
typename Vec3Gr
idT,
typename GradientT>
120 struct ComputeNeumannVelocityOp
122 using ValueT =
typename Vec3GridT::ValueType;
123 using VelocityAccessor =
typename Vec3GridT::ConstAccessor;
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;
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)) {
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;
213 template<
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());
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)
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();
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>
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 =
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,
337 const typename Vec3GridT::ValueType backgroundVelocity)
339 using Vec3T =
const typename Vec3GridT::ValueType;
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();
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) 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 acti...
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Apply an operator to an input grid to produce an output grid with the same active voxel topology but ...
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
Construct boolean mask grids from grids of arbitrary type.
Base class for interrupters.
Definition: NullInterrupter.h:25
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:1608
Definition: Exceptions.h:65
Implementation of morphological dilation and erosion.
SharedPtr< Grid > Ptr
Definition: Grid.h:573
Definition: Exceptions.h:13
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:484
Definition: Exceptions.h:64
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:683
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
#define OPENVDB_VEC3_TREE_INSTANTIATE(Function)
Definition: version.h.in:164
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218