diff options
Diffstat (limited to 'extern/mantaflow/preprocessed/plugin/flip.cpp')
-rw-r--r-- | extern/mantaflow/preprocessed/plugin/flip.cpp | 94 |
1 files changed, 43 insertions, 51 deletions
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; |