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];
}
}