AI时代：我写了一个有限元单元体积计算小软件
今天分享一个单元体积的计算程序。首先需要明白这个小程序的需求来源：
部分水利部门在做水压致裂或者一些流固耦合的计算时，需要提取出每一个单元的体积，并根据对应的公式得到每个单元的渗透力。但是MIDAS GTS只能得到网格组的体积，并不能得到每一个单元的体积。对比其他软件，我们可以知道，ansys是可以直接提取这个体积的。但是在理论上，这个功能的原理并不难，就是得到这个模型的节点坐标，根据每一个单元对应的节点坐标再次计算体积即可。其中节点坐标，单元对应的节点编号均可以直接在MIDAS中导出。
节点坐标
单元对应的节点号
方法与思路：采用字节的Trae辅助编程工具，做一个单元体积计算的小程序。
1.1 简化计算尝试
先在GTS中绘制几个最简单的体积，包含各种类型的单元的，我自己手动就可以计算出体积
参考Deepseek或者豆包提供的思路，采用将非四面体拆分成四面体的计算方法，开始做第一波的测试。为了测试出来各个单元的体积是不是正确。这里我收集的资料里面有多个体积计算方法：
后面我发现拆分为四面体和雅克比矩阵积分在计算结果上是一致的。
def hexahedron_volume_by_integration(vertices):"""使用数值积分计算六面体体积适用于任意六面体，包括扭曲的情况"""# 定义高斯积分点和权重（2点高斯积分）gauss_points = [-0.5773502691896257, 0.5773502691896257]gauss_weights = [1.0, 1.0]volume = 0.0# 遍历所有高斯积分点for i in range(2):for j in range(2):for k in range(2):xi, eta, zeta = gauss_points[i], gauss_points[j], gauss_points[k]wi, wj, wk = gauss_weights[i], gauss_weights[j], gauss_weights[k]# 计算雅可比行列式J = compute_jacobian(vertices, xi, eta, zeta)# 累加体积volume += wi * wj * wk * abs(J)return volumedef compute_jacobian(vertices, xi, eta, zeta):"""计算六面体在局部坐标(xi, eta, zeta)处的雅可比行列式"""# 形函数导数dN_dxi = np.array([-0.125*(1-eta)*(1-zeta), 0.125*(1-eta)*(1-zeta),0.125*(1+eta)*(1-zeta), -0.125*(1+eta)*(1-zeta),-0.125*(1-eta)*(1+zeta), 0.125*(1-eta)*(1+zeta),0.125*(1+eta)*(1+zeta), -0.125*(1+eta)*(1+zeta)])dN_deta = np.array([-0.125*(1-xi)*(1-zeta), -0.125*(1+xi)*(1-zeta),0.125*(1+xi)*(1-zeta), 0.125*(1-xi)*(1-zeta),-0.125*(1-xi)*(1+zeta), -0.125*(1+xi)*(1+zeta),0.125*(1+xi)*(1+zeta), 0.125*(1-xi)*(1+zeta)])dN_dzeta = np.array([-0.125*(1-xi)*(1-eta), -0.125*(1+xi)*(1-eta),-0.125*(1+xi)*(1+eta), -0.125*(1-xi)*(1+eta),0.125*(1-xi)*(1-eta), 0.125*(1+xi)*(1-eta),0.125*(1+xi)*(1+eta), 0.125*(1-xi)*(1+eta)])# 计算雅可比矩阵J = np.zeros((3, 3))for n in range(8):J[0, 0] += dN_dxi[n] * vertices[n, 0]J[0, 1] += dN_dxi[n] * vertices[n, 1]J[0, 2] += dN_dxi[n] * vertices[n, 2]J[1, 0] += dN_deta[n] * vertices[n, 0]J[1, 1] += dN_deta[n] * vertices[n, 1]J[1, 2] += dN_deta[n] * vertices[n, 2]J[2, 0] += dN_dzeta[n] * vertices[n, 0]J[2, 1] += dN_dzeta[n] * vertices[n, 1]J[2, 2] += dN_dzeta[n] * vertices[n, 2]# 返回雅可比行列式return np.linalg.det(J)
这个是拆分四面体计算方法
import numpy as npdef general_hexahedron_volume(vertices):"""计算任意六面体体积（通过8个顶点）使用分解为5个四面体的方法vertices: 8x3数组，顶点坐标"""# 确保顶点顺序正确（标准六面体顺序）# 计算一个中心点（顶点平均值）center = np.mean(vertices, axis=0)# 定义六面体的6个面（四边形）# 每个面由4个顶点组成faces = [[0, 1, 2, 3], # 底面[4, 5, 6, 7], # 顶面[0, 1, 5, 4], # 前面[1, 2, 6, 5], # 右面[2, 3, 7, 6], # 后面[3, 0, 4, 7] # 左面]# 方法1：将每个四边形面分解为2个三角形，与中心点形成四面体total_volume = 0for face in faces:# 四边形分解为两个三角形tri1 = [vertices[face[0]], vertices[face[1]], vertices[face[2]]]tri2 = [vertices[face[0]], vertices[face[2]], vertices[face[3]]]# 计算两个四面体的体积并相加volume1 = tetrahedron_volume(center, tri1[0], tri1[1], tri1[2])volume2 = tetrahedron_volume(center, tri2[0], tri2[1], tri2[2])total_volume += volume1 + volume2return total_volumedef tetrahedron_volume(a, b, c, d):"""计算四面体体积: V = |(AB · (AC × AD))| / 6"""ab = b - aac = c - aad = d - a# 三重积（混合积）triple_product = np.dot(ab, np.cross(ac, ad))return abs(triple_product) / 6.0# 示例：计算一个任意六面体vertices = np.array([[0, 0, 0], # 0[2, 0, 0], # 1[2, 3, 0], # 2[0, 3, 0], # 3[0, 0, 4], # 4[2, 0, 4], # 5[2, 3, 4], # 6[0, 3, 4] # 7])volume = general_hexahedron_volume(vertices)print(f"六面体体积: {volume}") # 24.0
具体做法是，我把这个代码直接给AI看，让他解释一下这个什么意思，然后再让他加到代码里面。
先让AI写一个初始版本的程序，计算一下这个简单模型的体积。在简单的模型体积计算成功了以后，再去测试大模型。很显然，小模型计算出来的体积和手算一致。
3.2 复杂模型体现
简单的模型计算出来以后，我们就可以尝试计算复杂的模型体积。比如我有一个总体积为108000的体积网格。里面包含各种类型的单元，然后我导出单元与节点的Excel。然后让他写一个UI界面，这个界面会比较简单，
只有两个输入窗口，一个输出窗口。
单元与节点文件，都是从GTS导出的。点击开始计算就可以在日志文件查看初始结果，点击保存即可。但是后面测试较大模型的时候，计算会稍微有些慢，于是我让他加速计算，开始用多线程一起处理。图片的测试案例的精确体积是108000，计算出来是108054，误差为0.05%，在误差范围内。
后面找了一下误差的原因，发现纯四面体的误差为0，混合网格的误差在于多面体的拆分过程产生的，如果多面体的形状较好，网格质量高，那么误差就几乎为0。部分较为扭曲的非四面体，进行拆分后的体积与实际体积会有一定差异，当然这个误差在可控范围内，是可以接受的。
3.3 封装
本次学习了一个新的技能，就是把这个py文件封装为exe文件，还是比较简单
在cmd的命令流中
其中-w为不显示黑色的命令流框，logo.ico为图标文件，inp2fpn.py为你要封装的文件，但是直接运行py文件和封装的是会有一点差异的，尤其是调用多线程的时候，所以这个时候要检查一下，然后要AI帮忙优化一下
原理如下
这个具体原理不太懂，大概解决方案如下，然后我让AI参考这个代码帮我整理了方案就好了。
