10 #ifndef OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED 11 #define OPENVDB_TOOLS_POTENTIAL_FLOW_HAS_BEEN_INCLUDED 29 template<
typename VecGr
idT>
32 typename VecGridT::template ValueConverter<typename VecGridT::ValueType::value_type>::Type;
33 using Ptr =
typename Type::Ptr;
43 template<typename GridT, typename MaskT = typename GridT::template ValueConverter<ValueMask>::Type>
57 template<
typename Vec3T,
typename Gr
idT,
typename MaskT>
58 typename GridT::template ValueConverter<Vec3T>::Type::Ptr
60 const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
61 const Vec3T& backgroundVelocity);
75 template<
typename Vec3Gr
idT,
typename MaskT,
typename InterrupterT = util::NullInterrupter>
78 InterrupterT* interrupter =
nullptr);
87 template<
typename Vec3Gr
idT>
88 typename Vec3GridT::Ptr
90 const Vec3GridT& neumann,
91 const typename Vec3GridT::ValueType backgroundVelocity =
92 zeroVal<typename Vec3GridT::TreeType::ValueType>());
99 namespace potential_flow_internal {
104 template<
typename Gr
idT>
105 typename GridT::TreeType::template ValueConverter<ValueMask>::Type::Ptr
106 extractOuterVoxelMask(GridT& inGrid)
108 using MaskTreeT =
typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
110 typename MaskTreeT::Ptr boundaryMask(
new MaskTreeT(inGrid.tree(),
false,
TopologyCopy()));
120 template<
typename Vec3Gr
idT,
typename GradientT>
121 struct ComputeNeumannVelocityOp
123 using ValueT =
typename Vec3GridT::ValueType;
124 using VelocityAccessor =
typename Vec3GridT::ConstAccessor;
126 typename Vec3GridT::ConstAccessor,
BoxSampler>;
127 using GradientValueT =
typename GradientT::TreeType::ValueType;
129 ComputeNeumannVelocityOp(
const GradientT&
gradient,
130 const Vec3GridT& velocity,
131 const ValueT& backgroundVelocity)
132 : mGradient(gradient)
133 , mVelocity(&velocity)
134 , mBackgroundVelocity(backgroundVelocity) { }
136 ComputeNeumannVelocityOp(
const GradientT&
gradient,
137 const ValueT& backgroundVelocity)
138 : mGradient(gradient)
139 , mBackgroundVelocity(backgroundVelocity) { }
141 void operator()(
typename Vec3GridT::TreeType::LeafNodeType& leaf,
size_t)
const {
142 auto gradientAccessor = mGradient.getConstAccessor();
144 std::unique_ptr<VelocityAccessor> velocityAccessor;
145 std::unique_ptr<VelocitySamplerT> velocitySampler;
147 velocityAccessor.reset(
new VelocityAccessor(mVelocity->getConstAccessor()));
148 velocitySampler.reset(
new VelocitySamplerT(*velocityAccessor, mVelocity->transform()));
151 for (
auto it = leaf.beginValueOn(); it; ++it) {
152 Coord ijk = it.getCoord();
153 auto gradient = gradientAccessor.getValue(ijk);
155 const Vec3d xyz = mGradient.transform().indexToWorld(ijk);
156 const ValueT sampledVelocity = velocitySampler ?
157 velocitySampler->wsSample(xyz) : zeroVal<ValueT>();
158 auto velocity = sampledVelocity + mBackgroundVelocity;
169 const GradientT& mGradient;
170 const Vec3GridT* mVelocity =
nullptr;
171 const ValueT& mBackgroundVelocity;
176 template<
typename Vec3Gr
idT,
typename MaskT>
177 struct SolveBoundaryOp
179 SolveBoundaryOp(
const Vec3GridT& velGrid,
const MaskT& domainGrid)
180 : mVoxelSize(domainGrid.voxelSize()[0])
182 , mDomainGrid(domainGrid)
185 void operator()(
const Coord& ijk,
const Coord& neighbor,
186 double& source,
double& diagonal)
const {
188 typename Vec3GridT::ConstAccessor velGridAccessor = mVelGrid.getAccessor();
189 const Coord diff = (ijk - neighbor);
191 if (velGridAccessor.isValueOn(ijk)) {
192 const typename Vec3GridT::ValueType& sampleVel = velGridAccessor.getValue(ijk);
193 source += mVoxelSize*diff[0]*sampleVel[0];
194 source += mVoxelSize*diff[1]*sampleVel[1];
195 source += mVoxelSize*diff[2]*sampleVel[2];
202 const double mVoxelSize;
203 const Vec3GridT& mVelGrid;
204 const MaskT& mDomainGrid;
214 template<
typename Gr
idT,
typename MaskT>
218 using MaskTreeT =
typename MaskT::TreeType;
220 if (!grid.hasUniformVoxels()) {
228 typename MaskTreeT::Ptr maskTree(
new MaskTreeT(interior->tree(),
false,
TopologyCopy()));
229 typename MaskT::Ptr mask = MaskT::create(maskTree);
230 mask->setTransform(grid.transform().copy());
235 mask->tree().topologyDifference(interior->tree());
241 template<
typename Vec3T,
typename Gr
idT,
typename MaskT>
243 const GridT& collider,
245 const typename GridT::template ValueConverter<Vec3T>::Type::ConstPtr boundaryVelocity,
246 const Vec3T& backgroundVelocity)
248 using Vec3GridT =
typename GridT::template ValueConverter<Vec3T>::Type;
249 using TreeT =
typename Vec3GridT::TreeType;
250 using ValueT =
typename TreeT::ValueType;
253 using potential_flow_internal::ComputeNeumannVelocityOp;
263 if (backgroundVelocity == zeroVal<Vec3T>() &&
264 (!boundaryVelocity || boundaryVelocity->empty())) {
265 auto neumann = Vec3GridT::create();
266 neumann->setTransform(collider.transform().copy());
271 using MaskTreeT =
typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
272 typename MaskTreeT::Ptr boundary(
new MaskTreeT(domain.tree(),
false,
TopologyCopy()));
273 boundary->topologyIntersection(collider.tree());
275 typename TreeT::Ptr neumannTree(
new TreeT(*boundary, zeroVal<ValueT>(),
TopologyCopy()));
276 neumannTree->voxelizeActiveTiles();
283 if (boundaryVelocity && !boundaryVelocity->empty()) {
284 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
285 neumannOp(*gradient, *boundaryVelocity, backgroundVelocity);
286 leafManager.
foreach(neumannOp,
false);
289 ComputeNeumannVelocityOp<Vec3GridT, GradientT>
290 neumannOp(*gradient, backgroundVelocity);
291 leafManager.
foreach(neumannOp,
false);
297 typename Vec3GridT::Ptr neumann(Vec3GridT::create(neumannTree));
298 neumann->setTransform(collider.transform().copy());
304 template<
typename Vec3Gr
idT,
typename MaskT,
typename InterrupterT>
309 using ScalarT =
typename Vec3GridT::ValueType::value_type;
310 using ScalarTreeT =
typename Vec3GridT::TreeType::template ValueConverter<ScalarT>::Type;
311 using ScalarGridT =
typename Vec3GridT::template ValueConverter<ScalarT>::Type;
313 using potential_flow_internal::SolveBoundaryOp;
316 ScalarTreeT solveTree(domain.tree(), zeroVal<ScalarT>(),
TopologyCopy());
317 solveTree.voxelizeActiveTiles();
320 if (!interrupter) interrupter = &nullInterrupt;
323 SolveBoundaryOp<Vec3GridT, MaskT>
solve(neumann, domain);
324 typename ScalarTreeT::Ptr potentialTree =
327 auto potential = ScalarGridT::create(potentialTree);
328 potential->setTransform(domain.transform().copy());
334 template<
typename Vec3Gr
idT>
335 typename Vec3GridT::Ptr
337 const Vec3GridT& neumann,
338 const typename Vec3GridT::ValueType backgroundVelocity)
340 using Vec3T =
const typename Vec3GridT::ValueType;
341 using potential_flow_internal::extractOuterVoxelMask;
356 auto applyNeumann = [&
gradient, &neumann] (
357 const MaskGrid::TreeType::LeafNodeType& leaf, size_t)
359 typename Vec3GridT::Accessor gradientAccessor =
gradient->getAccessor();
360 typename Vec3GridT::ConstAccessor neumannAccessor = neumann.getAccessor();
361 for (
auto it = leaf.beginValueOn(); it; ++it) {
362 const Coord ijk = it.getCoord();
363 typename Vec3GridT::ValueType
value;
364 if (neumannAccessor.probeValue(ijk, value)) {
365 gradientAccessor.setValue(ijk, value);
372 leafManager.
foreach(applyNeumann);
376 if (backgroundVelocity != zeroVal<Vec3T>()) {
377 auto applyBackgroundVelocity = [&backgroundVelocity] (
378 typename Vec3GridT::TreeType::LeafNodeType& leaf, size_t)
380 for (
auto it = leaf.beginValueOn(); it; ++it) {
381 it.setValue(it.getValue() - backgroundVelocity);
386 leafManager2.
foreach(applyBackgroundVelocity);
398 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION 400 #ifdef OPENVDB_INSTANTIATE_POTENTIALFLOW 404 #define _FUNCTION(TreeT) \ 405 Grid<TreeT>::Ptr createPotentialFlowNeumannVelocities(const FloatGrid&, const MaskGrid&, \ 406 const Grid<TreeT>::ConstPtr, const TreeT::ValueType&) 410 #define _FUNCTION(TreeT) \ 411 VectorToScalarGrid<Grid<TreeT>>::Ptr computeScalarPotential(const MaskGrid&, const Grid<TreeT>&, \ 412 math::pcg::State&, util::NullInterrupter*) 416 #define _FUNCTION(TreeT) \ 417 Grid<TreeT>::Ptr computePotentialFlow( \ 418 const VectorToScalarGrid<Grid<TreeT>>::Type&, const Grid<TreeT>&, const TreeT::ValueType) 422 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION 429 #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:1612
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Base class for interrupters.
Definition: NullInterrupter.h:25
Information about the state of a conjugate gradient solution.
Definition: ConjGradient.h:46
Apply an operator to an input grid to produce an output grid with the same active voxel topology but ...
Definition: Exceptions.h:64
Definition: Exceptions.h:65
Construct boolean mask grids from grids of arbitrary type.
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:84
Implementation of morphological dilation and erosion.
Definition: Exceptions.h:13
ValueT value
Definition: GridBuilder.h:1287
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
SharedPtr< Grid > Ptr
Definition: Grid.h:575
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:644
#define OPENVDB_VEC3_TREE_INSTANTIATE(Function)
Definition: version.h.in:149
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:202