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.
This commit is contained in:
Sebastián Barschkis 2021-06-18 12:18:21 +02:00
parent 7c68147709
commit adefdbc9df
3 changed files with 56 additions and 51 deletions

View File

@ -1,3 +1,3 @@
#define MANTA_GIT_VERSION "commit 9c505cd22e289b98c9aa717efba8ef3201c7e458"
#define MANTA_GIT_VERSION "commit 8fbebe02459b7f72575872c20961f7cb757db408"

View File

@ -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;

View File

@ -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();
// 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);
}
}
// 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;