| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // Copyright Contributors to the OpenVDB Project | ||
| 2 | // SPDX-License-Identifier: MPL-2.0 | ||
| 3 | |||
| 4 | /// @file unittest/TestPoissonSolver.cc | ||
| 5 | /// @authors D.J. Hill, Peter Cucka | ||
| 6 | |||
| 7 | #include <openvdb/openvdb.h> | ||
| 8 | #include <openvdb/Types.h> | ||
| 9 | #include <openvdb/math/ConjGradient.h> // for JacobiPreconditioner | ||
| 10 | #include <openvdb/tools/Composite.h> // for csgDifference/Union/Intersection | ||
| 11 | #include <openvdb/tools/LevelSetSphere.h> // for tools::createLevelSetSphere() | ||
| 12 | #include <openvdb/tools/LevelSetUtil.h> // for tools::sdfToFogVolume() | ||
| 13 | #include <openvdb/tools/MeshToVolume.h> // for createLevelSetBox() | ||
| 14 | #include <openvdb/tools/PoissonSolver.h> | ||
| 15 | |||
| 16 | #include <gtest/gtest.h> | ||
| 17 | |||
| 18 | #include <cmath> | ||
| 19 | |||
| 20 | |||
| 21 | 6 | class TestPoissonSolver: public ::testing::Test | |
| 22 | { | ||
| 23 | }; | ||
| 24 | |||
| 25 | |||
| 26 | //////////////////////////////////////// | ||
| 27 | |||
| 28 | |||
| 29 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | TEST_F(TestPoissonSolver, testIndexTree) |
| 30 | { | ||
| 31 | using namespace openvdb; | ||
| 32 | using tools::poisson::VIndex; | ||
| 33 | |||
| 34 | using VIdxTree = FloatTree::ValueConverter<VIndex>::Type; | ||
| 35 | using LeafNodeType = VIdxTree::LeafNodeType; | ||
| 36 | |||
| 37 | 2 | VIdxTree tree; | |
| 38 | /// @todo populate tree | ||
| 39 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | tree::LeafManager<const VIdxTree> leafManager(tree); |
| 40 | |||
| 41 | 1 | VIndex testOffset = 0; | |
| 42 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | for (size_t n = 0, N = leafManager.leafCount(); n < N; ++n) { |
| 43 | const LeafNodeType& leaf = leafManager.leaf(n); | ||
| 44 | ✗ | for (LeafNodeType::ValueOnCIter it = leaf.cbeginValueOn(); it; ++it, testOffset++) { | |
| 45 | ✗ | EXPECT_EQ(testOffset, *it); | |
| 46 | } | ||
| 47 | } | ||
| 48 | |||
| 49 | //if (testOffset != VIndex(tree.activeVoxelCount())) { | ||
| 50 | // std::cout << "--Testing offsetmap - " | ||
| 51 | // << testOffset<<" != " | ||
| 52 | // << tree.activeVoxelCount() | ||
| 53 | // << " has active tile count " | ||
| 54 | // << tree.activeTileCount()<<std::endl; | ||
| 55 | //} | ||
| 56 | |||
| 57 |
2/16✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
1 | EXPECT_EQ(VIndex(tree.activeVoxelCount()), testOffset); |
| 58 | 1 | } | |
| 59 | |||
| 60 | |||
| 61 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | TEST_F(TestPoissonSolver, testTreeToVectorToTree) |
| 62 | { | ||
| 63 | using namespace openvdb; | ||
| 64 | using tools::poisson::VIndex; | ||
| 65 | |||
| 66 | using VIdxTree = FloatTree::ValueConverter<VIndex>::Type; | ||
| 67 | |||
| 68 | FloatGrid::Ptr sphere = tools::createLevelSetSphere<FloatGrid>( | ||
| 69 |
1/2✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
1 | /*radius=*/10.f, /*center=*/Vec3f(0.f), /*voxelSize=*/0.25f); |
| 70 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | tools::sdfToFogVolume(*sphere); |
| 71 | FloatTree& inputTree = sphere->tree(); | ||
| 72 | |||
| 73 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | const Index64 numVoxels = inputTree.activeVoxelCount(); |
| 74 | |||
| 75 | // Generate an index tree. | ||
| 76 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | VIdxTree::Ptr indexTree = tools::poisson::createIndexTree(inputTree); |
| 77 |
1/16✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
1 | EXPECT_TRUE(bool(indexTree)); |
| 78 | |||
| 79 | // Copy the values of the active voxels of the tree into a vector. | ||
| 80 | math::pcg::VectorS::Ptr vec = | ||
| 81 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | tools::poisson::createVectorFromTree<float>(inputTree, *indexTree); |
| 82 |
2/16✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
1 | EXPECT_EQ(math::pcg::SizeType(numVoxels), vec->size()); |
| 83 | |||
| 84 | { | ||
| 85 | // Convert the vector back to a tree. | ||
| 86 | FloatTree::Ptr inputTreeCopy = tools::poisson::createTreeFromVector( | ||
| 87 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1 | *vec, *indexTree, /*bg=*/0.f); |
| 88 | |||
| 89 | // Check that voxel values were preserved. | ||
| 90 | FloatGrid::ConstAccessor inputAcc = sphere->getConstAccessor(); | ||
| 91 |
2/2✓ Branch 0 taken 267731 times.
✓ Branch 1 taken 1 times.
|
267732 | for (FloatTree::ValueOnCIter it = inputTreeCopy->cbeginValueOn(); it; ++it) { |
| 92 |
1/2✓ Branch 1 taken 267731 times.
✗ Branch 2 not taken.
|
267731 | const Coord ijk = it.getCoord(); |
| 93 | //if (!math::isApproxEqual(*it, inputTree.getValue(ijk))) { | ||
| 94 | // std::cout << " value error " << *it << " " | ||
| 95 | // << inputTree.getValue(ijk) << std::endl; | ||
| 96 | //} | ||
| 97 |
3/18✓ Branch 1 taken 267731 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 267731 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 267731 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
267731 | EXPECT_NEAR(inputAcc.getValue(ijk), *it, /*tolerance=*/1.0e-6); |
| 98 | } | ||
| 99 | } | ||
| 100 | 1 | } | |
| 101 | |||
| 102 | |||
| 103 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | TEST_F(TestPoissonSolver, testLaplacian) |
| 104 | { | ||
| 105 | using namespace openvdb; | ||
| 106 | using tools::poisson::VIndex; | ||
| 107 | |||
| 108 | using VIdxTree = FloatTree::ValueConverter<VIndex>::Type; | ||
| 109 | |||
| 110 | // For two different problem sizes, N = 8 and N = 20... | ||
| 111 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
|
3 | for (int N = 8; N <= 20; N += 12) { |
| 112 | // Construct an N x N x N volume in which the value of voxel (i, j, k) | ||
| 113 | // is sin(i) * sin(j) * sin(k), using a voxel spacing of pi / N. | ||
| 114 | 2 | const double delta = openvdb::math::pi<double>() / N; | |
| 115 | 4 | FloatTree inputTree(/*background=*/0.f); | |
| 116 | Coord ijk(0); | ||
| 117 | Int32 &i = ijk[0], &j = ijk[1], &k = ijk[2]; | ||
| 118 |
2/2✓ Branch 0 taken 26 times.
✓ Branch 1 taken 2 times.
|
28 | for (i = 1; i < N; ++i) { |
| 119 |
2/2✓ Branch 0 taken 410 times.
✓ Branch 1 taken 26 times.
|
436 | for (j = 1; j < N; ++j) { |
| 120 |
2/2✓ Branch 0 taken 7202 times.
✓ Branch 1 taken 410 times.
|
7612 | for (k = 1; k < N; ++k) { |
| 121 | 7202 | inputTree.setValue(ijk, static_cast<float>( | |
| 122 |
1/2✓ Branch 1 taken 7202 times.
✗ Branch 2 not taken.
|
7202 | std::sin(delta * i) * std::sin(delta * j) * std::sin(delta * k))); |
| 123 | } | ||
| 124 | } | ||
| 125 | } | ||
| 126 | const Index64 numVoxels = inputTree.activeVoxelCount(); | ||
| 127 | |||
| 128 | // Generate an index tree. | ||
| 129 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | VIdxTree::Ptr indexTree = tools::poisson::createIndexTree(inputTree); |
| 130 |
1/16✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
2 | EXPECT_TRUE(bool(indexTree)); |
| 131 | |||
| 132 | // Copy the values of the active voxels of the tree into a vector. | ||
| 133 | math::pcg::VectorS::Ptr source = | ||
| 134 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | tools::poisson::createVectorFromTree<float>(inputTree, *indexTree); |
| 135 |
2/18✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
2 | EXPECT_EQ(math::pcg::SizeType(numVoxels), source->size()); |
| 136 | |||
| 137 | // Create a mask of the interior voxels of the source tree. | ||
| 138 | 4 | BoolTree interiorMask(/*background=*/false); | |
| 139 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
|
4 | interiorMask.fill(CoordBBox(Coord(2), Coord(N-2)), /*value=*/true, /*active=*/true); |
| 140 | |||
| 141 | // Compute the Laplacian of the source: | ||
| 142 | // D^2 sin(i) * sin(j) * sin(k) = -3 sin(i) * sin(j) * sin(k) | ||
| 143 | tools::poisson::LaplacianMatrix::Ptr laplacian = | ||
| 144 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | tools::poisson::createISLaplacian(*indexTree, interiorMask, /*staggered=*/true); |
| 145 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | laplacian->scale(1.0 / (delta * delta)); // account for voxel spacing |
| 146 |
2/18✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
2 | EXPECT_EQ(math::pcg::SizeType(numVoxels), laplacian->size()); |
| 147 | |||
| 148 | math::pcg::VectorS result(source->size()); | ||
| 149 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | laplacian->vectorMultiply(*source, result); |
| 150 | |||
| 151 | // Dividing the result by the source should produce a vector of uniform value -3. | ||
| 152 | // Due to finite differencing, the actual ratio will be somewhat different, though. | ||
| 153 | const math::pcg::VectorS& src = *source; | ||
| 154 | const float expected = // compute the expected ratio using one of the corner voxels | ||
| 155 | 2 | float((3.0 * src[1] - 6.0 * src[0]) / (delta * delta * src[0])); | |
| 156 |
2/2✓ Branch 0 taken 7202 times.
✓ Branch 1 taken 2 times.
|
7204 | for (math::pcg::SizeType n = 0; n < result.size(); ++n) { |
| 157 |
1/2✓ Branch 1 taken 7202 times.
✗ Branch 2 not taken.
|
7202 | result[n] /= src[n]; |
| 158 |
2/16✓ Branch 1 taken 7202 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 7202 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
7202 | EXPECT_NEAR(expected, result[n], /*tolerance=*/1.0e-4); |
| 159 | } | ||
| 160 | } | ||
| 161 | 1 | } | |
| 162 | |||
| 163 | |||
| 164 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | TEST_F(TestPoissonSolver, testSolve) |
| 165 | { | ||
| 166 | using namespace openvdb; | ||
| 167 | |||
| 168 | FloatGrid::Ptr sphere = tools::createLevelSetSphere<FloatGrid>( | ||
| 169 |
1/2✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
1 | /*radius=*/10.f, /*center=*/Vec3f(0.f), /*voxelSize=*/0.25f); |
| 170 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | tools::sdfToFogVolume(*sphere); |
| 171 | |||
| 172 | math::pcg::State result = math::pcg::terminationDefaults<float>(); | ||
| 173 | 1 | result.iterations = 100; | |
| 174 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | result.relativeError = result.absoluteError = 1.0e-4; |
| 175 | |||
| 176 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | FloatTree::Ptr outTree = tools::poisson::solve(sphere->tree(), result); |
| 177 | |||
| 178 |
1/16✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
1 | EXPECT_TRUE(result.success); |
| 179 |
1/16✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
1 | EXPECT_TRUE(result.iterations < 60); |
| 180 | 1 | } | |
| 181 | |||
| 182 | |||
| 183 | //////////////////////////////////////// | ||
| 184 | |||
| 185 | |||
| 186 | namespace { | ||
| 187 | |||
| 188 | struct BoundaryOp { | ||
| 189 | void operator()(const openvdb::Coord& ijk, const openvdb::Coord& neighbor, | ||
| 190 | double& source, double& diagonal) const | ||
| 191 | { | ||
| 192 | ✗ | if (neighbor.x() == ijk.x() && neighbor.z() == ijk.z()) { | |
| 193 | // Workaround for spurious GCC 4.8 -Wstrict-overflow warning: | ||
| 194 | ✗ | const openvdb::Coord::ValueType dy = (ijk.y() - neighbor.y()); | |
| 195 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
200 | if (dy > 0) source -= 1.0; |
| 196 | 200 | else diagonal -= 1.0; | |
| 197 | } | ||
| 198 | } | ||
| 199 | }; | ||
| 200 | |||
| 201 | |||
| 202 | template<typename TreeType> | ||
| 203 | void | ||
| 204 | 4 | doTestSolveWithBoundaryConditions() | |
| 205 | { | ||
| 206 | using namespace openvdb; | ||
| 207 | |||
| 208 | using ValueType = typename TreeType::ValueType; | ||
| 209 | |||
| 210 | // Solve for the pressure in a cubic tank of liquid that is open at the top. | ||
| 211 | // Boundary conditions are P = 0 at the top, dP/dy = -1 at the bottom | ||
| 212 | // and dP/dx = 0 at the sides. | ||
| 213 | // | ||
| 214 | // P = 0 | ||
| 215 | // +------+ (N,-1,N) | ||
| 216 | // /| /| | ||
| 217 | // (0,-1,0) +------+ | | ||
| 218 | // | | | | dP/dx = 0 | ||
| 219 | // dP/dx = 0 | +----|-+ | ||
| 220 | // |/ |/ | ||
| 221 | // (0,-N-1,0) +------+ (N,-N-1,0) | ||
| 222 | // dP/dy = -1 | ||
| 223 | |||
| 224 | const int N = 9; | ||
| 225 | 4 | const ValueType zero = zeroVal<ValueType>(); | |
| 226 | const double epsilon = math::Delta<ValueType>::value(); | ||
| 227 | |||
| 228 | 8 | TreeType source(/*background=*/zero); | |
| 229 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
4 | source.fill(CoordBBox(Coord(0, -N-1, 0), Coord(N, -1, N)), /*value=*/zero); |
| 230 | |||
| 231 | math::pcg::State state = math::pcg::terminationDefaults<ValueType>(); | ||
| 232 | 4 | state.iterations = 100; | |
| 233 | 4 | state.relativeError = state.absoluteError = epsilon; | |
| 234 | |||
| 235 | 4 | util::NullInterrupter interrupter; | |
| 236 | |||
| 237 |
2/4✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
|
4 | typename TreeType::Ptr solution = tools::poisson::solveWithBoundaryConditions( |
| 238 | source, BoundaryOp(), state, interrupter, /*staggered=*/true); | ||
| 239 | |||
| 240 |
1/16✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
4 | EXPECT_TRUE(state.success); |
| 241 |
1/16✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
4 | EXPECT_TRUE(state.iterations < 60); |
| 242 | |||
| 243 | // Verify that P = -y throughout the solution space. | ||
| 244 |
2/2✓ Branch 0 taken 2000 times.
✓ Branch 1 taken 2 times.
|
4004 | for (typename TreeType::ValueOnCIter it = solution->cbeginValueOn(); it; ++it) { |
| 245 |
4/22✓ Branch 1 taken 2000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2000 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2000 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 2000 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
|
4000 | EXPECT_NEAR( |
| 246 | double(-it.getCoord().y()), double(*it), /*tolerance=*/10.0 * epsilon); | ||
| 247 | } | ||
| 248 | 4 | } | |
| 249 | |||
| 250 | } // unnamed namespace | ||
| 251 | |||
| 252 | |||
| 253 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | TEST_F(TestPoissonSolver, testSolveWithBoundaryConditions) |
| 254 | { | ||
| 255 | 1 | doTestSolveWithBoundaryConditions<openvdb::FloatTree>(); | |
| 256 | 1 | doTestSolveWithBoundaryConditions<openvdb::DoubleTree>(); | |
| 257 | 1 | } | |
| 258 | |||
| 259 | |||
| 260 | namespace { | ||
| 261 | |||
| 262 | openvdb::FloatGrid::Ptr | ||
| 263 | 3 | newCubeLS( | |
| 264 | const int outerLength, // in voxels | ||
| 265 | const int innerLength, // in voxels | ||
| 266 | const openvdb::Vec3I& centerIS, // in index space | ||
| 267 | const float dx, // grid spacing | ||
| 268 | bool openTop) | ||
| 269 | { | ||
| 270 | using namespace openvdb; | ||
| 271 | |||
| 272 | using BBox = math::BBox<Vec3f>; | ||
| 273 | |||
| 274 | // World space dimensions and center for this box | ||
| 275 | 3 | const float outerWS = dx * float(outerLength); | |
| 276 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
|
3 | const float innerWS = dx * float(innerLength); |
| 277 | Vec3f centerWS(centerIS); | ||
| 278 | centerWS *= dx; | ||
| 279 | |||
| 280 | // Construct world space bounding boxes | ||
| 281 | BBox outerBBox( | ||
| 282 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
|
3 | Vec3f(-outerWS / 2, -outerWS / 2, -outerWS / 2), |
| 283 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
|
3 | Vec3f( outerWS / 2, outerWS / 2, outerWS / 2)); |
| 284 | BBox innerBBox; | ||
| 285 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
|
3 | if (openTop) { |
| 286 | 1 | innerBBox = BBox( | |
| 287 | 1 | Vec3f(-innerWS / 2, -innerWS / 2, -innerWS / 2), | |
| 288 | 1 | Vec3f( innerWS / 2, innerWS / 2, outerWS)); | |
| 289 | } else { | ||
| 290 | 2 | innerBBox = BBox( | |
| 291 | 2 | Vec3f(-innerWS / 2, -innerWS / 2, -innerWS / 2), | |
| 292 | 2 | Vec3f( innerWS / 2, innerWS / 2, innerWS / 2)); | |
| 293 | } | ||
| 294 | outerBBox.translate(centerWS); | ||
| 295 | innerBBox.translate(centerWS); | ||
| 296 | |||
| 297 | 3 | math::Transform::Ptr xform = math::Transform::createLinearTransform(dx); | |
| 298 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | FloatGrid::Ptr cubeLS = tools::createLevelSetBox<FloatGrid>(outerBBox, *xform); |
| 299 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | FloatGrid::Ptr inside = tools::createLevelSetBox<FloatGrid>(innerBBox, *xform); |
| 300 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | tools::csgDifference(*cubeLS, *inside); |
| 301 | |||
| 302 | 3 | return cubeLS; | |
| 303 | } | ||
| 304 | |||
| 305 | |||
| 306 | class LSBoundaryOp | ||
| 307 | { | ||
| 308 | public: | ||
| 309 | 1 | LSBoundaryOp(const openvdb::FloatTree& lsTree): mLS(&lsTree) {} | |
| 310 |
2/8✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
|
5 | LSBoundaryOp(const LSBoundaryOp& other): mLS(other.mLS) {} |
| 311 | |||
| 312 |
2/2✓ Branch 0 taken 6758 times.
✓ Branch 1 taken 2914 times.
|
9672 | void operator()(const openvdb::Coord& ijk, const openvdb::Coord& neighbor, |
| 313 | double& source, double& diagonal) const | ||
| 314 | { | ||
| 315 | // Doing nothing is equivalent to imposing dP/dn = 0 boundary condition | ||
| 316 | |||
| 317 |
4/4✓ Branch 0 taken 6758 times.
✓ Branch 1 taken 2914 times.
✓ Branch 2 taken 3844 times.
✓ Branch 3 taken 2914 times.
|
9672 | if (neighbor.x() == ijk.x() && neighbor.y() == ijk.y()) { // on top or bottom |
| 318 |
2/2✓ Branch 1 taken 1922 times.
✓ Branch 2 taken 1922 times.
|
3844 | if (mLS->getValue(neighbor) <= 0.f) { |
| 319 | // closed boundary | ||
| 320 | 1922 | source -= 1.0; | |
| 321 | } else { | ||
| 322 | // open boundary | ||
| 323 | 1922 | diagonal -= 1.0; | |
| 324 | } | ||
| 325 | } | ||
| 326 | 9672 | } | |
| 327 | |||
| 328 | private: | ||
| 329 | const openvdb::FloatTree* mLS; | ||
| 330 | }; | ||
| 331 | |||
| 332 | } // unnamed namespace | ||
| 333 | |||
| 334 | |||
| 335 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | TEST_F(TestPoissonSolver, testSolveWithSegmentedDomain) |
| 336 | { | ||
| 337 | // In fluid simulations, incompressibility is enforced by the pressure, which is | ||
| 338 | // computed as a solution of a Poisson equation. Often, procedural animation | ||
| 339 | // of objects (e.g., characters) interacting with liquid will result in boundary | ||
| 340 | // conditions that describe multiple disjoint regions: regions of free surface flow | ||
| 341 | // and regions of trapped fluid. It is this second type of region for which | ||
| 342 | // there may be no consistent pressure (e.g., a shrinking watertight region | ||
| 343 | // filled with incompressible liquid). | ||
| 344 | // | ||
| 345 | // This unit test demonstrates how to use a level set and topological tools | ||
| 346 | // to separate the well-posed problem of a liquid with a free surface | ||
| 347 | // from the possibly ill-posed problem of fully enclosed liquid regions. | ||
| 348 | // | ||
| 349 | // For simplicity's sake, the physical boundaries are idealized as three | ||
| 350 | // non-overlapping cubes, one with an open top and two that are fully closed. | ||
| 351 | // All three contain incompressible liquid (x), and one of the closed cubes | ||
| 352 | // will be partially filled so that two of the liquid regions have a free surface | ||
| 353 | // (Dirichlet boundary condition on one side) while the totally filled cube | ||
| 354 | // would have no free surface (Neumann boundary conditions on all sides). | ||
| 355 | // ________________ ________________ | ||
| 356 | // __ __ | __________ | | __________ | | ||
| 357 | // | |x x x x x | | | | | | | |x x x x x | | | ||
| 358 | // | |x x x x x | | | |x x x x x | | | |x x x x x | | | ||
| 359 | // | |x x x x x | | | |x x x x x | | | |x x x x x | | | ||
| 360 | // | —————————— | | —————————— | | —————————— | | ||
| 361 | // |________________| |________________| |________________| | ||
| 362 | // | ||
| 363 | // The first two regions are clearly well-posed, while the third region | ||
| 364 | // may have no solution (or multiple solutions). | ||
| 365 | // -D.J.Hill | ||
| 366 | |||
| 367 | using namespace openvdb; | ||
| 368 | |||
| 369 | using PreconditionerType = | ||
| 370 | math::pcg::IncompleteCholeskyPreconditioner<tools::poisson::LaplacianMatrix>; | ||
| 371 | |||
| 372 | // Grid spacing | ||
| 373 | const float dx = 0.05f; | ||
| 374 | |||
| 375 | // Construct the solid boundaries in a single grid. | ||
| 376 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | FloatGrid::Ptr solidBoundary; |
| 377 | { | ||
| 378 | // Create three non-overlapping cubes. | ||
| 379 | const int outerDim = 41; | ||
| 380 | const int innerDim = 31; | ||
| 381 | FloatGrid::Ptr | ||
| 382 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1 | openDomain = newCubeLS(outerDim, innerDim, /*ctr=*/Vec3I(0, 0, 0), dx, /*open=*/true), |
| 383 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1 | closedDomain0 = newCubeLS(outerDim, innerDim, /*ctr=*/Vec3I(60, 0, 0), dx, false), |
| 384 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
1 | closedDomain1 = newCubeLS(outerDim, innerDim, /*ctr=*/Vec3I(120, 0, 0), dx, false); |
| 385 | |||
| 386 | // Union all three cubes into one grid. | ||
| 387 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | tools::csgUnion(*openDomain, *closedDomain0); |
| 388 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | tools::csgUnion(*openDomain, *closedDomain1); |
| 389 | |||
| 390 | // Strictly speaking the solidBoundary level set should be rebuilt | ||
| 391 | // (with tools::levelSetRebuild()) after the csgUnions to insure a proper | ||
| 392 | // signed distance field, but we will forgo the rebuild in this example. | ||
| 393 | solidBoundary = openDomain; | ||
| 394 | } | ||
| 395 | |||
| 396 | // Generate the source for the Poisson solver. | ||
| 397 | // For a liquid simulation this will be the divergence of the velocity field | ||
| 398 | // and will coincide with the liquid location. | ||
| 399 | // | ||
| 400 | // We activate by hand cells in distinct solution regions. | ||
| 401 | |||
| 402 |
1/2✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
2 | FloatTree source(/*background=*/0.f); |
| 403 | |||
| 404 | // The source is active in the union of the following "liquid" regions: | ||
| 405 | |||
| 406 | // Fill the open box. | ||
| 407 | const int N = 15; | ||
| 408 | 1 | CoordBBox liquidInOpenDomain(Coord(-N, -N, -N), Coord(N, N, N)); | |
| 409 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | source.fill(liquidInOpenDomain, 0.f); |
| 410 | |||
| 411 | // Totally fill closed box 0. | ||
| 412 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | CoordBBox liquidInClosedDomain0(Coord(-N, -N, -N), Coord(N, N, N)); |
| 413 | liquidInClosedDomain0.translate(Coord(60, 0, 0)); | ||
| 414 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | source.fill(liquidInClosedDomain0, 0.f); |
| 415 | |||
| 416 | // Half fill closed box 1. | ||
| 417 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | CoordBBox liquidInClosedDomain1(Coord(-N, -N, -N), Coord(N, N, 0)); |
| 418 | liquidInClosedDomain1.translate(Coord(120, 0, 0)); | ||
| 419 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | source.fill(liquidInClosedDomain1, 0.f); |
| 420 | |||
| 421 | // Compute the number of voxels in the well-posed region of the source. | ||
| 422 | const Index64 expectedWellPosedVolume = | ||
| 423 | 1 | liquidInOpenDomain.volume() + liquidInClosedDomain1.volume(); | |
| 424 | |||
| 425 | // Generate a mask that defines the solution domain. | ||
| 426 | // Inactive values of the source map to false and active values map to true. | ||
| 427 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
2 | const BoolTree totalSourceDomain(source, /*inactive=*/false, /*active=*/true, TopologyCopy()); |
| 428 | |||
| 429 | // Extract the "interior regions" from the solid boundary. | ||
| 430 | // The result will correspond to the the walls of the boxes unioned with inside of the full box. | ||
| 431 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | const BoolTree::ConstPtr interiorMask = tools::extractEnclosedRegion( |
| 432 | solidBoundary->tree(), /*isovalue=*/float(0), &totalSourceDomain); | ||
| 433 | |||
| 434 | // Identify the well-posed part of the problem. | ||
| 435 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
2 | BoolTree wellPosedDomain(source, /*inactive=*/false, /*active=*/true, TopologyCopy()); |
| 436 | wellPosedDomain.topologyDifference(*interiorMask); | ||
| 437 |
2/16✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
|
1 | EXPECT_EQ(expectedWellPosedVolume, wellPosedDomain.activeVoxelCount()); |
| 438 | |||
| 439 | // Solve the well-posed Poisson equation. | ||
| 440 | |||
| 441 | const double epsilon = math::Delta<float>::value(); | ||
| 442 | math::pcg::State state = math::pcg::terminationDefaults<float>(); | ||
| 443 | 1 | state.iterations = 200; | |
| 444 | 1 | state.relativeError = state.absoluteError = epsilon; | |
| 445 | |||
| 446 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | util::NullInterrupter interrupter; |
| 447 | |||
| 448 | // Define boundary conditions that are consistent with solution = 0 | ||
| 449 | // at the liquid/air boundary and with a linear response with depth. | ||
| 450 | LSBoundaryOp boundaryOp(solidBoundary->tree()); | ||
| 451 | |||
| 452 | // Compute the solution | ||
| 453 | FloatTree::Ptr wellPosedSolutionP = | ||
| 454 | tools::poisson::solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>( | ||
| 455 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | source, wellPosedDomain, boundaryOp, state, interrupter, /*staggered=*/true); |
| 456 | |||
| 457 |
3/20✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
|
1 | EXPECT_EQ(expectedWellPosedVolume, wellPosedSolutionP->activeVoxelCount()); |
| 458 |
1/16✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
1 | EXPECT_TRUE(state.success); |
| 459 |
1/16✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
|
1 | EXPECT_TRUE(state.iterations < 68); |
| 460 | |||
| 461 | // Verify that the solution is linear with depth. | ||
| 462 |
2/2✓ Branch 0 taken 45167 times.
✓ Branch 1 taken 1 times.
|
45168 | for (FloatTree::ValueOnCIter it = wellPosedSolutionP->cbeginValueOn(); it; ++it) { |
| 463 | Index32 depth; | ||
| 464 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 45167 times.
|
60543 | if (liquidInOpenDomain.isInside(it.getCoord())) { |
| 465 |
1/2✓ Branch 1 taken 29791 times.
✗ Branch 2 not taken.
|
29791 | depth = 1 + liquidInOpenDomain.max().z() - it.getCoord().z(); |
| 466 | } else { | ||
| 467 |
1/2✓ Branch 1 taken 15376 times.
✗ Branch 2 not taken.
|
15376 | depth = 1 + liquidInClosedDomain1.max().z() - it.getCoord().z(); |
| 468 | } | ||
| 469 |
3/18✓ Branch 1 taken 45167 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45167 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 45167 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
45167 | EXPECT_NEAR(double(depth), double(*it), /*tolerance=*/10.0 * epsilon); |
| 470 | } | ||
| 471 | |||
| 472 | #if 0 | ||
| 473 | // Optionally, one could attempt to compute the solution in the enclosed regions. | ||
| 474 | { | ||
| 475 | // Identify the potentially ill-posed part of the problem. | ||
| 476 | BoolTree illPosedDomain(source, /*inactive=*/false, /*active=*/true, TopologyCopy()); | ||
| 477 | illPosedDomain.topologyIntersection(source); | ||
| 478 | |||
| 479 | // Solve the Poisson equation in the two unconnected regions. | ||
| 480 | FloatTree::Ptr illPosedSoln = | ||
| 481 | tools::poisson::solveWithBoundaryConditionsAndPreconditioner<PreconditionerType>( | ||
| 482 | source, illPosedDomain, LSBoundaryOp(*solidBoundary->tree()), | ||
| 483 | state, interrupter, /*staggered=*/true); | ||
| 484 | } | ||
| 485 | #endif | ||
| 486 | 1 | } | |
| 487 |