使用Python计算二维截面的惯性矩和面积，中间带洞的截面也可计算，代码

因为最近到土木行业了，所以要计算二维截面的惯性矩和面积，用迈达斯太low了，所以想自己写python代码计算。

前言

主要思路是三角网划分，然后计算每个小三角的惯性矩和面积，求和得到最终截面的惯性矩和面积。
python中可用于三角划分的库有scipy.spatial.Delaunay,pytriangle,triangle,gmsh,pygmsh,shapely等，所有的通过尝试后，建议直接使用triangle，不要再用其他的了，其他的属实不行。这个库的链接以及示例代码和库函数api的链接如下：

https://rufat.be/triangle/examples.html

https://rufat.be/triangle/API.html

代码部分

下面就是最关注的代码部分：

import triangle as tr
import numpy as np
import matplotlib.pyplot as plt

def 面积和惯性矩(a, b):
    inertia_sum = 0
    areas = 0
    for i in range(len(b)):
        idx1, idx2, idx3 = b[i]
        point1, point2, point3 = a[idx1], a[idx2], a[idx3]
        x1, y1 = point1
        x2, y2 = point2
        x3, y3 = point3
        area = 0.5 * abs((x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)))
        inertia = (area * (x1**2 + x2**2 + x3**2 + x1*x2 + x2*x3 + x1*x3 + y1**2 + y2**2 + y3**2 + y1*y2 + y2*y3 + y1*y3)) / 6
        areas += area
        inertia_sum += inertia
    return areas,inertia_sum

def 分割截止线(matrix):
    num_edges = len(matrix)
    index_matrix = np.zeros((num_edges, 2), dtype=int)
    for i in range(num_edges):
        index_matrix[i] = [i, (i + 1) % num_edges]
    return index_matrix
# 下面的matrix代表的是二维矩阵，里面包含组成截面的所有顶点坐标，下面matrix1是一个示例
matrix1 = [[6390.000, -1036.050], [6390.000, -1696.033], [3382.499, -3700.000], [2795.000, -3700.000], [2495.000, -3400.000], [2495.000, -961.090], [3398.119, -678.052], [5495.699, -720.004]]
matrix = np.array(matrix1, dtype=float)
A= dict(vertices=matrix)
B= tr.triangulate(A,'a100000qlen',)
area, inertia = 面积和惯性矩(B['vertices'], B['triangles'])
print(area,inertia)
tr.compare(plt, A, B)
plt.show()

除了返回的面积和惯性矩，画出来的图应该是这样的：

左边图片是顶点位置可视化，右边是划分的三角形。
带洞的话就加个hole和segments两个参数，就可以了，具体可以私信我。

分割线
——————————————————————————————————
更新了gpu版本，计算速度可以从60s降低到4ms，原始的惯性矩计算稍有问题，主要是计算方法有问题，需要修改，不过聪明的同志从上面代码应该自己就可以修改出来了，我就不放代码了。我的最终计算结果和midas的没有区别。

带洞的三角网建立效果，使用的GPU平台为torch，注意在使用的时候要用矢量方法，类似于这种：

centroid_x = torch.sum(triangle_centroids[:, 0] * triangle_areas)
centroid_y = torch.sum(triangle_centroids[:, 1] * triangle_areas)