Loading...
Searching...
No Matches
PoissonSolver.h
Go to the documentation of this file.
// Copyright Contributors to the OpenVDB Project
// SPDX-License-Identifier: MPL-2.0
//
/// @file PoissonSolver.h
///
/// @authors D.J. Hill, Peter Cucka
///
/// @brief Solve Poisson's equation &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i>
/// for <i>x</i>, where @e b is a vector comprising the values of all of the active voxels
/// in a grid.
///
/// @par Example:
/// Solve for the pressure in a cubic tank of liquid, assuming uniform boundary conditions:
/// @code
/// FloatTree source(/*background=*/0.0f);
/// // Activate voxels to indicate that they contain liquid.
/// source.fill(CoordBBox(Coord(0, -10, 0), Coord(10, 0, 10)), /*value=*/0.0f);
///
/// math::pcg::State state = math::pcg::terminationDefaults<float>();
/// FloatTree::Ptr solution = tools::poisson::solve(source, state);
/// @endcode
///
/// @par Example:
/// Solve for the pressure, <i>P</i>, in a cubic tank of liquid that is open at the top.
/// Boundary conditions are <i>P</i>&nbsp;=&nbsp;0 at the top,
/// &part;<i>P</i>/&part;<i>y</i>&nbsp;=&nbsp;&minus;1 at the bottom
/// and &part;<i>P</i>/&part;<i>x</i>&nbsp;=&nbsp;0 at the sides:
/// <pre>
/// P = 0
/// +--------+ (N,0,N)
/// /| /|
/// (0,0,0) +--------+ |
/// | | | | dP/dx = 0
/// dP/dx = 0 | +------|-+
/// |/ |/
/// (0,-N,0) +--------+ (N,-N,0)
/// dP/dy = -1
/// </pre>
/// @code
/// const int N = 10;
/// DoubleTree source(/*background=*/0.0);
/// // Activate voxels to indicate that they contain liquid.
/// source.fill(CoordBBox(Coord(0, -N, 0), Coord(N, 0, N)), /*value=*/0.0);
///
/// auto boundary = [](const openvdb::Coord& ijk, const openvdb::Coord& neighbor,
/// double& source, double& diagonal)
/// {
/// if (neighbor.x() == ijk.x() && neighbor.z() == ijk.z()) {
/// if (neighbor.y() < ijk.y()) source -= 1.0;
/// else diagonal -= 1.0;
/// }
/// };
///
/// math::pcg::State state = math::pcg::terminationDefaults<double>();
/// util::NullInterrupter interrupter;
///
/// DoubleTree::Ptr solution = tools::poisson::solveWithBoundaryConditions(
/// source, boundary, state, interrupter);
/// @endcode
#ifndef OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
#define OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
#include <openvdb/Types.h>
#include <openvdb/math/ConjGradient.h>
#include <openvdb/tree/LeafManager.h>
#include <openvdb/tree/Tree.h>
#include <openvdb/util/NullInterrupter.h>
#include "Morphology.h" // for erodeActiveValues
#include <openvdb/openvdb.h>
namespace openvdb {
OPENVDB_USE_VERSION_NAMESPACE
namespace OPENVDB_VERSION_NAME {
namespace tools {
namespace poisson {
// This type should be at least as wide as math::pcg::SizeType.
using VIndex = Int32;
/// The type of a matrix used to represent a three-dimensional %Laplacian operator
using LaplacianMatrix = math::pcg::SparseStencilMatrix<double, 7>;
//@{
/// @brief Solve &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i> for <i>x</i>,
/// where @e b is a vector comprising the values of all of the active voxels
/// in the input tree.
/// @return a new tree, with the same active voxel topology as the input tree,
/// whose voxel values are the elements of the solution vector <i>x</i>.
/// @details On input, the State object should specify convergence criteria
/// (minimum error and maximum number of iterations); on output, it gives
/// the actual termination conditions.
/// @details The solution is computed using the conjugate gradient method
/// with (where possible) incomplete Cholesky preconditioning, falling back
/// to Jacobi preconditioning.
/// @sa solveWithBoundaryConditions
template<typename TreeType>
typename TreeType::Ptr
solve(const TreeType&, math::pcg::State&, bool staggered = false);
template<typename TreeType, typename Interrupter>
typename TreeType::Ptr
solve(const TreeType&, math::pcg::State&, Interrupter&, bool staggered = false);
//@}
//@{
/// @brief Solve &nabla;<sup><small>2</small></sup><i>x</i> = <i>b</i> for <i>x</i>
/// with user-specified boundary conditions, where @e b is a vector comprising
/// the values of all of the active voxels in the input tree or domain mask if provided
/// @return a new tree, with the same active voxel topology as the input tree,
/// whose voxel values are the elements of the solution vector <i>x</i>.
/// @details On input, the State object should specify convergence criteria
/// (minimum error and maximum number of iterations); on output, it gives
/// the actual termination conditions.
/// @details The solution is computed using the conjugate gradient method with
/// the specified type of preconditioner (default: incomplete Cholesky),
/// falling back to Jacobi preconditioning if necessary.
/// @details Each thread gets its own copy of the BoundaryOp, which should be
/// a functor of the form
/// @code
/// struct BoundaryOp {
/// using ValueType = LaplacianMatrix::ValueType;
/// void operator()(
/// const Coord& ijk, // coordinates of a boundary voxel
/// const Coord& ijkNeighbor, // coordinates of an exterior neighbor of ijk
/// ValueType& source, // element of b corresponding to ijk
/// ValueType& diagonal // element of Laplacian matrix corresponding to ijk
/// ) const;
/// };
/// @endcode
/// The functor is called for each of the exterior neighbors of each boundary voxel @ijk,
/// and it must specify a boundary condition for @ijk by modifying one or both of two
/// provided values: the entry in the source vector @e b corresponding to @ijk and
/// the weighting coefficient for @ijk in the Laplacian operator matrix.
///
/// @sa solve
template<typename TreeType, typename BoundaryOp, typename Interrupter>
typename TreeType::Ptr
solveWithBoundaryConditions(
const TreeType&,
const BoundaryOp&,
math::pcg::State&,
Interrupter&,
bool staggered = false);
template<
typename PreconditionerType,
typename TreeType,
typename BoundaryOp,
typename Interrupter>
typename TreeType::Ptr
solveWithBoundaryConditionsAndPreconditioner(
const TreeType&,
const BoundaryOp&,
math::pcg::State&,
Interrupter&,
bool staggered = false);
template<
typename PreconditionerType,
typename TreeType,
typename DomainTreeType,
typename BoundaryOp,
typename Interrupter>
typename TreeType::Ptr
solveWithBoundaryConditionsAndPreconditioner(
const TreeType&,
const DomainTreeType&,
const BoundaryOp&,
math::pcg::State&,
Interrupter&,
bool staggered = false);
//@}
/// @name Low-level functions
//@{
// The following are low-level routines that can be used to assemble custom solvers.
/// @brief Overwrite each active voxel in the given scalar tree
/// with a sequential index, starting from zero.
template<typename VIndexTreeType>
void populateIndexTree(VIndexTreeType&);
/// @brief Iterate over the active voxels of the input tree and for each one
/// assign its index in the iteration sequence to the corresponding voxel
/// of an integer-valued output tree.
template<typename TreeType>
typename TreeType::template ValueConverter<VIndex>::Type::Ptr
createIndexTree(const TreeType&);
/// @brief Return a vector of the active voxel values of the scalar-valued @a source tree.
/// @details The <i>n</i>th element of the vector corresponds to the voxel whose value
/// in the @a index tree is @e n.
/// @param source a tree with a scalar value type
/// @param index a tree of the same configuration as @a source but with
/// value type VIndex that maps voxels to elements of the output vector
template<typename VectorValueType, typename SourceTreeType>
typename math::pcg::Vector<VectorValueType>::Ptr
createVectorFromTree(
const SourceTreeType& source,
const typename SourceTreeType::template ValueConverter<VIndex>::Type& index);
/// @brief Return a tree with the same active voxel topology as the @a index tree
/// but whose voxel values are taken from the the given vector.
/// @details The voxel whose value in the @a index tree is @e n gets assigned
/// the <i>n</i>th element of the vector.
/// @param index a tree with value type VIndex that maps voxels to elements of @a values
/// @param values a vector of values with which to populate the active voxels of the output tree
/// @param background the value for the inactive voxels of the output tree
template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
createTreeFromVector(
const math::pcg::Vector<VectorValueType>& values,
const VIndexTreeType& index,
const TreeValueType& background);
/// @brief Generate a sparse matrix of the index-space (&Delta;<i>x</i> = 1) %Laplacian operator
/// using second-order finite differences.
/// @details This construction assumes homogeneous Dirichlet boundary conditions
/// (exterior grid points are zero).
template<typename BoolTreeType>
LaplacianMatrix::Ptr
createISLaplacian(
const typename BoolTreeType::template ValueConverter<VIndex>::Type& vectorIndexTree,
const BoolTreeType& interiorMask,
bool staggered = false);
/// @brief Generate a sparse matrix of the index-space (&Delta;<i>x</i> = 1) %Laplacian operator
/// with user-specified boundary conditions using second-order finite differences.
/// @details Each thread gets its own copy of @a boundaryOp, which should be
/// a functor of the form
/// @code
/// struct BoundaryOp {
/// using ValueType = LaplacianMatrix::ValueType;
/// void operator()(
/// const Coord& ijk, // coordinates of a boundary voxel
/// const Coord& ijkNeighbor, // coordinates of an exterior neighbor of ijk
/// ValueType& source, // element of source vector corresponding to ijk
/// ValueType& diagonal // element of Laplacian matrix corresponding to ijk
/// ) const;
/// };
/// @endcode
/// The functor is called for each of the exterior neighbors of each boundary voxel @ijk,
/// and it must specify a boundary condition for @ijk by modifying one or both of two
/// provided values: an entry in the given @a source vector corresponding to @ijk and
/// the weighting coefficient for @ijk in the %Laplacian matrix.
template<typename BoolTreeType, typename BoundaryOp>
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);
/// @brief Dirichlet boundary condition functor
/// @details This is useful in describing fluid/air interfaces, where the pressure
/// of the air is assumed to be zero.
template<typename ValueType>
struct DirichletBoundaryOp {
inline void operator()(const Coord&, const Coord&, ValueType&, ValueType& diag) const {
// Exterior neighbors are empty, so decrement the weighting coefficient
// as for interior neighbors but leave the source vector unchanged.
diag -= 1;
}
};
//@}
////////////////////////////////////////
/// @cond OPENVDB_DOCS_INTERNAL
namespace internal {
/// @brief Functor for use with LeafManager::foreach() to populate an array
/// with per-leaf active voxel counts
template<typename LeafType>
struct LeafCountOp
{
VIndex* count;
LeafCountOp(VIndex* count_): count(count_) {}
void operator()(const LeafType& leaf, size_t leafIdx) const {
count[leafIdx] = static_cast<VIndex>(leaf.onVoxelCount());
}
};
/// @brief Functor for use with LeafManager::foreach() to populate
/// active leaf voxels with sequential indices
template<typename LeafType>
struct LeafIndexOp
{
const VIndex* count;
LeafIndexOp(const VIndex* count_): count(count_) {}
void operator()(LeafType& leaf, size_t leafIdx) const {
VIndex idx = (leafIdx == 0) ? 0 : count[leafIdx - 1];
for (typename LeafType::ValueOnIter it = leaf.beginValueOn(); it; ++it) {
it.setValue(idx++);
}
}
};
} // namespace internal
/// @endcond
template<typename VIndexTreeType>
inline void
populateIndexTree(VIndexTreeType& result)
{
using LeafT = typename VIndexTreeType::LeafNodeType;
using LeafMgrT = typename tree::LeafManager<VIndexTreeType>;
// Linearize the tree.
LeafMgrT leafManager(result);
const size_t leafCount = leafManager.leafCount();
if (leafCount == 0) return;
// Count the number of active voxels in each leaf node.
std::unique_ptr<VIndex[]> perLeafCount(new VIndex[leafCount]);
VIndex* perLeafCountPtr = perLeafCount.get();
leafManager.foreach(internal::LeafCountOp<LeafT>(perLeafCountPtr));
// The starting index for each leaf node is the total number
// of active voxels in all preceding leaf nodes.
for (size_t i = 1; i < leafCount; ++i) {
perLeafCount[i] += perLeafCount[i - 1];
}
// The last accumulated value should be the total of all active voxels.
assert(Index64(perLeafCount[leafCount-1]) == result.activeVoxelCount());
// Parallelize over the leaf nodes of the tree, storing a unique index
// in each active voxel.
leafManager.foreach(internal::LeafIndexOp<LeafT>(perLeafCountPtr));
}
template<typename TreeType>
inline typename TreeType::template ValueConverter<VIndex>::Type::Ptr
createIndexTree(const TreeType& tree)
{
using VIdxTreeT = typename TreeType::template ValueConverter<VIndex>::Type;
// Construct an output tree with the same active voxel topology as the input tree.
const VIndex invalidIdx = -1;
typename VIdxTreeT::Ptr result(
new VIdxTreeT(tree, /*background=*/invalidIdx, TopologyCopy()));
// All active voxels are degrees of freedom, including voxels contained in active tiles.
result->voxelizeActiveTiles();
populateIndexTree(*result);
return result;
}
////////////////////////////////////////
/// @cond OPENVDB_DOCS_INTERNAL
namespace internal {
/// @brief Functor for use with LeafManager::foreach() to populate a vector
/// with the values of a tree's active voxels
template<typename VectorValueType, typename SourceTreeType>
struct CopyToVecOp
{
using VIdxTreeT = typename SourceTreeType::template ValueConverter<VIndex>::Type;
using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
using LeafT = typename SourceTreeType::LeafNodeType;
using TreeValueT = typename SourceTreeType::ValueType;
using VectorT = typename math::pcg::Vector<VectorValueType>;
const SourceTreeType* tree;
VectorT* vector;
CopyToVecOp(const SourceTreeType& t, VectorT& v): tree(&t), vector(&v) {}
void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
{
VectorT& vec = *vector;
if (const LeafT* leaf = tree->probeLeaf(idxLeaf.origin())) {
// If a corresponding leaf node exists in the source tree,
// copy voxel values from the source node to the output vector.
for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
vec[*it] = leaf->getValue(it.pos());
}
} else {
// If no corresponding leaf exists in the source tree,
// fill the vector with a uniform value.
const TreeValueT& value = tree->getValue(idxLeaf.origin());
for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
vec[*it] = value;
}
}
}
};
} // namespace internal
/// @endcond
template<typename VectorValueType, typename SourceTreeType>
inline typename math::pcg::Vector<VectorValueType>::Ptr
createVectorFromTree(const SourceTreeType& tree,
const typename SourceTreeType::template ValueConverter<VIndex>::Type& idxTree)
{
using VIdxTreeT = typename SourceTreeType::template ValueConverter<VIndex>::Type;
using VIdxLeafMgrT = tree::LeafManager<const VIdxTreeT>;
using VectorT = typename math::pcg::Vector<VectorValueType>;
// Allocate the vector.
const size_t numVoxels = idxTree.activeVoxelCount();
typename VectorT::Ptr result(new VectorT(static_cast<math::pcg::SizeType>(numVoxels)));
// Parallelize over the leaf nodes of the index tree, filling the output vector
// with values from corresponding voxels of the source tree.
VIdxLeafMgrT leafManager(idxTree);
leafManager.foreach(internal::CopyToVecOp<VectorValueType, SourceTreeType>(tree, *result));
return result;
}
////////////////////////////////////////
/// @cond OPENVDB_DOCS_INTERNAL
namespace internal {
/// @brief Functor for use with LeafManager::foreach() to populate a tree
/// with values from a vector
template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
struct CopyFromVecOp
{
using OutTreeT = typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
using OutLeafT = typename OutTreeT::LeafNodeType;
using VIdxLeafT = typename VIndexTreeType::LeafNodeType;
using VectorT = typename math::pcg::Vector<VectorValueType>;
const VectorT* vector;
OutTreeT* tree;
CopyFromVecOp(const VectorT& v, OutTreeT& t): vector(&v), tree(&t) {}
void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
{
const VectorT& vec = *vector;
OutLeafT* leaf = tree->probeLeaf(idxLeaf.origin());
assert(leaf != nullptr);
for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
leaf->setValueOnly(it.pos(), static_cast<TreeValueType>(vec[*it]));
}
}
};
} // namespace internal
/// @endcond
template<typename TreeValueType, typename VIndexTreeType, typename VectorValueType>
inline typename VIndexTreeType::template ValueConverter<TreeValueType>::Type::Ptr
createTreeFromVector(
const math::pcg::Vector<VectorValueType>& vector,
const VIndexTreeType& idxTree,
const TreeValueType& background)
{
using OutTreeT = typename VIndexTreeType::template ValueConverter<TreeValueType>::Type;
using VIdxLeafMgrT = typename tree::LeafManager<const VIndexTreeType>;
// Construct an output tree with the same active voxel topology as the index tree.
typename OutTreeT::Ptr result(new OutTreeT(idxTree, background, TopologyCopy()));
OutTreeT& tree = *result;
// Parallelize over the leaf nodes of the index tree, populating voxels
// of the output tree with values from the input vector.
VIdxLeafMgrT leafManager(idxTree);
leafManager.foreach(
internal::CopyFromVecOp<TreeValueType, VIndexTreeType, VectorValueType>(vector, tree));
return result;
}
////////////////////////////////////////
/// @cond OPENVDB_DOCS_INTERNAL
namespace internal {
/// Functor for use with LeafManager::foreach() to populate a sparse %Laplacian matrix
template<typename BoolTreeType, typename BoundaryOp>
struct ISStaggeredLaplacianOp
{
using VIdxTreeT = typename BoolTreeType::template ValueConverter<VIndex>::Type;
using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
using ValueT = LaplacianMatrix::ValueType;
using VectorT = typename math::pcg::Vector<ValueT>;
LaplacianMatrix* laplacian;
const VIdxTreeT* idxTree;
const BoolTreeType* interiorMask;
const BoundaryOp boundaryOp;
VectorT* source;
ISStaggeredLaplacianOp(LaplacianMatrix& m, const VIdxTreeT& idx,
const BoolTreeType& mask, const BoundaryOp& op, VectorT& src):
laplacian(&m), idxTree(&idx), interiorMask(&mask), boundaryOp(op), source(&src) {}
void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
{
// Local accessors
typename tree::ValueAccessor<const BoolTreeType> interior(*interiorMask);
typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
Coord ijk;
VIndex column;
const ValueT diagonal = -6.f, offDiagonal = 1.f;
// Loop over active voxels in this leaf.
for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
assert(it.getValue() > -1);
const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
LaplacianMatrix::RowEditor row = laplacian->getRowEditor(rowNum);
ijk = it.getCoord();
if (interior.isValueOn(ijk)) {
// The current voxel is an interior voxel.
// All of its neighbors are in the solution domain.
// -x direction
row.setValue(vectorIdx.getValue(ijk.offsetBy(-1, 0, 0)), offDiagonal);
// -y direction
row.setValue(vectorIdx.getValue(ijk.offsetBy(0, -1, 0)), offDiagonal);
// -z direction
row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, -1)), offDiagonal);
// diagonal
row.setValue(rowNum, diagonal);
// +z direction
row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 0, 1)), offDiagonal);
// +y direction
row.setValue(vectorIdx.getValue(ijk.offsetBy(0, 1, 0)), offDiagonal);
// +x direction
row.setValue(vectorIdx.getValue(ijk.offsetBy(1, 0, 0)), offDiagonal);
} else {
// The current voxel is a boundary voxel.
// At least one of its neighbors is outside the solution domain.
ValueT modifiedDiagonal = 0.f;
// -x direction
if (vectorIdx.probeValue(ijk.offsetBy(-1, 0, 0), column)) {
row.setValue(column, offDiagonal);
modifiedDiagonal -= 1;
} else {
boundaryOp(ijk, ijk.offsetBy(-1, 0, 0), source->at(rowNum), modifiedDiagonal);
}
// -y direction
if (vectorIdx.probeValue(ijk.offsetBy(0, -1, 0), column)) {
row.setValue(column, offDiagonal);
modifiedDiagonal -= 1;
} else {
boundaryOp(ijk, ijk.offsetBy(0, -1, 0), source->at(rowNum), modifiedDiagonal);
}
// -z direction
if (vectorIdx.probeValue(ijk.offsetBy(0, 0, -1), column)) {
row.setValue(column, offDiagonal);
modifiedDiagonal -= 1;
} else {
boundaryOp(ijk, ijk.offsetBy(0, 0, -1), source->at(rowNum), modifiedDiagonal);
}
// +z direction
if (vectorIdx.probeValue(ijk.offsetBy(0, 0, 1), column)) {
row.setValue(column, offDiagonal);
modifiedDiagonal -= 1;
} else {
boundaryOp(ijk, ijk.offsetBy(0, 0, 1), source->at(rowNum), modifiedDiagonal);
}
// +y direction
if (vectorIdx.probeValue(ijk.offsetBy(0, 1, 0), column)) {
row.setValue(column, offDiagonal);
modifiedDiagonal -= 1;
} else {
boundaryOp(ijk, ijk.offsetBy(0, 1, 0), source->at(rowNum), modifiedDiagonal);
}
// +x direction
if (vectorIdx.probeValue(ijk.offsetBy(1, 0, 0), column)) {
row.setValue(column, offDiagonal);
modifiedDiagonal -= 1;
} else {
boundaryOp(ijk, ijk.offsetBy(1, 0, 0), source->at(rowNum), modifiedDiagonal);
}
// diagonal
row.setValue(rowNum, modifiedDiagonal);
}
} // end loop over voxels
}
};
// Stencil 1 is the correct stencil, but stencil 2 requires
// half as many comparisons and produces smoother results at boundaries.
//#define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 1
#define OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL 2
/// Functor for use with LeafManager::foreach() to populate a sparse %Laplacian matrix
template<typename VIdxTreeT, typename BoundaryOp>
struct ISLaplacianOp
{
using VIdxLeafT = typename VIdxTreeT::LeafNodeType;
using ValueT = LaplacianMatrix::ValueType;
using VectorT = typename math::pcg::Vector<ValueT>;
LaplacianMatrix* laplacian;
const VIdxTreeT* idxTree;
const BoundaryOp boundaryOp;
VectorT* source;
ISLaplacianOp(LaplacianMatrix& m, const VIdxTreeT& idx, const BoundaryOp& op, VectorT& src):
laplacian(&m), idxTree(&idx), boundaryOp(op), source(&src) {}
void operator()(const VIdxLeafT& idxLeaf, size_t /*leafIdx*/) const
{
typename tree::ValueAccessor<const VIdxTreeT> vectorIdx(*idxTree);
const int kNumOffsets = 6;
const Coord ijkOffset[kNumOffsets] = {
#if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
Coord(-1,0,0), Coord(1,0,0), Coord(0,-1,0), Coord(0,1,0), Coord(0,0,-1), Coord(0,0,1)
#else
Coord(-2,0,0), Coord(2,0,0), Coord(0,-2,0), Coord(0,2,0), Coord(0,0,-2), Coord(0,0,2)
#endif
};
// For each active voxel in this leaf...
for (typename VIdxLeafT::ValueOnCIter it = idxLeaf.cbeginValueOn(); it; ++it) {
assert(it.getValue() > -1);
const Coord ijk = it.getCoord();
const math::pcg::SizeType rowNum = static_cast<math::pcg::SizeType>(it.getValue());
LaplacianMatrix::RowEditor row = laplacian->getRowEditor(rowNum);
ValueT modifiedDiagonal = 0.f;
// For each of the neighbors of the voxel at (i,j,k)...
for (int dir = 0; dir < kNumOffsets; ++dir) {
const Coord neighbor = ijk + ijkOffset[dir];
VIndex column;
// For collocated vector grids, the central differencing stencil requires
// access to neighbors at a distance of two voxels in each direction
// (-x, +x, -y, +y, -z, +z).
#if OPENVDB_TOOLS_POISSON_LAPLACIAN_STENCIL == 1
const bool ijkIsInterior = (vectorIdx.probeValue(neighbor + ijkOffset[dir], column)
&& vectorIdx.isValueOn(neighbor));
#else
const bool ijkIsInterior = vectorIdx.probeValue(neighbor, column);
#endif
if (ijkIsInterior) {
// If (i,j,k) is sufficiently far away from the exterior,
// set its weight to one and adjust the center weight accordingly.
row.setValue(column, 1.f);
modifiedDiagonal -= 1.f;
} else {
// If (i,j,k) is adjacent to or one voxel away from the exterior,
// invoke the boundary condition functor.
boundaryOp(ijk, neighbor, source->at(rowNum), modifiedDiagonal);
}
}
// Set the (possibly modified) weight for the voxel at (i,j,k).
row.setValue(rowNum, modifiedDiagonal);
}
}
};
} // namespace internal
/// @endcond
template<typename BoolTreeType>
inline LaplacianMatrix::Ptr
createISLaplacian(const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
const BoolTreeType& interiorMask, bool staggered)
{
using ValueT = LaplacianMatrix::ValueType;
math::pcg::Vector<ValueT> unused(
static_cast<math::pcg::SizeType>(idxTree.activeVoxelCount()));
DirichletBoundaryOp<ValueT> op;
return createISLaplacianWithBoundaryConditions(idxTree, interiorMask, op, unused, staggered);
}
template<typename BoolTreeType, typename BoundaryOp>
inline LaplacianMatrix::Ptr
createISLaplacianWithBoundaryConditions(
const typename BoolTreeType::template ValueConverter<VIndex>::Type& idxTree,
const BoolTreeType& interiorMask,
const BoundaryOp& boundaryOp,
typename math::pcg::Vector<LaplacianMatrix::ValueType>& source,
bool staggered)
{
using VIdxTreeT = typename BoolTreeType::template ValueConverter<VIndex>::Type;
using VIdxLeafMgrT = typename tree::LeafManager<const VIdxTreeT>;
// The number of active voxels is the number of degrees of freedom.
const Index64 numDoF = idxTree.activeVoxelCount();
// Construct the matrix.
LaplacianMatrix::Ptr laplacianPtr(
new LaplacianMatrix(static_cast<math::pcg::SizeType>(numDoF)));
LaplacianMatrix& laplacian = *laplacianPtr;
// Populate the matrix using a second-order, 7-point CD stencil.
VIdxLeafMgrT idxLeafManager(idxTree);
if (staggered) {
idxLeafManager.foreach(internal::ISStaggeredLaplacianOp<BoolTreeType, BoundaryOp>(
laplacian, idxTree, interiorMask, boundaryOp, source));
} else {
idxLeafManager.foreach(internal::ISLaplacianOp<VIdxTreeT, BoundaryOp>(
laplacian, idxTree, boundaryOp, source));
}
return laplacianPtr;
}
////////////////////////////////////////
template<typename TreeType>
typename TreeType::Ptr
solve(const TreeType& inTree, math::pcg::State& state, bool staggered)
{
util::NullInterrupter interrupter;
return solve(inTree, state, interrupter, staggered);
}
template<typename TreeType, typename Interrupter>
typename TreeType::Ptr
solve(const TreeType& inTree, math::pcg::State& state, Interrupter& interrupter, bool staggered)
{
DirichletBoundaryOp<LaplacianMatrix::ValueType> boundaryOp;
return solveWithBoundaryConditions(inTree, boundaryOp, state, interrupter, staggered);
}
template<typename TreeType, typename BoundaryOp, typename Interrupter>
typename TreeType::Ptr
solveWithBoundaryConditions(const TreeType& inTree, const BoundaryOp& boundaryOp,
math::pcg::State& state, Interrupter& interrupter, bool staggered)
{
using DefaultPrecondT = math::pcg::IncompleteCholeskyPreconditioner<LaplacianMatrix>;
return solveWithBoundaryConditionsAndPreconditioner<DefaultPrecondT>(
inTree, boundaryOp, state, interrupter, staggered);
}
template<
typename PreconditionerType,
typename TreeType,
typename BoundaryOp,
typename Interrupter>
typename TreeType::Ptr
solveWithBoundaryConditionsAndPreconditioner(
const TreeType& inTree,
const BoundaryOp& boundaryOp,
math::pcg::State& state,
Interrupter& interrupter,
bool staggered)
{
return solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>(
/*source=*/inTree, /*domain mask=*/inTree, boundaryOp, state, interrupter, staggered);
}
template<
typename PreconditionerType,
typename TreeType,
typename DomainTreeType,
typename BoundaryOp,
typename Interrupter>
typename TreeType::Ptr
solveWithBoundaryConditionsAndPreconditioner(
const TreeType& inTree,
const DomainTreeType& domainMask,
const BoundaryOp& boundaryOp,
math::pcg::State& state,
Interrupter& interrupter,
bool staggered)
{
using TreeValueT = typename TreeType::ValueType;
using VecValueT = LaplacianMatrix::ValueType;
using VectorT = typename math::pcg::Vector<VecValueT>;
using VIdxTreeT = typename TreeType::template ValueConverter<VIndex>::Type;
using MaskTreeT = typename TreeType::template ValueConverter<bool>::Type;
// 1. Create a mapping from active voxels of the input tree to elements of a vector.
typename VIdxTreeT::ConstPtr idxTree = createIndexTree(domainMask);
// 2. Populate a vector with values from the input tree.
typename VectorT::Ptr b = createVectorFromTree<VecValueT>(inTree, *idxTree);
// 3. Create a mask of the interior voxels of the input tree (from the densified index tree).
/// @todo Is this really needed?
typename MaskTreeT::Ptr interiorMask(
new MaskTreeT(*idxTree, /*background=*/false, TopologyCopy()));
tools::erodeActiveValues(*interiorMask, /*iterations=*/1, tools::NN_FACE, tools::IGNORE_TILES);
// 4. Create the Laplacian matrix.
LaplacianMatrix::Ptr laplacian = createISLaplacianWithBoundaryConditions(
*idxTree, *interiorMask, boundaryOp, *b, staggered);
// 5. Solve the Poisson equation.
laplacian->scale(-1.0); // matrix is negative-definite; solve -M x = -b
b->scale(-1.0);
typename VectorT::Ptr x(new VectorT(b->size(), zeroVal<VecValueT>()));
typename math::pcg::Preconditioner<VecValueT>::Ptr precond(
new PreconditionerType(*laplacian));
if (!precond->isValid()) {
precond.reset(new math::pcg::JacobiPreconditioner<LaplacianMatrix>(*laplacian));
}
state = math::pcg::solve(*laplacian, *b, *x, *precond, interrupter, state);
// 6. Populate the output tree with values from the solution vector.
/// @todo if (state.success) ... ?
return createTreeFromVector<TreeValueT>(*x, *idxTree, /*background=*/zeroVal<TreeValueT>());
}
////////////////////////////////////////
// Explicit Template Instantiation
#ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
#ifdef OPENVDB_INSTANTIATE_POISSONSOLVER
#include <openvdb/util/ExplicitInstantiation.h>
#endif
#define _FUNCTION(TreeT) \
TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
math::pcg::IncompleteCholeskyPreconditioner<LaplacianMatrix>>( \
const TreeT&, const TreeT&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
math::pcg::State&, util::NullInterrupter&, bool)
OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
#undef _FUNCTION
#define _FUNCTION(TreeT) \
TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
math::pcg::IncompleteCholeskyPreconditioner<LaplacianMatrix>>( \
const TreeT&, const BoolTree&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
math::pcg::State&, util::NullInterrupter&, bool)
OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
#undef _FUNCTION
#define _FUNCTION(TreeT) \
TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
math::pcg::JacobiPreconditioner<LaplacianMatrix>>( \
const TreeT&, const TreeT&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
math::pcg::State&, util::NullInterrupter&, bool)
OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
#undef _FUNCTION
#define _FUNCTION(TreeT) \
TreeT::Ptr solveWithBoundaryConditionsAndPreconditioner< \
math::pcg::JacobiPreconditioner<LaplacianMatrix>>( \
const TreeT&, const BoolTree&, const DirichletBoundaryOp<LaplacianMatrix::ValueType>&, \
math::pcg::State&, util::NullInterrupter&, bool)
OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
#undef _FUNCTION
#endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
} // namespace poisson
} // namespace tools
} // namespace OPENVDB_VERSION_NAME
} // namespace openvdb
#endif // OPENVDB_TOOLS_POISSONSOLVER_HAS_BEEN_INCLUDED
ConjGradient.h
Preconditioned conjugate gradient solver (solves Ax = b using the conjugate gradient method with one ...
ExplicitInstantiation.h
value
ValueT value
Definition GridBuilder.h:1290
LeafManager.h
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
Morphology.h
Implementation of morphological dilation and erosion.
NullInterrupter.h
Tree.h
openvdb::v10_0::TopologyCopy
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition Types.h:644
openvdb::v10_0::math::Coord
Signed (x, y, z) 32-bit integer coordinates.
Definition Coord.h:25
openvdb::v10_0::math::Coord::offsetBy
Coord offsetBy(Int32 dx, Int32 dy, Int32 dz) const
Definition Coord.h:91
openvdb::v10_0::math::pcg::IncompleteCholeskyPreconditioner
Preconditioner using incomplete Cholesky factorization.
Definition ConjGradient.h:1347
openvdb::v10_0::math::pcg::JacobiPreconditioner
Diagonal preconditioner.
Definition ConjGradient.h:1281
openvdb::v10_0::math::pcg::Preconditioner::isValid
virtual bool isValid() const
Definition ConjGradient.h:473
openvdb::v10_0::math::pcg::Preconditioner::Ptr
SharedPtr< Preconditioner > Ptr
Definition ConjGradient.h:468
openvdb::v10_0::math::pcg::SparseStencilMatrix
Sparse, square matrix representing a 3D stencil operator of size STENCIL_SIZE.
Definition ConjGradient.h:238
openvdb::v10_0::math::pcg::SparseStencilMatrix< double, 7 >::ValueType
double ValueType
Definition ConjGradient.h:240
openvdb::v10_0::math::pcg::SparseStencilMatrix< double, 7 >::Ptr
SharedPtr< SparseStencilMatrix > Ptr
Definition ConjGradient.h:242
openvdb::v10_0::math::pcg::Vector
Lightweight, variable-length vector.
Definition ConjGradient.h:139
openvdb::v10_0::math::pcg::Vector::Ptr
SharedPtr< Vector > Ptr
Definition ConjGradient.h:142
openvdb::v10_0::tree::LeafManager
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition LeafManager.h:85
openvdb::v10_0::tree::ValueAccessor
Definition ValueAccessor.h:191
openvdb::v10_0::internal
Definition openvdb.h:114
openvdb::v10_0::math::pcg::SizeType
Index32 SizeType
Definition ConjGradient.h:33
openvdb::v10_0::math::pcg::solve
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
openvdb::v10_0::tools::poisson
Definition PoissonSolver.h:77
openvdb::v10_0::tools::poisson::createISLaplacianWithBoundaryConditions
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
openvdb::v10_0::tools::poisson::solveWithBoundaryConditions
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
openvdb::v10_0::tools::poisson::populateIndexTree
void populateIndexTree(VIndexTreeType &)
Overwrite each active voxel in the given scalar tree with a sequential index, starting from zero.
Definition PoissonSolver.h:320
openvdb::v10_0::tools::poisson::createISLaplacian
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
openvdb::v10_0::tools::poisson::VIndex
Int32 VIndex
Definition PoissonSolver.h:80
openvdb::v10_0::tools::poisson::createIndexTree
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
openvdb::v10_0::tools::poisson::LaplacianMatrix
math::pcg::SparseStencilMatrix< double, 7 > LaplacianMatrix
The type of a matrix used to represent a three-dimensional Laplacian operator.
Definition PoissonSolver.h:83
openvdb::v10_0::tools::poisson::createVectorFromTree
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
openvdb::v10_0::tools::poisson::createTreeFromVector
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
openvdb::v10_0::tools::poisson::solve
TreeType::Ptr solve(const TreeType &, math::pcg::State &, 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:749
openvdb::v10_0::tools::poisson::solveWithBoundaryConditionsAndPreconditioner
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
openvdb::v10_0::tools::NN_FACE
@ NN_FACE
Definition Morphology.h:59
openvdb::v10_0::tools::interiorMask
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
openvdb::v10_0::tools::erodeActiveValues
void erodeActiveValues(TreeOrLeafManagerT &tree, const int iterations=1, const NearestNeighbors nn=NN_FACE, const TilePolicy mode=PRESERVE_TILES, const bool threaded=true)
Topologically erode all active values (i.e. both voxels and tiles) in a tree using one of three neare...
Definition Morphology.h:1133
openvdb::v10_0::tools::IGNORE_TILES
@ IGNORE_TILES
Definition Morphology.h:81
openvdb::v10_0::tools::laplacian
GridType::Ptr laplacian(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the Laplacian of the given scalar grid.
Definition GridOperators.h:1018
openvdb::v10_0::tree
Definition PointDataGrid.h:169
openvdb::v10_0::zeroVal
T zeroVal()
Return the value of type T that corresponds to zero.
Definition Math.h:70
openvdb::v10_0::Int32
int32_t Int32
Definition Types.h:56
openvdb::v10_0::Index64
uint64_t Index64
Definition Types.h:53
openvdb
Definition Exceptions.h:13
Types.h
openvdb.h
openvdb::v10_0::math::pcg::State
Information about the state of a conjugate gradient solution.
Definition ConjGradient.h:46
openvdb::v10_0::tools::poisson::DirichletBoundaryOp
Dirichlet boundary condition functor.
Definition PoissonSolver.h:269
openvdb::v10_0::tools::poisson::DirichletBoundaryOp::operator()
void operator()(const Coord &, const Coord &, ValueType &, ValueType &diag) const
Definition PoissonSolver.h:270
openvdb::v10_0::util::NullInterrupter
Base class for interrupters.
Definition NullInterrupter.h:26
OPENVDB_VERSION_NAME
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
OPENVDB_USE_VERSION_NAMESPACE
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212
OPENVDB_REAL_TREE_INSTANTIATE
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition version.h.in:157