如何使用matplotlib绘制2D有限元法结果？

有限元法广泛应用于各个领域的工程和科学计算中，作为一种常用的数值计算方法，有限元法的结果可以通过绘图的形式展示给用户，方便用户理解分析。在有限元的数值计算过程中，使用Python的matplotlib库可以快速绘制出2D有限元结果，本文将介绍如何使用matplotlib绘制2D有限元法结果。

网格剖分

有限元分析的第一步就是对所分析区域进行网格剖分，将区域划分为多个小的单元，在每个单元内进行局部近似计算。在Python中使用的有限元分析包是FEniCS，可以使用meshio库生成网格，代码如下：

import meshio
import numpy as np

# 定义矩形的四个点
p1 = [0.0, 0.0]
p2 = [1.0, 0.0]
p3 = [1.0, 1.0]
p4 = [0.0, 1.0]

# 定义矩形的四个节点
points = np.array([p1, p2, p3, p4])

# 定义矩形的两个三角形
cells = {"triangle": np.array([[0, 1, 2], [0, 2, 3]])}

# 将节点和三角形组成的信息进行合并
mesh = meshio.Mesh(points=points, cells=cells)

# 将生成的网格写入到文件中
meshio.write("rectangle.vtk", mesh)

此时，该段代码将生成一个名为“rectangle.vtk”的文件，里面存放着网格信息，可以在ParaView中打开查看该网格，总共有两个三角形单元。

载入有限元方法数据

准备好网格之后，可以在Python中执行有限元方法，计算每个单元内的物理量，例如坐标、位移、应力、应变等。在这里我们以简单的弹性问题为例，使用FEniCS计算线弹性平面应变问题的位移和应力，并将结果输出到文件中。

from fenics import *

# 定义杨氏模量与泊松比
E, nu = 70e9, 0.3
# 定义拉梅常数
mu, lmbda = E / (2 * (1 + nu)), E * nu / ((1 + nu) * (1 - 2 * nu))

# 定义网格文件名
mesh_file = "rectangle.vtk"
# 从网格文件中读取网格信息
mesh = Mesh(mesh_file)

# 定义有限元函数空间
V = VectorFunctionSpace(mesh, "CG", 1)

# 定义试验函数和变分
u = TrialFunction(V)
v = TestFunction(V)
d = u.geometric_dimension()

# 定义应力张量和弹性能量密度
sigma = lambda u: 2 * mu * sym(grad(u)) + lmbda * tr(sym(grad(u))) * Identity(d)
W = lambda u: mu * tr(sym(grad(u)) * sym(grad(u))) + 0.5 * lmbda * (tr(sym(grad(u))) ** 2)

# 定义右端力和边界条件
f = Constant((0, 0))
bc = DirichletBC(V, Constant((0, 0)), "on_boundary")

# 定义弹性问题
u = Function(V)
problem = LinearVariationalProblem(
    lhs=(W(u) + dot(f, u)) * dx,
    rhs=dot(v, -sigma(u)) * dx,
    u=u,
    bcs=[bc],
)

# 解决问题并进行后处理
solver = LinearVariationalSolver(problem)
solver.parameters["linear_solver"] = "default"
solver.parameters["preconditioner"] = "default"
solver.solve()

# 将位移和应力输出到文件中
File("displacement.pvd") << u
File("stress.pvd") << project(sigma(u), TensorFunctionSpace(mesh, "CG", 1))

执行完该段代码后，会在当前目录下生成两个文件：displacement.pvd和stress.pvd。这两个文件包含有限元法计算得到的位移和应力场数据。

可视化位移和应力

位移和应力场数据都是定义在每个单元内部的，因此需要利用插值的方式将这些物理量映射到单元的节点上。在Python中，可以使用FEniCS的project函数将位移和应力场数据从单元中心插值到每个单元的节点上，并使用matplotlib将其可视化出来。

import matplotlib.pyplot as plt
import numpy as np
from fenics import *
from mpl_toolkits.mplot3d import Axes3D

# 定义位移场和应力场数据文件名
displacement_file = "displacement.pvd"
stress_file = "stress.pvd"

# 读取场数据
displacement = Function(VectorFunctionSpace(mesh, "CG", 1))
displacement.vector()[:] = np.load(displacement_file.replace(".pvd", ".npy"))
stress = Function(TensorFunctionSpace(mesh, "DG", 0))
stress.vector()[:] = np.load(stress_file.replace(".pvd", ".npy"))

# 将位移和应力场数据插值到每个单元的节点上
displacement_projection = project(displacement, VectorFunctionSpace(mesh, "CG", 1))
stress_projection = project(stress, TensorFunctionSpace(mesh, "CG", 1))

# 从插值后的场数据中提取每个节点的坐标和物理量的数值
x, y = mesh.coordinates()[:, 0], mesh.coordinates()[:, 1]
u, v = displacement_projection.split(deepcopy=True)
sxx, sxy, syy = stress_projection.split(deepcopy=True)

# 绘制位移
fig = plt.figure()
ax = fig.add_subplot(111)
ax.tricontourf(x, y, u.vector().get_local(), 100)
ax.set_aspect("equal")
plt.colorbar()
plt.title("Displacement")
plt.savefig("displacement.png")

# 绘制应力
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot_trisurf(x, y, sxx.vector().get_local(), cmap=plt.cm.jet, linewidth=0.2)
plt.title("xx Stress")
plt.savefig("stress_xx.png")

该段代码将读取之前生成的位移和应力场数据，进行插值并可视化出来。绘制位移图像使用了三角形插值法，绘制应力图像使用了三维曲面绘制方法。执行完该段代码后，会在当前目录下生成三张图片：displacement.png、stress_xx.png和stress_xy.png。这些图片展示了位移和应力场数据在矩形区域内的分布情况。