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;
124 typename Vec3GridT::ConstAccessor,
BoxSampler>;
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);
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:560
Definition: openvdb/Exceptions.h:64
Definition: openvdb/Exceptions.h:65
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:26
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
@ GRID_LEVEL_SET
Definition: openvdb/Types.h:333
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
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.in:116
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:178