----
#contents
----

***Poly6 kernel[1] [#hcde5047]
密度計算など.rは2乗項しか含まれていないので,平方根を計算しなくてよいのがメリット.
#ref(eq_poly6.gif)
勾配は,
#ref(eq_poly6_gradient.gif)
ラプラシアンは,
#ref(eq_poly6_laplacian.gif)

グラフにすると,
#ref(kernel_poly6.png,,50%)
ラプラシアンの値は第2軸(右側)参照

コード例
#code(C){{
/*!
 * Poly6カーネル関数値の計算
 * @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カーネル関数勾配値の計算
 * @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カーネル関数ラプラシアンの計算
 * @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] [#yb8b47cb]
圧力項計算用.Poly6カーネルがr=0で勾配が0になり,パーティクルがクラスターを構成してしまう問題を解決する.
#ref(eq_spiky.gif)
勾配は,
#ref(eq_spiky_gradient.gif)
ラプラシアンは,
#ref(eq_spiky_laplacian.gif)
[2]での定義は,
#ref(eq_spiky_desbrun.gif)

グラフにすると,
#ref(kernel_spiky.png,,50%)
ラプラシアンの値は第2軸(右側)参照

コード例
#code(C){{
/*!
 * Spikyカーネル関数値の計算
 * @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カーネル関数勾配値の計算
 * @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カーネル関数ラプラシアンの計算
 * @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] [#a85e448a]
粘性拡散項計算用.Laplacianが線型な関数になる.
#ref(eq_visc.gif)
勾配は,
#ref(eq_visc_gradient.gif)
ラプラシアンは,
#ref(eq_visc_laplacian.gif)

グラフにすると,
#ref(kernel_viscosity.png,,50%)
ラプラシアンの値は第2軸(右側)参照

コード例
#code(C){{
/*!
 * Viscosityカーネル関数値の計算
 * @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カーネル関数勾配値の計算
 * @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カーネル関数ラプラシアンの計算
 * @param[in] r 距離
 * @param[in] h 有効半径
 * @return ラプラシアンの値
 */
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] [#o437ed20]
#ref(sph_kernel_spline.eq1.gif,nolink,70%)
ここで,
#ref(sph_kernel_spline.eq2.gif,nolink,70%)
勾配は,
#ref(sph_kernel_spline.eq3.gif,nolink,70%)
ここで,
#ref(sph_kernel_spline.eq4.gif,nolink,70%)
ラプラシアンは,
#ref(sph_kernel_spline.eq5.gif,nolink,70%)
ここで,
#ref(sph_kernel_spline.eq6.gif,nolink,70%)
グラフにすると,
#ref(kernel_spline.png,,50%)
ラプラシアンの値は第2軸(右側)参照



***Super Gaussian kernel[3] [#s1074f09]
#ref(eq_gaussian2.gif)
勾配は,
#ref(eq_gaussian2_gradient.gif)
ラプラシアンは,
#ref(eq_gaussian2_laplacian.gif)

グラフにすると,
#ref(kernel_gaussian2.png,,50%)
ラプラシアンの値は第2軸(右側)参照

コード例
#code(C){{
/*!
 * Gaussianカーネル関数値の計算
 * @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カーネル関数勾配値の計算
 * @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カーネル関数ラプラシアンの計算
 * @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 [#u853dbcc]
#ref(eq_fourth.gif)
勾配は,
#ref(eq_fourth_gradient.gif)
ラプラシアンは,
#ref(eq_fourth_laplacian.gif)

グラフにすると,
#ref(kernel_fourth.png,,50%)
ラプラシアンの値は第2軸(右側)参照

コード例
#code(C){{
/*!
 * 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が定数ならば,上式のα部分はグローバル変数にでもあらかじめ計算して入れておいて使った方がよい.



***参考文献 [#zeacb0e2]
>> [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.


トップ   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS