简介
有限元方法, 把物体分割成一个个有体积的单元来模拟.
线性有限元方法在二维空间中把物体分割成三角形/四边形, 在三维空间中把物体分割成四面体/六面体.
有限元方法由能量对位置求导得到力,有限体积方法对 traction vector 积分得到力。
TIPS
STVK模型, PPT中写错了, 正确的在 https://www.bilibili.com/read/cv14754926 上有讲解
\[W=\frac{s_{0}}{8}(I-3)^{2}+\frac{s_{1}}{4}(I I-2 I+3)
\]
核心参考PPT
SVTK实现
FVM实现, 说实话 我也不知道 下面实现的是有限元还是有限体积, 逻辑上应该是有限体积
关于Laplacian 速度平滑
就是V_SUM 周围顶点速度的和
V_NUM 周围顶点的个数
附录: 师弟代码
code
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
using System;
using System.IO;
public class FVM : MonoBehaviour
{
float dt = 0.003f;
float mass = 1;
float stiffness_0 = 7500.0f;
float stiffness_1 = 5000.0f;
float damp = 0.999f;
float Un = 0.5f; // for collision
float Ut = 0.5f; // for collision
Vector3 P = new Vector3(0, -3, 0);
Vector3 N = new Vector3(0, 1, 0);
int[] Tet; //四面体数据
int tet_number; //The number of tetrahedra 四面体数目
Vector3[] Force; //力
Vector3[] V; //速度
Vector3[] X; //点数据
int number; //The number of vertices
Matrix4x4[] inv_Dm;
//For Laplacian smoothing.
Vector3[] V_sum;
int[] V_num;
SVD svd = new SVD();
// Start is called before the first frame update
void Start()
{
// FILO IO: Read the house model from files.
// The model is from Jonathan Schewchuk's Stellar lib.
{
string fileContent = File.ReadAllText("Assets/house2.ele");
string[] Strings = fileContent.Split(new char[] { ' ', '\t', '\r', '\n' }, StringSplitOptions.RemoveEmptyEntries);
tet_number = int.Parse(Strings[0]);
Tet = new int[tet_number * 4];
for (int tet = 0; tet < tet_number; tet++)
{
Tet[tet * 4 + 0] = int.Parse(Strings[tet * 5 + 4]) - 1;
Tet[tet * 4 + 1] = int.Parse(Strings[tet * 5 + 5]) - 1;
Tet[tet * 4 + 2] = int.Parse(Strings[tet * 5 + 6]) - 1;
Tet[tet * 4 + 3] = int.Parse(Strings[tet * 5 + 7]) - 1;
}
}
{
string fileContent = File.ReadAllText("Assets/house2.node");
string[] Strings = fileContent.Split(new char[] { ' ', '\t', '\r', '\n' }, StringSplitOptions.RemoveEmptyEntries);
number = int.Parse(Strings[0]);
X = new Vector3[number];
for (int i = 0; i < number; i++)
{
X[i].x = float.Parse(Strings[i * 5 + 5]) * 0.4f;
X[i].y = float.Parse(Strings[i * 5 + 6]) * 0.4f;
X[i].z = float.Parse(Strings[i * 5 + 7]) * 0.4f;
}
//Centralize the model.
Vector3 center = Vector3.zero;
for (int i = 0; i < number; i++) center += X[i];
center = center / number;
for (int i = 0; i < number; i++)
{
X[i] -= center;
float temp = X[i].y;
X[i].y = X[i].z;
X[i].z = temp;
}
}
/*tet_number=1;
Tet = new int[tet_number*4];
Tet[0]=0;
Tet[1]=1;
Tet[2]=2;
Tet[3]=3;
number=4;
X = new Vector3[number];
V = new Vector3[number];
Force = new Vector3[number];
X[0]= new Vector3(0, 0, 0);
X[1]= new Vector3(1, 0, 0);
X[2]= new Vector3(0, 1, 0);
X[3]= new Vector3(0, 0, 1);*/
//Create triangle mesh.
Vector3[] vertices = new Vector3[tet_number * 12];
int vertex_number = 0;
for (int tet = 0; tet < tet_number; tet++)
{
vertices[vertex_number++] = X[Tet[tet * 4 + 0]];
vertices[vertex_number++] = X[Tet[tet * 4 + 2]];
vertices[vertex_number++] = X[Tet[tet * 4 + 1]];
vertices[vertex_number++] = X[Tet[tet * 4 + 0]];
vertices[vertex_number++] = X[Tet[tet * 4 + 3]];
vertices[vertex_number++] = X[Tet[tet * 4 + 2]];
vertices[vertex_number++] = X[Tet[tet * 4 + 0]];
vertices[vertex_number++] = X[Tet[tet * 4 + 1]];
vertices[vertex_number++] = X[Tet[tet * 4 + 3]];
vertices[vertex_number++] = X[Tet[tet * 4 + 1]];
vertices[vertex_number++] = X[Tet[tet * 4 + 2]];
vertices[vertex_number++] = X[Tet[tet * 4 + 3]];
}
int[] triangles = new int[tet_number * 12];
for (int t = 0; t < tet_number * 4; t++)
{
triangles[t * 3 + 0] = t * 3 + 0;
triangles[t * 3 + 1] = t * 3 + 1;
triangles[t * 3 + 2] = t * 3 + 2;
}
Mesh mesh = GetComponent<MeshFilter>().mesh;
mesh.vertices = vertices;
mesh.triangles = triangles;
mesh.RecalculateNormals();
V = new Vector3[number];
Force = new Vector3[number];
V_sum = new Vector3[number];
V_num = new int[number];
for (int tet = 0; tet < tet_number; tet++)
{
V_num[Tet[tet * 4]] += 4;
V_num[Tet[tet * 4 + 1]] += 4;
V_num[Tet[tet * 4 + 2]] += 4;
V_num[Tet[tet * 4 + 3]] += 4;
}
//TODO: Need to allocate and assign inv_Dm
inv_Dm = new Matrix4x4[tet_number];
for (int tet = 0; tet < tet_number; tet++)
{
inv_Dm[tet] = Build_Edge_Matrix(tet).inverse;
//Debug.Log("inv_Dm:" + inv_Dm[tet]);
}
}
Matrix4x4 Build_Edge_Matrix(int tet)
{
Matrix4x4 ret = Matrix4x4.zero;
//TODO: Need to build edge matrix here.
ret[0, 0] = X[Tet[tet * 4]][0] - X[Tet[tet * 4 + 1]][0];
ret[1, 0] = X[Tet[tet * 4]][1] - X[Tet[tet * 4 + 1]][1];
ret[2, 0] = X[Tet[tet * 4]][2] - X[Tet[tet * 4 + 1]][2];
ret[0, 1] = X[Tet[tet * 4]][0] - X[Tet[tet * 4 + 2]][0];
ret[1, 1] = X[Tet[tet * 4]][1] - X[Tet[tet * 4 + 2]][1];
ret[2, 1] = X[Tet[tet * 4]][2] - X[Tet[tet * 4 + 2]][2];
ret[0, 2] = X[Tet[tet * 4]][0] - X[Tet[tet * 4 + 3]][0];
ret[1, 2] = X[Tet[tet * 4]][1] - X[Tet[tet * 4 + 3]][1];
ret[2, 2] = X[Tet[tet * 4]][2] - X[Tet[tet * 4 + 3]][2];
ret[3, 3] = 1;
return ret;
}
Matrix4x4 M_Division(Matrix4x4 m, float scale)
{
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
m[i, j] /= scale;
}
}
return m;
}
Matrix4x4 M_Multipy(Matrix4x4 m, float scale)
{
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
m[i, j] *= scale;
}
}
return m;
}
Matrix4x4 M_addtion(Matrix4x4 m1, Matrix4x4 m2)
{
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
m1[i, j] += m2[i, j];
}
}
return m1;
}
Matrix4x4 M_Dec(Matrix4x4 m1, Matrix4x4 m2)
{
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
m1[i, j] -= m2[i, j];
}
}
return m1;
}
float TR(Matrix4x4 m)
{
float res = 0f;
res += m.m00;
res += m.m11;
res += m.m22;
return res;
}
void Update_Collision_Point(int idx, float distance, Vector3 N)
{
X[idx] -= distance * N;
float viN = Vector3.Dot(V[idx], N);
if (viN >= 0) return;
Vector3 v_Ni = viN * N; // (vi * N)N 反弹面的分量
Vector3 v_Ti = V[idx] - v_Ni; // 沿着面的分量
//定义摩擦系数a
float a = Mathf.Max(1.0f - Ut * (1.0f + Un) * v_Ni.magnitude / v_Ti.magnitude, 0.0f);
//计算新速度
Vector3 v_Ni_new = -Un * v_Ni;
Vector3 v_Ti_new = a * v_Ti;
V[idx] = v_Ni_new + v_Ti_new;
}
Matrix4x4 SVtk(Matrix4x4 F)
{
Matrix4x4 P = Matrix4x4.zero;
Matrix4x4 U = Matrix4x4.zero;
Matrix4x4 S = Matrix4x4.zero;
Matrix4x4 V = Matrix4x4.zero;
Matrix4x4 Plamda = Matrix4x4.zero;
svd.svd(F, ref U, ref S, ref V);
float lamda0 = S[0, 0];
float lamda1 = S[1, 1];
float lamda2 = S[2, 2];
float Ic = lamda0 * lamda0 + lamda1 * lamda1 +
lamda2 * lamda2;
float Ic_l0 = 2f * lamda0;
float Ic_l1 = 2f * lamda1;
float Ic_l2 = 2f * lamda2;
//float IIc_l0 = 2f * lamda0 * (lamda1 * lamda1 + lamda2 * lamda2);
//float IIc_l1 = 2f * lamda1 * (lamda0 * lamda0 + lamda2 * lamda2);
//float IIc_l2 = 2f * lamda2 * (lamda1 * lamda1 + lamda0 * lamda0);
float IIc_l0 = 4f * lamda0 * lamda0 * lamda0;
float IIc_l1 = 4f * lamda1 * lamda1 * lamda1;
float IIc_l2 = 4f * lamda2 * lamda2 * lamda2;
//float IIc_l0 = 4f * lamda0 * lamda0 * lamda0 * lamda1 * lamda1 * lamda1 * lamda1 * lamda2 * lamda2 * lamda2 * lamda2;
//float IIc_l1 = 4f * lamda1 * lamda1 * lamda1 * lamda0 * lamda0 * lamda0 * lamda0 * lamda2 * lamda2 * lamda2 * lamda2;
//float IIc_l2 = 4f * lamda2 * lamda2 * lamda2 * lamda1 * lamda1 * lamda1 * lamda1 * lamda0 * lamda0 * lamda0 * lamda0;
Plamda[0, 0] = stiffness_0 * (Ic - 3f) * Ic_l0 /4f +
stiffness_1 * (IIc_l0 - 2 * Ic_l0) / 4f;
Plamda[1, 1] = stiffness_0 * (Ic - 3f) * Ic_l1 /4f +
stiffness_1 * (IIc_l1 - 2 * Ic_l1) / 4f;
Plamda[2, 2] = stiffness_0 * (Ic - 3f) * Ic_l2 /4f +
stiffness_1 * (IIc_l2 - 2 * Ic_l2) / 4f;
Plamda[3, 3] = 1f;
P = U * Plamda * V.transpose;
return P;
}
void _Update()
{
// Jump up.
if (Input.GetKeyDown(KeyCode.Space))
{
for (int i = 0; i < number; i++)
V[i].y += 0.2f;
}
for (int i = 0; i < number; i++)
{
//TODO: Add gravity to Force.
Force[i] = new Vector3(0, -9.8f, 0);
}
for (int tet = 0; tet < tet_number; tet++)
{
Matrix4x4 I = Matrix4x4.identity;
//TODO: Deformation Gradient
Matrix4x4 F = Build_Edge_Matrix(tet);
F = F * inv_Dm[tet];
//TODO: Green Strain
Matrix4x4 G = M_Division(M_Dec(F.transpose * F, I)
, 2f);
//TODO: Second PK Stress
Matrix4x4 S = M_addtion(M_Multipy(G,
stiffness_1 * 2f), M_Multipy(I, TR(G) *
stiffness_0));
Matrix4x4 P = F * S;
//SVTK
P = SVtk(F);
//TODO: Elastic Force
Matrix4x4 EF = M_Multipy(P * (inv_Dm[tet].transpose)
, -1f / (6f * inv_Dm[tet].determinant));
//Debug.Log("EF:" + EF);
Force[Tet[tet * 4 + 1]][0] += EF[0, 0];
Force[Tet[tet * 4 + 1]][1] += EF[1, 0];
Force[Tet[tet * 4 + 1]][2] += EF[2, 0];
Force[Tet[tet * 4 + 2]][0] += EF[0, 1];
Force[Tet[tet * 4 + 2]][1] += EF[1, 1];
Force[Tet[tet * 4 + 2]][2] += EF[2, 1];
Force[Tet[tet * 4 + 3]][0] += EF[0, 2];
Force[Tet[tet * 4 + 3]][1] += EF[1, 2];
Force[Tet[tet * 4 + 3]][2] += EF[2, 2];
Force[Tet[tet * 4]][0] += -EF[0, 0] - EF[0, 1] - EF[0, 2];
Force[Tet[tet * 4]][1] += -EF[1, 0] - EF[1, 1] - EF[1, 2];
Force[Tet[tet * 4]][2] += -EF[2, 0] - EF[2, 1] - EF[2, 2];
}
for (int i = 0; i < number; i++)
{
//TODO: Update X and V here.
V[i] += Force[i] * dt;
V[i] *= damp;
X[i] += V[i] * dt;
//TODO: (Particle) collision with floor.
float distance = Vector3.Dot(X[i] - P, N);
if (distance < 0)
{
Update_Collision_Point(i, distance, N);
}
V_sum[i] = Vector3.zero;
}
//laplacian
for (int tet = 0; tet < tet_number; tet++)
{
Vector3 v0 = V[Tet[tet * 4]];
Vector3 v1 = V[Tet[tet * 4 + 1]];
Vector3 v2 = V[Tet[tet * 4 + 2]];
Vector3 v3 = V[Tet[tet * 4 + 3]];
Vector3 vs = v0 + v1 + v2 + v3;
V_sum[Tet[tet * 4]] += vs;
V_sum[Tet[tet * 4 + 1]] += vs;
V_sum[Tet[tet * 4 + 2]] += vs;
V_sum[Tet[tet * 4 + 3]] += vs;
}
for (int i = 0; i < number; i++)
{
V[i] = V_sum[i] / V_num[i];
}
Un = Mathf.Max(0.0f, Un - 0.005f);
}
// Update is called once per frame
void Update()
{
for (int l = 0; l < 10; l++)
_Update();
// Dump the vertex array for rendering.
Vector3[] vertices = new Vector3[tet_number * 12];
int vertex_number = 0;
for (int tet = 0; tet < tet_number; tet++)
{
vertices[vertex_number++] = X[Tet[tet * 4 + 0]];
vertices[vertex_number++] = X[Tet[tet * 4 + 2]];
vertices[vertex_number++] = X[Tet[tet * 4 + 1]];
vertices[vertex_number++] = X[Tet[tet * 4 + 0]];
vertices[vertex_number++] = X[Tet[tet * 4 + 3]];
vertices[vertex_number++] = X[Tet[tet * 4 + 2]];
vertices[vertex_number++] = X[Tet[tet * 4 + 0]];
vertices[vertex_number++] = X[Tet[tet * 4 + 1]];
vertices[vertex_number++] = X[Tet[tet * 4 + 3]];
vertices[vertex_number++] = X[Tet[tet * 4 + 1]];
vertices[vertex_number++] = X[Tet[tet * 4 + 2]];
vertices[vertex_number++] = X[Tet[tet * 4 + 3]];
}
Mesh mesh = GetComponent<MeshFilter>().mesh;
mesh.vertices = vertices;
mesh.RecalculateNormals();
}
}
image
参考链接
https://zhuanlan.zhihu.com/p/447263178
https://github.com/Teafox-Yang/GAMES103_HW_TEAFOX
https://www.bilibili.com/read/cv14754926
---------------------------我的天空里没有太阳,总是黑夜,但并不暗,因为有东西代替了太阳。虽然没有太阳那么明亮,但对我来说已经足够。凭借着这份光,我便能把黑夜当成白天。我从来就没有太阳,所以不怕失去。
--------《白夜行》