diff options
author | Sebastián Barschkis <sebbas@sebbas.org> | 2021-06-18 13:18:21 +0300 |
---|---|---|
committer | Sebastián Barschkis <sebbas@sebbas.org> | 2021-06-18 13:18:21 +0300 |
commit | adefdbc9dfa34eed505820fdd8a1569a19eb0cb7 (patch) | |
tree | 49dc3d174290ab143b8de6fb7a4ac4b9c7c0a699 /extern | |
parent | 7c681477094e16dd386201c364337dc068ec595a (diff) |
Fluid: Optimization for FLIP neighbor search radius
Contributed by @erik85 in D11400. The idea from this patch was placed in
a more generic context: A new FOR macro has been added that loops
over the neighbors of a cell within a given radius.
Diffstat (limited to 'extern')
-rw-r--r-- | extern/mantaflow/preprocessed/gitinfo.h | 2 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/kernel.h | 13 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/plugin/flip.cpp | 94 |
3 files changed, 57 insertions, 52 deletions
diff --git a/extern/mantaflow/preprocessed/gitinfo.h b/extern/mantaflow/preprocessed/gitinfo.h index 03fd0112095..6bc92278a33 100644 --- a/extern/mantaflow/preprocessed/gitinfo.h +++ b/extern/mantaflow/preprocessed/gitinfo.h @@ -1,3 +1,3 @@ -#define MANTA_GIT_VERSION "commit 9c505cd22e289b98c9aa717efba8ef3201c7e458" +#define MANTA_GIT_VERSION "commit 8fbebe02459b7f72575872c20961f7cb757db408" diff --git a/extern/mantaflow/preprocessed/kernel.h b/extern/mantaflow/preprocessed/kernel.h index 90e30cd21e1..dbcc2342a11 100644 --- a/extern/mantaflow/preprocessed/kernel.h +++ b/extern/mantaflow/preprocessed/kernel.h @@ -71,6 +71,19 @@ class ParticleBase; for (int j = bnd; j < (grid).getSizeY() - bnd; ++j) \ for (int i = bnd; i < (grid).getSizeX() - bnd; ++i) +#define FOR_NEIGHBORS_BND(grid, radius, bnd) \ + for (int zj = ((grid).is3D() ? std::max(bnd, k - radius) : 0); \ + zj <= ((grid).is3D() ? std::min(k + radius, (grid).getSizeZ() - 1 - bnd) : 0); \ + zj++) \ + for (int yj = std::max(bnd, j - radius); \ + yj <= std::min(j + radius, (grid).getSizeY() - 1 - bnd); \ + yj++) \ + for (int xj = std::max(bnd, i - radius); \ + xj <= std::min(i + radius, (grid).getSizeX() - 1 - bnd); \ + xj++) + +#define FOR_NEIGHBORS(grid, radius) FOR_NEIGHBORS_BND(grid, radius, 0) + //! Basic data structure for kernel data, initialized based on kernel type (e.g. single, idx, etc). struct KernelBase { int maxX, maxY, maxZ, minZ, maxT, minT; diff --git a/extern/mantaflow/preprocessed/plugin/flip.cpp b/extern/mantaflow/preprocessed/plugin/flip.cpp index d6fd3437959..8757958d4b0 100644 --- a/extern/mantaflow/preprocessed/plugin/flip.cpp +++ b/extern/mantaflow/preprocessed/plugin/flip.cpp @@ -822,33 +822,29 @@ struct ComputeUnionLevelsetPindex : public KernelBase { { const Vec3 gridPos = Vec3(i, j, k) + Vec3(0.5); // shifted by half cell Real phiv = radius * 1.0; // outside + const int r = int(radius) + 1; - int r = int(radius) + 1; - int rZ = phi.is3D() ? r : 0; - for (int zj = k - rZ; zj <= k + rZ; zj++) - for (int yj = j - r; yj <= j + r; yj++) - for (int xj = i - r; xj <= i + r; xj++) { - if (!phi.isInBounds(Vec3i(xj, yj, zj))) - continue; + FOR_NEIGHBORS(phi, r) + { - // note, for the particle indices in indexSys the access is periodic (ie, dont skip for - // eg inBounds(sx,10,10) - IndexInt isysIdxS = index.index(xj, yj, zj); - IndexInt pStart = index(isysIdxS), pEnd = 0; - if (phi.isInBounds(isysIdxS + 1)) - pEnd = index(isysIdxS + 1); - else - pEnd = indexSys.size(); - - // now loop over particles in cell - for (IndexInt p = pStart; p < pEnd; ++p) { - const int psrc = indexSys[p].sourceIndex; - if (ptype && ((*ptype)[psrc] & exclude)) - continue; - const Vec3 pos = parts[psrc].pos; - phiv = std::min(phiv, fabs(norm(gridPos - pos)) - radius); - } - } + // note, for the particle indices in indexSys the access is periodic (ie, dont skip for eg + // inBounds(sx,10,10) + IndexInt isysIdxS = index.index(xj, yj, zj); + IndexInt pStart = index(isysIdxS), pEnd = 0; + if (phi.isInBounds(isysIdxS + 1)) + pEnd = index(isysIdxS + 1); + else + pEnd = indexSys.size(); + + // now loop over particles in cell + for (IndexInt p = pStart; p < pEnd; ++p) { + const int psrc = indexSys[p].sourceIndex; + if (ptype && ((*ptype)[psrc] & exclude)) + continue; + const Vec3 pos = parts[psrc].pos; + phiv = std::min(phiv, fabs(norm(gridPos - pos)) - radius); + } + } phi(i, j, k) = phiv; } inline const Grid<int> &getArg0() @@ -1026,39 +1022,35 @@ struct ComputeAveragedLevelsetWeight : public KernelBase { // loop over neighborhood, similar to ComputeUnionLevelsetPindex const Real sradiusInv = 1. / (4. * radius * radius); - int r = int(1. * radius) + 1; - int rZ = phi.is3D() ? r : 0; + const int r = int(radius) + 1; // accumulators Real wacc = 0.; Vec3 pacc = Vec3(0.); Real racc = 0.; - for (int zj = k - rZ; zj <= k + rZ; zj++) - for (int yj = j - r; yj <= j + r; yj++) - for (int xj = i - r; xj <= i + r; xj++) { - if (!phi.isInBounds(Vec3i(xj, yj, zj))) - continue; + FOR_NEIGHBORS(phi, r) + { - IndexInt isysIdxS = index.index(xj, yj, zj); - IndexInt pStart = index(isysIdxS), pEnd = 0; - if (phi.isInBounds(isysIdxS + 1)) - pEnd = index(isysIdxS + 1); - else - pEnd = indexSys.size(); - for (IndexInt p = pStart; p < pEnd; ++p) { - IndexInt psrc = indexSys[p].sourceIndex; - if (ptype && ((*ptype)[psrc] & exclude)) - continue; + IndexInt isysIdxS = index.index(xj, yj, zj); + IndexInt pStart = index(isysIdxS), pEnd = 0; + if (phi.isInBounds(isysIdxS + 1)) + pEnd = index(isysIdxS + 1); + else + pEnd = indexSys.size(); + for (IndexInt p = pStart; p < pEnd; ++p) { + IndexInt psrc = indexSys[p].sourceIndex; + if (ptype && ((*ptype)[psrc] & exclude)) + continue; - Vec3 pos = parts[psrc].pos; - Real s = normSquare(gridPos - pos) * sradiusInv; - // Real w = std::max(0., cubed(1.-s) ); - Real w = std::max(0., (1. - s)); // a bit smoother - wacc += w; - racc += radius * w; - pacc += pos * w; - } - } + Vec3 pos = parts[psrc].pos; + Real s = normSquare(gridPos - pos) * sradiusInv; + // Real w = std::max(0., cubed(1.-s) ); + Real w = std::max(0., (1. - s)); // a bit smoother + wacc += w; + racc += radius * w; + pacc += pos * w; + } + } if (wacc > VECTOR_EPSILON) { racc /= wacc; |