КТМИАД 2
Пронумеровать Лагранжевы базисные функции второго порядка и грани треугольников в сетке.
Сетка
Загрузка
# mesh = meshio.read("/work/ctmda1.msh")
# mesh = meshio.read("/work/Test1.msh")
mesh = meshio.read("/work/ctmda1-ugly.msh")
Визуализация
boundary_points = np.unique(np.concatenate(mesh.points[mesh.get_cells_type('line')]), axis=0)
fig = plt.figure(figsize=(16, 9))
ax = fig.add_subplot(111)
ax.triplot(mesh.points[:, 0], mesh.points[:, 1], 'yo--', ms=18, triangles=mesh.get_cells_type('triangle'))
plt.plot(boundary_points[:, 0], boundary_points[:, 1], "co", ms=18)
plt.tight_layout()
plt.show()
Нумерация
mesh
tris = np.sort(mesh.get_cells_type("triangle"), axis=1, kind='mergesort')
edges = np.array([
np.concatenate([tris[:, 0, None], tris[:, 1, None]], axis=1),
np.concatenate([tris[:, 1, None], tris[:, 2, None]], axis=1),
np.concatenate([tris[:, 0, None], tris[:, 2, None]], axis=1)
])
unique_edges = np.unique(edges.reshape(-1, 2), axis=0)
matches = np.all(unique_edges.reshape(-1, 1, 2) == edges.reshape(1, -1, 2), axis=2).T
edges_ids = np.argwhere(matches)[:, 1].reshape(3, -1).T
functions_ids = np.concatenate([tris, mesh.points.shape[0] + edges_ids], axis=1)
funcs_ids_edges_pos = (np.sum(mesh.points[unique_edges], axis=1) / 2)[:, :2]
funcs_ids_edges_pos = np.concatenate([mesh.points[:, :2], funcs_ids_edges_pos], axis=0)
funcs_ids_edges_vals = np.arange(unique_edges.shape[0] + mesh.points.shape[0])
boundary_points = np.unique(np.concatenate(mesh.points[mesh.get_cells_type('line')]), axis=0)
fig = plt.figure(figsize=(16, 9))
ax = fig.add_subplot(111)
ax.triplot(mesh.points[:, 0], mesh.points[:, 1], 'yo--', ms=18, triangles=mesh.get_cells_type('triangle'))
plt.plot(boundary_points[:, 0], boundary_points[:, 1], "co", ms=18)
for coord, ind in zip(funcs_ids_edges_pos, funcs_ids_edges_vals):
xy = (coord[0], coord[1])
plt.annotate(str(ind), xy, ha='center', va='center', fontsize=14)
plt.tight_layout()
plt.show()