如何在Delaunay三角剖分中识别删除零内角单元（保留节点）

我用Python代码生成三角形网格，但生成的网格文件里存在含零内角的单元，这类单元无法用于后续计算。请问如何在已创建的Delaunay三角剖分中识别并删除这些单元，同时保留原有节点？

CSV文件结构：每行存储节点的X、Y坐标，部分文件首行可能是无效的[0,0]值。
问题网格情况：存在三点共线、内角为0的三角形单元。

原代码：

import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv
import os
from scipy.spatial import Delaunay

path = r"C:\Users\\"
filename = "old_mesh.csv"
datos = np.genfromtxt(path + filename, delimiter=',')

# delete X,Y first line
if np.array_equal(datos[0], np.array([0, 0])):
    datos = np.delete(datos, 0, axis=0)

#verifiy if there are NaN values
if np.isnan(datos).any():
    print("El archivo CSV contiene valores faltantes o NaN")
    # elimina las filas con valores faltantes o NaN
    datos = datos[~np.isnan(datos).any(axis=1)]

# creation of the mesh based in the nodes

tri = Delaunay(datos)

# vtu file creation 
with open(path + os.path.basename(filename)[0:-4] + '_outputmesh.vtu', "w") as f:
    f.write('<?xml version="1.0"?>\n')
    f.write('<VTKFile type="UnstructuredGrid" version="1" byte_order="LittleEndian" header_type="UInt64">\n')
    f.write('\t<UnstructuredGrid>\n')
    f.write('\t	<Piece NumberOfPoints="{}" NumberOfCells="{}">\n'.format(len(datos), len(tri.simplices)))
    f.write('\t		<Points>\n')
    f.write('\t			<DataArray type="Float64" Name="Points" NumberOfComponents="3" format="ascii">\n')
    for i in range(len(datos)):
        f.write('\t				{} {} {}\n'.format(datos[i][0], datos[i][1], 0))
    f.write('\t			</DataArray>\n')
    f.write('\t		</Points>\n')
    f.write('\t		<Cells>\n')
    f.write('\t			<DataArray type="Int64" Name="connectivity" NumberOfComponents="1" format="ascii">\n')
    for i in range(len(tri.simplices)):
        f.write('\t				{} {} {}\n'.format(tri.simplices[i][0], tri.simplices[i][1], tri.simplices[i][2]))
    f.write('\t			</DataArray>\n')
    f.write('\t			<DataArray type="Int64" Name="offsets" NumberOfComponents="1" format="ascii">\n')
    for i in range(len(tri.simplices)):
        f.write('\t				{}\n'.format((i+1)*3))
    f.write('\t			</DataArray>\n')
    f.write('\t			<DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">\n')
    for i in range(len(tri.simplices)):
        f.write('\t				5\n')
    f.write('\t			</DataArray>\n')
    f.write('\t		</Cells>\n')
    f.write('\t	</Piece>\n')
    f.write('\t</UnstructuredGrid>\n')
    f.write('</VTKFile>')

核心思路是计算每个三角形的三个内角，筛选出所有内角都大于极小阈值（避免浮点精度误差）的单元，剔除含零或接近零内角的无效单元，同时保留全部原始节点。

1. 实现逻辑 #

• 编写函数计算单个三角形的三个内角，利用向量点积公式转换为角度值
• 遍历所有Delaunay剖分生成的单元，筛选出内角均大于阈值的有效单元
• 修改VTU文件生成逻辑，仅写入有效单元

2. 修改后的完整代码 #

import numpy as np
import os
from scipy.spatial import Delaunay

def calculate_triangle_angles(points, simplex):
    """计算三角形三个内角（单位：度）"""
    # 获取三个顶点坐标
    p0, p1, p2 = points[simplex]
    
    # 计算边向量
    v0 = p1 - p0
    v1 = p2 - p0
    v2 = p0 - p1
    v3 = p2 - p1
    
    # 计算向量长度，避免三点共线时除以零
    len_v0 = np.linalg.norm(v0)
    len_v1 = np.linalg.norm(v1)
    len_v2 = np.linalg.norm(v2)
    len_v3 = np.linalg.norm(v3)
    
    if len_v0 == 0 or len_v1 == 0 or len_v2 == 0 or len_v3 == 0:
        return [0.0, 0.0, 0.0]
    
    # 计算余弦值并转换为角度，np.clip避免浮点误差导致超出[-1,1]范围
    angle0 = np.arccos(np.clip(np.dot(v0, v1) / (len_v0 * len_v1), -1.0, 1.0)) * 180 / np.pi
    angle1 = np.arccos(np.clip(np.dot(v2, v3) / (len_v2 * len_v3), -1.0, 1.0)) * 180 / np.pi
    angle2 = 180 - angle0 - angle1
    
    return [angle0, angle1, angle2]

path = r"C:\Users\\"
filename = "old_mesh.csv"
datos = np.genfromtxt(path + filename, delimiter=',')

# 删除首行无效的[0,0]
if np.array_equal(datos[0], np.array([0, 0])):
    datos = np.delete(datos, 0, axis=0)

# 清理含NaN的行
if np.isnan(datos).any():
    print("CSV文件包含缺失值或NaN，已自动清理")
    datos = datos[~np.isnan(datos).any(axis=1)]

# 生成Delaunay三角剖分
tri = Delaunay(datos)

# 筛选有效单元：内角都大于0.0001度（可根据精度需求调整阈值）
valid_simplices = []
threshold = 1e-4
for simplex in tri.simplices:
    angles = calculate_triangle_angles(datos, simplex)
    if all(angle > threshold for angle in angles):
        valid_simplices.append(simplex)

valid_simplices = np.array(valid_simplices)
print(f"原始单元数：{len(tri.simplices)}，有效单元数：{len(valid_simplices)}")

# 生成VTU文件
output_filename = path + os.path.basename(filename)[0:-4] + '_outputmesh.vtu'
with open(output_filename, "w") as f:
    f.write('<?xml version="1.0"?>\n')
    f.write('<VTKFile type="UnstructuredGrid" version="1" byte_order="LittleEndian" header_type="UInt64">\n')
    f.write('\t<UnstructuredGrid>\n')
    # 更新有效单元数量
    f.write('\t	<Piece NumberOfPoints="{}" NumberOfCells="{}">\n'.format(len(datos), len(valid_simplices)))
    f.write('\t		<Points>\n')
    f.write('\t			<DataArray type="Float64" Name="Points" NumberOfComponents="3" format="ascii">\n')
    for i in range(len(datos)):
        f.write('\t				{} {} {}\n'.format(datos[i][0], datos[i][1], 0))
    f.write('\t			</DataArray>\n')
    f.write('\t		</Points>\n')
    f.write('\t		<Cells>\n')
    f.write('\t			<DataArray type="Int64" Name="connectivity" NumberOfComponents="1" format="ascii">\n')
    for i in range(len(valid_simplices)):
        f.write('\t				{} {} {}\n'.format(valid_simplices[i][0], valid_simplices[i][1], valid_simplices[i][2]))
    f.write('\t			</DataArray>\n')
    f.write('\t			<DataArray type="Int64" Name="offsets" NumberOfComponents="1" format="ascii">\n')
    for i in range(len(valid_simplices)):
        f.write('\t				{}\n'.format((i+1)*3))
    f.write('\t			</DataArray>\n')
    f.write('\t			<DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">\n')
    for i in range(len(valid_simplices)):
        f.write('\t				5\n')
    f.write('\t			</DataArray>\n')
    f.write('\t		</Cells>\n')
    f.write('\t	</Piece>\n')
    f.write('\t</UnstructuredGrid>\n')
    f.write('</VTKFile>')

关键说明 #

• 角度阈值：设置threshold = 1e-4是为了规避浮点计算的精度误差，实际可根据网格精度需求调整。如果三点严格共线，函数会直接返回三个0度角，这类单元会被自动剔除。
• 节点保留：整个过程仅删除无效单元，所有原始节点都会保留在输出的VTU文件中，完全满足需求。
• 效率优化：如果节点数量极大，可将循环计算改为向量化运算，进一步提升处理速度。

内容的提问来源于Stack Exchange，提问作者Luis Camilo