【CG物理模拟系列】流体模拟--粒子法之SPH法的加权函数计算
来源:互联网 发布:软件中心 编辑:程序博客网 时间:2024/03/29 10:19
SPH法中的加权函数
Poly6 kernel[1]
用于密度计算等。由于r只存在于2次项中,可以省去平方根的计算。
梯度为,
拉普拉斯算子为,
如图所示,
右侧坐标为拉普拉斯值
Code Sample
/*!
* Poly6 kernel函数值的计算
* @param[in] r 距离
* @param[in] h 有效半径
* @return 函数值
*/inline double RxKernelPoly6(const double &r, const double &h)
{
if(r >= 0.0 && r < h){
double q = h*h-r*r;
return 315.0/(64.0*RX_PI*pow(h, 9))*q*q*q;
}
else{
return 0.0;
}
}
/*!
* Poly6 kernel梯度值的计算
* @param[in] r 距离
* @param[in] rij 相对位置
* @param[in] h 有效半径
* @return 梯度值
*/inline Vec3 RxKernelPoly6G(const double &r, const Vec3 &rij, const double &h)
{
if(r >= 0.0 && r < h){
double q = h*h-r*r;
return -945.0/(32.0*RX_PI*pow(h, 9))*q*q*rij;
}
else{
return Vec3(0.0);
}
}
/*!
* Poly6 kernel函数拉普拉斯算子的计算
* @param[in] r 距离
* @param[in] h 有效半径
* @return 拉普拉斯值
*/inline double RxKernelPoly6L(const double &r, const double &h)
{
if(r >= 0.0 && r < h){
double q = h*h-r*r;
return -945.0/(32.0*RX_PI*pow(h, 9))*(3*q*q-4*r*r*q);
}
else{
return 0.0;
}
}
h为定值的时候,上式的α部分作为全局变量计算即可。
↑
Spiky kernel[1,2]
用于计算圧力项。当Poly6 kernel的r=0时,其梯度也为0,把粒子当作是由cluster构成的问题求解即可。
梯度为,
拉普拉斯值为,
[2]的定义为,
如图所示,
右侧轴为拉普拉斯值
Code Sample
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
/*! * Spiky kernel函数值的计算 * @param[in] r 距离 * @param[in] h 有效半径 * @return 函数值 */inline double RxKernelSpiky(const double &r, const double &h){ if(r >= 0.0 && r < h){ double q = h-r; return 15/(RX_PI*pow(h, 6))*q*q*q; } else{ return 0.0; }} /*! * Spiky kernel函数梯度值的计算 * @param[in] r 距离 * @param[in] rij 相对位置 * @param[in] h 有效半径 * @return 梯度值 */inline Vec3 RxKernelSpikyG(const double &r, const Vec3 &rij, const double &h){ if(r >= 0.0 && r < h){ double q = h-r; return -45.0/(RX_PI*pow(h, 6))*q*q*rij/r; } else{ return Vec3(0.0); }} /*! * Spiky kernel函数拉普拉斯值的计算 * @param[in] r 距离 * @param[in] h 有效半径 * @return 拉普拉斯值 */inline double RxKernelSpikyL(const double &r, const double &h){ if(r >= 0.0 && r < h){ double q = h-r; return -90.0/(RX_PI*pow(h, 6))*(q*q/r-q); } else{ return 0.0; }}
h为定值的时候,上式的α部分当作全局变量计算即可。
↑
Viscosity kernel[1]
用于计算粘性扩散项。Laplacian为线性函数。
梯度为,
拉普拉斯算子为,
如图所示,
右侧轴为拉普拉斯值。
Code Sample
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
/*! * Viscosity kernel函数值的计算 * @param[in] r 距离 * @param[in] h 有效半径 * @return 函数值 */inline double RxKernelViscosity(const double &r, const double &h){ if(r >= 0.0 && r < h){ double q = r/h; return 15.0/(2.0*RX_PI*pow(h, 3))*(-0.5*q*q*q+q*q+0.5/q-1); } else{ return 0.0; }} /*! * Viscosity kernel函数梯度的计算 * @param[in] r 距离 * @param[in] rij 相对位置 * @param[in] h 有效半径 * @return 梯度值 */inline Vec3 RxKernelViscosityG(const double &r, const Vec3 &rij, const double &h){ if(r >= 0.0 && r < h){ double q = r/h; return 15.0/(2.0*RX_PI*pow(h, 5))*(-1.5*q+2-0.5/(q*q*q))*rij; } else{ return Vec3(0.0); }} /*! * Viscosity kernel函数拉普拉斯算子的计算 * @param[in] r 距离 * @param[in] h 有效半径 * @return Laplacian值 */inline double RxKernelViscosityL(const double &r, const double &h){ if(r >= 0.0 && r < h){ double q = h-r; return 45.0/(RX_PI*pow(h, 6))*q; } else{ return 0.0; }}
h为定值的时候,上式的α部分作为全局变量计算即可。
↑
Spline kernel[3,4]
这里,
剃度为,
这里,
拉普拉斯算子为,
这里,
如图所示,
右侧轴为拉普拉斯值。
↑
Super Gaussian kernel[3]
剃度为,
拉普拉斯算子为i,
如图所示,
右侧轴为拉普拉斯算子的值
コード例
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
/*! * Gaussian kernel函数值的计算 * @param[in] r 距离 * @param[in] h 有效半径 * @return 函数值 */inline double RxKernelSuperGaussian(const double &r, const double &h){ double q = r/h; return 1.0/(pow(RX_PI, 1.5)*pow(h, 3))*(2.5-r*r)*exp(-q*q); //return 1.0/(sqrt(RX_PI)*h)*exp(-q*q);} /*! * Gaussian kernel函数梯度值的计算 * @param[in] r 距离 * @param[in] rij 相对距离 * @param[in] h 有效半径 * @return 梯度值 */inline Vec3 RxKernelSuperGaussianG(const double &r, const Vec3 &rij, const double &h){ double q = r/h; return -2.0/(pow(RX_PI, 1.5)*pow(h, 3))*(1.0+2.5/(h*h)-q*q)*exp(-q*q)*rij;} /*! * Gaussian kernel函数拉普拉斯算子的计算 * @param[in] r 距离 * @param[in] h 有效半径 * @return 拉普拉斯算子的值 */inline double RxKernelSuperGaussianL(const double &r, const double &h){ double q = r/h; return 2.0/(pow(RX_PI, 1.5)*pow(h, 3))*(2.0*q*q-(1.0+2.5/(h*h)-q*q)*(1.0-2.0*q*q))*exp(-q*q);}
h为定值的时候,上式的α部分作为全局变量计算即可。
↑
Fourth-order weighting function
梯度函数为,
拉普拉斯算子为,
如图所示,
右侧轴为拉普拉斯值。
コード例
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
/*! * Fourth-order 加权函数值的计算 * @param[in] r 距离 * @param[in] h 有效半径 * @return 函数值 */inline double RxKernelFourth(const double &r, const double &h){ double q = 3*r/h; if(q >= 0.0 && q < 1.0){ double q3 = (3-q)*(3-q)*(3-q)*(3-q)*(3-q); double q2 = (2-q)*(2-q)*(2-q)*(2-q)*(2-q); double q1 = (1-q)*(1-q)*(1-q)*(1-q)*(1-q); return 9.0/(40.0*RX_PI*pow(h, 3))*(q3-6.0*q2+15.0*q1); } else if(q >= 1.0 && q < 2.0){ double q3 = (3-q)*(3-q)*(3-q)*(3-q)*(3-q); double q2 = (2-q)*(2-q)*(2-q)*(2-q)*(2-q); return 9.0/(40.0*RX_PI*pow(h, 3))*(q3-6.0*q2); } else if(q >= 2.0 && q < 3.0){ double q3 = (3-q)*(3-q)*(3-q)*(3-q)*(3-q); return 9.0/(40.0*RX_PI*pow(h, 3))*(q3); } else{ return 0.0; }} /*! * Fourth-order加权函数梯度值的计算 * @param[in] r 距离 * @param[in] rij 相对位置 * @param[in] h 有效半径 * @return 梯度值 */inline Vec3 RxKernelFourthG(const double &r, const Vec3 &rij, const double &h){ double q = 3*r/h; if(q >= 0.0 && q < 1.0){ double q3 = (3-q)*(3-q)*(3-q)*(3-q); double q2 = (2-q)*(2-q)*(2-q)*(2-q); double q1 = (1-q)*(1-q)*(1-q)*(1-q); return -27.0/(8.0*RX_PI*pow(h, 4))*(q3-6.0*q2+15.0*q1)*rij/r; } else if(q >= 1.0 && q < 2.0){ double q3 = (3-q)*(3-q)*(3-q)*(3-q); double q2 = (2-q)*(2-q)*(2-q)*(2-q); return -27.0/(8.0*RX_PI*pow(h, 4))*(q3-6.0*q2)*rij/r; } else if(q >= 2.0 && q < 3.0){ double q3 = (3-q)*(3-q)*(3-q)*(3-q); return -27.0/(8.0*RX_PI*pow(h, 4))*(q3)*rij/r; } else{ return Vec3(0.0); }} /*! * Fourth-order加权函数拉普拉斯算子的计算 * @param[in] r 距离 * @param[in] h 有效半径 * @return 拉普拉斯算子的值 */inline double RxKernelFourthL(const double &r, const double &h){ double q = 3*r/h; if(q >= 0.0 && q < 1.0){ double q3 = (3-q)*(3-q)*(3-q); double q2 = (2-q)*(2-q)*(2-q); double q1 = (1-q)*(1-q)*(1-q); return 81.0/(2.0*RX_PI*pow(h, 5))*(q3-6.0*q2+15.0*q1); } else if(q >= 1.0 && q < 2.0){ double q3 = (3-q)*(3-q)*(3-q); double q2 = (2-q)*(2-q)*(2-q); return 81.0/(2.0*RX_PI*pow(h, 5))*(q3-6.0*q2); } else if(q >= 2.0 && q < 3.0){ double q3 = (3-q)*(3-q)*(3-q); return 81.0/(2.0*RX_PI*pow(h, 5))*(q3); } else{ return 0.0; }}
h为定值的时候,上式的α部分作为全局变量计算即可。
↑
Kernel函数系数的计算
在SPH法的kernel函数W中,我们假设以下的条件成立。
由此,便可以计算出各kernel函数的系数。以Poly6 kernel函数()为例。
这里,为所求系数。
在平面空间中,使用有效半径为h的圆积分。 如下图所示,用微小范围积分。
使其值为1,得出。
2维条件下
下面讨论3维空间的计算。如下图左,以球中心到表面延伸出的锥形内的微小体积做计算。 各边的长为, , (图右), kernel函数的积分计算如下。
带入Poly6 kernel中计算。
的条件下,得出。
3维条件下
↑
参考文献
[1] Muller et al., Particle-Based Fluid Simulation for Interactive Applications, SCA2003.
[2] Desbrun and Cani, Smoothed Particles: A new paradigm for animating highly deformable bodies, Eurographics Workshop on Computer Animation and Simulation (EGCAS), 1996.
[3] Monaghan, Smoothed particle hydrodynamics, Annual Review of Astronomy and Astrophysics, 30, pp543-574, 1992.
[4] Becker and Teschner, Weakly compressible SPH for free surface flows, SCA2007.
0 0
- 【CG物理模拟系列】流体模拟--粒子法之SPH法的加权函数计算
- 【CG物理模拟系列】流体模拟--粒子法之SPH(理论)
- 【CG物理模拟系列】流体模拟--粒子法之SPH(实现)
- 【CG物理模拟系列】流体模拟--粒子法之SPH(代码讲解)
- 【CG物理模拟系列】流体模拟--粒子法之MPS法(理论)
- 【CG物理模拟系列】流体模拟--粒子法之Position Based Fluids
- 【CG物理模拟系列】粒子法--表面生成手法(上)
- 【CG物理模拟系列】粒子法--表面生成手法(下)
- 【CG物理模拟系列】弹性体模拟--Position-based法之Shape Matching
- 【CG物理模拟系列】弹性体模拟--Position-based法之Shape Matching(代码实现)
- 基于SPH的流体模拟实践和一些技巧总结
- 基于SPH的流体模拟实践和一些技巧总结
- 基于SPH的流体模拟实践和一些技巧总结
- 【CG物理模拟】风筝模拟
- 【CG物理模拟系列】开篇:介绍(上)
- 【CG物理模拟系列】开篇:介绍(下)
- 【CG物理模拟系列】弹性体模拟--Mass Spring及绳子模拟
- 粒子物理蒙特卡罗模拟库Geant4之能谱制作
- Generator 函数的含义与用法
- SQL模型类-QSqlTableModel模型
- POJ-1426-Find The Multiple(BFS DFS)
- poj 3678
- 移植Busybox
- 【CG物理模拟系列】流体模拟--粒子法之SPH法的加权函数计算
- 线程池的简单使用
- [dsu] codeforces 375D. Tree and Queries
- js获取上传文件的url
- 每天专篇技术贴——浅谈C++中内存分配方式
- ImportError: numpy.core.multiarray failed to import解决方法
- centos中添加sudo用户和日志
- 龙书D3D11章节习题答案(第十章)
- 树状数组