Go to the documentation of this file.
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>
43 inline typename MaskT::Ptr
56 template<
typename Vec3T,
typename Gr
idT,
typename MaskT>
57 inline 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 inline typename Vec3GridT::Ptr
89 const Vec3GridT& neumann,
90 const typename Vec3GridT::ValueType backgroundVelocity =
91 zeroVal<typename Vec3GridT::TreeType::ValueType>());
97 namespace potential_flow_internal {
102 template<
typename Gr
idT>
103 inline typename GridT::TreeType::template ValueConverter<ValueMask>::Type::Ptr
104 extractOuterVoxelMask(GridT& inGrid)
106 using MaskTreeT =
typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
108 typename MaskTreeT::Ptr boundaryMask(
new MaskTreeT(inGrid.tree(),
false,
TopologyCopy()));
118 template<
typename Vec3Gr
idT,
typename GradientT>
121 using ValueT =
typename Vec3GridT::ValueType;
128 const Vec3GridT& velocity,
129 const ValueT& backgroundVelocity)
131 , mVelocity(&velocity)
132 , mBackgroundVelocity(backgroundVelocity) { }
135 const ValueT& backgroundVelocity)
137 , mBackgroundVelocity(backgroundVelocity) { }
139 void operator()(
typename Vec3GridT::TreeType::LeafNodeType& leaf,
size_t)
const {
140 auto gradientAccessor = mGradient.getConstAccessor();
142 std::unique_ptr<VelocityAccessor> velocityAccessor;
143 std::unique_ptr<VelocitySamplerT> velocitySampler;
145 velocityAccessor.reset(
new VelocityAccessor(mVelocity->getConstAccessor()));
146 velocitySampler.reset(
new VelocitySamplerT(*velocityAccessor, mVelocity->transform()));
149 for (
auto it = leaf.beginValueOn(); it; ++it) {
150 Coord ijk = it.getCoord();
151 auto gradient = gradientAccessor.getValue(ijk);
153 const Vec3d xyz = mGradient.transform().indexToWorld(ijk);
154 const ValueT sampledVelocity = velocitySampler ?
155 velocitySampler->wsSample(xyz) : zeroVal<ValueT>();
156 auto velocity = sampledVelocity + mBackgroundVelocity;
167 const GradientT& mGradient;
168 const Vec3GridT* mVelocity =
nullptr;
169 const ValueT& mBackgroundVelocity;
174 template<
typename Vec3Gr
idT,
typename MaskT>
178 : mVoxelSize(domainGrid.voxelSize()[0])
180 , mDomainGrid(domainGrid)
184 double& source,
double& diagonal)
const {
186 typename Vec3GridT::ConstAccessor velGridAccessor = mVelGrid.getAccessor();
187 const Coord diff = (ijk - neighbor);
189 if (velGridAccessor.isValueOn(ijk)) {
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];
211 template<
typename Gr
idT,
typename MaskT>
212 inline typename MaskT::Ptr
215 using MaskTreeT =
typename MaskT::TreeType;
217 if (!grid.hasUniformVoxels()) {
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());
232 mask->tree().topologyDifference(interior->tree());
238 template<
typename Vec3T,
typename Gr
idT,
typename MaskT>
240 const GridT& collider,
242 const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
243 const Vec3T& backgroundVelocity)
245 using Vec3GridT =
typename GridT::template ValueConverter<Vec3T>::Type;
246 using TreeT =
typename Vec3GridT::TreeType;
247 using ValueT =
typename TreeT::ValueType;
255 !std::is_floating_point<typename GridT::TreeType::ValueType>::value) {
260 if (backgroundVelocity == zeroVal<Vec3T>() &&
261 (!boundaryVelocity || boundaryVelocity->empty())) {
262 auto neumann = Vec3GridT::create();
263 neumann->setTransform(collider.transform().copy());
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());
272 typename TreeT::Ptr neumannTree(
new TreeT(*boundary, zeroVal<ValueT>(),
TopologyCopy()));
273 neumannTree->voxelizeActiveTiles();
280 if (boundaryVelocity && !boundaryVelocity->empty()) {
281 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
282 neumannOp(*
gradient, *boundaryVelocity, backgroundVelocity);
283 leafManager.
foreach(neumannOp,
false);
286 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
287 neumannOp(*
gradient, backgroundVelocity);
288 leafManager.
foreach(neumannOp,
false);
294 typename Vec3GridT::Ptr neumann(Vec3GridT::create(neumannTree));
295 neumann->setTransform(collider.transform().copy());
301 template<
typename Vec3Gr
idT,
typename MaskT,
typename InterrupterT>
302 inline typename VectorToScalarGrid<Vec3GridT>::Ptr
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;
313 ScalarTreeT solveTree(domain.tree(), zeroVal<ScalarT>(),
TopologyCopy());
314 solveTree.voxelizeActiveTiles();
317 if (!interrupter) interrupter = &nullInterrupt;
320 SolveBoundaryOp<Vec3GridT, MaskT>
solve(neumann, domain);
321 typename ScalarTreeT::Ptr potentialTree =
324 auto potential = ScalarGridT::create(potentialTree);
325 potential->setTransform(domain.transform().copy());
331 template<
typename Vec3Gr
idT>
332 inline typename Vec3GridT::Ptr
334 const Vec3GridT& neumann,
335 const typename Vec3GridT::ValueType backgroundVelocity)
337 using Vec3T =
const typename Vec3GridT::ValueType;
338 using potential_flow_internal::extractOuterVoxelMask;
353 auto applyNeumann = [&
gradient, &neumann] (
354 const MaskGrid::TreeType::LeafNodeType& leaf, size_t)
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);
369 leafManager.
foreach(applyNeumann);
373 if (backgroundVelocity != zeroVal<Vec3T>()) {
374 auto applyBackgroundVelocity = [&backgroundVelocity] (
375 typename Vec3GridT::TreeType::LeafNodeType& leaf, size_t)
377 for (
auto it = leaf.beginValueOn(); it; ++it) {
378 it.setValue(it.getValue() - backgroundVelocity);
383 leafManager2.
foreach(applyBackgroundVelocity);
396 #endif // OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED
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
#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
Solve Poisson's equation ∇2x = b for x, where b is a vector comprising the values of all of the activ...
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:84
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:25
#define OPENVDB_THROW(exception, message)
Definition: openvdb/Exceptions.h:74
Definition: openvdb/Exceptions.h:65
Apply an operator to an input grid to produce an output grid with the same active voxel topology but ...
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
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
Construct boolean mask grids from grids of arbitrary type.
@ GRID_LEVEL_SET
Definition: openvdb/Types.h:333
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: openvdb/Types.h:560
Vec3< double > Vec3d
Definition: Vec3.h:668
Definition: openvdb/Exceptions.h:64
Implementation of morphological dilation and erosion.
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
Definition: openvdb/Exceptions.h:13
SharedPtr< Grid > Ptr
Definition: Grid.h:579