0 前言
宣告：此篇部落格僅用於個人學習記錄之用，並非是分享程式碼。Hornor the Code
本來想仔細學一下有限元再進行LAB3的製作，之前也是這麼做的。問題是咱是做遊戲的，不是做cae的，養家餬口之後真沒有時間去研究那些東西。所以後面的像 張量分析 之類的也就沒仔細去看。只是做了矩陣微積分的那一部分。
後面的流體模擬，也是一個深坑。淚目。
這次作業的Mesh可以看作是沒有索引的頂點陣列。
索引陣列降低了記憶體佔用，但是帶來的負面效果就是，沒法做頂點的刪減增加。這個缺點在三角形細分和三角形裁剪的時候尤為恐怖。
感謝 Games103王華民老師的《Intro to Physics-Based Animation!》。
1 拉格朗日模擬
先從能量說起。能量分為動能和勢能。物理學有一大定律，名為能量守恆。我們知道動能和勢能是可以轉換的。動能和勢能的轉換依靠的就是力。
力總是希望快速的降低勢能。想象一個靜止的彈簧，抓著兩個端點狠狠的一拉，不放手。這個時候動能為0，但是你的手已經為彈簧做了功。彈簧中充滿了潛在的能量，勢能。
這個時候彈簧就會自然的產生力，來減少這個勢能。
牛頓定律三條，說了一堆。核心就是第三條，力是哪來的？兩個物體。最後兩個物體的合力是多少，零。但是對於每一個單獨的物體，卻是受到力的。
萬物都儘量想保持原狀，但是總會有外部的物體來給與你互動施加能量，由此萬物開始運動。只不過因為大家碰到了，只不過因為大家想保持原狀。這就是宇宙。宇宙的合力為0。
質點彈簧系統，一根彈簧，兩個質點作用。中間都是空的，彈簧由能量定義。能量由位置定義。
有限元系統，一個三角形，三個質點作用。中間都是空的，面積由能量定義。能量由位置定義。
有限體積法， 一個四面體，四個質點的作用。中間都是空的，體積由能量定義。能量由位置定義。
僅此而已。
2 四面體資料格式
2.1 house2.ele
/*
First line: <# of tetrahedra> <nodes per tetrahedron> <# of attributes>
Remaining lines list of # of tetrahedra:
<tetrahedron #> <node> <node> <node> <node> ... [attributes]
*/
1389 4 0
0 1 75 175 52
1 1 75 52 360
0 1 11 46 174
0 1 11 174 57
https://wias-berlin.de/software/tetgen/fformats.ele.html
2.2 house2.node
/*
First line: <# of points> <dimension (must be 3)> <# of attributes> <# of boundary markers (0 or 1)>
Remaining lines list # of points:
<point #> <x> <y> <z> [attributes] [boundary marker]
*/
400 3 0 1
1 6.5 1 7.5430812757201648 1
2 4.4456243310934527 5.3728065092415962 8.3129834982130344 0
3 4.1326492294372636 8.1243529245273844 4.2828387696230079 0
4 4.75 7.7548808372563904 7.6928224992291403 0
https://wias-berlin.de/software/tetgen/fformats.node.html
3. Impl
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 = 20000.0f;
float stiffness_1 = 5000.0f;
float damp = 0.999f;
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_number = 1;
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];
//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;
}
}
Matrix4x4 Build_Edge_Matrix(int tet)
{
Matrix4x4 ret = Matrix4x4.zero;
Vector3 X_10 = X[Tet[tet * 4 + 1]] - X[Tet[tet * 4 + 0]];
Vector3 X_20 = X[Tet[tet * 4 + 2]] - X[Tet[tet * 4 + 0]];
Vector3 X_30 = X[Tet[tet * 4 + 3]] - X[Tet[tet * 4 + 0]];
//TODO: Need to build edge matrix here.
ret.SetColumn(0, new Vector4(X_10.x, X_10.y, X_10.z, 0.0f));
ret.SetColumn(1, new Vector4(X_20.x, X_20.y, X_20.z, 0.0f));
ret.SetColumn(2, new Vector4(X_30.x, X_30.y, X_30.z, 0.0f));
ret.SetColumn(3, new Vector4(0.0f, 0.0f, 0.0f, 1.0f));
return ret;
}
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) * mass;
}
for (int tet = 0; tet < tet_number; tet++)
{
//TODO: Deformation Gradient
Vector3 x_10 = X[Tet[tet * 4 + 1]] - X[Tet[tet * 4 + 0]];
Vector3 x_20 = X[Tet[tet * 4 + 2]] - X[Tet[tet * 4 + 0]];
Vector3 x_30 = X[Tet[tet * 4 + 3]] - X[Tet[tet * 4 + 0]];
Matrix4x4 Dm = new Matrix4x4();
Dm.SetColumn(0, new Vector4(x_10.x, x_10.y, x_10.z, 0.0f));
Dm.SetColumn(1, new Vector4(x_20.x, x_20.y, x_20.z, 0.0f));
Dm.SetColumn(2, new Vector4(x_30.x, x_30.y, x_30.z, 0.0f));
Dm.SetColumn(3, new Vector4(0.0f, 0.0f, 0.0f, 1.0f));
Matrix4x4 F = Dm * inv_Dm[tet];
//TODO: Green Strain
Matrix4x4 G = M_times_scalar(0.5f, M_minus_M(F.transpose * F, Matrix4x4.identity));
//TODO: Second PK Stress
Matrix4x4 SPS = M_plus_M(M_times_scalar(2 * stiffness_1, G), M_times_scalar(stiffness_0 * trace_M(G), Matrix4x4.identity));
Matrix4x4 P = F * SPS;
//TODO: Elastic Force
Matrix4x4 EForce = M_times_scalar(-1.0f / (6.0f * inv_Dm[tet].determinant), P * inv_Dm[tet].transpose);
Vector3 f1 = new Vector3(EForce.GetColumn(0).x, EForce.GetColumn(0).y, EForce.GetColumn(0).z);
Vector3 f2 = new Vector3(EForce.GetColumn(1).x, EForce.GetColumn(1).y, EForce.GetColumn(1).z);
Vector3 f3 = new Vector3(EForce.GetColumn(2).x, EForce.GetColumn(2).y, EForce.GetColumn(2).z);
Vector3 f0 = -(f1 + f2 + f3);
Force[Tet[tet * 4 + 0]] += f0;
Force[Tet[tet * 4 + 1]] += f1;
Force[Tet[tet * 4 + 2]] += f2;
Force[Tet[tet * 4 + 3]] += f3;
}
for (int i = 0; i < number; i++)
{
//TODO: Update X and V here.
V[i] = V[i] + Force[i] * dt / mass;
V[i] = damp * V[i];
X[i] = X[i] + V[i] * dt;
//TODO: (Particle) collision with floor.
if (X[i].y < -3)
{
X[i].y = -3;
V[i] = V[i] * 0.5f;
}
}
}
// 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();
}
private Matrix4x4 M_minus_M(Matrix4x4 lhs, Matrix4x4 rhs)
{
Matrix4x4 A = Matrix4x4.zero;
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
A[i, j] = lhs[i, j] - rhs[i, j];
}
}
// Metion here, if A[3,3] == 0 then the matrix is sigular, which doesn't have a inverse,
// however, in the singular case, the inverse is Matrix.zero;
A[3, 3] = 1;
return A;
}
private Matrix4x4 M_plus_M(Matrix4x4 lhs, Matrix4x4 rhs)
{
Matrix4x4 A = Matrix4x4.zero;
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
A[i, j] = lhs[i, j] + rhs[i, j];
}
}
// Metion here, if A[3,3] == 0 then the matrix is sigular, which doesn't have a inverse,
// however, in the singular case, the inverse is Matrix.zero;
A[3, 3] = 1;
return A;
}
private Matrix4x4 M_times_scalar(float scalar, Matrix4x4 rhs)
{
Matrix4x4 A = Matrix4x4.zero;
for (int i = 0; i < 4; i++)
{
for (int j = 0; j < 4; j++)
{
A[i, j] = scalar * rhs[i, j];
}
}
// Metion here, if A[3,3] == 0 then the matrix is sigular, which doesn't have a inverse,
// however, in the singular case, the inverse is Matrix.zero;
A[3, 3] = 1;
return A;
}
private float trace_M(Matrix4x4 lhs)
{
return lhs[0, 0] + lhs[1, 1] + lhs[2, 2];
}
}