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

一个选项通过对所有元素进行三角测量然后使用 matplotlib tripcolor() 绘图来工作我现在发现的功能。两个有用的演示是 here和 here .

我的全局域的自动三角剖分可能会有问题，但单个四边形的 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()