因为最近到土木行业了,所以要计算二维截面的惯性矩和面积,用迈达斯太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两个参数,就可以了,具体可以私信我。