有限元仿真 有限体积仿真

0 前言

声明:此篇博客仅用于个人学习记录之用,并非是分享代码。Hornor the Code

本来想仔细学一下有限元再进行LAB3的制作,之前也是这么做的。问题是咱是做游戏的,不是做cae的,养家糊口之后真没有时间去研究那些东西。所以后面的像 张量分析 之类的也就没仔细去看。只是做了矩阵微积分的那一部分。

后面的流体模拟,也是一个深坑。泪目。

这次作业的Mesh可以看作是没有索引的顶点数组。

索引数组降低了内存占用,但是带来的负面效果就是,没法做顶点的删减增加。这个缺点在三角形细分和三角形裁剪的时候尤为恐怖。

感谢 Games103王华民老师的《Intro to Physics-Based Animation!》

1 拉格朗日模拟

先从能量说起。能量分为动能和势能。物理学有一大定律,名为能量守恒。我们知道动能和势能是可以转换的。动能和势能的转换依靠的就是力。

力总是希望快速的降低势能。想象一个静止的弹簧,抓着两个端点狠狠的一拉,不放手。这个时候动能为0,但是你的手已经为弹簧做了功。弹簧中充满了潜在的能量,势能。

这个时候弹簧就会自然的产生力,来减少这个势能。

牛顿定律三条,说了一堆。核心就是第三条,力是哪来的?两个物体。最后两个物体的合力是多少,零。但是对于每一个单独的物体,却是受到力的。

万物都尽量想保持原状,但是总会有外部的物体来给与你交互施加能量,由此万物开始运动。只不过因为大家碰到了,只不过因为大家想保持原状。这就是宇宙。宇宙的合力为0。

质点弹簧系统,一根弹簧,两个质点作用。中间都是空的,弹簧由能量定义。能量由位置定义。

有限元系统,一个三角形,三个质点作用。中间都是空的,面积由能量定义。能量由位置定义。

有限体积法, 一个四面体,四个质点的作用。中间都是空的,体积由能量定义。能量由位置定义。

仅此而已。
image

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

}

posted @ 2024-04-05 12:45  Dba_sys  阅读(13)  评论(0编辑  收藏  举报