# MPM 雪、海绵

# 简介

在上一篇文章写完了 MLS-MPM 对于水的模型后，想着记录下其他材质的模拟方法。其实这一篇对现在的项目来说没什么卵用，但打算还是记录下，因为我知道未来在游戏肯定会用上实时模拟，现在所学的没用的东西都是为了未来做铺垫。那么这一篇就写其他弹性材质的其他数学模型。这就是用 MPM 的好处。就是对于其他弹性材料我们只需要换一个数学模型就可以模拟出来并且还可以把他们整合在一起产生互相交互，并且只需要在 P2G 上面做材料判断就可以实现，非常简单，不需要水用 FLIP，软体用 FEM 这样重新写一遍。

Elastic Object （塑形物体）：就是用 Neo Hookean 和 Corotated 模型，这次我们用 Corotated 模型，上面的是 Neo Hookean 模型，下面是 Corotated 模型

Elastoplastic objects（弹塑性物体）：例如雪、沙子这些就是软绵绵的又是分散的物体，有很多个模型，例如 ad-hoc boxing、Cam-clay、Drucker-prager，在实现方面，ad-hoc boxing 是最简单最不物理的，所以我们就用这个模型，后面比较复杂所以我也没看

# 雪

雪的解算需要定义弹塑性能量密度方程，就是这个方程将会作用于后面解算的力，那么方程是这样的 $\Psi(F)$ , 根据流动塑性理论我们可以把它分解为 elastic deformation gradient $F_E$ 和 plastic deformation gradient $F_P$ , 那么方程就变成了 $\Psi(F_E,F_P)$

, 而 elastic deformation gradient 部分我们用 Corotated 模型 [Stomakhin et al. 2012]，那么方程就变成了

其中 $J_E = det F_E,J_P = det F_P$ $J$ 就是体积的变化率而 $F_E= R_ES_E$ 就是使奇异值分解，因为 HLSL 没有自带的奇异值分解，我们就翻译宾夕法尼亚大学开源库 ziran 中的用于实现 MPM 的 SVD 分解，地址在

顺便说一下 SVD 是什么，在二维里面一个形变矩阵可以分解为

U・Sigma・V^t，就是先是旋转，缩放，然后再旋转，就得到最终的变换。而上面的 $R_E$ 就是旋转矩阵部分，就是 mul（U，transpos（V））。那么为什么要分解形变矩阵 F 呢，是因为我们要对 F 旋转部分给剔除掉，只要他的施加弹性力形变的大小

说回函数，Corotated 模型里面定义了 $\mu和\lambda$ 两个参数

然后我们还要限制一个弹性力的临界压缩，因为我们经过 SVD 分解得到 sigma 这个值，就是弹性力大小，我们对他进行限制，然后把力分给塑形力，那么为什么这么做的呢？想象一下你手抓一个雪球用力挤压他，这时候就是对他施加弹性力，现在他的大型改变了，但他还是一个球，但超过一个临界值，他爆开了，那么多出这一部分就是他塑形力，所以我们

$F_E$ 限制在 $[1 - \theta_c,1+\theta_s]$ ，然后把这部分的值补充回塑形形变的 $J$ 上，的在论文中推荐了几个值

这时候一个 grid 总弹性势能可以写成

对于更新形变矩阵 F，我们还是用回上一篇文章的 F 更新

${F}_p^{n+1}=(\textbf I+\Delta t \textbf C_p^n)\textbf{F}_p^{n} \\$

那么现在总弹性势能就是 stress force，替换上一篇文章的公式就可以写成

$(m\textbf v)_i^{n}=\sum_pw_{i,p} \left\{ m_p \textbf v_p^n+\left[m_p\textbf C_p^n-\frac{4\Delta t}{\Delta x^2}\sum_p V_p^0\textbf P(\textbf F_p^{n+1})(\textbf F_p^{n+1})^T \right](\textbf x_i-\textbf x_p^n)\right\}$

那么这样就是一个雪模拟公式，后面的还是跟上一篇文章一样，在 P2G 上修改一下。

总结一下上面流程

1、第一步首先更新形变矩阵 F
2、然后对他进行 SVD 分解
3、约束 SVD 解出来的 sigma
4、重建原来 F 矩阵

# Niagara 实现

雪 P2G 部分
``

P2G

Functions Fu;
int Cellx;
int Celly;
int Cellz;
RasterizationGrid3D.GetNumCells(Cellx, Celly, Cellz);
float dx = 1/ (float)Cellx;//((float)Cellx / WorldSize.x);
float DeltaDx = dx*dx;

float3 CellPosition =  Position * Cellx ;

float3 BasePosition = floor(CellPosition - 0.5);
float3 FPosition = CellPosition - BasePosition;

float3 Weight[3];
float IGNORE ;
Weight[0] = 0.5 * (1.5 - FPosition) * (1.5 - FPosition);
Weight[1] = 0.75 - (FPosition - 1) * (FPosition - 1);
Weight[2] = 0.5 * (FPosition - 0.5) * (FPosition - 0.5);

float pvol = pow(dx*0.5 ,3);
float4x4 Identity = float4x4(1,0,0,0,
0,1,0,0,
0,0,1,0,
0,0,0,1);

F = mul(Identity + DeltaTime * C,F);

float h = exp(10 * (1-J));
float mu = (E / (2 * (1 + nu))) * h;
float la = (E * nu / ((1 + nu) * (1 - 2 * nu))) * h;
TestVector = float3(mu,la,0) / h;
float3x3 u = 0.f;
float3 s = 0.f;
float3x3 v = 0.f;
float3x3 TempF;
Fu.ConverMatrix(F,TempF);
Fu.svd(TempF,u,s,v);

float JE = 1.f;

for(int i = 0;i<3;i++)
{
	s[i] = max(s[i],1e-6);
	float MaxValue = 1.0075;//1 + 4.5e-3;//1.0075;//
	float MinValue =0.975;//1 - 2.5e-2;//0.975;//
	float NewSig = s[i];
		NewSig= clamp(s[i],MinValue,MaxValue);
	J *= (s[i] / NewSig);
	s[i] = NewSig;
	JE *= s[i];
}

float4x4 NewU;
Fu.ConverMatrix(u,NewU);

float4x4 NewV;
Fu.ConverMatrix(v,NewV);
TestM = NewU;
float4x4 SMatrix = float4x4(s.x,0,0,0,
	0,s.y,0,0,
	0,0,s.z,0,
	0,0,0,0);
F = mul(mul(NewU,SMatrix),transpose(NewV));
OutF = F;
OutJ = J;

float4x4 Stress = 2 * mu * mul(transpose(F), F -mul( NewU,transpose(NewV)) ) + (Identity * la * JE * (JE - 1));
Stress = (-DeltaTime * pvol * 4 /DeltaDx) * Stress;
float4x4 Affine =  Stress + Mass * C;

for(int x = 0;x<3;x++)
{
    for(int y = 0; y < 3; y++)
    {
        for(int z = 0; z<3; z++)
        {
            float3 offset = float3(x,y,z);
            float3 Dpos = (offset - FPosition) * dx;
            int IndexX = BasePosition.x + offset.x;
            int IndexY = BasePosition.y + offset.y;
            int IndexZ = BasePosition.z + offset.z;
            if(IndexX < Cellx&&IndexY < Celly&&IndexZ < Cellz&&
                IndexX > 0&&IndexY > 0&&IndexZ > 0)
            {
                float CellWeight = Weight[x].x*Weight[y].y*Weight[z].z;
                float4 GridData = 0;
               // Grid3D.GetPreviousVector4ValueAtIndex(IndexX, IndexY, IndexZ, 0, GridData);
                float4 GridVelocity ;//= GridData.xyz ;
                float4 Mpos = float4(Dpos,0);
                GridVelocity = CellWeight*(Mass * float4(Velocity,0) + mul(Affine,Mpos));
               
                float GridMass ;//= GridData.w;
                GridMass =  CellWeight * Mass;
                RasterizationGrid3D.InterlockedAddFloatGridValue(IndexX, IndexY, IndexZ, 0, GridVelocity.x,IGNORE);
                RasterizationGrid3D.InterlockedAddFloatGridValue(IndexX, IndexY, IndexZ, 1, GridVelocity.y,IGNORE);
                RasterizationGrid3D.InterlockedAddFloatGridValue(IndexX, IndexY, IndexZ, 2, GridVelocity.z,IGNORE);
                RasterizationGrid3D.InterlockedAddFloatGridValue(IndexX, IndexY, IndexZ, 3, GridMass,IGNORE);
                
                //Grid3D.SetVector4ValueAtIndex(IndexX, IndexY, IndexZ, 0, float4(GridVelocity,GridMass));

            }
            

        }
    }
}

# 弹性体

因为上面 Corotated 模型本来是个弹塑性物体的模拟模型，我们把雪关于弹性

$F_E$ 和塑形 $F_P$

的一部分删除了就可以实现 softbody 模拟，顺便说一下这对于 FEM 来说最大的优势就在自碰撞，我文章写过一个 FEM，但其中可以看出，他没有自碰撞，只是单纯求 F。最后 MPM 的 SoftBody 代码是这样的，这个 Softbody 还可以断裂的

G2P

Functions Fu;
int Cellx;
int Celly;
int Cellz;
RasterizationGrid3D.GetNumCells(Cellx, Celly, Cellz);
float dx = 1/ (float)Cellx;//((float)Cellx / WorldSize.x);
float DeltaDx = dx*dx;

float3 CellPosition =  Position * Cellx ;

float3 BasePosition = floor(CellPosition - 0.5);
float3 FPosition = CellPosition - BasePosition;

float3 Weight[3];
float IGNORE ;
Weight[0] = 0.5 * (1.5 - FPosition) * (1.5 - FPosition);
Weight[1] = 0.75 - (FPosition - 1) * (FPosition - 1);
Weight[2] = 0.5 * (FPosition - 0.5) * (FPosition - 0.5);

float pvol = pow(dx*0.5 ,3);
float4x4 Identity = float4x4(1,0,0,0,
0,1,0,0,
0,0,1,0,
0,0,0,1);

F = mul(Identity + DeltaTime * C,F);

float h = 0.3;//exp(10 * (1-J));
float mu = (E / (2 * (1 + nu))) * h;
float la = (E * nu / ((1 + nu) * (1 - 2 * nu))) * h;
TestVector = float3(mu,la,0) / h;
float3x3 u = 0.f;
float3 s = 0.f;
float3x3 v = 0.f;
float3x3 TempF;
Fu.ConverMatrix(F,TempF);
Fu.svd(TempF,u,s,v);

float JT = 1.f;

float4x4 NewU;
Fu.ConverMatrix(u,NewU);

float4x4 NewV;
Fu.ConverMatrix(v,NewV);
TestM = NewU;

OutF = F;
OutJ = J;

float4x4 Stress = 2 * mu * mul(transpose(F), F -mul( NewU,transpose(NewV)) ) + (Identity * la * JT * (JT - 1));
Stress = (-DeltaTime * pvol * 4 /DeltaDx) * Stress;
float4x4 Affine =  Stress + Mass * C;

for(int x = 0;x<3;x++)
{
    for(int y = 0; y < 3; y++)
    {
        for(int z = 0; z<3; z++)
        {
            float3 offset = float3(x,y,z);
            float3 Dpos = (offset - FPosition) * dx;
            int IndexX = BasePosition.x + offset.x;
            int IndexY = BasePosition.y + offset.y;
            int IndexZ = BasePosition.z + offset.z;
            if(IndexX < Cellx&&IndexY < Celly&&IndexZ < Cellz&&
                IndexX > 0&&IndexY > 0&&IndexZ > 0)
            {
                float CellWeight = Weight[x].x*Weight[y].y*Weight[z].z;
                float4 GridData = 0;
               // Grid3D.GetPreviousVector4ValueAtIndex(IndexX, IndexY, IndexZ, 0, GridData);
                float4 GridVelocity ;//= GridData.xyz ;
                float4 Mpos = float4(Dpos,0);
                GridVelocity = CellWeight*(Mass * float4(Velocity,0) + mul(Affine,Mpos));
               
                float GridMass ;//= GridData.w;
                GridMass =  CellWeight * Mass;
                RasterizationGrid3D.InterlockedAddFloatGridValue(IndexX, IndexY, IndexZ, 0, GridVelocity.x,IGNORE);
                RasterizationGrid3D.InterlockedAddFloatGridValue(IndexX, IndexY, IndexZ, 1, GridVelocity.y,IGNORE);
                RasterizationGrid3D.InterlockedAddFloatGridValue(IndexX, IndexY, IndexZ, 2, GridVelocity.z,IGNORE);
                RasterizationGrid3D.InterlockedAddFloatGridValue(IndexX, IndexY, IndexZ, 3, GridMass,IGNORE);
                
                //Grid3D.SetVector4ValueAtIndex(IndexX, IndexY, IndexZ, 0, float4(GridVelocity,GridMass));

            }
            

        }
    }
}

# SVD 代码

SVD

struct Functions
{
#define EPSILON 1e-10

	// Function to create a 2x2 rotation matrix
	float2x2 G2(float c, float s)
	{
		return float2x2(c, s, -s, c);
	}

	// Function to compute Givens rotation for constrained case
	float2 G2_Con(float a, float b)
	{
		float d = a * a + b * b;
		float c = 1;
		float s = 0;

		if (abs(d) > 0)
		{
			float t = 1.0 / sqrt(d);
			c = a * t;
			s = -b * t;
		}

		return float2(c, s);
	}

	// Function to compute Givens rotation for unconstrained case
	float2 G2_unCon(float a, float b)
	{
		float d = a * a + b * b;
		float c = 1;
		float s = 0;

		if (abs(d) > 0)
		{
			float t = 1.0 / sqrt(d);
			c = a * t;
			s = b * t;
		}

		return float2(c, s);
	}

	// Function to create a 3x3 matrix for Givens rotation in the 1-2 plane
	float3x3 G3_12(float c, float s, bool con = true)
	{
		float2 cs = con ? G2_Con(c, s) : G2_unCon(c, s);
		return float3x3(float3(cs.x, cs.y, 0), float3(-cs.y, cs.x, 0), float3(0, 0, 1));
	}

	float3x3 G3_12_Direct(float c, float s)
	{
		return float3x3(float3(c, s, 0), float3(-s, c, 0), float3(0, 0, 1));
	}

	// Function to create a 3x3 matrix for Givens rotation in the 2-3 plane
	float3x3 G3_23(float c, float s, bool con = true)
	{
		float2 cs = con ? G2_Con(c, s) : G2_unCon(c, s);
		return float3x3(float3(1, 0, 0), float3(0, cs.x, cs.y), float3(0, -cs.y, cs.x));
	}

	float3x3 G3_23_Direct(float c, float s)
	{
		return float3x3(float3(1, 0, 0), float3(0, c, s), float3(0, -s, c));
	}

	// Function to create a 3x3 matrix for Givens rotation in the 1-3 plane
	float3x3 G3_13(float c, float s, bool con = true)
	{
		float2 cs = con ? G2_Con(c, s) : G2_unCon(c, s);
		return float3x3(float3(cs.x, 0, cs.y), float3(0, 1, 0), float3(-cs.y, 0, cs.x));
	}

	float3x3 G3_13_Direct(float c, float s)
	{
		return float3x3(float3(c, 0, s), float3(0, 1, 0), float3(-s, 0, c));
	}

	// Function for polar decomposition
	void PolarDecompostion(float2x2 A, out float2x2 R, out float2x2 S)
	{
		float x = A[0].x + A[1].y;
		float y = A[1].x - A[0].y;
		float d = sqrt(x * x + y * y);

		R = float2x2(1, 0, 0, 1);

		if (abs(d) > EPSILON)
		{
			d = 1 / d;
			R = G2(x * d, -y * d);
		}

		S = mul(transpose(R), A);
	}

	// Function for singular value decomposition of a 2x2 matrix
	void SVD2x2(float2x2 A, out float2x2 U, out float2 D, out float2x2 V)
	{
		U = float2x2(0, 0, 0, 0);
		D = float2(0, 0);
		V = float2x2(0, 0, 0, 0);

		float2x2 R, S;
		PolarDecompostion(A, R, S);

		float c = 1;
		float s = 0;

		if (abs(S[0].y) < EPSILON)
		{
			D.x = S[0][0];
			D.y = S[1][1];
		}
		else
		{
			float taw = 0.5 * (S[0][0] - S[1][1]);
			float w = sqrt(taw * taw + S[0][1] * S[0][1]);

			float t = 0;

			if (taw > 0)
				t = S[0].y / (taw + w);
			else
				t = S[0].y / (taw - w);

			c = 1.0 / sqrt(t * t + 1);
			s = -t * c;

			D.x = (c * c * S[0][0]) - (2 * c * s * S[0][1]) + (s * s * S[1][1]);
			D.y = (s * s * S[0][0]) + (2 * c * s * S[0][1]) + (c * c * S[1][1]);
		}

		if (D.x < D.y)
		{
			float temp = D.x;
			D.x = D.y;
			D.y = temp;
			V = G2(-s, c);
		}
		else
		{
			V = G2(c, s);
		}

		U = mul(R, V);
	}

	// Other functions and main SVD3x3 function go here...
	// Function to perform zero chasing
	void ZeroChasing(inout float3x3 U, inout float3x3 A, inout float3x3 V)
	{
		float3x3 G = G3_12(A[0][0], A[1][0]);
		A = mul(transpose(G), A);
		U = mul(U, G);

		float c = A[0][1];
		float s = A[0][2];

		if (abs(A[1][0]) > EPSILON)
		{
			c = A[0].x * A[0][1] + A[1][0] * A[1][1];
			s = A[0].x * A[0][2] + A[1][0] * A[1][2];
		}

		G = G3_23(c, s);
		A = mul(A, G);
		V = mul(V, G);

		G = G3_23(A[1][1], A[2][1]);
		A = mul(transpose(G), A);
		U = mul(U, G);
	}

	// Function to perform bidiagonalization
	void Bidiag(inout float3x3 U, inout float3x3 A, inout float3x3 V)
	{
		float3x3 G = G3_23(A[1][0], A[2][0]);
		A = mul(transpose(G), A);
		U = mul(U, G);
		ZeroChasing(U, A, V);
	}

	// Function to compute Frobenius norm of a 3x3 matrix
	float FrobeniusNorm(float3x3 B)
	{
		float ret = 0;
		for(int i =0;i<3;i++ )
		{
			for(int j =0;j<3;j++ )
			{
				ret += B[i][j]*B[i][j];
			}
		}
		return ret;//dot(B[0], B[0]) + dot(B[1], B[1]) + dot(B[2], B[2]);
	}

	// Function to flip signs in a matrix and sigma vector
	void FlipSign(int idx, inout float3x3 mat, inout float3 sigma)
	{
		mat[0][idx] = -mat[0][idx];
		mat[1][idx] = -mat[1][idx];
		mat[2][idx] = -mat[2][idx];
		sigma[idx] = -sigma[idx];
	}

	// Function to flip signs in a specific column of a matrix
	void FlipSignColumn(inout float3x3 mat, int idx)
	{
		mat[0][idx] = -mat[0][idx];
		mat[1][idx] = -mat[1][idx];
		mat[2][idx] = -mat[2][idx];
	}

	// Function to swap columns in a matrix
	void SwapColumn(inout float3x3 a, int col_a, int col_b)
	{
		float3 acopy = a[col_a];
		a[col_a] = a[col_b];
		a[col_b] = acopy;
	}

	// Function to sort with top-left singular value
	void SortWithTopLeft(inout float3x3 U, inout float3 sigma, inout float3x3 V)
	{
		if (abs(sigma[1]) >= abs(sigma[2]))
		{
			if (sigma[1] < 0)
			{
				FlipSign(1, U, sigma);
				FlipSign(2, U, sigma);
			}

			return;
		}

		if (sigma[2] < 0)
		{
			FlipSign(1, U, sigma);
			FlipSign(2, U, sigma);
		}

		float temp = sigma[1];
		sigma[1] = sigma[2];
		sigma[2] = temp;

		SwapColumn(U, 1, 2);
		SwapColumn(V, 1, 2);

		if (sigma[1] > sigma[0])
		{
			temp = sigma[0];
			sigma[0] = sigma[1];
			sigma[1] = temp;

			SwapColumn(U, 0, 1);
			SwapColumn(V, 0, 1);
		}
		else
		{
			FlipSignColumn(U, 2);
			FlipSignColumn(V, 2);
		}
	}

	// Function to sort with bottom-right singular value
	void SortWithBotRight(inout float3x3 U, inout float3 sigma, inout float3x3 V)
	{
		if (abs(sigma[0]) >= abs(sigma[1]))
		{
			if (sigma[0] < 0)
			{
				FlipSign(0, U, sigma);
				FlipSign(2, U, sigma);
			}

			return;
		}

		float temp = sigma[0];
		sigma[0] = sigma[1];
		sigma[1] = temp;

		SwapColumn(U, 0, 1);
		SwapColumn(V, 0, 1);

		if (abs(sigma[1]) < abs(sigma[2]))
		{
			temp = sigma[1];
			sigma[1] = sigma[2];
			sigma[2] = temp;

			SwapColumn(U, 1, 2);
			SwapColumn(V, 1, 2);
		}
		else
		{
			FlipSignColumn(U, 2);
			FlipSignColumn(V, 2);
		}

		if (sigma[1] < 0)
		{
			FlipSign(1, U, sigma);
			FlipSign(2, U, sigma);
		}
	}

	// Function to solve for reduced top-left matrix
	void SolveReducedTopLeft(inout float3x3 B, inout float3x3 U, inout float3x3 V, inout float3 sigma)
	{
		float s3 = B[2][2];
		float2x2 top_left = float2x2(B[0].xy, B[1].xy);
		float2x2 U2, V2;
		float2 D2;
		SVD2x2(top_left, U2, D2, V2);

		float3x3 u3 = G3_12_Direct(U2[0].x, U2[0].y);
		float3x3 v3 = G3_12_Direct(V2[0].x, V2[0].y);

		U = mul(U, u3);
		V = mul(V, v3);

		sigma = float3(D2[0], D2[1], s3);
	}

	// Function to solve for reduced bottom-right matrix
	void SolveReducedBotRight(inout float3x3 B, inout float3x3 U, inout float3x3 V, inout float3 sigma)
	{
		float s1 = B[0][0];
		float2x2 botRight = float2x2(B[1].yz, B[2].yz);
		float2x2 U2, V2;
		float2 D2;
		SVD2x2(botRight, U2, D2, V2);

		float3x3 u3 = G3_23_Direct(U2[0].x, U2[0].y);
		float3x3 v3 = G3_23_Direct(V2[0].x, V2[0].y);

		U = mul(U, u3);
		V = mul(V, v3);

		sigma = float3(s1, D2[0], D2[1]);
	}
	// Function for post-processing
	void PostProcess(float3x3 B, inout float3x3 U,out float3 sigma, inout float3x3 V, float alpha1, float alpha2, float alpha3, float beta1, float beta2, float gamma1, float gamma2, float tao)
	{

		if (abs(beta2) <= tao)
		{
			SolveReducedTopLeft(B, U, V, sigma);
			SortWithTopLeft(U, sigma, V);
		}
		else if (abs(beta1) <= tao)
		{
			SolveReducedBotRight(B, U, V, sigma);
			SortWithBotRight(U, sigma, V);
		}
		else if (abs(alpha2) <= tao)
		{
			float3x3 G = G3_23(B[1].z, B[2].z, false);
			B = mul(transpose(G), B);
			U = mul(U, G);
			SolveReducedTopLeft(B, U, V, sigma);
			SortWithTopLeft(U, sigma, V);
		}
		else if (abs(alpha3) <= tao)
		{
			float3x3 G = G3_23(B[1].y, B[1].z);
			B = mul(B, G);
			V = mul(V, G);
			G = G3_13(B[0].x, B[0].z);
			B = mul(B, G);
			V = mul(V, G);
			SolveReducedTopLeft(B, U, V, sigma);
			SortWithTopLeft(U, sigma, V);
		}
		else if (abs(alpha1) <= tao)
		{
			float3x3 G = G3_12(B[0].y, B[1].y, false);
			B = mul(transpose(G), B);
			U = mul(U, G);
			G = G3_13(B[0].z, B[2].z, false);
			B = mul(transpose(G), B);
			U = mul(U, G);
			SolveReducedBotRight(B, U, V, sigma);
			SortWithBotRight(U, sigma, V);
		}
	}

	// Function for SVD of a 3x3 matrix
	void svd(float3x3 A, out float3x3 U, out float3 Sigma, out float3x3 V)
	{
		U = float3x3(1, 0, 0, 0, 1, 0, 0, 0, 1);
		V = float3x3(1, 0, 0, 0, 1, 0, 0, 0, 1);

		Bidiag(U, A, V);
		float3x3 B = A;

		float alpha1 = B[0].x;
		float alpha2 = B[1].y;
		float alpha3 = B[2].z;
		float beta1 = B[0].y;
		float beta2 = B[1].z;
		float gamma1 = alpha1 * beta1;
		float gamma2 = alpha2 * beta2;

		float tol = 128 * EPSILON;
		float tao = tol * max(0.5 * FrobeniusNorm(B), 1);
		int Num = 0;
		while (abs(beta1) > tao && abs(beta2) > tao &&
			   abs(alpha1) > tao && abs(alpha2) > tao && abs(alpha3) > tao)
		{
            if(Num>50)
            {
                break;
            }
			float a1 = alpha2 * alpha2 + beta1 * beta1;
			float a2 = alpha3 * alpha3 + beta1 * beta1;
			float b1 = gamma1;
			float d = (a1 - a2) * 0.5;
			float mu = (b1 * b1) / (abs(d) + sqrt(d * d + b1 * b1));

			if (d < 0)
				mu = -mu;

			mu = a2 - mu;
			float3x3 G = G3_12(alpha1 * alpha1 - mu, gamma1);
			B = mul(B, G);
			V = mul(V, G);
			ZeroChasing(U, B, V);

			alpha1 = B[0].x;
			alpha2 = B[1].y;
			alpha3 = B[2].z;
			beta1 = B[0].y;
			beta2 = B[1].z;
			gamma1 = alpha1 * beta1;
			gamma2 = alpha2 * beta2;
            Num++;

		}

		PostProcess(B, U,Sigma, V, alpha1, alpha2, alpha3, beta1, beta2, gamma1, gamma2, tao);

	}
	void ConverMatrix(float3x3 Inmat, out float4x4 OutMat)
	{
		OutMat =  float4x4(Inmat[0].x,Inmat[0].y,Inmat[0].z,0,
			Inmat[1].x,Inmat[1].y,Inmat[1].z,0,
			Inmat[2].x,Inmat[2].y,Inmat[2].z,0,
			0,0,0,1);

	}
	void ConverMatrix(float4x4 Inmat,out float3x3 OutMat)
	{
		OutMat =  float3x3(Inmat[0].x,Inmat[0].y,Inmat[0].z,
			Inmat[1].x,Inmat[1].y,Inmat[1].z,
			Inmat[2].x,Inmat[2].y,Inmat[2].z);
	}
};

模拟文章

# MPM 雪、海绵

# 简介

# 雪

# Niagara 实现

# 弹性体

# SVD 代码

MPM水

奥高公式