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()));
117 template<
typename Vec3Gr
idT,
typename GradientT>
120 using ValueT =
typename Vec3GridT::ValueType;
127 const Vec3GridT& velocity,
128 const ValueT& backgroundVelocity)
130 , mVelocity(&velocity)
131 , mBackgroundVelocity(backgroundVelocity) { }
134 const ValueT& backgroundVelocity)
136 , mBackgroundVelocity(backgroundVelocity) { }
138 void operator()(
typename Vec3GridT::TreeType::LeafNodeType& leaf,
size_t)
const {
139 auto gradientAccessor = mGradient.getConstAccessor();
141 std::unique_ptr<VelocityAccessor> velocityAccessor;
142 std::unique_ptr<VelocitySamplerT> velocitySampler;
144 velocityAccessor.reset(
new VelocityAccessor(mVelocity->getConstAccessor()));
145 velocitySampler.reset(
new VelocitySamplerT(*velocityAccessor, mVelocity->transform()));
148 for (
auto it = leaf.beginValueOn(); it; ++it) {
149 Coord ijk = it.getCoord();
150 auto gradient = gradientAccessor.getValue(ijk);
152 const Vec3d xyz = mGradient.transform().indexToWorld(ijk);
153 const ValueT sampledVelocity = velocitySampler ?
154 velocitySampler->wsSample(xyz) : zeroVal<ValueT>();
155 auto velocity = sampledVelocity + mBackgroundVelocity;
166 const GradientT& mGradient;
167 const Vec3GridT* mVelocity =
nullptr;
168 const ValueT& mBackgroundVelocity;
173 template<
typename Vec3Gr
idT,
typename MaskT>
177 : mVoxelSize(domainGrid.voxelSize()[0])
179 , mDomainGrid(domainGrid)
183 double& source,
double& diagonal)
const {
185 typename Vec3GridT::ConstAccessor velGridAccessor = mVelGrid.getAccessor();
186 const Coord diff = (ijk - neighbor);
188 if (velGridAccessor.isValueOn(ijk)) {
189 const typename Vec3GridT::ValueType& sampleVel = velGridAccessor.getValue(ijk);
190 source += mVoxelSize*diff[0]*sampleVel[0];
191 source += mVoxelSize*diff[1]*sampleVel[1];
192 source += mVoxelSize*diff[2]*sampleVel[2];
210 template<
typename Gr
idT,
typename MaskT>
211 inline typename MaskT::Ptr
214 using MaskTreeT =
typename MaskT::TreeType;
216 if (!grid.hasUniformVoxels()) {
224 typename MaskTreeT::Ptr maskTree(
new MaskTreeT(interior->tree(),
false,
TopologyCopy()));
225 typename MaskT::Ptr mask = MaskT::create(maskTree);
226 mask->setTransform(grid.transform().copy());
231 mask->tree().topologyDifference(interior->tree());
237 template<
typename Vec3T,
typename Gr
idT,
typename MaskT>
239 const GridT& collider,
241 const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
242 const Vec3T& backgroundVelocity)
244 using Vec3GridT =
typename GridT::template ValueConverter<Vec3T>::Type;
245 using TreeT =
typename Vec3GridT::TreeType;
246 using ValueT =
typename TreeT::ValueType;
254 !std::is_floating_point<typename GridT::TreeType::ValueType>::value) {
259 if (backgroundVelocity == zeroVal<Vec3T>() &&
260 (!boundaryVelocity || boundaryVelocity->empty())) {
261 auto neumann = Vec3GridT::create();
262 neumann->setTransform(collider.transform().copy());
267 using MaskTreeT =
typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
268 typename MaskTreeT::Ptr boundary(
new MaskTreeT(domain.tree(),
false,
TopologyCopy()));
269 boundary->topologyIntersection(collider.tree());
271 typename TreeT::Ptr neumannTree(
new TreeT(*boundary, zeroVal<ValueT>(),
TopologyCopy()));
272 neumannTree->voxelizeActiveTiles();
279 if (boundaryVelocity && !boundaryVelocity->empty()) {
280 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
281 neumannOp(*
gradient, *boundaryVelocity, backgroundVelocity);
282 leafManager.
foreach(neumannOp,
false);
285 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
286 neumannOp(*
gradient, backgroundVelocity);
287 leafManager.
foreach(neumannOp,
false);
293 typename Vec3GridT::Ptr neumann(Vec3GridT::create(neumannTree));
294 neumann->setTransform(collider.transform().copy());
300 template<
typename Vec3Gr
idT,
typename MaskT,
typename InterrupterT>
301 inline typename VectorToScalarGrid<Vec3GridT>::Ptr
305 using ScalarT =
typename Vec3GridT::ValueType::value_type;
306 using ScalarTreeT =
typename Vec3GridT::TreeType::template ValueConverter<ScalarT>::Type;
307 using ScalarGridT =
typename Vec3GridT::template ValueConverter<ScalarT>::Type;
312 ScalarTreeT solveTree(domain.tree(), zeroVal<ScalarT>(),
TopologyCopy());
313 solveTree.voxelizeActiveTiles();
316 if (!interrupter) interrupter = &nullInterrupt;
319 SolveBoundaryOp<Vec3GridT, MaskT>
solve(neumann, domain);
320 typename ScalarTreeT::Ptr potentialTree =
323 auto potential = ScalarGridT::create(potentialTree);
324 potential->setTransform(domain.transform().copy());
330 template<
typename Vec3Gr
idT>
331 inline typename Vec3GridT::Ptr
333 const Vec3GridT& neumann,
334 const typename Vec3GridT::ValueType backgroundVelocity)
336 using Vec3T =
const typename Vec3GridT::ValueType;
337 using potential_flow_internal::extractOuterVoxelMask;
352 auto applyNeumann = [&
gradient, &neumann] (
353 const MaskGrid::TreeType::LeafNodeType& leaf, size_t)
355 typename Vec3GridT::Accessor gradientAccessor =
gradient->getAccessor();
356 typename Vec3GridT::ConstAccessor neumannAccessor = neumann.getAccessor();
357 for (
auto it = leaf.beginValueOn(); it; ++it) {
358 const Coord ijk = it.getCoord();
359 typename Vec3GridT::ValueType value;
360 if (neumannAccessor.probeValue(ijk, value)) {
361 gradientAccessor.setValue(ijk, value);
368 leafManager.
foreach(applyNeumann);
372 if (backgroundVelocity != zeroVal<Vec3T>()) {
373 auto applyBackgroundVelocity = [&backgroundVelocity] (
374 typename Vec3GridT::TreeType::LeafNodeType& leaf, size_t)
376 for (
auto it = leaf.beginValueOn(); it; ++it) {
377 it.setValue(it.getValue() - backgroundVelocity);
382 leafManager2.
foreach(applyBackgroundVelocity);
395 #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...
Vec3< double > Vec3d
Definition: Vec3.h:662
SharedPtr< Grid > Ptr
Definition: Grid.h:574
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:82
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.
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
Definition: Exceptions.h:65
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
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:13
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.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:82
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:681
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:496
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
Definition: Exceptions.h:64