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.

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;