请问关于mpm88的fx和base到底代表什么?

各位好
关于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]?

1 Like

您好,首先非常欢迎您加入Taichi社区。

关于您的问题,这涉及到MPM的具体数学公式。

  1. 关于base和fx的求解,可以参考here 的公式121.
  2. 其次形函数w的使用,可以参考上面文档的Page 33.
1 Like

您好
感谢您的回答

我认真读过您所发的 蒋老师的 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

那么 px 轴上,根据距离差可以算出和对应网格顶点的权重
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] 的三个函数来表示。依然用图片给你展示一下:

  1. N(x) \in [0.5,1.5] = f(x) \in [0.5,1.5]

  2. N(x) \in [-0.5,0.5] = g(x) \in [0.5,1.5]

  3. 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)

5 Likes

过这么久了,自己动手写了一下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个单位的时候,就假设它穿过了边界,然后处理。

2 Likes

感谢解答,对于理解代码很有用。

当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 Like

系统梳理了一遍原理~如有不严谨的地方请多指教 :grin:

/(ㄒoㄒ)/~~不让俺上传多张图片,那就只能传到B站上了

https://www.bilibili.com/read/cv16719283

1 Like

原理与程序部分:

例子:
https://www.bilibili.com/read/cv16719283

1 Like

同学,谢谢你。我仔细看看!

谢谢指出错误,已经更改了。