二维稳态问题三角单元有限元法结果的可视化

利用Python中的matplotlib库实现低阶三角形单元有限元法计算结果的可视化，二维稳态问题中，有限元方法的解答是一个列向量，元素为每一个节点值，当然因为是采用的三角剖分，所以相对应的节点坐标，以及单元的节点组成矩阵都应当作为输入的参数。如果以一个二维Poisson方程的解答为例。方程为

$$ \frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} = f(x,y)\\ 其中f(x,y)表达为 f(x,y) = 2(x + y - x^2 - y^2) $$

稳态问题属于边值问题，边界条件对解的分布影响十分重要，模型尺寸为正方形，长宽均为1，当然选择正方形只是为了图个方便，有限元法本身对于几何形状并不敏感，四个边均为第一类齐次边界条件。

利用Python实现一个后处理类，不同类型的云图绑定为实例的方法，从而实现结构的后处理。具体代码为：

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation, LinearTriInterpolator
from mpl_toolkits.mplot3d import Axes3D


class View2DMode(object):
    def __init__(self, node, element, sol):
        self.node = node
        self.element = element
        self.sol = sol
        self.triang = Triangulation(node[:, 0], node[:, 1], element)

    def view_contourf(self, showlines=False, unitlabel=None, cmap="coolwarm",
                      levels=25, save=False, name=None, fontsize=8):
        fig, ax = plt.subplots()
        cf = ax.tricontourf(self.triang, self.sol, levels=levels, cmap=cmap)
        cb = fig.colorbar(cf, ax=ax)
        cb.set_label(unitlabel)
        if showlines:
            cs = ax.tricontour(self.triang, self.sol, levels=levels,
                               colors='k', linewidth=0.8)
            cs.clabel(inline=True, fontsize=fontsize, fmt="%.3f")
        if save:
            fig.savefig(name, dpi=600)
        plt.show()

    def view_colormap(self, cmap='jet', shading='flat', unitlabel=None,
                      mesh=False, lw=0.7, linecolor="k", save=False, name=None):
        fig, ax = plt.subplots()
        cp = ax.tripcolor(self.triang, self.sol, cmap=cmap, shading=shading)
        cb = fig.colorbar(cp, ax=ax)
        cb.set_label(unitlabel)
        ax.set_xlabel("x / m")
        ax.set_ylabel("y / m")
        ax.set_aspect("equal")
        if mesh:
            ax.triplot(self.triang, linewidth=lw, color=linecolor,
                       linestyle="-.")
        if save:
            fig.colorbar(name, dpi=600)
        plt.show()

    def view_3d(self, cmap="RdBu", unitlabel=None, save=False, name=None):
        fig = plt.figure()
        ax = fig.add_subplot(111, projection="3d")
        ax.plot_trisurf(self.triang, self.sol, cmap=cmap, antialiased=False,
                        shade=False)
        ax.set_xlabel("x / m")
        ax.set_ylabel("y / m")
        ax.set_zlabel(unitlabel)
        if save:
            fig.savefig(name, dpi=600)
        plt.show()

    def view_quiver(self, cmap='viridis'):
        ltri = LinearTriInterpolator(self.triang, -self.sol)
        sol_x, sol_y = ltri.gradient(self.triang.x, self.triang.y)
        sol_norm = np.sqrt(sol_x**2 + sol_y**2)
        fig, ax = plt.subplots()
        ax.quiver(self.triang.x, self.triang.y, sol_x/sol_norm,
                  sol_y/sol_norm, sol_norm, cmap=cmap)
        ax.set_aspect("equal")
        ax.set_xlabel("x/ m")
        ax.set_ylabel("y / m")
        plt.show()

    def view_intepolation_quiver(self, xn, yn, cmap='viridis', scatter=False,
                                 s=None, sco='k'):
        xmax = self.node[:, 0].max()
        xmin = self.node[:, 0].min()
        ymax = self.node[:, 1].max()
        ymin = self.node[:, 1].min()
        x = np.linspace(xmin, xmax, xn)
        y = np.linspace(ymin, ymax, yn)
        xx, yy = np.meshgrid(x, y)
        ltri = LinearTriInterpolator(self.triang,-self.sol)
        sol_x, sol_y = ltri.gradient(xx, yy)
        sol_norm = np.sqrt(sol_x**2 + sol_y**2)
        fig, ax = plt.subplots()
        ax.quiver(xx, yy, sol_x/sol_norm, sol_y/sol_norm, sol_norm, cmap=cmap)
        if scatter:
            ax.scatter(xx, yy, s=s, color=sco)
        ax.set_aspect("equal")
        ax.set_xlabel("x/ m")
        ax.set_ylabel("y / m")
        plt.show()

    def view_stream(self, xn, yn):
        xmax = self.node[:, 0].max()
        xmin = self.node[:, 0].min()
        ymax = self.node[:, 1].max()
        ymin = self.node[:, 1].min()
        x = np.linspace(xmin, xmax, xn)
        y = np.linspace(ymin, ymax, yn)
        xx, yy = np.meshgrid(x, y)
        ltri = LinearTriInterpolator(self.triang, -self.sol)
        sol_x, sol_y = ltri.gradient(xx, yy)
        sol_norm = np.sqrt(sol_x ** 2 + sol_y ** 2)
        fig, ax = plt.subplots()
        ax.streamplot(xx, yy, sol_x/sol_norm, sol_y/sol_norm, color=sol_norm,
                      cmap='coolwarm', linewidth=1)
        ax.set_aspect("equal")
        ax.set_xlabel("x/ m")
        ax.set_ylabel("y / m")
        plt.show()

    def view_contours(self, unitlabel=None, levels=25, cmap='RdBu', lw=1.0, fontsize=8):
        fig, ax = plt.subplots()
        cs = ax.tricontour(self.triang, self.sol, cmap=cmap, levels=levels,
                           linewidths=lw)
        cs.clabel(inline=True, fontsize=fontsize, colors='k')
        ax.set_aspect("equal")
        ax.set_xlabel("x / m")
        ax.set_ylabel("y / m")
        cb = fig.colorbar(cs)
        cb.set_label(unitlabel)
        plt.show()

    def view_3dscatter(self):
        fig = plt.figure()
        ax = fig.add_subplot(111, projection='3d')
        ax.scatter(self.node[:, 0], self.node[:, 1], self.sol, c=self.sol)
        ax.set_xlabel("x / m")
        ax.set_ylabel("y / m")
        plt.show()

网格划分后得到node矩阵和element矩阵，node矩阵为2列，第一列为x坐标，第二列为y坐标，由于是三节点三角形单元，因此一个单元由3个点组成，element矩阵列数为3，每一行表示一个单元的三个节点。有限元计算后的结果为sol，是一个列向量，由于本问题只是一个自由度，因此，节点个数就是sol列向量的长度。利用上述的三个参数，生成一个View2DModel对象的实例。然后直接调用实例方法即可得到云图结果。这个类中绑定了多个类型的云图，如填充的等值线图，色彩云图，三维图，矢量图、流线图，等值线图。由于矢量图和流线图通常是依据采样点来绘制的，如果网格太密，则可视化效果不好，为了避免这个问题，重新采用插值映射，进行了重采样，只需要在方法中输入间隔数即可。

本站文章如非特别说明，均为原创。未经本人同意，请勿转载！转载请务必注明出处！