Python绘制FEA应力图：如何消除三角剖分边缘多余三角形？

问题根源 #

matplotlib默认的Triangulation采用Delaunay三角剖分，仅根据点的分布生成三角形，不会识别FEA网格的实际边界，因此会在边缘区域生成超出网格范围的三角形，导致彩色应力区域溢出。

解决方案 #

方案1：使用FEA原始单元数据（最优） #

如果你的FEA求解器能导出单元节点索引信息（比如每个三角形/四边形单元对应的3个或4个节点索引，四边形可拆分为两个三角形），直接用这些单元构建三角剖分，完全避免多余三角形：

# 假设已从FEA导出单元节点索引数组elements，每个元素为单元的3个节点索引
triang = tri.Triangulation(x, y, triangles=elements)

这种方式生成的三角剖分完全匹配原始FEA网格，不会出现边界外的多余三角形。

方案2：过滤超出边界的三角形（无单元数据时） #

若无法获取单元数据，可通过两种方式过滤异常三角形：

方法A：基于三角形尺寸过滤

利用TriAnalyzer计算三角形外接圆半径，过滤掉远大于平均尺寸的异常三角形：

triang = tri.Triangulation(x, y)
# 分析三角剖分，获取所有三角形的外接圆半径
tri_analyzer = tri.TriAnalyzer(triang)
# 设置半径阈值（示例为平均半径的2倍，可根据实际调整）
radius_threshold = np.mean(tri_analyzer.circumcircle_radii()) * 2
# 生成异常三角形掩码
mask = tri_analyzer.get_flat_tri_mask(radius_threshold)
# 应用掩码，过滤多余三角形
triang.set_mask(mask)

方法B：基于凸包边界过滤

通过计算点集的凸包，判断三角形重心是否在凸包内部，过滤掉凸包外的三角形：

from scipy.spatial import ConvexHull

triang = tri.Triangulation(x, y)
# 计算点集的凸包
hull = ConvexHull(np.vstack((x, y)).T)
A = hull.equations[:, :-1]
b = hull.equations[:, -1]

# 生成掩码：标记重心在凸包外的三角形
mask = []
for tri_nodes in triang.triangles:
    # 计算三角形重心
    centroid_x = (x[tri_nodes[0]] + x[tri_nodes[1]] + x[tri_nodes[2]]) / 3
    centroid_y = (y[tri_nodes[0]] + y[tri_nodes[1]] + y[tri_nodes[2]]) / 3
    # 判断重心是否在凸包内
    in_hull = np.all(A @ np.array([centroid_x, centroid_y]) + b <= 1e-6)
    mask.append(not in_hull)

triang.set_mask(mask)

修改后的完整代码示例（以方案2的方法A为例） #

import matplotlib as mpl
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri

for orient in ['top', 'bot', 'side']: 
    x = []
    y = []
    z = []
    stress = []
    with open('data.txt') as file:
        for line in file:
            cur_line = line.split('\t')
            cur_x_old = cur_line[0]
            cur_y_old = cur_line[1]
            cur_z_old = cur_line[2]
            cur_s_old = cur_line[3]
            
            if cur_x_old == 'X Location (mm)':
                continue
            # 转换数值格式并添加到列表
            cur_x = float(cur_x_old.replace(",", "."))
            cur_y = float(cur_y_old.replace(",", "."))
            cur_z = float(cur_z_old.replace(",", "."))
            cur_s = float(cur_s_old.replace(",", "."))
            
            x.append(cur_x)
            y.append(cur_y)
            z.append(cur_z)
            stress.append(cur_s)
                
    stress = np.array(stress)
    x = np.array(x)
    y = np.array(y)
    z = np.array(z)
    
    levels = np.linspace(stress.min(), stress.max(), num=100)
    
    if orient == 'side':
        plt.figure(figsize=(max(x)/50, abs(min(y))/50))
        # 构建三角剖分并过滤多余三角形
        triang = tri.Triangulation(x, y)
        tri_analyzer = tri.TriAnalyzer(triang)
        radius_threshold = np.mean(tri_analyzer.circumcircle_radii()) * 2
        mask = tri_analyzer.get_flat_tri_mask(radius_threshold)
        triang.set_mask(mask)
        
        plt.tricontourf(triang, stress, cmap='jet', norm=mpl.colors.Normalize(0, 100), levels=levels, extend='max')
        plt.scatter(x, y, color='k')
    else:
        plt.figure(figsize=(max(x)/50, max(z)*2/50))
        # 对x-z平面做同样的三角剖分与过滤
        triang = tri.Triangulation(x, z)
        tri_analyzer = tri.TriAnalyzer(triang)
        radius_threshold = np.mean(tri_analyzer.circumcircle_radii()) * 2
        mask = tri_analyzer.get_flat_tri_mask(radius_threshold)
        triang.set_mask(mask)
        
        plt.tricontourf(triang, stress, cmap='jet', norm=mpl.colors.Normalize(0, 100), levels=levels)

plt.show()

注意事项 #

• 方案1精度最高，建议优先从FEA求解器导出单元数据
• 方案2的阈值需根据点集密度调整，若2倍平均半径效果不佳，可尝试1.5或3倍
• 对于非凸的复杂网格，凸包方法可能不够准确，此时优先使用方案1，或借助shapely库创建自定义边界多边形来判断三角形位置

内容的提问来源于stack exchange，提问作者Ian Venter