OpenVDB 10.0.1
Loading...
Searching...
No Matches
PoissonSolver.h
Go to the documentation of this file.
1// Copyright Contributors to the OpenVDB Project
2// SPDX-License-Identifier: MPL-2.0
3//
4/// @file PoissonSolver.h
5///
6/// @authors D.J. Hill, Peter Cucka
7///
8/// @brief Solve Poisson's equation &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i>
9/// for <i>x</i>, where @e b is a vector comprising the values of all of the active voxels
10/// in a grid.
11///
12/// @par Example:
13/// Solve for the pressure in a cubic tank of liquid, assuming uniform boundary conditions:
14/// @code
15/// FloatTree source(/*background=*/0.0f);
16/// // Activate voxels to indicate that they contain liquid.
17/// source.fill(CoordBBox(Coord(0, -10, 0), Coord(10, 0, 10)), /*value=*/0.0f);
18///
19/// math::pcg::State state = math::pcg::terminationDefaults<float>();
20/// FloatTree::Ptr solution = tools::poisson::solve(source, state);
21/// @endcode
22///
23/// @par Example:
24/// Solve for the pressure, <i>P</i>, in a cubic tank of liquid that is open at the top.
25/// Boundary conditions are <i>P</i>&nbsp;=&nbsp;0 at the top,
26/// &part;<i>P</i>/&part;<i>y</i>&nbsp;=&nbsp;&minus;1 at the bottom
27/// and &part;<i>P</i>/&part;<i>x</i>&nbsp;=&nbsp;0 at the sides:
28/// <pre>
29/// P = 0
30/// +--------+ (N,0,N)
31/// /| /|
32/// (0,0,0) +--------+ |
33/// | | | | dP/dx = 0
34/// dP/dx = 0 | +------|-+
35/// |/ |/
36/// (0,-N,0) +--------+ (N,-N,0)
37/// dP/dy = -1
38/// </pre>
39/// @code
40/// const int N = 10;
41/// DoubleTree source(/*background=*/0.0);
42/// // Activate voxels to indicate that they contain liquid.
43/// source.fill(CoordBBox(Coord(0, -N, 0), Coord(N, 0, N)), /*value=*/0.0);
44///
45/// auto boundary = [](const openvdb::Coord& ijk, const openvdb::Coord& neighbor,
46/// double& source, double& diagonal)
47/// {
48/// if (neighbor.x() == ijk.x() && neighbor.z() == ijk.z()) {
49/// if (neighbor.y() < ijk.y()) source -= 1.0;
50/// else diagonal -= 1.0;
51/// }
52/// };
53///
54/// math::pcg::State state = math::pcg::terminationDefaults<double>();
55/// util::NullInterrupter interrupter;
56///
57/// DoubleTree::Ptr solution = tools::poisson::solveWithBoundaryConditions(
58/// source, boundary, state, interrupter);
59/// @endcode
60
61#ifndef OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
62#define OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
63
64#include <openvdb/Types.h>
67#include <openvdb/tree/Tree.h>
69#include "Morphology.h" // for erodeActiveValues
70#include <openvdb/openvdb.h>
71
72
73namespace openvdb {
75namespace OPENVDB_VERSION_NAME {
76namespace tools {
77namespace poisson {
78
79// This type should be at least as wide as math::pcg::SizeType.
80using VIndex = Int32;
81
82/// The type of a matrix used to represent a three-dimensional %Laplacian operator
84
85
86//@{
87/// @brief Solve &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i> for <i>x</i>,
88/// where @e b is a vector comprising the values of all of the active voxels
89/// in the input tree.
90/// @return a new tree, with the same active voxel topology as the input tree,
91/// whose voxel values are the elements of the solution vector <i>x</i>.
92/// @details On input, the State object should specify convergence criteria
93/// (minimum error and maximum number of iterations); on output, it gives
94/// the actual termination conditions.
95/// @details The solution is computed using the conjugate gradient method
96/// with (where possible) incomplete Cholesky preconditioning, falling back
97/// to Jacobi preconditioning.
98/// @sa solveWithBoundaryConditions
99template<typename TreeType>
100typename TreeType::Ptr
101solve(const TreeType&, math::pcg::State&, bool staggered = false);
102
103template<typename TreeType, typename Interrupter>
104typename TreeType::Ptr
105solve(const TreeType&, math::pcg::State&, Interrupter&, bool staggered = false);
106//@}
107
108
109//@{
110/// @brief Solve &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i> for <i>x</i>
111/// with user-specified boundary conditions, where @e b is a vector comprising
112/// the values of all of the active voxels in the input tree or domain mask if provided
113/// @return a new tree, with the same active voxel topology as the input tree,
114/// whose voxel values are the elements of the solution vector <i>x</i>.
115/// @details On input, the State object should specify convergence criteria
116/// (minimum error and maximum number of iterations); on output, it gives
117/// the actual termination conditions.
118/// @details The solution is computed using the conjugate gradient method with
119/// the specified type of preconditioner (default: incomplete Cholesky),
120/// falling back to Jacobi preconditioning if necessary.
121/// @details Each thread gets its own copy of the BoundaryOp, which should be
122/// a functor of the form
123/// @code
124/// struct BoundaryOp {
125/// using ValueType = LaplacianMatrix::ValueType;
126/// void operator()(
127/// const Coord& ijk, // coordinates of a boundary voxel
128/// const Coord& ijkNeighbor, // coordinates of an exterior neighbor of ijk
129/// ValueType& source, // element of b corresponding to ijk
130/// ValueType& diagonal // element of Laplacian matrix corresponding to ijk
131/// ) const;
132/// };
133/// @endcode
134/// The functor is called for each of the exterior neighbors of each boundary voxel @ijk,
135/// and it must specify a boundary condition for @ijk by modifying one or both of two
136/// provided values: the entry in the source vector @e b corresponding to @ijk and
137/// the weighting coefficient for @ijk in the Laplacian operator matrix.
138///
139/// @sa solve
140template<typename TreeType, typename BoundaryOp, typename Interrupter>
141typename TreeType::Ptr
143 const TreeType&,
144 const BoundaryOp&,
146 Interrupter&,
147 bool staggered = false);
148
149template<
150 typename PreconditionerType,
151 typename TreeType,
152 typename BoundaryOp,
153 typename Interrupter>
154typename TreeType::Ptr
156 const TreeType&,
157 const BoundaryOp&,
159 Interrupter&,
160 bool staggered = false);
161
162template<
163 typename PreconditionerType,
164 typename TreeType,
165 typename DomainTreeType,
166 typename BoundaryOp,
167 typename Interrupter>
168typename TreeType::Ptr
170 const TreeType&,
171 const DomainTreeType&,
172 const BoundaryOp&,
174 Interrupter&,
175 bool staggered = false);
176//@}
177
178
179/// @name Low-level functions
180//@{
181// The following are low-level routines that can be used to assemble custom solvers.
182
183/// @brief Overwrite each active voxel in the given scalar tree
184/// with a sequential index, starting from zero.
185template<typename VIndexTreeType>
186void populateIndexTree(VIndexTreeType&);
187
188/// @brief Iterate over the active voxels of the input tree and for each one
189/// assign its index in the iteration sequence to the corresponding voxel
190/// of an integer-valued output tree.
191template<typename TreeType>
192typename TreeType::template ValueConverter<VIndex>::Type::Ptr
193createIndexTree(const TreeType&);
194
195
196/// @brief Return a vector of the active voxel values of the scalar-valued @a source tree.
197/// @details The <i>n</i>th element of the vector corresponds to the voxel whose value
198/// in the @a index tree is @e n.
199/// @param source a tree with a scalar value type
200/// @param index a tree of the same configuration as @a source but with
201/// value type VIndex that maps voxels to elements of the output vector
202template<typename VectorValueType, typename SourceTreeType>
205 const SourceTreeType& source,
206 const typename SourceTreeType::template ValueConverter<VIndex>::Type& index);
207
208
209/// @brief Return a tree with the same active voxel topology as the @a index tree
210/// but whose voxel values are taken from the the given vector.
211/// @details The voxel whose value in the @a index tree is @e n gets assigned
212/// the <i>n</i>th element of the vector.
213/// @param index a tree with value type VIndex that maps voxels to elements of @a values
214/// @param values a vector of values with which to populate the active voxels of the output tree
215/// @param background the value for the inactive voxels of the output tree
216template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
217typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
220 const VIndexTreeType& index,
221 const TreeValueType& background);
222
223
224/// @brief Generate a sparse matrix of the index-space (&Delta;<i>x</i> = 1) %Laplacian operator
225/// using second-order finite differences.
226/// @details This construction assumes homogeneous Dirichlet boundary conditions
227/// (exterior grid points are zero).
228template<typename BoolTreeType>
231 const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
232 const BoolTreeType& interiorMask,
233 bool staggered = false);
234
235
236/// @brief Generate a sparse matrix of the index-space (&Delta;<i>x</i> = 1) %Laplacian operator
237/// with user-specified boundary conditions using second-order finite differences.
238/// @details Each thread gets its own copy of @a boundaryOp, which should be
239/// a functor of the form
240/// @code
241/// struct BoundaryOp {
242/// using ValueType = LaplacianMatrix::ValueType;
243/// void operator()(
244/// const Coord& ijk, // coordinates of a boundary voxel
245/// const Coord& ijkNeighbor, // coordinates of an exterior neighbor of ijk
246/// ValueType& source, // element of source vector corresponding to ijk
247/// ValueType& diagonal // element of Laplacian matrix corresponding to ijk
248/// ) const;
249/// };
250/// @endcode
251/// The functor is called for each of the exterior neighbors of each boundary voxel @ijk,
252/// and it must specify a boundary condition for @ijk by modifying one or both of two
253/// provided values: an entry in the given @a source vector corresponding to @ijk and
254/// the weighting coefficient for @ijk in the %Laplacian matrix.
255template<typename BoolTreeType, typename BoundaryOp>
258 const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
259 const BoolTreeType& interiorMask,
260 const BoundaryOp& boundaryOp,
262 bool staggered = false);
263
264
265/// @brief Dirichlet boundary condition functor
266/// @details This is useful in describing fluid/air interfaces, where the pressure
267/// of the air is assumed to be zero.
268template<typename ValueType>
270 inline void operator()(const Coord&, const Coord&, ValueType&, ValueType& diag) const {
271 // Exterior neighbors are empty, so decrement the weighting coefficient
272 // as for interior neighbors but leave the source vector unchanged.
273 diag -= 1;
274 }
275};
276
277//@}
278
279
280////////////////////////////////////////
281
282/// @cond OPENVDB_DOCS_INTERNAL
283
284namespace internal {
285
286/// @brief Functor for use with LeafManager::foreach() to populate an array
287/// with per-leaf active voxel counts
288template<typename LeafType>
289struct LeafCountOp
290{
291 VIndex* count;
292 LeafCountOp(VIndex* count_): count(count_) {}
293 void operator()(const LeafType& leaf, size_t leafIdx) const {
294 count[leafIdx] = static_cast<VIndex>(leaf.onVoxelCount());
295 }
296};
297
298
299/// @brief Functor for use with LeafManager::foreach() to populate
300/// active leaf voxels with sequential indices
301template<typename LeafType>
302struct LeafIndexOp
303{
304 const VIndex* count;
305 LeafIndexOp(const VIndex* count_): count(count_) {}
306 void operator()(LeafType& leaf, size_t leafIdx) const {
307 VIndex idx = (leafIdx == 0) ? 0 : count[leafIdx - 1];
308 for (typename LeafType::ValueOnIter it = leaf.beginValueOn(); it; ++it) {
309 it.setValue(idx++);
310 }
311 }
312};
313
314} // namespace internal
315
316/// @endcond
317
318template<typename VIndexTreeType>
319inline void
320populateIndexTree(VIndexTreeType& result)
321{
322 using LeafT = typename VIndexTreeType::LeafNodeType;
323 using LeafMgrT = typename tree::LeafManager<VIndexTreeType>;
324
325 // Linearize the tree.
326 LeafMgrT leafManager(result);
327 const size_t leafCount = leafManager.leafCount();
328
329 if (leafCount == 0) return;
330
331 // Count the number of active voxels in each leaf node.
332 std::unique_ptr<VIndex[]> perLeafCount(new VIndex[leafCount]);
333 VIndex* perLeafCountPtr = perLeafCount.get();
334 leafManager.foreach(internal::LeafCountOp<LeafT>(perLeafCountPtr));
335
336 // The starting index for each leaf node is the total number
337 // of active voxels in all preceding leaf nodes.
338 for (size_t i = 1; i < leafCount; ++i) {
339 perLeafCount[i] += perLeafCount[i - 1];
340 }
341
342 // The last accumulated value should be the total of all active voxels.
343 assert(Index64(perLeafCount[leafCount-1]) == result.activeVoxelCount());
344
345 // Parallelize over the leaf nodes of the tree, storing a unique index
346 // in each active voxel.
347 leafManager.foreach(internal::LeafIndexOp<LeafT>(perLeafCountPtr));
348}
349
350
351template<typename TreeType>
352inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
353createIndexTree(const TreeType& tree)
354{
355 using VIdxTreeT = typename TreeType::template ValueConverter<VIndex>::Type;
356
357 // Construct an output tree with the same active voxel topology as the input tree.
358 const VIndex invalidIdx = -1;
359 typename VIdxTreeT::Ptr result(
360 new VIdxTreeT(tree, /*background=*/invalidIdx, TopologyCopy()));
361
362 // All active voxels are degrees of freedom, including voxels contained in active tiles.
363 result->voxelizeActiveTiles();
364
365 populateIndexTree(*result);
366
367 return result;
368}
369
370
371////////////////////////////////////////
372
373/// @cond OPENVDB_DOCS_INTERNAL
374
375namespace internal {
376
377/// @brief Functor for use with LeafManager::foreach() to populate a vector
378/// with the values of a tree's active voxels
379template<typename VectorValueType, typename SourceTreeType>
380struct CopyToVecOp
381{
382 using VIdxTreeT = typename SourceTreeType::template ValueConverter<VIndex>::Type;
383 using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
384 using LeafT = typename SourceTreeType::LeafNodeType;
385 using TreeValueT = typename SourceTreeType::ValueType;
386 using VectorT = typename math::pcg::Vector<VectorValueType>;
387
388 const SourceTreeType* tree;
389 VectorT* vector;
390
391 CopyToVecOp(const SourceTreeType& t, VectorT& v): tree(&t), vector(&v) {}
392
393 void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
394 {
395 VectorT& vec = *vector;
396 if (const LeafT* leaf = tree->probeLeaf(idxLeaf.origin())) {
397 // If a corresponding leaf node exists in the source tree,
398 // copy voxel values from the source node to the output vector.
399 for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
400 vec[*it] = leaf->getValue(it.pos());
401 }
402 } else {
403 // If no corresponding leaf exists in the source tree,
404 // fill the vector with a uniform value.
405 const TreeValueT& value = tree->getValue(idxLeaf.origin());
406 for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
407 vec[*it] = value;
408 }
409 }
410 }
411};
412
413} // namespace internal
414
415/// @endcond
416
417template<typename VectorValueType, typename SourceTreeType>
418inline typename math::pcg::Vector<VectorValueType>::Ptr
419createVectorFromTree(const SourceTreeType& tree,
420 const typename SourceTreeType::template ValueConverter<VIndex>::Type& idxTree)
421{
422 using VIdxTreeT = typename SourceTreeType::template ValueConverter<VIndex>::Type;
423 using VIdxLeafMgrT = tree::LeafManager<const VIdxTreeT>;
424 using VectorT = typename math::pcg::Vector<VectorValueType>;
425
426 // Allocate the vector.
427 const size_t numVoxels = idxTree.activeVoxelCount();
428 typename VectorT::Ptr result(new VectorT(static_cast<math::pcg::SizeType>(numVoxels)));
429
430 // Parallelize over the leaf nodes of the index tree, filling the output vector
431 // with values from corresponding voxels of the source tree.
432 VIdxLeafMgrT leafManager(idxTree);
433 leafManager.foreach(internal::CopyToVecOp<VectorValueType, SourceTreeType>(tree, *result));
434
435 return result;
436}
437
438
439////////////////////////////////////////
440
441/// @cond OPENVDB_DOCS_INTERNAL
442
443namespace internal {
444
445/// @brief Functor for use with LeafManager::foreach() to populate a tree
446/// with values from a vector
447template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
448struct CopyFromVecOp
449{
450 using OutTreeT = typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
451 using OutLeafT = typename OutTreeT::LeafNodeType;
452 using VIdxLeafT = typename VIndexTreeType::LeafNodeType;
453 using VectorT = typename math::pcg::Vector<VectorValueType>;
454
455 const VectorT* vector;
456 OutTreeT* tree;
457
458 CopyFromVecOp(const VectorT& v, OutTreeT& t): vector(&v), tree(&t) {}
459
460 void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
461 {
462 const VectorT& vec = *vector;
463 OutLeafT* leaf = tree->probeLeaf(idxLeaf.origin());
464 assert(leaf != nullptr);
465 for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
466 leaf->setValueOnly(it.pos(), static_cast<TreeValueType>(vec[*it]));
467 }
468 }
469};
470
471} // namespace internal
472
473/// @endcond
474
475template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
476inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
479 const VIndexTreeType& idxTree,
480 const TreeValueType& background)
481{
482 using OutTreeT = typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
483 using VIdxLeafMgrT = typename tree::LeafManager<const VIndexTreeType>;
484
485 // Construct an output tree with the same active voxel topology as the index tree.
486 typename OutTreeT::Ptr result(new OutTreeT(idxTree, background, TopologyCopy()));
487 OutTreeT& tree = *result;
488
489 // Parallelize over the leaf nodes of the index tree, populating voxels
490 // of the output tree with values from the input vector.
491 VIdxLeafMgrT leafManager(idxTree);
492 leafManager.foreach(
493 internal::CopyFromVecOp<TreeValueType, VIndexTreeType, VectorValueType>(vector, tree));
494
495 return result;
496}
497
498
499////////////////////////////////////////
500
501/// @cond OPENVDB_DOCS_INTERNAL
502
503namespace internal {
504
505/// Functor for use with LeafManager::foreach() to populate a sparse %Laplacian matrix
506template<typename BoolTreeType, typename BoundaryOp>
507struct ISStaggeredLaplacianOp
508{
509 using VIdxTreeT = typename BoolTreeType::template ValueConverter<VIndex>::Type;
510 using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
511 using ValueT = LaplacianMatrix::ValueType;
512 using VectorT = typename math::pcg::Vector<ValueT>;
513
514 LaplacianMatrix* laplacian;
515 const VIdxTreeT* idxTree;
516 const BoolTreeType* interiorMask;
517 const BoundaryOp boundaryOp;
518 VectorT* source;
519
520 ISStaggeredLaplacianOp(LaplacianMatrix& m, const VIdxTreeT& idx,
521 const BoolTreeType& mask, const BoundaryOp& op, VectorT& src):
522 laplacian(&m), idxTree(&idx), interiorMask(&mask), boundaryOp(op), source(&src) {}
523
524 void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
525 {
526 // Local accessors
527 typename tree::ValueAccessor<const BoolTreeType> interior(*interiorMask);
528 typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
529
530 Coord ijk;
531 VIndex column;
532 const ValueT diagonal = -6.f, offDiagonal = 1.f;
533
534 // Loop over active voxels in this leaf.
535 for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
536 assert(it.getValue() > -1);
537 const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
538
539 LaplacianMatrix::RowEditor row = laplacian->getRowEditor(rowNum);
540
541 ijk = it.getCoord();
542 if (interior.isValueOn(ijk)) {
543 // The current voxel is an interior voxel.
544 // All of its neighbors are in the solution domain.
545
546 // -x direction
547 row.setValue(vectorIdx.getValue(ijk.offsetBy(-1, 0, 0)), offDiagonal);
548 // -y direction
549 row.setValue(vectorIdx.getValue(ijk.offsetBy(0, -1, 0)), offDiagonal);
550 // -z direction
551 row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, -1)), offDiagonal);
552 // diagonal
553 row.setValue(rowNum, diagonal);
554 // +z direction
555 row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, 1)), offDiagonal);
556 // +y direction
557 row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 1, 0)), offDiagonal);
558 // +x direction
559 row.setValue(vectorIdx.getValue(ijk.offsetBy(1, 0, 0)), offDiagonal);
560
561 } else {
562 // The current voxel is a boundary voxel.
563 // At least one of its neighbors is outside the solution domain.
564
565 ValueT modifiedDiagonal = 0.f;
566
567 // -x direction
568 if (vectorIdx.probeValue(ijk.offsetBy(-1, 0, 0), column)) {
569 row.setValue(column, offDiagonal);
570 modifiedDiagonal -= 1;
571 } else {
572 boundaryOp(ijk, ijk.offsetBy(-1, 0, 0), source->at(rowNum), modifiedDiagonal);
573 }
574 // -y direction
575 if (vectorIdx.probeValue(ijk.offsetBy(0, -1, 0), column)) {
576 row.setValue(column, offDiagonal);
577 modifiedDiagonal -= 1;
578 } else {
579 boundaryOp(ijk, ijk.offsetBy(0, -1, 0), source->at(rowNum), modifiedDiagonal);
580 }
581 // -z direction
582 if (vectorIdx.probeValue(ijk.offsetBy(0, 0, -1), column)) {
583 row.setValue(column, offDiagonal);
584 modifiedDiagonal -= 1;
585 } else {
586 boundaryOp(ijk, ijk.offsetBy(0, 0, -1), source->at(rowNum), modifiedDiagonal);
587 }
588 // +z direction
589 if (vectorIdx.probeValue(ijk.offsetBy(0, 0, 1), column)) {
590 row.setValue(column, offDiagonal);
591 modifiedDiagonal -= 1;
592 } else {
593 boundaryOp(ijk, ijk.offsetBy(0, 0, 1), source->at(rowNum), modifiedDiagonal);
594 }
595 // +y direction
596 if (vectorIdx.probeValue(ijk.offsetBy(0, 1, 0), column)) {
597 row.setValue(column, offDiagonal);
598 modifiedDiagonal -= 1;
599 } else {
600 boundaryOp(ijk, ijk.offsetBy(0, 1, 0), source->at(rowNum), modifiedDiagonal);
601 }
602 // +x direction
603 if (vectorIdx.probeValue(ijk.offsetBy(1, 0, 0), column)) {
604 row.setValue(column, offDiagonal);
605 modifiedDiagonal -= 1;
606 } else {
607 boundaryOp(ijk, ijk.offsetBy(1, 0, 0), source->at(rowNum), modifiedDiagonal);
608 }
609 // diagonal
610 row.setValue(rowNum, modifiedDiagonal);
611 }
612 } // end loop over voxels
613 }
614};
615
616
617// Stencil 1 is the correct stencil, but stencil 2 requires
618// half as many comparisons and produces smoother results at boundaries.
619//#define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 1
620#define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 2
621
622/// Functor for use with LeafManager::foreach() to populate a sparse %Laplacian matrix
623template<typename VIdxTreeT, typename BoundaryOp>
624struct ISLaplacianOp
625{
626 using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
627 using ValueT = LaplacianMatrix::ValueType;
628 using VectorT = typename math::pcg::Vector<ValueT>;
629
631 const VIdxTreeT* idxTree;
632 const BoundaryOp boundaryOp;
633 VectorT* source;
634
635 ISLaplacianOp(LaplacianMatrix& m, const VIdxTreeT& idx, const BoundaryOp& op, VectorT& src):
636 laplacian(&m), idxTree(&idx), boundaryOp(op), source(&src) {}
637
638 void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
639 {
640 typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
641
642 const int kNumOffsets = 6;
643 const Coord ijkOffset[kNumOffsets] = {
644#if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
645 Coord(-1,0,0), Coord(1,0,0), Coord(0,-1,0), Coord(0,1,0), Coord(0,0,-1), Coord(0,0,1)
646#else
647 Coord(-2,0,0), Coord(2,0,0), Coord(0,-2,0), Coord(0,2,0), Coord(0,0,-2), Coord(0,0,2)
648#endif
649 };
650
651 // For each active voxel in this leaf...
652 for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
653 assert(it.getValue() > -1);
654
655 const Coord ijk = it.getCoord();
656 const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
657
658 LaplacianMatrix::RowEditor row = laplacian->getRowEditor(rowNum);
659
660 ValueT modifiedDiagonal = 0.f;
661
662 // For each of the neighbors of the voxel at (i,j,k)...
663 for (int dir = 0; dir < kNumOffsets; ++dir) {
664 const Coord neighbor = ijk + ijkOffset[dir];
665 VIndex column;
666 // For collocated vector grids, the central differencing stencil requires
667 // access to neighbors at a distance of two voxels in each direction
668 // (-x, +x, -y, +y, -z, +z).
669#if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
670 const bool ijkIsInterior = (vectorIdx.probeValue(neighbor + ijkOffset[dir], column)
671 && vectorIdx.isValueOn(neighbor));
672#else
673 const bool ijkIsInterior = vectorIdx.probeValue(neighbor, column);
674#endif
675 if (ijkIsInterior) {
676 // If (i,j,k) is sufficiently far away from the exterior,
677 // set its weight to one and adjust the center weight accordingly.
678 row.setValue(column, 1.f);
679 modifiedDiagonal -= 1.f;
680 } else {
681 // If (i,j,k) is adjacent to or one voxel away from the exterior,
682 // invoke the boundary condition functor.
683 boundaryOp(ijk, neighbor, source->at(rowNum), modifiedDiagonal);
684 }
685 }
686 // Set the (possibly modified) weight for the voxel at (i,j,k).
687 row.setValue(rowNum, modifiedDiagonal);
688 }
689 }
690};
691
692} // namespace internal
693
694/// @endcond
695
696
697template<typename BoolTreeType>
698inline LaplacianMatrix::Ptr
699createISLaplacian(const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
700 const BoolTreeType& interiorMask, bool staggered)
701{
702 using ValueT = LaplacianMatrix::ValueType;
704 static_cast<math::pcg::SizeType>(idxTree.activeVoxelCount()));
706 return createISLaplacianWithBoundaryConditions(idxTree, interiorMask, op, unused, staggered);
707}
708
709
710template<typename BoolTreeType, typename BoundaryOp>
711inline LaplacianMatrix::Ptr
713 const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
714 const BoolTreeType& interiorMask,
715 const BoundaryOp& boundaryOp,
717 bool staggered)
718{
719 using VIdxTreeT = typename BoolTreeType::template ValueConverter<VIndex>::Type;
720 using VIdxLeafMgrT = typename tree::LeafManager<const VIdxTreeT>;
721
722 // The number of active voxels is the number of degrees of freedom.
723 const Index64 numDoF = idxTree.activeVoxelCount();
724
725 // Construct the matrix.
726 LaplacianMatrix::Ptr laplacianPtr(
727 new LaplacianMatrix(static_cast<math::pcg::SizeType>(numDoF)));
728 LaplacianMatrix& laplacian = *laplacianPtr;
729
730 // Populate the matrix using a second-order, 7-point CD stencil.
731 VIdxLeafMgrT idxLeafManager(idxTree);
732 if (staggered) {
733 idxLeafManager.foreach(internal::ISStaggeredLaplacianOp<BoolTreeType, BoundaryOp>(
734 laplacian, idxTree, interiorMask, boundaryOp, source));
735 } else {
736 idxLeafManager.foreach(internal::ISLaplacianOp<VIdxTreeT, BoundaryOp>(
737 laplacian, idxTree, boundaryOp, source));
738 }
739
740 return laplacianPtr;
741}
742
743
744////////////////////////////////////////
745
746
747template<typename TreeType>
748typename TreeType::Ptr
749solve(const TreeType& inTree, math::pcg::State& state, bool staggered)
750{
751 util::NullInterrupter interrupter;
752 return solve(inTree, state, interrupter, staggered);
753}
754
755
756template<typename TreeType, typename Interrupter>
757typename TreeType::Ptr
758solve(const TreeType& inTree, math::pcg::State& state, Interrupter& interrupter, bool staggered)
759{
761 return solveWithBoundaryConditions(inTree, boundaryOp, state, interrupter, staggered);
762}
763
764
765template<typename TreeType, typename BoundaryOp, typename Interrupter>
766typename TreeType::Ptr
767solveWithBoundaryConditions(const TreeType& inTree, const BoundaryOp& boundaryOp,
768 math::pcg::State& state, Interrupter& interrupter, bool staggered)
769{
771 return solveWithBoundaryConditionsAndPreconditioner<DefaultPrecondT>(
772 inTree, boundaryOp, state, interrupter, staggered);
773}
774
775
776template<
777 typename PreconditionerType,
778 typename TreeType,
779 typename BoundaryOp,
780 typename Interrupter>
781typename TreeType::Ptr
783 const TreeType& inTree,
784 const BoundaryOp& boundaryOp,
785 math::pcg::State& state,
786 Interrupter& interrupter,
787 bool staggered)
788{
789 return solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>(
790 /*source=*/inTree, /*domain mask=*/inTree, boundaryOp, state, interrupter, staggered);
791}
792
793template<
794 typename PreconditionerType,
795 typename TreeType,
796 typename DomainTreeType,
797 typename BoundaryOp,
798 typename Interrupter>
799typename TreeType::Ptr
801 const TreeType& inTree,
802 const DomainTreeType& domainMask,
803 const BoundaryOp& boundaryOp,
804 math::pcg::State& state,
805 Interrupter& interrupter,
806 bool staggered)
807{
808 using TreeValueT = typename TreeType::ValueType;
809 using VecValueT = LaplacianMatrix::ValueType;
810 using VectorT = typename math::pcg::Vector<VecValueT>;
811 using VIdxTreeT = typename TreeType::template ValueConverter<VIndex>::Type;
812 using MaskTreeT = typename TreeType::template ValueConverter<bool>::Type;
813
814 // 1. Create a mapping from active voxels of the input tree to elements of a vector.
815 typename VIdxTreeT::ConstPtr idxTree = createIndexTree(domainMask);
816
817 // 2. Populate a vector with values from the input tree.
818 typename VectorT::Ptr b = createVectorFromTree<VecValueT>(inTree, *idxTree);
819
820 // 3. Create a mask of the interior voxels of the input tree (from the densified index tree).
821 /// @todo Is this really needed?
822 typename MaskTreeT::Ptr interiorMask(
823 new MaskTreeT(*idxTree, /*background=*/false, TopologyCopy()));
824 tools::erodeActiveValues(*interiorMask, /*iterations=*/1, tools::NN_FACE, tools::IGNORE_TILES);
825
826 // 4. Create the Laplacian matrix.
828 *idxTree, *interiorMask, boundaryOp, *b, staggered);
829
830 // 5. Solve the Poisson equation.
831 laplacian->scale(-1.0); // matrix is negative-definite; solve -M x = -b
832 b->scale(-1.0);
833 typename VectorT::Ptr x(new VectorT(b->size(), zeroVal<VecValueT>()));
835 new PreconditionerType(*laplacian));
836 if (!precond->isValid()) {
838 }
839
840 state = math::pcg::solve(*laplacian, *b, *x, *precond, interrupter, state);
841
842 // 6. Populate the output tree with values from the solution vector.
843 /// @todo if (state.success) ... ?
844 return createTreeFromVector<TreeValueT>(*x, *idxTree, /*background=*/zeroVal<TreeValueT>());
845}
846
847
848////////////////////////////////////////
849
850
851// Explicit Template Instantiation
852
853#ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
854
855#ifdef OPENVDB_INSTANTIATE_POISSONSOLVER
857#endif
858
859#define _FUNCTION(TreeT) \
860 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
861 math::pcg::IncompleteCholeskyPreconditioner<LaplacianMatrix>>( \
862 const TreeT&, const TreeT&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
863 math::pcg::State&, util::NullInterrupter&, bool)
865#undef _FUNCTION
866
867#define _FUNCTION(TreeT) \
868 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
869 math::pcg::IncompleteCholeskyPreconditioner<LaplacianMatrix>>( \
870 const TreeT&, const BoolTree&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
871 math::pcg::State&, util::NullInterrupter&, bool)
873#undef _FUNCTION
874
875#define _FUNCTION(TreeT) \
876 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
877 math::pcg::JacobiPreconditioner<LaplacianMatrix>>( \
878 const TreeT&, const TreeT&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
879 math::pcg::State&, util::NullInterrupter&, bool)
881#undef _FUNCTION
882
883#define _FUNCTION(TreeT) \
884 TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
885 math::pcg::JacobiPreconditioner<LaplacianMatrix>>( \
886 const TreeT&, const BoolTree&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
887 math::pcg::State&, util::NullInterrupter&, bool)
889#undef _FUNCTION
890
891#endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
892
893
894} // namespace poisson
895} // namespace tools
896} // namespace OPENVDB_VERSION_NAME
897} // namespace openvdb
898
899#endif // OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
Preconditioned conjugate gradient solver (solves Ax = b using the conjugate gradient method with one ...
ValueT value
Definition GridBuilder.h:1290
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
Implementation of morphological dilation and erosion.
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition Types.h:644
Signed (x, y, z) 32-bit integer coordinates.
Definition Coord.h:25
Coord offsetBy(Int32 dx, Int32 dy, Int32 dz) const
Definition Coord.h:91
Preconditioner using incomplete Cholesky factorization.
Definition ConjGradient.h:1347
Diagonal preconditioner.
Definition ConjGradient.h:1281
virtual bool isValid() const
Definition ConjGradient.h:473
SharedPtr< Preconditioner > Ptr
Definition ConjGradient.h:468
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition ConjGradient.h:238
ValueType_ ValueType
Definition ConjGradient.h:240
SharedPtr< SparseStencilMatrix > Ptr
Definition ConjGradient.h:242
Lightweight, variable-length vector.
Definition ConjGradient.h:139
SharedPtr< Vector > Ptr
Definition ConjGradient.h:142
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition LeafManager.h:85
Definition ValueAccessor.h:191
Index32 SizeType
Definition ConjGradient.h:33
LaplacianMatrix::Ptr createISLaplacianWithBoundaryConditions(const typename BoolTreeType::template ValueConverter< VIndex >::Type &vectorIndexTree, const BoolTreeType &interiorMask, const BoundaryOp &boundaryOp, typename math::pcg::Vector< LaplacianMatrix::ValueType > &source, bool staggered=false)
Generate a sparse matrix of the index-space (Δx = 1) Laplacian operator with user-specified boundary ...
Definition PoissonSolver.h:712
TreeType::Ptr solveWithBoundaryConditions(const TreeType &, const BoundaryOp &, math::pcg::State &, Interrupter &, bool staggered=false)
Solve ∇2x = b for x with user-specified boundary conditions, where b is a vector comprising the value...
Definition PoissonSolver.h:767
void populateIndexTree(VIndexTreeType &)
Overwrite each active voxel in the given scalar tree with a sequential index, starting from zero.
Definition PoissonSolver.h:320
LaplacianMatrix::Ptr createISLaplacian(const typename BoolTreeType::template ValueConverter< VIndex >::Type &vectorIndexTree, const BoolTreeType &interiorMask, bool staggered=false)
Generate a sparse matrix of the index-space (Δx = 1) Laplacian operator using second-order finite dif...
Definition PoissonSolver.h:699
Int32 VIndex
Definition PoissonSolver.h:80
TreeType::template ValueConverter< VIndex >::Type::Ptr createIndexTree(const TreeType &)
Iterate over the active voxels of the input tree and for each one assign its index in the iteration s...
Definition PoissonSolver.h:353
math::pcg::SparseStencilMatrix< double, 7 > LaplacianMatrix
The type of a matrix used to represent a three-dimensional Laplacian operator.
Definition PoissonSolver.h:83
math::pcg::Vector< VectorValueType >::Ptr createVectorFromTree(const SourceTreeType &source, const typename SourceTreeType::template ValueConverter< VIndex >::Type &index)
Return a vector of the active voxel values of the scalar-valued source tree.
Definition PoissonSolver.h:419
VIndexTreeType::template ValueConverter< TreeValueType >::Type::Ptr createTreeFromVector(const math::pcg::Vector< VectorValueType > &values, const VIndexTreeType &index, const TreeValueType &background)
Return a tree with the same active voxel topology as the index tree but whose voxel values are taken ...
Definition PoissonSolver.h:477
TreeType::Ptr solveWithBoundaryConditionsAndPreconditioner(const TreeType &, const BoundaryOp &, math::pcg::State &, Interrupter &, bool staggered=false)
Solve ∇2x = b for x, where b is a vector comprising the values of all of the active voxels in the inp...
Definition PoissonSolver.h:782
GridType::template ValueConverter< bool >::Type::Ptr interiorMask(const GridType &grid, const double isovalue=0.0)
Given an input grid of any type, return a new, boolean grid whose active voxel topology matches the i...
Definition Mask.h:105
GridType::Ptr laplacian(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the Laplacian of the given scalar grid.
Definition GridOperators.h:1018
int32_t Int32
Definition Types.h:56
uint64_t Index64
Definition Types.h:53
Definition Exceptions.h:13
Information about the state of a conjugate gradient solution.
Definition ConjGradient.h:46
Dirichlet boundary condition functor.
Definition PoissonSolver.h:269
void operator()(const Coord &, const Coord &, ValueType &, ValueType &diag) const
Definition PoissonSolver.h:270
Base class for interrupters.
Definition NullInterrupter.h:26
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition version.h.in:157