如何使用 matplotlib 绘制二维 FEM 结果？

Question

我正在开发一个二维平面有限元工具。其中一项功能是能够可视化特定对象上的应力。

此工具使用以下数据创建四边形网格：

• 节点：numpy 数组[[x1 y1], [x2 y2], etc]->x和y网格中每个节点的坐标

• 元素：numpy 数组[[1 2 3 4], [2 3 5 6]]->数组的每一行对应于网格的一个特定元素的 4 个点。

我能够实现一种绘制网格的方法：

import matplotlib.pyplot as plt
import matplotlib.collections
import matplotlib.cm as cm

import numpy as np


def showMeshPlot(nodes, elements):

    y = nodes[:,0]
    z = nodes[:,1]

    #/sf/ask/3474821801/
    def quatplot(y,z, quatrangles, ax=None, **kwargs):

        if not ax: ax=plt.gca()
        yz = np.c_[y,z]
        verts= yz[quatrangles]
        pc = matplotlib.collections.PolyCollection(verts, **kwargs)
        ax.add_collection(pc)
        ax.autoscale()

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

    quatplot(y,z, np.asarray(elements), ax=None, color="crimson", facecolor="None")
    if nodes:            
        plt.plot(y,z, marker="o", ls="", color="crimson")

    plt.title('This is the plot for: quad')
    plt.xlabel('Y Axis')
    plt.ylabel('Z Axis')


    plt.show()

nodes = np.array([[0,0], [0,0.5],[0,1],[0.5,0], [0.5,0.5], [0.5,1], [1,0], 
                  [1,0.5],[1,1]])
elements = np.array([[0,3,4,1],[1,4,5,2],[3,6,7,4],[4,7,8,5]])
stresses = np.array([1,2,3,4])

showMeshPlot(nodes, elements)

这产生了这样的情节：

现在，我有一个一维数组，其中包含对象上的应力，其长度与元素数组相同。

我的问题是如何使用 matplotlib 可视化这些压力（使用标量条）？我查看了 pcolormesh，但我无法理解它如何处理我的数据。这是我想要实现的示例（归功于 robbievanleeuwen）：

注意：我无法复制上面的示例，因为他使用三角形网格而不是四边形。

提前致谢！

Answer 1

经过一段时间的思考，以下代码是使用 绘制有限元网格（具有节点标量场）的最简单方法之一matplotlib。

此解决方案基于matplotlib.pyplot.tricontourf(). 不幸的是，matplotlib如果您的有限元网格中有四边形或更高阶的元素，则没有一种简单的方法来绘制填充轮廓。为了绘制轮廓，必须首先将所有元素“切割”成三角形，例如，一个四边形可以被分割或切割成2个三角形，依此类推……

还必须使用用于绘制网格线的自定义方法，因为matplotlib.pyplot.tricontourf()仅适用于三角形网格/网格。为此，matplotlib.pyplot.fill()使用了。

这是带有简单示例的完整代码：

import matplotlib.pyplot as plt
import matplotlib.tri as tri

# converts quad elements into tri elements
def quads_to_tris(quads):
    tris = [[None for j in range(3)] for i in range(2*len(quads))]
    for i in range(len(quads)):
        j = 2*i
        n0 = quads[i][0]
        n1 = quads[i][1]
        n2 = quads[i][2]
        n3 = quads[i][3]
        tris[j][0] = n0
        tris[j][1] = n1
        tris[j][2] = n2
        tris[j + 1][0] = n2
        tris[j + 1][1] = n3
        tris[j + 1][2] = n0
    return tris

# plots a finite element mesh
def plot_fem_mesh(nodes_x, nodes_y, elements):
    for element in elements:
        x = [nodes_x[element[i]] for i in range(len(element))]
        y = [nodes_y[element[i]] for i in range(len(element))]
        plt.fill(x, y, edgecolor='black', fill=False)

# FEM data
nodes_x = [0.0, 1.0, 2.0, 0.0, 1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 3.0]
nodes_y = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0]
nodal_values = [1.0, 0.9, 1.1, 0.9, 2.1, 2.1, 0.9, 1.0, 1.0, 0.9, 0.8]
elements_tris = [[2, 6, 5], [5, 6, 10], [10, 9, 5]]
elements_quads = [[0, 1, 4, 3], [1, 2, 5, 4], [3, 4, 8, 7], [4, 5, 9, 8]]
elements = elements_tris + elements_quads

# convert all elements into triangles
elements_all_tris = elements_tris + quads_to_tris(elements_quads)

# create an unstructured triangular grid instance
triangulation = tri.Triangulation(nodes_x, nodes_y, elements_all_tris)

# plot the finite element mesh
plot_fem_mesh(nodes_x, nodes_y, elements)

# plot the contours
plt.tricontourf(triangulation, nodal_values)

# show
plt.colorbar()
plt.axis('equal')
plt.show()

哪些输出：

只需更改 FEM 数据（节点、节点值、元素），上述代码就可以用于更复杂的网格，但是，该代码仅准备处理包含三角形和四边形的网格：

您可能会注意到，对于大网格，matplotlib会变慢。同样，matplotlib 您无法可视化 3D 元素。因此，为了获得更高的效率和更多的功能，请考虑使用 VTK，例如。

Answer 2

PolyCollection 是一个ScalarMappable. 它可以有一个值数组、一个颜色图和一个标准化集。在这里，您可以将stresses数组提供给 PolyCollection 并选择一些要使用的颜色图。剩下的就是稍微重新排列函数，以便可以将附加数据作为输入并创建颜色条。

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


def showMeshPlot(nodes, elements, values):

    y = nodes[:,0]
    z = nodes[:,1]

    def quatplot(y,z, quatrangles, values, ax=None, **kwargs):

        if not ax: ax=plt.gca()
        yz = np.c_[y,z]
        verts= yz[quatrangles]
        pc = matplotlib.collections.PolyCollection(verts, **kwargs)
        pc.set_array(values)
        ax.add_collection(pc)
        ax.autoscale()
        return pc

    fig, ax = plt.subplots()
    ax.set_aspect('equal')

    pc = quatplot(y,z, np.asarray(elements), values, ax=ax, 
             edgecolor="crimson", cmap="rainbow")
    fig.colorbar(pc, ax=ax)        
    ax.plot(y,z, marker="o", ls="", color="crimson")

    ax.set(title='This is the plot for: quad', xlabel='Y Axis', ylabel='Z Axis')

    plt.show()

nodes = np.array([[0,0], [0,0.5],[0,1],[0.5,0], [0.5,0.5], [0.5,1], [1,0], 
                  [1,0.5],[1,1]])
elements = np.array([[0,3,4,1],[1,4,5,2],[3,6,7,4],[4,7,8,5]])
stresses = np.array([1,2,3,4])

showMeshPlot(nodes, elements, stresses)