python - python：如何绘制以节点为中心的二维不连续数据？

Question

我有一个二维数据和二维四边形网格，描述了一个细分为补丁的域。数据在每个网格节点上定义。数据中的不连续性存在于补丁边界处，即数据在同一位置被多次定义。

如何使用 Python 在节点之间使用线性插值绘制这些数据，并正确表示每个补丁面的不连续值？

下面是三个示例元素或补丁，每个元素或补丁都有六个节点值。

节点位置和值数据可能存储在[Kx3x2]数组中，其中 K 是元素的数量。例如，

x = np.array( [
[ [0.0, 1.0], [0.0, 1.0], [0.0, 1.0]  ],  #element 0
[ [1.0, 2.0], [1.0, 2.0], [1.0, 2.0]  ],  #element 1
[ [2.0, 3.0], [2.0, 3.0], [2.0, 3.0]  ],  #element 2
] )

y = np.array( [
[ [0.0, 0.0], [0.5, 0.5], [1.0, 1.0]  ],  #element 0
[ [0.0, 1.0], [0.5, 1.5], [1.0, 2.0]  ],  #element 1
[ [1.0, 1.0], [1.5, 1.5], [2.0, 2.0]  ],  #element 2
] )

z = np.array( [
[ [0.0, 0.5], [0.0, 0.8], [0.0, 1.0]  ],  #element 0
[ [0.3, 1.0], [0.6, 1.2], [0.8, 1.3]  ],  #element 1
[ [1.2, 1.5], [1.3, 1.4], [1.5, 1.7]  ],  #element 2
] )

我考虑过pyplot.imshow()。这不能同时考虑整个域并且仍然代表多值不连续节点。imshow()为每个补丁单独调用可能会起作用。但是，我将如何在同一轴上绘制每个补丁图像？ imshow()对于非矩形补丁也是有问题的，这是我的一般情况。

我考虑过pyplot.pcolormesh()，但它似乎只适用于以细胞为中心的数据。

Answer 1

一种选择是通过对所有元素进行三角测量，然后使用tripcolor()我现在发现的 matplotlib 函数进行绘图。这里和这里有两个有用的演示。

我的全局域的自动三角剖分可能会出现问题，但单个四边形的 Delaunay 三角剖分效果很好：

我通过附加每个元素的三角剖分来创建全局三角剖分。这意味着共享节点实际上在位置数组和值数组中是重复的。这允许元素面上的不连续数据。

可以使用tripcolor()函数实现所需的线性插值和不连续性绘图，提供节点位置和每个节点的值。

我有点担心等高线绘图可能如何工作，因为元素面不再在逻辑上连接。 tricontour()仍然按预期工作。（此处显示叠加了三角测量）

使用以下代码复制：

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

x = np.array( [
[ [0.0, 1.0], [0.0, 1.0], [0.0, 1.0]  ],  #element 0
[ [1.0, 2.0], [1.0, 2.0], [1.0, 2.0]  ],  #element 1
[ [2.0, 3.0], [2.0, 3.0], [2.0, 3.0]  ],  #element 2
] )

y = np.array( [
[ [0.0, 0.0], [0.5, 0.5], [1.0, 1.0]  ],  #element 0
[ [0.0, 1.0], [0.5, 1.5], [1.0, 2.0]  ],  #element 1
[ [1.0, 1.0], [1.5, 1.5], [2.0, 2.0]  ],  #element 2
] )

z = np.array( [
[ [0.0, 0.5], [0.0, 0.8], [0.0, 1.0]  ],  #element 0
[ [0.3, 1.0], [0.6, 1.2], [0.8, 1.3]  ],  #element 1
[ [1.2, 1.5], [1.3, 1.4], [1.5, 1.7]  ],  #element 2
] )



global_num_pts =  z.size
global_x = np.zeros( global_num_pts )
global_y = np.zeros( global_num_pts )
global_z = np.zeros( global_num_pts )
global_triang_list = list()

offset = 0;
num_triangles = 0;

#process triangulation element-by-element
for k in range(z.shape[0]):
    points_x = x[k,...].flatten()
    points_y = y[k,...].flatten()
    z_element = z[k,...].flatten()
    num_points_this_element = points_x.size

    #auto-generate Delauny triangulation for the element, which should be flawless due to quadrilateral element shape
    triang = tri.Triangulation(points_x, points_y)
    global_triang_list.append( triang.triangles + offset ) #offseting triangle indices by start index of this element

    #store results for this element in global triangulation arrays
    global_x[offset:(offset+num_points_this_element)] = points_x
    global_y[offset:(offset+num_points_this_element)] = points_y
    global_z[offset:(offset+num_points_this_element)] = z_element

    num_triangles += triang.triangles.shape[0]
    offset += num_points_this_element


#go back and turn all of the triangle indices into one global triangle array
offset = 0
global_triang = np.zeros( (num_triangles, 3) )
for t in global_triang_list:
    global_triang[ offset:(offset+t.shape[0] )] = t
    offset += t.shape[0]

plt.figure()
plt.gca().set_aspect('equal')

plt.tripcolor(global_x, global_y, global_triang, global_z, shading='gouraud' )
#plt.tricontour(global_x, global_y, global_triang, global_z )
#plt.triplot(global_x, global_y, global_triang, 'go-') #plot just the triangle mesh

plt.xlim((-0.25, 3.25))
plt.ylim((-0.25, 2.25))
plt.show()