计算机图形学games103——作业lab3——有限元弹性体

首先看一看lab3的在unity中给我们的文件打开后是这样

首先可以看到hierachy窗口,除去本身有的camera外,给了三个游戏对象,分别是一个光源,一个地板和房子,其中房子和地板都使用了unity自身带有的网格组件,没有rigidbody组件是因为碰撞将由我们自己写的代码实现,房子和地板相比多了一个用于渲染的材质组件和一个我们自己写的脚本FVM,这个FVM就是我们要实现软体模拟的脚本组件

FVM写法

整个文件可以看做几个大步骤

一是声明变量部分

其实,开始写的时候,不一定会将变量全部声明

时间步长,质量,刚度系数(就是弹簧的K值),阻尼系数,几乎一定需要声明的

恢复系数,摩擦系数是用于碰撞检测的

接下来是对网格处理的变量声明

由于软体几乎都是用的四面体网格tet接下来需要声明的就是四面体顶点数组,以及用于储存每个四面体索引的数组

同样的在这里碰撞检测是和地面平面的碰撞检测,所以还需要声明用于求解每一个点与平面距离投影所需的变量,分别是地面上的一个点和地面的法向量,这个涉及到的计算后面会说明

最后是用于更新每一个顶点需要的数组,分别是顶点的速度,位置和力

为了实现速度平滑,还需要声明用于拉普拉斯平滑的变量(后面说明)

    float dt = 0.003f; //dt: 时间步长,用于控制模拟的更新频率。较小的时间步长可以提高模拟的精度,但会增加计算量。
    float mass = 1;//mass: 质量,每个顶点的质量。在物理计算中,用于计算力和加速度
    //float stiffness_0 = 20000.0f;
    float stiffness_0 = 5000.0f; //刚度系数,用于材料的弹性计算。这些参数决定了材料在受到外力作用时的变形程度。stiffness_0 通常表示线性刚度,stiffness_1 表示二次非线性刚度
    float stiffness_1 = 5000.0f;
    float damp = 0.999f;//阻尼系数,用于减少系统的能量,以防止振荡过大。值越接近1,阻尼效果越小
    float restitution = 0.5f;//恢复系数,表示碰撞后的恢复力。值为0表示完全非弹性碰撞,值为1表示完全弹性碰撞
    float friction = 0.5f;// 摩擦系数,用于表示接触表面之间的摩擦力


    int[] Tet;
    int tet_number;			//The number of tetrahedra

    Vector3 gravity = new Vector3(0.0f, -9.8f, 0.0f);//重力向量,通常指向负y轴,大小为9.8米每二次方秒
    Vector3 P = new Vector3(0.0f, -3.0f, 0.0f);//碰撞平面上的一点,用于检测和处理与地面的碰撞
    Vector3 N = new Vector3(0.0f, 1.0f, 0.0f);//碰撞平面的法线向量,通常指向正y轴

    Vector3[] Force; // 存储每个顶点所受的力的数组
    Vector3[] V;//存储每个顶点速度的数组
    Vector3[] X;//存储每个顶点位置的数组

    int number;             //The number of vertices 顶点的数量

    Matrix4x4[] inv_Dm;//存储每个四面体的逆变形梯度矩阵的数组,用于计算变形梯度

    //For Laplacian smoothing.
    Vector3[] V_sum;// 用于Laplacian平滑的顶点速度总和数组
    int[] V_num;// 用于Laplacian平滑的顶点邻居数量数组

    SVD svd = new SVD();//svd方法,在使用sdvt应力应变转换方式时不需要用到

第二部分就是初始化,初始化也涉及很多的内容

首先是在第一帧之前需要实现的内容,也就是void.start()

这里主要需要实现的就是获取house模型,这个模型是一个四面体网格,主要要用到的就是ele文件和node模型

首先是对ele文件的处理,读取ele文件储存在字符串后需要处理,使用空格、制表符、回车和换行符作为分割符,将其存入一个字符串数组之中

ele文件存储的是四面体文件每个顶点的索引信息,ele文件的格式例子如下

                             

其中第一个数字代表了四面体的数量,后面每一列是每一个四面体的顶点索引,比如第0个四面体是由4429,1344,4917,2344这几个索引的四面体构成的

因此基于此文件的格式,我们可以取出的第一个string就是四面体的数量,那么一共有这么多的四面体,同样也就代表最多会有四倍于四面体的顶点数量(实际上当然会有很多重叠的),我们声明一个整数int数组存储每一个四面体的顶点索引

接下来遍历每一个四面体,实际上遍历的是所有的顶点,将string转为int存入Tet数组中,这一块代码是根据ele文件的格式来的,因为我们不需要四面体索引,而且-1的原因是因为ele文件顶点的索引是从1开始的

处理完ele文件后就是node文件,和ele文件类似,但node文件存储的是网格顶点的位置向量,node文件格式例子如下

 同样的第一个数字是顶点的数量,这里就可以发现顶点的索引是从1开始的而四面体是从0开始的,同样的我们声明一个向量数组存储每一个顶点的位置

遍历所有顶点,将位置存入数组中,这里*0.4是做了一个缩放,不用在意

此时,已经做完了这个模型的导入,但实际上还有一步,是将模型中心化,就是将几何中心移动到原点,这样易于后续的处理

设置中心为0,遍历所有顶点的位置,将其相加,再除以顶点的数量,得到模型的中心点

将模型的每一个顶点位置减去中心店位置吗,得到中心化的模型

最后互换yz轴是因为unity符合的是左手定律,和我们平时用的很多软件,yz轴相反

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");//读取 house2.ele 文件的全部内容,并将其存储在 fileContent 字符串中
            string[] Strings = fileContent.Split(new char[] { ' ', '\t', '\r', '\n' }, StringSplitOptions.RemoveEmptyEntries);//使用空格、制表符、回车和换行符作为分割符,将文件内容分割成一个字符串数组 Strings

            tet_number = int.Parse(Strings[0]);//文件的第一个条目表示四面体的数量,转换为整数并赋值给 tet_number//parse转换为整数并赋值
            Tet = new int[tet_number * 4];//根据四面体数量初始化一个大小为 tet_number * 4 的整型数组,字面体顶点数组 Tet

            // index of vertices in tetrahedral
            for (int tet = 0; tet < tet_number; tet++)//使用一个循环遍历所有四面体,并解析每个四面体的顶点索引。每个四面体有4个顶点,因此每次循环读取4个索引并减1(因为文件中的索引是从1开始的,而数组索引是从0开始的
            {
                Tet[tet * 4 + 0] = int.Parse(Strings[tet * 5 + 4]) - 1;//这一块是根据ele文件来的
                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");//读取 house2.node 文件的全部内容,并将其存储在 fileContent 字符串中
            string[] Strings = fileContent.Split(new char[] { ' ', '\t', '\r', '\n' }, StringSplitOptions.RemoveEmptyEntries);//使用空格、制表符、回车和换行符作为分割符,将文件内容分割成一个字符串数组 Strings
            number = int.Parse(Strings[0]);//文件的第一个条目表示顶点的数量,转换为整数并赋值给 number
            X = new Vector3[number];//根据顶点数量初始化一个大小为 number 的 Vector3 数组 X 初始化顶点坐标数组
            for (int i = 0; i < number; i++)//使用一个循环遍历所有顶点,并解析每个顶点的坐标。每个顶点有3个坐标,因此每次循环读取3个坐标,并将它们乘以0.4(用于缩放)
            {
                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;
                // swap y and z 交换每个顶点的 y 和 z 坐标,因为使用的是unity,法则是通常使用左手坐标系
                float temp = X[i].y;
                X[i].y = X[i].z;
                X[i].z = temp;
            }
        }

做到这一步,仅仅只是将西面体网格读进来了,接下来是如何根据这些顶点来创建需要渲染的网格

因此接下来一步就是创建三角形网格

第一步 创建一个3维向量的数组,这个数组存储是三角形的顶点位置信息,每一个四面体有四个三角面,每一个三角面有3各顶点,所以这个数组一共是12倍的四面体数量,再定义一个整型用以追踪已经添加的顶点数量

遍历所有顶点将其位置信息存入三角形顶点数组,存入的逻辑是先从Tet中取出顶点的索引,再根据索引从X中取出顶点的位置

完成后再定义一个三角形的索引数组,每一个三角形有三个顶点麻,遍历全部三角形,为其生成索引即可

//Create triangle mesh.创建一个用于渲染的三角形网格
        Vector3[] vertices = new Vector3[tet_number * 12];//创建一个Vector3数组vertices,长度为 tet_number * 12。每个四面体有四个面,每个面有三个顶点,因此一个四面体需要12个顶点
        int vertex_number = 0;//定义一个计数器 vertex_number,用于跟踪已经添加的顶点数量

        // one tetrahedral mesh has four triangle mesh 遍历所有四面体,生成三角形顶点,每个四面体都有四个三角形组成
        for (int tet = 0; tet < tet_number; tet++)//通过 vertices[vertex_number++] 将每个三角形的顶点坐标添加到 vertices 数组中
        {
            vertices[vertex_number++] = X[Tet[tet * 4 + 0]];//第一个三角形由顶点 0, 2, 1 组成
            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];//创建一个 int 数组 triangles,长度为 tet_number * 12
        for (int t = 0; t < tet_number * 4; t++)//通过 triangles[t * 3 + 0] = t * 3 + 0 等操作,将顶点索引按照每个三角形的顺序添加到triangles数组中
        {
            triangles[t * 3 + 0] = t * 3 + 0;
            triangles[t * 3 + 1] = t * 3 + 1;
            triangles[t * 3 + 2] = t * 3 + 2;
        }

接下来,创建unity中的网格,这些都是unity中的基础操作,首先获取网格组件,将我们的三角形顶点位置和索引赋予网格,完成后重新计算法线确保渲染效果

到这个地步,start中关于网格的初始化就完成了,因此接下来需要初始化物理量,分别是速度和力

V_sum和V_num是关于拉普拉斯平滑用到的数组

最后一步是定义和计算逆边矩阵,这个矩阵是后续物理变形计算的变量之一,这个矩阵是一个4x4的矩阵,我们知道对于一个四面体来说,顶点有四个,那么一个顶点到另外三个顶点的向量就有三个,那么边矩阵应该是一个3x3的矩阵,为什么会是4X4,在后面定义边矩阵计算的时候说明

 //创建和更新 Unity 网格
        Mesh mesh = GetComponent<MeshFilter>().mesh; //通过 GetComponent<MeshFilter>().mesh 获取当前对象的 MeshFilter 组件,并访问其 mesh 属性
        mesh.vertices = vertices;//将生成的顶点数组赋值给网格的顶点属性
        mesh.triangles = triangles;//将生成的三角形索引数组赋值给网格的三角形属性
        mesh.RecalculateNormals();//重新计算网格的法线,以确保光照和渲染效果正确

        //初始化物理模拟所需的数组,初始化速度、力和Laplacian平滑相关的数组
        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];//通过 Build_Edge_Matrix(tet) 计算每个四面体的边矩阵,并求逆
        for (int tet = 0; tet < tet_number; tet++)
        {
            inv_Dm[tet] = Build_Edge_Matrix(tet).inverse;
        }

 下一步就是我们定义的一些用于矩阵计算的方程

第一个是如何计算边矩阵,首先遍历每一个四面体,矩阵的前三列是第一个顶点到另外三个顶点的向量,这里计算的逻辑注意不要弄错,最后一列填充一个(0,0,0,1)的列向量是将这个矩阵转换为齐次坐标,因为3x3的矩阵仅仅能表达顶点旋转和缩放,变为四维后将位移也添加了进来,在计算机图形学中,这是一个很重要的简化计算操作

剩下的一些矩阵计算就是因为本身unity不支持这些计算,都比较简单

//构建四面体的边矩阵
    Matrix4x4 Build_Edge_Matrix(int tet)
    {
        Matrix4x4 ret = Matrix4x4.zero;//初始化边矩阵
        //TODO: Need to build edge matrix here.

        for (int i = 0; i < 3; i++)//使用 tet 四面体的顶点坐标构建边矩阵。矩阵的前三列存储了三个边向量,这些向量是从四面体的第一个顶点到其他三个顶点的向量
        {
            ret[0, i] = X[Tet[tet * 4 + i + 1]].x - X[Tet[tet * 4 + 0]].x;
            ret[1, i] = X[Tet[tet * 4 + i + 1]].y - X[Tet[tet * 4 + 0]].y;
            ret[2, i] = X[Tet[tet * 4 + i + 1]].z - X[Tet[tet * 4 + 0]].z;
        }
        ret[3, 3] = 1.0f;//ret[3, 3] = 1.0f 用于将矩阵转换为齐次坐标

        return ret;
    }
    //写出矩阵减法,因为unity本身不支持
    Matrix4x4 Matrix4x4_minus_Matrix4x4(Matrix4x4 m1, Matrix4x4 m2)
    {
        Matrix4x4 res = Matrix4x4.zero;
        for (int i = 0; i < 4; i++)
            for (int j = 0; j < 4; j++)
            {
                res[i, j] = m1[i, j] - m2[i, j];
            }
        return res;
    }
    //写出矩阵加法,因为unity本身不支持
    Matrix4x4 Matrix4x4_add_Matrix4x4(Matrix4x4 m1, Matrix4x4 m2)
    {
        Matrix4x4 res = Matrix4x4.zero;
        for (int i = 0; i < 4; i++)
            for (int j = 0; j < 4; j++)
            {
                res[i, j] = m1[i, j] + m2[i, j];
            }
        return res;
    }
    //写出矩阵和标量的乘法,因为unity本身不支持
    Matrix4x4 Matrix4x4_mul_float(Matrix4x4 m, float x)
    {
        Matrix4x4 res = Matrix4x4.zero;
        for (int i = 0; i < 4; i++)
            for (int j = 0; j < 4; j++)
            {
                res[i, j] = m[i, j] * x;
            }
        return res;
    }

 接下来是拉普拉斯平滑的定义

这个算法是怎么做的,首先遍历所有四面体,将四个定点的速度相加存在sum中,对于每一个顶点我们计算它所有邻居的速度之和以及所有的邻居顶点的数量,然后使用加权的计算公式计算新的速度,新速度由自身速度和邻居速度的平均决定,这里的权重w就是平滑权重

 void Smooth_V(float w)
    {
        for (int i = 0; i < number; i++) //将 V_sum 数组中的所有元素初始化为零向量,将 V_num 数组中的所有元素初始化为0
        {
            V_sum[i] = Vector3.zero;
            V_num[i] = 0;
        }

        for (int tet = 0; tet < tet_number; tet++)//计算每个顶点的累加速度和邻居数量
        {
            Vector3 sum = V[Tet[tet * 4 + 0]] + V[Tet[tet * 4 + 1]] + V[Tet[tet * 4 + 2]] + V[Tet[tet * 4 + 3]];//遍历每个四面体,计算四个顶点速度的和 sum
            //对每个顶点,累加其邻居的速度到 V_sum,并增加邻居数量 V_num,
            //对于每个顶点,将其他三个顶点的速度相加,并将结果累加到当前顶点的 V_sum 中,同时将 V_num 加3(因为每个顶点有三个邻居
            V_sum[Tet[tet * 4 + 0]] += sum - V[Tet[tet * 4 + 0]];
            V_sum[Tet[tet * 4 + 1]] += sum - V[Tet[tet * 4 + 1]];
            V_sum[Tet[tet * 4 + 2]] += sum - V[Tet[tet * 4 + 2]];
            V_sum[Tet[tet * 4 + 3]] += sum - V[Tet[tet * 4 + 3]];
            V_num[Tet[tet * 4 + 0]] += 3;
            V_num[Tet[tet * 4 + 1]] += 3;
            V_num[Tet[tet * 4 + 2]] += 3;
            V_num[Tet[tet * 4 + 3]] += 3;
        }

        for (int i = 0; i < number; i++)//根据权重 w,计算每个顶点的平滑速度
        {
            使用加权平均计算每个顶点的新速度:w * V[i] + (1.0f - w) * V_sum[i] / V_num[i]。其中,w 是当前速度的权重,(1.0f - w) 是邻居速度的权重
            V[i] = w * V[i] + (1.0f - w) * V_sum[i] / V_num[i];
        }
    }

接下来就是核心的物理部分,这部分的代码写在_Update方法中,注意我们没有写在update中中,因为后续我们每一帧执行多次物理计算

在这个方法中首先实现重力和我们交互施加的力(其实就是给一个初始的仿真状态)

每次我们按下空格会给每一个顶点一个向上的速度

同时在每一个顶点上施加重力

 void _Update()//Update 方法实现了物理模拟的核心部分,包括处理用户输入、施加重力、计算变形梯度、应变、应力和弹性力,以及应用拉普拉斯平滑
    {
        // Jump up.用户输入部分
        if (Input.GetKeyDown(KeyCode.Space))
        {
            for (int i = 0; i < number; i++)
                V[i].y += 0.2f; //如果按下,遍历所有顶点,将每个顶点的 y 方向速度增加 0.2f
        }
        //施加重力
        for (int i = 0; i < number; i++)
        {
            //TODO: Add gravity to Force.
            Force[i] = gravity * mass;//将重力加速度 gravity 乘以质量 mass,得到重力作用力,并赋值给 Force 数组中的相应元素
        }

接下来就是计算的核心,我们如何实现线性有限元的方法模拟软体组织,原理在之前更新的有限元方法中已经介绍了,物理更新的伪代码如下

遍历所有的四边形

第一步 首先是计算D_m也就是边矩阵,这个直接由之前定义的函数Build_Edge_martix()就可以或得,且我们在初始化的使用已经计算量D_m的逆矩阵并存入了inv_D_m的数组当中

第二步 计算形变梯度,按照公式运算即可,最后是一个4x4矩阵的乘法,unity本身有实现

第三步 计算格林应变,这里涉及到矩阵和标量的乘法以及矩阵的减法,需要用到之前写好的公式

第四步 计算第一应力计算方法如下stvk方法,用形变梯度与其相乘即可

第五步计算力,线性有限元和体积有限元是等价的,我们这实际用的是体积有限元,我们使用determined方法求出体积缩放因子这是个标量,再求其和我们求出的应力与逆边矩阵的转置的乘积即可

求出力后将其施加在每一个顶点上,这里施加在每一个四面体前三个顶点上

最后一个顶点的力是根据力守恒的原则计算出来的

//FEM / FVM
        //计算变形梯度、应变、应力和弹性力
        for (int tet = 0; tet < tet_number; tet++)
        {
            //TODO: Deformation Gradient这个基本上根据公式顺下去就OK了,前面已经设置好了
            Matrix4x4 F = Build_Edge_Matrix(tet) * inv_Dm[tet];
            //TODO: Green Strain
            Matrix4x4 G = Matrix4x4_mul_float(Matrix4x4_minus_Matrix4x4(F.transpose * F, Matrix4x4.identity), 0.5f);
            //TODO: Second PK Stress
            Matrix4x4 S = Matrix4x4.zero;
            float trace = G[0, 0] + G[1, 1] + G[2, 2];
            S = Matrix4x4_add_Matrix4x4(Matrix4x4_mul_float(G, 2.0f * stiffness_1), Matrix4x4_mul_float(Matrix4x4.identity, stiffness_0 * trace));
            //TODO: Elastic Force
            float volume = 1.0f / (inv_Dm[tet].determinant * 6);
            Matrix4x4 Elastic_force = Matrix4x4_mul_float(F * S * inv_Dm[tet].transpose, -volume);
            for (int k = 0; k < 3; k++)//将弹性力应用到顶点
            {
                Force[Tet[tet * 4 + k + 1]].x += Elastic_force[0, k];
                Force[Tet[tet * 4 + k + 1]].y += Elastic_force[1, k];
                Force[Tet[tet * 4 + k + 1]].z += Elastic_force[2, k];
            }

            // Elastic_force0 = -Elastic_force1-Elastic_force2-Elastic_force3 
            //将弹性力添加到四面体的顶点上。前三个顶点的力加上 Elastic_force,第一个顶点的力减去所有其他顶点的力,以保证力的平衡
            Force[Tet[tet * 4 + 0]].x -= Elastic_force[0, 0] + Elastic_force[0, 1] + Elastic_force[0, 2];
            Force[Tet[tet * 4 + 0]].y -= Elastic_force[1, 0] + Elastic_force[1, 1] + Elastic_force[1, 2];
            Force[Tet[tet * 4 + 0]].z -= Elastic_force[2, 0] + Elastic_force[2, 1] + Elastic_force[2, 2];
        }

接着,我们施加一个速度平衡

再就是依据力更新每个顶点的位置和速度,先根据力求解速度,再由速度求解位置

做完这一步,再进行一下和地面的碰撞检测就OK了

这里的碰撞检测比较容易,我们计算每一个顶点和地面的距离,这个距离通过向量的点乘,也就是空间几何中的知识,物体上每一个顶点和地面某一个点形成的向量与地面法向量的乘积可以计算出距离的大小,因为某一向量和目标向量点乘后除以目标向量模的平方便可以得出模的大小,这里用的法向量模为1

接下来进行判断,若是距离小于0,说明碰撞发生了,碰撞发生后我们计算法向的速度分量,再做判断,若是法向的速度也是小于0.就说明物体继续往下,此时我们计算法向的速度向量和切向的速度向量,法线速度的变化由恢复系数决定,切向的权重则由恢复系数,摩擦系数决定,最终还要做一个比较,都计算完后相加赋值新的速度

Smooth_V(0.75f);
        //
        for (int i = 0; i < number; i++)
        {
            //TODO: Update X and V here.接下来更新顶点位置和速度
            V[i] = (V[i] + dt * Force[i] / mass) * damp;//Force[i] / mass 计算每个顶点所受的加速度,dt * Force[i] / mass 计算时间步长 dt 内的速度变化量
            X[i] = X[i] + V[i] * dt;//X[i] + V[i] * dt 使用更新后的速度计算新位置

            //TODO: (Particle) collision with floor.处理与地面的碰撞
            float d = Vector3.Dot(X[i] - P, N);//计算顶点与地面的距离 Vector3.Dot(X[i] - P, N) 计算顶点到地面某一点的向量在地面法线 N 上的投影,得到顶点到地面的垂直距离 d
            if (d < 0.0f) // collision occur
            {
                //计算法向速度分量 v_N_size
                float v_N_size = Vector3.Dot(V[i], N);//Vector3.Dot(V[i], N) 计算速度 V[i] 在法线方向上的分量 v_N_size
                // check velocity检查法向速度方向
                if (v_N_size < 0.0f)//v_N_size < 0.0f,表示顶点在向下运动,需要处理碰撞
                {
                    Vector3 v_N = v_N_size * N;//Vector3 v_N = v_N_size * N 计算法向速度分量 v_N
                    Vector3 v_T = V[i] - v_N;//Vector3 v_T = V[i] - v_N 计算切向速度分量 v_T
                    Vector3 v_N_new = -1.0f * restitution * v_N;//Vector3 v_N_new = -1.0f * restitution * v_N 计算碰撞后的法向速度,应用恢复系数 restitution 进行反弹
                    float a = Math.Max(1.0f - friction * (1.0f + restitution) * v_N.magnitude / v_T.magnitude, 0.0f);//计算切向速度的衰减系数,考虑摩擦系数 friction 和恢复系数 restitution
                    Vector3 v_T_new = a * v_T;//计算碰撞后的切向速度
                    V[i] = v_N_new + v_T_new;//更新顶点速度,包含新的法向速度和切向速度
                }
            }

        }
    }

最后,写一个每一帧的更新,每一迭代10次物理计算

创建一个新的顶点位置数组,遍历这些数组,将新的位置更新进入顶点数组,获取网格组件,将新的顶点位置更新并计算法向量,整个模拟就完成了

void Update()
    {
        for (int l = 0; l < 10; l++)//在每一帧中调用 _Update 方法10次。这样可以增加物理模拟的稳定性,确保每一帧中的物理计算足够精确
            _Update();

        // Dump the vertex array for rendering.更新用于渲染的顶点数组,确保模型在屏幕上正确显示
        Vector3[] vertices = new Vector3[tet_number * 12];//创建一个新的 Vector3 数组 vertices,长度为 tet_number * 12,用于存储所有四面体的顶点
        int vertex_number = 0;//初始化一个计数器 vertex_number,用于跟踪已添加的顶点数量
        for (int tet = 0; tet < tet_number; tet++)//遍历每个四面体,按照特定顺序将顶点坐标添加到 vertices 数组中。每个四面体由四个顶点组成,可以形成12个顶点(三个三角形)
        {
            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;//将新的顶点数组 vertices 赋值给网格的 vertices 属性
        mesh.RecalculateNormals();//调用 mesh.RecalculateNormals() 方法,重新计算网格的法线,以确保光照和渲染效果正确
    }
}

  • 22
    点赞
  • 9
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值