显示以下代码的运行结果import numpy as np
import matplotlib.pyplot as plt
class Edge:
def __init__(self, n1, n2):
self.nodes = (n1, n2) # 边的两个端点坐标
self.active = True # 是否在波前中
class Triangle:
def __init__(self, n1, n2, n3):
self.nodes = (n1, n2, n3) # 三角形的三个端点坐标
class Mesh:
def __init__(self):
self.nodes = [] # 节点坐标列表 [(x1,y1), (x2,y2), ...]
self.edges = [] # 边对象列表
self.triangles = [] # 三角形对象列表
def initialize_front(mesh, boundary_points):
for i in range(len(boundary_points)):
n1 = boundary_points[i]
n2 = boundary_points[(i+1)%len(boundary_points)]
if n1 not in mesh.nodes:
mesh.nodes.append(n1)
if n2 not in mesh.nodes:
mesh.nodes.append(n2)
mesh.edges.append(Edge(n1, n2))
def compute_normal(p1, p2):
dx = p2[0] - p1[0]
dy = p2[1] - p1[1]
normal = np.array([dy, -dx]) # 法向量指向顺时针边界的内部
length = np.linalg.norm(normal)
return normal / length if length != 0 else np.zeros(2)
def advancing_front(mesh, target_length=0.1):
while True:
active_edges = [e for e in mesh.edges if e.active]
if not active_edges:
break
front_edge = active_edges[0]
n1, n2 = front_edge.nodes
# 计算新节点坐标
midpoint = (np.array(n1) + np.array(n2)) / 2
normal = compute_normal(n1, n2)
new_node = midpoint + target_length * normal
# 手动实现碰撞检测
if len(mesh.nodes) > 0:
nodes_array = np.array(mesh.nodes)
new_node_array = np.array(new_node)
# 计算所有节点到新节点的距离
distances = np.linalg.norm(nodes_array - new_node_array, axis=1)
candidates = []
# 排除当前边端点并收集候选节点
for i, (node, d) in enumerate(zip(mesh.nodes, distances)):
node_arr = np.array(node)
if not np.allclose(node_arr, np.array(n1)) and \
not np.allclose(node_arr, np.array(n2)):
candidates.append((d, i))
# 按距离排序并检查前5个候选节点
candidates.sort()
too_close = any(d < 0.8*target_length for d, _ in candidates[:5])
if too_close:
front_edge.active = False
continue
# 添加新节点和三角形
new_node_tuple = tuple(new_node.round(6)) # 限制浮点精度避免重复
mesh.nodes.append(new_node_tuple)
mesh.triangles.append(Triangle(n1, n2, new_node_tuple))
# 更新波前边
front_edge.active = False
mesh.edges.append(Edge(n1, new_node_tuple))
mesh.edges.append(Edge(new_node_tuple, n2))
# 清理不活跃边
mesh.edges = [e for e in mesh.edges if e.active]
def plot_mesh(mesh):
plt.figure(figsize=(8, 8))
# 收集所有唯一边
edges = set()
for tri in mesh.triangles:
nodes = tri.nodes
for i in range(3):
edge = tuple(sorted([nodes[i], nodes[(i+1)%3]]))
edges.add(edge)
# 绘制边
for edge in edges:
x = [edge[0][0], edge[1][0]]
y = [edge[0][1], edge[1][1]]
plt.plot(x, y, 'b-', linewidth=0.5)
# 绘制节点
nodes = np.array(mesh.nodes)
plt.scatter(nodes[:,0], nodes[:,1], c='r', s=10)
plt.axis('equal')
plt.title(f"生成网格:共{len(mesh.nodes)}个节点,{len(mesh.triangles)}个三角形")
plt.show()
# 定义初始边界(顺时针矩形)
boundary = [(0,0), (1,0), (1,1), (0,1)]
# 创建网格并初始化
mesh = Mesh()
initialize_front(mesh, boundary)
# 执行推进波前算法
advancing_front(mesh, target_length=0.15)
# 可视化结果
plot_mesh(mesh)
最新发布