各位好
关于mpm88里面，总是先有这么一行
base=int(x[p]/dx-0.5)
fx=x[p]/dx-base
请教一下这两个变量到底是什么含义啊？
为什么不直接用x而要非得减去0.5呢？
而且后面的形函数w，为什么要带入X=fx, fx-1, fx-2啊？
后面又为什么非得用w[0]*w[1]?
3 个赞
您好，首先非常欢迎您加入Taichi社区。
关于您的问题，这涉及到MPM的具体数学公式。
- 关于base和fx的求解，可以参考here 的公式121.
- 其次形函数
w的使用，可以参考上面文档的Page 33.
1 个赞
您好
感谢您的回答
我认真读过您所发的 蒋老师的 MPM course
然而
我所疑问的是
为什么该插值在代码中是这样处理的。关于这一点，我思考了很多天都没有理解。
尤其是base 和fx的具体几何含义。如果可以，可否用图片标注fx和base的具体所指？
以及C++注释版本中，提到把fx fx-1 fx-2分别带入公式的x中，我实在不能理解为什么这样做。
祝好
嗯嗯，这是一个很好的问题，确实很容易让你迷惑。
首先，减去0.5的原因是我们一般认为粒子在网格的中间。胡老师在Games 201, Lec 7, page 16，讲过为什么减去0.5，可以在B站具体听一下讲解here。
MPM中P2G阶段，需要将particle上的信息按照权重 w ，分散到周围的网格点上，这个例子里面就是 3\times3 的网格大小。
base 就是 这个网格的左下角网格点的位置。fx就是particle 到base网格点的差值。
我们可以以 p=[1.2,0.6], dx = 1.0 为例，研究一下权重 w 的求解。
我们这个例子中使用的是Quadratic kernels
那么 p 在 x 轴上，根据距离差可以算出和对应网格顶点的权重
对base的影响（权重 w ）应该是 N(1.2)
对B的影响（权重 w ）应该是 N(0.2)
对C的影响（权重 w ）应该是 N(-0.8)
你现在应该都还是很清楚的，但是你会看到这个公式和代码中的公式是不一样的。这就要涉及到一些数学了。我们将核函数 N 根据定义域分为三个部分 [-1.5, -0.5], [-0.5, 0.5], [0.5,1.5] 。这三个部分都可以用 [0.5,1.5] 的三个函数来表示。依然用图片给你展示一下：
-
N(x) \in [0.5,1.5] = f(x) \in [0.5,1.5]
-
N(x) \in [-0.5,0.5] = g(x) \in [0.5,1.5]
-
N(x) \in [-1.5,-0.5] = h(x) \in [0.5,1.5]
综上，你可以看出来
对base的影响（权重 w ）应该是 N(1.2), 也就是 f(1.2)
对B的影响（权重 w ）应该是 N(0.2), 也就是 g(1.2)
对C的影响（权重 w ）应该是 N(-0.8) ， 也就是 h(1.2)
其中 1.2 = fx.x ，
那么 mpm88.py的Line 37
w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]
中的 w[0].x = f(1.2) = N(1.2), w[1].x = g(1.2) = N(0.2) , w[2].x = h(1.2) = N(-0.8) 。
Y轴同理 w[0].y = f(0.6)=N(0.6), w[1].y = g(0.6)=N(-0.4), w[2].y = h(0.6) = N(-1.4)
7 个赞
过这么久了，自己动手写了一下SPH，我总算是明白为什么会减去0.5了。
不知道我想的对不对，如果不对请指正。
这是因为粒子在右边界和上边界计算出的网格编号会溢出，减去0.5相当于把整张网格向左下方平移了0.5个单位。
为什么会溢出，是这样的：
假设整个计算区域大小是100x100， 网格尺寸是1x1
我们利用int()来向下取整，假如粒子在上边界，比如position=（50，100），那么，当粒子的y坐标100转换为网格ID的时候，由于是向下取整的，网格的y编号就会变为100，但是网格的y编号的范围是0~99(因为从0开始算100个网格)。这就导致计算的网格编号超出范围了！
有边界也是同理的。
假如不是恰好在边界，而是稍微超过边界，比如100.1，那也会造成溢出。
但是下边界和左边界就没事，因为-0.1取整之后还是0，没有超出范围。所以左下边界是本身就具有一个单位的裕量的。
解决方案就是把整个网格整体向左下平移0.5个单位。这样原本100.1也会变成99.6，取整之后变为99，不会超出范围。移动之后，上下左右四个边界都拥有了半个单位的裕量。
同时，还应该注意边界处理的时候，把超过边界的值挪动到离着边界还有一个裕量的位置，而不是恰好在边界上，比如100-eps，eps可以设置为0.1
同时，最好给边界多预留一些距离，比如MPM99里面处理边界的时候，就在离着边界的距离小于3个单位的时候，就假设它穿过了边界，然后处理。
3 个赞
感谢解答，对于理解代码很有用。
当0.5<=x<1.5, base node 为0
当1.5<=1.5<2.5, base node 为1
…
但是我有一个疑问, 当粒子 x<0.5时，base node 依然为0，但是weight的计算（w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]）并不适用，对于这些粒子，怎么处理？
Hi @xiangzi , 应该还是适用的。你可以跟着代码和公式过一遍。即使base node 还是0， 也还是可以以base node 为原点，遍历3x3的网格顶点，而对应的weight计算公式还是不变。
假设找到一个粒子，使得fx=0.1, base node 是0， 那么把fx代入w计算得到的weight是=[0.98, -0.06, 0.08], 这是不合理的，是我哪里理解错了吗
请问，计算出来的fx在x方向是fx.x=1.2，然后带入w；但为什么y方向的fx.y却是不一样的值f(0.6),g(1.2),h(1.2)；以x的f(x),g(x),h(x)来计算，在计算出来的y方向的不应该是fx.y=0.6,f(0.6),g(0.6),h(0.6)这样吗？
1 个赞
同学，谢谢你。我仔细看看！
谢谢指出错误，已经更改了。
大佬您好，我带了几个点检验了一下胡老师代码(pic_vs_apic.py)里计算base的方法，得到的结果似乎是B点（借用一下您的图）
当dx等于网格维度的倒数的时候，好像base是当前位置左下角的网格点。这里有点没搞明白。。。
于是我按照自己的理解写了一个计算base的方法。base是距离该位置最近的网格点，然后以base为中心计算9个点的权重，最后完成P2G和G2P操作。如图所示，红框是以D点为中心的网格，在红框范围内的粒子计算出来的base就是D点。
代码如下：
@ti.kernel
def my_PIC():
for p in x:
# 胡老师的算法
# base = (x[p] * inv_dx - 0.5).cast(int)
# fx = x[p] * inv_dx - base.cast(float)
# Quadratic B-spline
# w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2]
# 我的理解
base = ti.round((x[p] * inv_dx - 0.5), int) # 根据粒子位置计算最近网格位置,以该位置作为中心网格取3*3个网格点分别计算权值
fx = x[p] * inv_dx - base.cast(float) - 0.5 # 粒子位置-中心网格位置=粒子偏移中心网格的量
# 以x轴为例，左中右与粒子位置的距离分别为：左=1+fx.x , 中=fx.x , 右=1-fx.x
# 又因为-0.5<fx.x<0.5，所以0.5<1+fx.x<1.5，0.5<1-fx.x<1.5
# 所以w计算如下
w = [0.5 * (0.5 - fx) ** 2, 0.75 - fx ** 2, 0.5 * (0.5 + fx) ** 2]
for i in ti.static(range(3)):
for j in ti.static(range(3)):
# offset = ti.Vector([i, j]) # 胡老师的算法
offset = ti.Vector([i - 1, j - 1]) # 我的理解：因为base是中心点，所以他的左下邻居的偏移量应该是[-1,-1]
weight = w[i][0] * w[j][1]
grid_v[base + offset] += weight * v[p]
grid_m[base + offset] += weight
for i, j in grid_m:
if grid_m[i, j] > 0:
inv_m = 1 / grid_m[i, j]
grid_v[i, j] = inv_m * grid_v[i, j]
for p in x:
# base = (x[p] * inv_dx - 0.5).cast(int)
# fx = x[p] * inv_dx - base.cast(float)
# Quadratic B-spline
# w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2]
# 我的理解
base = ti.round((x[p] * inv_dx - 0.5), int) # 根据粒子位置计算最近网格位置,以该位置作为中心网格取3*3个网格点分别计算权值
fx = x[p] * inv_dx - base.cast(float) - 0.5 # 粒子位置-网格中心位置=粒子偏移中心网格的量
w = [0.5 * (0.5 - fx) ** 2, 0.75 - fx ** 2, 0.5 * (0.5 + fx) ** 2]
new_v = ti.Vector.zero(ti.f32, 2)
for i in ti.static(range(3)):
for j in ti.static(range(3)):
weight = w[i][0] * w[j][1]
# new_v += weight * grid_v[base + ti.Vector([i, j])] # 胡老师的算法
new_v += weight * grid_v[base + ti.Vector([i - 1, j - 1])] # 我的理解：因为base是中心点，所以他的左下邻居的偏移量应该是[-1,-1]
x[p] = clamp_pos(x[p] + v[p] * dt)
v[p] = new_v
效果看起来和胡老师的代码是一样的。
我有两个疑问：
1.我的理解到底对不对。
2.如果我的理解是对的，那么我的方法和胡老师的方法在效率上哪个更好一点。
感谢大佬解惑！
1 个赞
base 点的坐标，你可以在python里跑一下看看
base.x = int(1.2 - 0.5) = int(0.7) = 0
base.y = int(0.6 - 0.5) = int(0.1) = 0
非常详细的解释哈，确实如你所说，是在黑色格点上计算权重的。
您好 请问 在flip 中为何一般都是在x y z 三个方向上为何是用 4X3X3 而不是 3X3X3来算权重？
@ti.func
def interp_grid(base, frac, vp):
# Index on sides
#idx_side = ti.Vector([base-1, base, base+1, base+2])
idx_side = [base-1, base, base+1, base+2]
# Weight on sides
#w_side = ti.Vector([quadratic_kernel(1.0+frac), quadratic_kernel(frac), quadratic_kernel(1.0-frac), quadratic_kernel(2.0-frac)])
w_side = [quadratic_kernel(1.0+frac), quadratic_kernel(frac), quadratic_kernel(1.0-frac), quadratic_kernel(2.0-frac)]
# Index on center
#idx_center = ti.Vector([base-1, base, base+1])
idx_center = [base-1, base, base+1]
# weight on center
#w_center= ti.Vector([quadratic_kernel(0.5+frac), quadratic_kernel(ti.abs(0.5-frac)), quadratic_kernel(ti.abs(0.5-frac)), quadratic_kernel(1.5-frac)])
w_center= [quadratic_kernel(0.5+frac), quadratic_kernel(ti.abs(0.5-frac)), quadratic_kernel(1.5-frac)]
#for i, j, k in ti.ndrange(4, 3, 3):
for i in ti.static(range(4)):
for j in ti.static(range(3)):
for k in ti.static(range(3)):
w = w_side[i].x * w_center[j].y * w_center[k].z
idx = (idx_side[i].x, idx_center[j].y, idx_center[k].z)
grid_v_x[idx] += vp.x * w
grid_w_x[idx] += w
#for i, j, k in ti.ndrange(3, 4, 3):
for i in ti.static(range(3)):
for j in ti.static(range(4)):
for k in ti.static(range(3)):
w = w_center[i].x * w_side[j].y * w_center[k].z
idx = (idx_center[i].x, idx_side[j].y, idx_center[k].z)
grid_v_y[idx] += vp.y *w
grid_w_y[idx] += w
#for i, j, k in ti.ndrange(3, 3, 4):
for i in ti.static(range(3)):
for j in ti.static(range(3)):
for k in ti.static(range(4)):
w = w_center[i].x * w_center[j].y * w_side[k].z
idx = (idx_center[i].x, idx_center[j].y, idx_side[k].z)
grid_v_z[idx] += vp.z * w
grid_w_z[idx] += w
这个很简单，你只要在P2G和P2G这两个部分，根据particle位置计算grid index的时候考虑grid的平移就可以了。
其实[3,3,3]也是可以的吧。