[Bf-blender-cvs] [8bdf191461a] master: Fluid: Added APIC simulation method
Sebastián Barschkis
noreply at git.blender.org
Fri Oct 30 09:52:58 CET 2020
Commit: 8bdf191461a6a09c71258e69246d13fa3034e357
Author: Sebastián Barschkis
Date: Fri Oct 30 00:09:03 2020 +0100
Branches: master
https://developer.blender.org/rB8bdf191461a6a09c71258e69246d13fa3034e357
Fluid: Added APIC simulation method
Basic support for velocity updates with the APIC method.
This commit adds APIC to the already existing dropdown menu for the simulation method. The APIC plugin within Mantaflow has been updated to the latest version.
===================================================================
M extern/mantaflow/preprocessed/plugin/apic.cpp
M intern/mantaflow/intern/MANTA_main.cpp
M intern/mantaflow/intern/strings/liquid_script.h
M release/scripts/startup/bl_ui/properties_physics_fluid.py
M source/blender/blenkernel/intern/fluid.c
M source/blender/makesrna/intern/rna_fluid.c
===================================================================
diff --git a/extern/mantaflow/preprocessed/plugin/apic.cpp b/extern/mantaflow/preprocessed/plugin/apic.cpp
index 6ff893014c9..100b26a4091 100644
--- a/extern/mantaflow/preprocessed/plugin/apic.cpp
+++ b/extern/mantaflow/preprocessed/plugin/apic.cpp
@@ -6,7 +6,7 @@
// ----------------------------------------------------------------------------
//
// MantaFlow fluid solver framework
-// Copyright 2016-2017 Kiwon Um, Nils Thuerey
+// Copyright 2016-2020 Kiwon Um, Nils Thuerey
//
// This program is free software, distributed under the terms of the
// Apache License, Version 2.0
@@ -21,6 +21,54 @@
namespace Manta {
+#define FOR_INT_IJK(num) \
+ for (int i = 0; i < num; ++i) \
+ for (int j = 0; j < num; ++j) \
+ for (int k = 0; k < num; ++k)
+
+static inline IndexInt indexUFace(const Vec3 &pos, const MACGrid &ref)
+{
+ const Vec3i f = toVec3i(pos), c = toVec3i(pos - 0.5);
+ const IndexInt index = f.x * ref.getStrideX() + c.y * ref.getStrideY() + c.z * ref.getStrideZ();
+ assertDeb(ref.isInBounds(index),
+ "Grid index out of bounds for particle position [" << pos.x << ", " << pos.y << ", "
+ << pos.z << "]");
+ return (ref.isInBounds(index)) ? index : -1;
+}
+
+static inline IndexInt indexVFace(const Vec3 &pos, const MACGrid &ref)
+{
+ const Vec3i f = toVec3i(pos), c = toVec3i(pos - 0.5);
+ const IndexInt index = c.x * ref.getStrideX() + f.y * ref.getStrideY() + c.z * ref.getStrideZ();
+ assertDeb(ref.isInBounds(index),
+ "Grid index out of bounds for particle position [" << pos.x << ", " << pos.y << ", "
+ << pos.z << "]");
+ return (ref.isInBounds(index)) ? index : -1;
+}
+
+static inline IndexInt indexWFace(const Vec3 &pos, const MACGrid &ref)
+{
+ const Vec3i f = toVec3i(pos), c = toVec3i(pos - 0.5);
+ const IndexInt index = c.x * ref.getStrideX() + c.y * ref.getStrideY() + f.z * ref.getStrideZ();
+ assertDeb(ref.isInBounds(index),
+ "Grid index out of bounds for particle position [" << pos.x << ", " << pos.y << ", "
+ << pos.z << "]");
+ return (ref.isInBounds(index)) ? index : -1;
+}
+
+static inline IndexInt indexOffset(
+ const IndexInt gidx, const int i, const int j, const int k, const MACGrid &ref)
+{
+ const IndexInt dX[2] = {0, ref.getStrideX()};
+ const IndexInt dY[2] = {0, ref.getStrideY()};
+ const IndexInt dZ[2] = {0, ref.getStrideZ()};
+ const IndexInt index = gidx + dX[i] + dY[j] + dZ[k];
+ assertDeb(ref.isInBounds(index),
+ "Grid index out of bounds for particle position [" << pos.x << ", " << pos.y << ", "
+ << pos.z << "]");
+ return (ref.isInBounds(index)) ? index : -1;
+}
+
struct knApicMapLinearVec3ToMACGrid : public KernelBase {
knApicMapLinearVec3ToMACGrid(const BasicParticleSystem &p,
MACGrid &mg,
@@ -30,7 +78,8 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
const ParticleDataImpl<Vec3> &cpy,
const ParticleDataImpl<Vec3> &cpz,
const ParticleDataImpl<int> *ptype,
- const int exclude)
+ const int exclude,
+ const int boundaryWidth)
: KernelBase(p.size()),
p(p),
mg(mg),
@@ -40,7 +89,8 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
cpy(cpy),
cpz(cpz),
ptype(ptype),
- exclude(exclude)
+ exclude(exclude),
+ boundaryWidth(boundaryWidth)
{
runMessage();
run();
@@ -54,73 +104,91 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
const ParticleDataImpl<Vec3> &cpy,
const ParticleDataImpl<Vec3> &cpz,
const ParticleDataImpl<int> *ptype,
- const int exclude)
+ const int exclude,
+ const int boundaryWidth)
{
if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
return;
- const IndexInt dX[2] = {0, vg.getStrideX()};
- const IndexInt dY[2] = {0, vg.getStrideY()};
- const IndexInt dZ[2] = {0, vg.getStrideZ()};
-
- const Vec3 &pos = p[idx].pos, &vel = vp[idx];
- const IndexInt fi = static_cast<IndexInt>(pos.x), fj = static_cast<IndexInt>(pos.y),
- fk = static_cast<IndexInt>(pos.z);
- const IndexInt ci = static_cast<IndexInt>(pos.x - 0.5),
- cj = static_cast<IndexInt>(pos.y - 0.5),
- ck = static_cast<IndexInt>(pos.z - 0.5);
- const Real wfi = clamp(pos.x - fi, Real(0), Real(1)),
- wfj = clamp(pos.y - fj, Real(0), Real(1)),
- wfk = clamp(pos.z - fk, Real(0), Real(1));
- const Real wci = clamp(Real(pos.x - ci - 0.5), Real(0), Real(1)),
- wcj = clamp(Real(pos.y - cj - 0.5), Real(0), Real(1)),
- wck = clamp(Real(pos.z - ck - 0.5), Real(0), Real(1));
- // TODO: check index for safety
+ if (!vg.isInBounds(p.getPos(idx), boundaryWidth)) {
+ debMsg("Skipping particle at index " << idx
+ << ". Is out of bounds and cannot be applied to grid.",
+ 1);
+ return;
+ }
+
+ const Vec3 &pos = p.getPos(idx), &vel = vp[idx];
+ const Vec3i f = toVec3i(pos);
+ const Vec3i c = toVec3i(pos - 0.5);
+ const Vec3 wf = clamp(pos - toVec3(f), Vec3(0.), Vec3(1.));
+ const Vec3 wc = clamp(pos - toVec3(c) - 0.5, Vec3(0.), Vec3(1.));
+
{ // u-face
- const IndexInt gidx = fi * dX[1] + cj * dY[1] + ck * dZ[1];
- const Vec3 gpos(fi, cj + 0.5, ck + 0.5);
- const Real wi[2] = {Real(1) - wfi, wfi};
- const Real wj[2] = {Real(1) - wcj, wcj};
- const Real wk[2] = {Real(1) - wck, wck};
- for (int i = 0; i < 2; ++i)
- for (int j = 0; j < 2; ++j)
- for (int k = 0; k < 2; ++k) {
- const Real w = wi[i] * wj[j] * wk[k];
- mg[gidx + dX[i] + dY[j] + dZ[k]].x += w;
- vg[gidx + dX[i] + dY[j] + dZ[k]].x += w * vel.x;
- vg[gidx + dX[i] + dY[j] + dZ[k]].x += w * dot(cpx[idx], gpos + Vec3(i, j, k) - pos);
- }
+ const IndexInt gidx = indexUFace(pos, vg);
+ if (gidx < 0)
+ return; // debug will fail
+
+ const Vec3 gpos(f.x, c.y + 0.5, c.z + 0.5);
+ const Real wi[2] = {Real(1) - wf.x, wf.x};
+ const Real wj[2] = {Real(1) - wc.y, wc.y};
+ const Real wk[2] = {Real(1) - wc.z, wc.z};
+
+ FOR_INT_IJK(2)
+ {
+ const Real w = wi[i] * wj[j] * wk[k];
+ const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+ if (vidx < 0)
+ continue; // debug will fail
+
+ mg[vidx].x += w;
+ vg[vidx].x += w * vel.x;
+ vg[vidx].x += w * dot(cpx[idx], gpos + Vec3(i, j, k) - pos);
+ }
}
{ // v-face
- const IndexInt gidx = ci * dX[1] + fj * dY[1] + ck * dZ[1];
- const Vec3 gpos(ci + 0.5, fj, ck + 0.5);
- const Real wi[2] = {Real(1) - wci, wci};
- const Real wj[2] = {Real(1) - wfj, wfj};
- const Real wk[2] = {Real(1) - wck, wck};
- for (int i = 0; i < 2; ++i)
- for (int j = 0; j < 2; ++j)
- for (int k = 0; k < 2; ++k) {
- const Real w = wi[i] * wj[j] * wk[k];
- mg[gidx + dX[i] + dY[j] + dZ[k]].y += w;
- vg[gidx + dX[i] + dY[j] + dZ[k]].y += w * vel.y;
- vg[gidx + dX[i] + dY[j] + dZ[k]].y += w * dot(cpy[idx], gpos + Vec3(i, j, k) - pos);
- }
+ const IndexInt gidx = indexVFace(pos, vg);
+ if (gidx < 0)
+ return;
+
+ const Vec3 gpos(c.x + 0.5, f.y, c.z + 0.5);
+ const Real wi[2] = {Real(1) - wc.x, wc.x};
+ const Real wj[2] = {Real(1) - wf.y, wf.y};
+ const Real wk[2] = {Real(1) - wc.z, wc.z};
+
+ FOR_INT_IJK(2)
+ {
+ const Real w = wi[i] * wj[j] * wk[k];
+ const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+ if (vidx < 0)
+ continue;
+
+ mg[vidx].y += w;
+ vg[vidx].y += w * vel.y;
+ vg[vidx].y += w * dot(cpy[idx], gpos + Vec3(i, j, k) - pos);
+ }
}
if (!vg.is3D())
return;
{ // w-face
- const IndexInt gidx = ci * dX[1] + cj * dY[1] + fk * dZ[1];
- const Vec3 gpos(ci + 0.5, cj + 0.5, fk);
- const Real wi[2] = {Real(1) - wci, wci};
- const Real wj[2] = {Real(1) - wcj, wcj};
- const Real wk[2] = {Real(1) - wfk, wfk};
- for (int i = 0; i < 2; ++i)
- for (int j = 0; j < 2; ++j)
- for (int k = 0; k < 2; ++k) {
- const Real w = wi[i] * wj[j] * wk[k];
- mg[gidx + dX[i] + dY[j] + dZ[k]].z += w;
- vg[gidx + dX[i] + dY[j] + dZ[k]].z += w * vel.z;
- vg[gidx + dX[i] + dY[j] + dZ[k]].z += w * dot(cpz[idx], gpos + Vec3(i, j, k) - pos);
- }
+ const IndexInt gidx = indexWFace(pos, vg);
+ if (gidx < 0)
+ return;
+
+ const Vec3 gpos(c.x + 0.5, c.y + 0.5, f.z);
+ const Real wi[2] = {Real(1) - wc.x, wc.x};
+ const Real wj[2] = {Real(1) - wc.y, wc.y};
+ const Real wk[2] = {Real(1) - wf.z, wf.z};
+
+ FOR_INT_IJK(2)
+ {
+ const Real w = wi[i] * wj[j] * wk[k];
+ const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+ if (vidx < 0)
+ continue;
+
+ mg[vidx].z += w;
+ vg[vidx].z += w * vel.z;
+ vg[vidx].z += w * dot(cpz[idx], gpos + Vec3(i, j, k) - pos);
+ }
}
}
inline const BasicParticleSystem &getArg0()
@@ -168,6 +236,11 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
return exclude;
}
typedef int type8;
+ inline const int &getArg9()
+ {
+ return boundaryWidth;
+ }
+ typedef int type9;
void runMessage()
{
debMsg("Executing kernel knApicMapLinearVec3ToMACGrid ", 3);
@@ -179,7 +252,7 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
{
const IndexInt _sz = size;
for (IndexInt i = 0; i < _sz; i++)
- op(i, p, mg, vg, vp, cpx, cpy, cpz, ptype,
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list