37 #ifndef OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED 38 #define OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED 40 #include <openvdb/openvdb.h> 55 template<
typename VecGr
idT>
58 typename VecGridT::template ValueConverter<typename VecGridT::ValueType::value_type>::Type;
59 using Ptr =
typename Type::Ptr;
69 template<typename GridT, typename MaskT = typename GridT::template ValueConverter<ValueMask>::Type>
70 inline typename MaskT::Ptr
83 template<
typename Vec3T,
typename Gr
idT,
typename MaskT>
84 inline typename GridT::template ValueConverter<Vec3T>::Type::Ptr
86 const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
87 const Vec3T& backgroundVelocity);
101 template<
typename Vec3Gr
idT,
typename MaskT,
typename InterrupterT = util::NullInterrupter>
104 InterrupterT* interrupter =
nullptr);
113 template<
typename Vec3Gr
idT>
114 inline typename Vec3GridT::Ptr
116 const Vec3GridT& neumann,
117 const typename Vec3GridT::ValueType backgroundVelocity =
118 zeroVal<typename Vec3GridT::TreeType::ValueType>());
124 namespace potential_flow_internal {
129 template<
typename Gr
idT>
130 inline typename GridT::TreeType::template ValueConverter<ValueMask>::Type::Ptr
131 extractOuterVoxelMask(GridT& inGrid)
133 using MaskTreeT =
typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
135 typename MaskTreeT::Ptr boundaryMask(
new MaskTreeT(inGrid.tree(),
false,
TopologyCopy()));
138 boundaryMask->topologyDifference(*interiorMask);
144 template<
typename Vec3Gr
idT,
typename GradientT>
147 using ValueT =
typename Vec3GridT::ValueType;
154 const Vec3GridT& velocity,
155 const ValueT& backgroundVelocity)
156 : mGradient(gradient)
157 , mVelocity(&velocity)
158 , mBackgroundVelocity(backgroundVelocity) { }
161 const ValueT& backgroundVelocity)
162 : mGradient(gradient)
163 , mBackgroundVelocity(backgroundVelocity) { }
165 void operator()(
typename Vec3GridT::TreeType::LeafNodeType& leaf,
size_t)
const {
166 auto gradientAccessor = mGradient.getConstAccessor();
168 std::unique_ptr<VelocityAccessor> velocityAccessor;
169 std::unique_ptr<VelocitySamplerT> velocitySampler;
171 velocityAccessor.reset(
new VelocityAccessor(mVelocity->getConstAccessor()));
172 velocitySampler.reset(
new VelocitySamplerT(*velocityAccessor, mVelocity->transform()));
175 for (
auto it = leaf.beginValueOn(); it; ++it) {
176 Coord ijk = it.getCoord();
177 auto gradient = gradientAccessor.getValue(ijk);
179 const Vec3d xyz = mGradient.transform().indexToWorld(ijk);
180 const ValueT sampledVelocity = velocitySampler ?
181 velocitySampler->wsSample(xyz) : zeroVal<ValueT>();
182 auto velocity = sampledVelocity + mBackgroundVelocity;
193 const GradientT& mGradient;
194 const Vec3GridT* mVelocity =
nullptr;
195 const ValueT& mBackgroundVelocity;
200 template<
typename Vec3Gr
idT,
typename MaskT>
204 : mVoxelSize(domainGrid.voxelSize()[0])
206 , mDomainGrid(domainGrid)
210 double& source,
double& diagonal)
const {
212 typename Vec3GridT::ConstAccessor velGridAccessor = mVelGrid.getAccessor();
213 const Coord diff = (ijk - neighbor);
215 if (velGridAccessor.isValueOn(ijk)) {
216 const typename Vec3GridT::ValueType& sampleVel = velGridAccessor.getValue(ijk);
217 source += mVoxelSize*diff[0]*sampleVel[0];
218 source += mVoxelSize*diff[1]*sampleVel[1];
219 source += mVoxelSize*diff[2]*sampleVel[2];
237 template<
typename Gr
idT,
typename MaskT>
238 inline typename MaskT::Ptr
241 using MaskTreeT =
typename MaskT::TreeType;
243 if (!grid.hasUniformVoxels()) {
251 typename MaskTreeT::Ptr maskTree(
new MaskTreeT(interior->tree(),
false,
TopologyCopy()));
252 typename MaskT::Ptr mask = MaskT::create(maskTree);
253 mask->setTransform(grid.transform().copy());
258 mask->tree().topologyDifference(interior->tree());
264 template<
typename Vec3T,
typename Gr
idT,
typename MaskT>
266 const GridT& collider,
268 const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
269 const Vec3T& backgroundVelocity)
271 using Vec3GridT =
typename GridT::template ValueConverter<Vec3T>::Type;
272 using TreeT =
typename Vec3GridT::TreeType;
273 using ValueT =
typename TreeT::ValueType;
281 !std::is_floating_point<typename GridT::TreeType::ValueType>::value) {
286 if (backgroundVelocity == zeroVal<Vec3T>() &&
287 (!boundaryVelocity || boundaryVelocity->empty())) {
288 auto neumann = Vec3GridT::create();
289 neumann->setTransform(collider.transform().copy());
294 using MaskTreeT =
typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
295 typename MaskTreeT::Ptr boundary(
new MaskTreeT(domain.tree(),
false,
TopologyCopy()));
296 boundary->topologyIntersection(collider.tree());
298 typename TreeT::Ptr neumannTree(
new TreeT(*boundary, zeroVal<ValueT>(),
TopologyCopy()));
299 neumannTree->voxelizeActiveTiles();
306 if (boundaryVelocity && !boundaryVelocity->empty()) {
307 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
308 neumannOp(*gradient, *boundaryVelocity, backgroundVelocity);
309 leafManager.
foreach(neumannOp,
false);
312 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
313 neumannOp(*gradient, backgroundVelocity);
314 leafManager.
foreach(neumannOp,
false);
320 typename Vec3GridT::Ptr neumann(Vec3GridT::create(neumannTree));
321 neumann->setTransform(collider.transform().copy());
327 template<
typename Vec3Gr
idT,
typename MaskT,
typename InterrupterT>
332 using ScalarT =
typename Vec3GridT::ValueType::value_type;
333 using ScalarTreeT =
typename Vec3GridT::TreeType::template ValueConverter<ScalarT>::Type;
334 using ScalarGridT =
typename Vec3GridT::template ValueConverter<ScalarT>::Type;
339 ScalarTreeT solveTree(domain.tree(), zeroVal<ScalarT>(),
TopologyCopy());
340 solveTree.voxelizeActiveTiles();
343 if (!interrupter) interrupter = &nullInterrupt;
346 SolveBoundaryOp<Vec3GridT, MaskT>
solve(neumann, domain);
347 typename ScalarTreeT::Ptr potentialTree =
350 auto potential = ScalarGridT::create(potentialTree);
351 potential->setTransform(domain.transform().copy());
357 template<
typename Vec3Gr
idT>
358 inline typename Vec3GridT::Ptr
360 const Vec3GridT& neumann,
361 const typename Vec3GridT::ValueType backgroundVelocity)
363 using Vec3T =
const typename Vec3GridT::ValueType;
364 using potential_flow_internal::extractOuterVoxelMask;
379 auto applyNeumann = [&
gradient, &neumann] (
380 const MaskGrid::TreeType::LeafNodeType& leaf, size_t)
382 typename Vec3GridT::Accessor gradientAccessor =
gradient->getAccessor();
383 typename Vec3GridT::ConstAccessor neumannAccessor = neumann.getAccessor();
384 for (
auto it = leaf.beginValueOn(); it; ++it) {
385 const Coord ijk = it.getCoord();
386 typename Vec3GridT::ValueType value;
387 if (neumannAccessor.probeValue(ijk, value)) {
388 gradientAccessor.setValue(ijk, value);
395 leafManager.
foreach(applyNeumann);
399 if (backgroundVelocity != zeroVal<Vec3T>()) {
400 auto applyBackgroundVelocity = [&backgroundVelocity] (
401 typename Vec3GridT::TreeType::LeafNodeType& leaf, size_t)
403 for (
auto it = leaf.beginValueOn(); it; ++it) {
404 it.setValue(it.getValue() - backgroundVelocity);
409 leafManager2.
foreach(applyBackgroundVelocity);
422 #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...
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:1639
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:109
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:51
Apply an operator to an input grid to produce an output grid with the same active voxel topology but ...
Definition: Exceptions.h:91
Construct boolean mask grids from grids of arbitrary type.
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:73
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:524
Vec3< double > Vec3d
Definition: Vec3.h:679
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:136
Definition: Exceptions.h:92
SharedPtr< Grid > Ptr
Definition: Grid.h:502
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:40
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:110
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:188
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:518
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52