上一期，我们介绍了MULTI-2D的基本数据结构，并根据此基本数据结构使用Python编写了MULTI-2D的数据读取脚本，本期我们将根据读取到的MULTI-2D计算数据，利用Python绘制一些必要的可视化结果。
上一期传送门：
下面我们将以MULTI-2D自带的算例 D51_H_Implo 为例对数据读取脚本"multi2d.py"编程修改。
一、针对数据读取部分的修改说明
现在让我们将上一节中跑出来的 D51_H_Implo 的数据结果回顾一下：数据文件（形如 0 ）、.d类型文件（形如0.d),以及我们修改源代码输出的拓扑信息文件t.data.除此之外，还包括很多其他杂七杂八的文件，我们暂时用不到。以上这些文件都被存储在文件夹D51_H_Implo内。
现在让我们打开Python，将工作目录调整至包含文件夹D51_H_Implo的目录，开始我们的编程。
首先让我们打开.d类型文件（随便打开一个就可以，我打开的是1.01e-08.d），如图1：
图1 1.01e-08.d内容
文件里从第二行开始存储了输出的变量是啥，以及网格节点数，可以看到，一共有8个变量，按照顺序依次是：rho, x, y, vx, vy, P, TC, T.那么我们需要在读取数据的函数中按照这个这个顺序修改一下代码，如下：
import math
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import struct
import sys
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
#from scipy.fftpack import fft #用于傅里叶分析
from scipy.fftpack import fft,fftshift
from matplotlib.pylab import mpl
import os
import matplotlib.ticker as ticker
#from sympy.vector import CoordSys3D,Del #用于求场的梯度
#定义函数read_multi01，该函数作用是返回'file01.d'文件中的nrho,nfrac1,
#nfrac2,nfrac3,nx,ny,nvx,nvy,nP,nTC,nT,nmid,nrays,nDM,nDV,nTR,nlambda.
def read_multi01(filename):
file01 = open(filename,'r') #读取文件存入file01
line = file01.readline() #读取第一行存入line,返回'MF01\n'
line = file01.readline() #读取第二行存入line,例如，返回'rho 1 48362\n'
tmp = line.split() #将line按照空格拆分，例如，得到['rho', '1', '48362']列表记为tmp
nrho = int(tmp[2]) #第三个元素记为nrho
line = file01.readline()
tmp = line.split()
nx = int(tmp[2])
line = file01.readline()
tmp = line.split()
ny = int(tmp[2])
line = file01.readline()
tmp = line.split()
nvx = int(tmp[2])
line = file01.readline()
tmp = line.split()
nvy = int(tmp[2])
line = file01.readline()
tmp = line.split()
nP = int(tmp[2])
line = file01.readline()
tmp = line.split()
nTC = int(tmp[2])
line = file01.readline()
tmp = line.split()
nT = int(tmp[2])
file01.close()
return(nrho,nx,ny,nvx,nvy,nP,nTC,nT)
#定义函数read_multi02，读取file02(数据文件例如6e-10),根据数据文件读取出
#rho,frac1,frac2,frac3,x,y,vx,vy,P,TC,T,mid,rays,DM,DV,TR,lamb
def read_multi02(filename,size):
file02 = open(filename,'rb')
line = file02.read(4)
real = struct.unpack('4c',line) #对二进制数据进行解析
nrho = size[0]
nx = size[1]
ny = size[2]
nvx = size[3]
nvy = size[4]
nP = size[5]
nTC = size[6]
nT = size[7]
rho = struct.unpack('f'*size[0] ,file02.read(4*size[0]))
x = struct.unpack('f'*size[1] ,file02.read(4*size[1]))
y = struct.unpack('f'*size[2] ,file02.read(4*size[2]))
vx = struct.unpack('f'*size[3] ,file02.read(4*size[3]))
vy = struct.unpack('f'*size[4] ,file02.read(4*size[4]))
P = struct.unpack('f'*size[5] ,file02.read(4*size[5]))
TC = struct.unpack('f'*size[6] ,file02.read(4*size[6]))
T = struct.unpack('f'*size[7],file02.read(4*size[7]))
file02.close()
#return（0， 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
return(rho,x,y,vx,vy,P,TC, T)
#定义函数read_malla,读取t.data文件，获得网格信息，n_node指节点数，存在.d文件的frac3里
def read_malla(filename,n_node):
mallafile = open(filename,'r')
node_t = []
for i in range(n_node):
ia = int(mallafile.readline())-1
ib = int(mallafile.readline())-1
ic = int(mallafile.readline())-1
node_t.append([ia,ib,ic])
mallafile.close()
return(node_t)
上面三个函数分别读取了三种类型的文件，需要注意的是，函数read_multi01()和函数read_multi02()里面读取的行数都是根据我们刚刚查看的图1的内容决定的，这里一共八个变量，所以按照python的索引顺序我们就从0依次索引到7，将这些数据读取出来。
像上节一样，我们还是将一些常用函数定义下来，以备不时之需：
# get the position of each triangle element
def get_nodex(x,y,node_t,n_node):
xnode = []
ynode = []
for i in range(n_node):
xtmp = (x[node_t[i][0]] + x[node_t[i][1]] + x[node_t[i][2]])/3.0
ytmp = (y[node_t[i][0]] + y[node_t[i][1]] + y[node_t[i][2]])/3.0
xnode.append(xtmp)
ynode.append(ytmp)
return(xnode,ynode)
#get the volume of each triangle element
def get_volume(x,y,node_t,n_node):
volumelist = []
for i in range(n_node):
x1 = x[node_t[i][0]]
y1 = y[node_t[i][0]]
x2 = x[node_t[i][1]]
y2 = y[node_t[i][1]]
x3 = x[node_t[i][2]]
y3 = y[node_t[i][2]]
r = (y1+y2+y3)/3.0
# 计算三边长度
a = math.sqrt((x1 - x2)**2 + (y1 - y2)**2)
b = math.sqrt((x2 - x3)**2 + (y2 - y3)**2)
c = math.sqrt((x1 - x3)**2 + (y1 - y3)**2)
# 计算半周长
p = (a + b + c) / 2
# 计算面积
area = math.sqrt(p * (p - a) * (p - b) * (p - c))
#print("三角形面积为：", area)
#计算体积
volume = (3.1415926/3)*(y1 + y2 +y3 ) * ((x2 -x1)*(y3- y1) - (x3 - x1)*(y2 - y1))
volumelist.append(volume)
#print("三角形体积为："，volume)
return volumelist
# change cell data to point data
def cell2point(fdata,node_t,npoint,nnode):
pdata = [0]*npoint
pcount = [0]*npoint
for i in range(nnode):
pdata[node_t[i][0]] += fdata[i]
pcount[node_t[i][0]] += 1
pdata[node_t[i][1]] += fdata[i]
pcount[node_t[i][1]] += 1
pdata[node_t[i][2]] += fdata[i]
pcount[node_t[i][2]] += 1
for i in range(npoint):
#print("i=%d pcount=%d" % (i,pcount[i]))
if pcount[i] > 0 :
pdata[i] = pdata[i]/pcount[i]
return pdata
二、单张分布图绘图函数和主程序编写
我们首先编写用于绘制单张分布图的函数，为了使这个函数可以被反复使用，设置形参时要把必要信息都放进去，比如时刻、变量名、有无网格等等，如下：
def plot_distribution(quantity,size,data,nodet,title,output_path,grid = 0):
# 设置中文字体
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
if quantity == 'P' or 'Pressure':
quantity_id = 5
if quantity == 'rho':
quantity_id = 0
x = data[1];y = data[2]
nnode = size[0]; npoint = size[7]
quantity_data = cell2point(data[quantity_id],nodet,npoint,nnode)
plt.figure(dpi=120,figsize=(10,4))
triangulation = tri.Triangulation(x,y,nodet)
# pcm = plt.tripcolor(triangulation,quantity_data,
# norm=matplotlib.colors.LogNorm(),#vmin = 1e7,vmax = 1e9),
# cmap='rainbow', shading='flat')
# plt.colorbar(pcm, extend='max')#,format='%d'
if grid == 1:
# 绘制三角网格的网格线
print('Grid open!')
plt.triplot(x,y, color='b',linewidth=0.4)
# plt.xlim(-0.1,0.1)
# plt.ylim(0,0.1)
plt.title(title)
plt.xlabel('X /cm')
plt.ylabel('Y /cm')
# plt.xlim(min(x),max(x))
# plt.ylim(min(y),max(y))
if not os.path.exists(output_path+'\\'+quantity):
os.makedirs(output_path+'\\'+quantity)
plt.savefig(output_path+'\\'+quantity+'\\'+title+'.png')
在使用时，我们可以通过调节quantitity的内容绘制不同的图，通过设置title的传入内容指定图片的标题，size和data分别是我们的两个函数read_multi01()和read_multi02()的返回结果。
下面，我们编写主程序完成数据的读取和可视化，代码如下：
input_path ='D51_H_Implo'
time = '5e-09'
output_path = '绘图文件夹'
filename01 = input_path+'\\'+time+'.d'
size = read_multi01(filename01)
filename02 = input_path+'\\'+time
data = read_multi02(filename02,size)
nodet = read_malla(input_path+'\\'+'t.data',size[0])
quantity = ''
title = input_path+' '+'Var='+quantity + ' ' +'time='+time
title = '网格设置展示:5e-09'
#
plot_distribution(quantity,size,data,nodet,title,output_path,grid = 1)
上面的代码中有一些变量需要说明，1. input_path指定了要读取的文件夹，可以根据需求自行更改；2.time指定了要读取的时刻，可以根据需求自行更改; 3. output_path制定了要输出的文件夹名称，该文件夹需要提前创建； 4. title部分复写了两次，第一次是默认设置，会在标题部分标记出当前算例、物理量、时间，下面的那一次是根据具体需求写的，比如目前上面的代码直接执行将会绘制5e-09时刻时的网格分布图(grid=1)，如果不需要可以直接注释掉。
直接执行上述代码，我们将看到一张网格分布，如图2，这里在绘制网格分布的时候，我们采用了三角剖分算法，之所以能这样画，是因为我们拿到了MULTI-2D输出的网格和节点的位置信息。
图2 5e-09时的网格分布
除此之外，我们可以直接修改quantity的内容为P，并将网格关掉，把plot_distribution（）中pcm到if grid那一行之间的注释都取消掉，绘制压强分布图，修改后的代码如下，图见图3：
def plot_distribution(quantity,size,data,nodet,title,output_path,grid = 0):
# 设置中文字体
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
if quantity == 'P' or 'Pressure':
quantity_id = 5
if quantity == 'rho':
quantity_id = 0
x = data[1];y = data[2]
nnode = size[0]; npoint = size[7]
quantity_data = cell2point(data[quantity_id],nodet,npoint,nnode)
plt.figure(dpi=120,figsize=(10,4))
triangulation = tri.Triangulation(x,y,nodet)
pcm = plt.tripcolor(triangulation,quantity_data,
norm=matplotlib.colors.LogNorm(),#vmin = 1e7,vmax = 1e9),
cmap='rainbow', shading='flat')
plt.colorbar(pcm, extend='max')#,format='%d'
if grid == 1:
# 绘制三角网格的网格线
print('Grid open!')
plt.triplot(x,y, color='b',linewidth=0.4)
# plt.xlim(-0.1,0.1)
# plt.ylim(0,0.1)
plt.title(title)
plt.xlabel('X /cm')
plt.ylabel('Y /cm')
# plt.xlim(min(x),max(x))
# plt.ylim(min(y),max(y))
if not os.path.exists(output_path+'\\'+quantity):
os.makedirs(output_path+'\\'+quantity)
plt.savefig(output_path+'\\'+quantity+'\\'+title+'.png')
#以下为主程序
input_path ='D51_H_Implo'
time = '5e-09'
output_path = '绘图文件夹'
filename01 = input_path+'\\'+time+'.d'
size = read_multi01(filename01)
filename02 = input_path+'\\'+time
data = read_multi02(filename02,size)
nodet = read_malla(input_path+'\\'+'t.data',size[0])
quantity = 'P'
title = input_path+' '+'Var='+quantity + ' ' +'time='+time
#title = '网格设置展示:5e-09'
#
plot_distribution(quantity,size,data,nodet,title,output_path,grid = 0)
图3 5e-09时刻压强分布图
直接将quantity改为’rho’，绘制密度分布图如下：
图4 5e-09时刻密度分布图
至此，我们基本复现了MULTI-2D自带的可视化的软件的可视化功能，并且拥有更灵活的调整空间，下面我们针对MULTI-2D没有的功能进行完善。
三、组合子图的物理量分布函数和主程序编写
在论文中，我们往往要将物理量的分布图做成组合图，标记上子图(a)(b)©(d)以节省论文版面，这个功能无法用MULTI-2D的自带后处理软件直接实现，因此需要编程处理。为了将上面形成的分布图按照时刻的不同组合成一张新图，以体现物理量随时间的演化，我们编写了一个新的函数plot_Combination().该函数可以将四张不同的子图组合起来，并且你可以根据需要调整子图的格式（例如title，units，子图标签，刻度轴范围等），代码如下：
def plot_Combination(quantity,size1,data1,nodet1,size2,data2,nodet2,size3,data3,nodet3,
size4,data4,nodet4,title,output_path):
# 设置中文字体
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['axes.unicode_minus'] = False # 用来正常显示负号
if quantity == 'P' or 'Pressure':
quantity_id = 5
unit = 'Unit: dyn/cm^2'
if quantity == 'rho':
quantity_id = 0
unit = 'Unit: g/cm^3' #g/cm3
lo_legend = (0.2,0.2)
xlim = (0,0.22)
ylim = (0,0.22)
vmin = None
vmax = None
x1 = data1[1];y1 = data1[2]
nnode1 = size1[0]; npoint1 = size1[7]
quantity_data1 = cell2point(data1[quantity_id],nodet1,npoint1,nnode1)
x2 = data2[1];y2 = data2[2]
nnode2 = size2[0]; npoint2 = size2[7]
quantity_data2 = cell2point(data2[quantity_id],nodet2,npoint2,nnode2)
x3 = data3[1];y3 = data3[2]
nnode3 = size3[0]; npoint3 = size3[7]
quantity_data3 = cell2point(data3[quantity_id],nodet3,npoint3,nnode3)
x4 = data4[1];y4 = data4[2]
nnode4 = size4[0]; npoint4 = size4[7]
quantity_data4 = cell2point(data4[quantity_id],nodet4,npoint4,nnode4)
fig, axs = plt.subplots(2,2,dpi=150,subplot_kw={'aspect': 'equal'})
axs = axs.flatten()
triangulation = tri.Triangulation(x1,y1,nodet1)
pcm = axs[0].tripcolor(triangulation,quantity_data1,
norm=matplotlib.colors.LogNorm(vmin = vmin,vmax = vmax),
cmap='rainbow', shading='flat')
fig.colorbar(pcm,ax=axs[0], extend='max',shrink=0.7, pad=0.05,format='%.1f')
axs[0].set_title(title[0], fontsize=8)
axs[0].set_xlabel('X /cm')
axs[0].set_ylabel('Y /cm')
axs[0].text(lo_legend[0], lo_legend[1], '(a)', fontsize=12, ha='center', va='center')
# axs[0].set_xlim(min(x1), max(x1))
# axs[0].set_ylim(min(y1), max(y1))
axs[0].set_xlim(xlim[0], xlim[1])
axs[0].set_ylim(ylim[0], ylim[1])
triangulation = tri.Triangulation(x2,y2,nodet2)
pcm = axs[1].tripcolor(triangulation,quantity_data2,
norm=matplotlib.colors.LogNorm(vmin = vmin,vmax = vmax),#vmin = 1e7,vmax = 1e9),
cmap='rainbow', shading='flat')
fig.colorbar(pcm,ax=axs[1], extend='max',shrink=0.7, pad=0.05,format='%.1f')
axs[1].set_title(title[1],fontsize=8)
axs[1].set_xlabel('X /cm')
axs[1].set_ylabel('Y /cm')
axs[1].text(lo_legend[0], lo_legend[1], '(b)', fontsize=12, ha='center', va='center')
# axs[1].set_xlim(min(x2), max(x2))
# axs[1].set_ylim(min(y2), max(y2))
axs[1].set_xlim(xlim[0], xlim[1])
axs[1].set_ylim(ylim[0], ylim[1])
triangulation = tri.Triangulation(x3,y3,nodet3)
pcm = axs[2].tripcolor(triangulation,quantity_data3,
norm=matplotlib.colors.LogNorm(vmin = vmin,vmax = vmax),#vmin = 1e7,vmax = 1e9),
cmap='rainbow', shading='flat')
fig.colorbar(pcm,ax=axs[2], extend='max',shrink=0.7, pad=0.05,format='%.1f')
axs[2].set_title(title[2],fontsize=8)
axs[2].set_xlabel('X /cm')
axs[2].set_ylabel('Y /cm')
axs[2].text(lo_legend[0], lo_legend[1], '(c)', fontsize=12, ha='center', va='center')
# axs[2].set_xlim(min(x3), max(x3))
#axs[2].set_ylim(min(y3), max(y3))
axs[2].set_xlim(xlim[0], xlim[1])
axs[2].set_ylim(ylim[0], ylim[1])
triangulation = tri.Triangulation(x4,y4,nodet4)
pcm = axs[3].tripcolor(triangulation,quantity_data4,
norm=matplotlib.colors.LogNorm(vmin = vmin,vmax = vmax),#vmin = 1e7,vmax = 1e9),
cmap='rainbow', shading='flat')
fig.colorbar(pcm,ax=axs[3], extend='max',shrink=0.7, pad=0.05,format='%.2f')
axs[3].set_title(title[3],fontsize=8)
axs[3].set_xlabel('X /cm')
axs[3].set_ylabel('Y /cm')
axs[3].text(lo_legend[0],lo_legend[1], '(d)', fontsize=12, ha='center', va='center')
axs[3].set_xlim(xlim[0], xlim[1])
axs[3].set_ylim(ylim[0], ylim[1])
# 手动设置轴的刻度位置和标签
# for ax in axs.flatten():
# xticks = [-5, -2.5,0,2.5,5]
# yticks = [0,2.5,5]
# ax.set_xticks(xticks)
# ax.set_yticks(yticks)
# ax.set_xticklabels(['a', 'b', 'c', 'd', 'e'])
# ax.set_yticklabels(['zero', 'half', 'one'])
#
fig.text(0.8,0.95,unit,fontsize=8)
#plt.suptitle('')
# 调整子图之间的间距
plt.subplots_adjust(wspace=0.1, hspace=0.5)
# plt.tight_layout()
函数编写好之后，我们需要编写主程序执行，这里我们以密度分布为例，选取了四个时刻，绘制了四张子图，主程序如下：
size_li = []
data_li = []
nodet_li = []
title_li = []
time_li = ['9e-10','1.1e-09','8e-09','1.1e-08']
for i in time_li:
#input_path = 'h-Au{}1-f'.format(i)
#print(input_path)
input_path = 'D51_H_Implo'
time = i #'5.5e-09'
output_path = '绘图文件夹'
filename01 = input_path+'\\'+time+'.d'
size = read_multi01(filename01)
size_li.append(size)
filename02 = input_path+'\\'+time
data = read_multi02(filename02,size)
data_li.append(data)
nodet = read_malla(input_path+'\\'+'t.data',size[0])
nodet_li.append(nodet)
quantity = 'rho'
title = input_path+' '+'Var='+quantity + ' ' +'time='+time
title_li.append(title)
plot_Combination(quantity,size_li[0],data_li[0],nodet_li[0],
size_li[1],data_li[1],nodet_li[1],
size_li[2],data_li[2],nodet_li[2],
size_li[3],data_li[3],nodet_li[3],
title_li,output_path)
我们这里使用了一个for循环来设置要读取的信息，直接执行上面的代码，输出结果如图5所示：
图5 不同时刻的密度分布
需要注意的是，很多时候我们为了对比更明显，会设置四张子图的colorbar取值范围相同，这个可以直接在plot_Combination()函数前半部分的vmin和vmax处设置，只需要写一个能囊括所有图的corlorbar取值的数值范围就可以了。压强分布也是一样的，修改quantity的值为’P’,可以直接得到压强分布如下：
图6 不同时刻压强分布
我们注意到由于上图中压强数值很大，corlorbar以整数方式显示会很长，所以需要修改一下代码中的colorbar显示方式，具体来说，将plot_Combination()函数中的四个
fig.colorbar(pcm,ax=axs[3], extend='max',shrink=0.7, pad=0.05,format='%.2f')
中的format='%.2f’删掉，这样重新执行代码，就会发现colorbar都以科学计数法显示了：
图7 不同时刻压强分布（corlorbar显示为科学计数法）
以上就是“MULTI-2D：压强、密度、网格分布图的Python实现”的全部内容，希望对你有所帮助。
关于Python技术储备
学好 Python 不论是就业还是做副业赚钱都不错，但要学会 Python 还是要有一个学习规划。最后大家分享一份全套的 Python 学习资料，给那些想学习 Python 的小伙伴们一点帮助！
一、Python所有方向的学习路线
Python所有方向的技术点做的整理，形成各个领域的知识点汇总，它的用处就在于，你可以按照上面的知识点去找对应的学习资源，保证自己学得较为全面。
二、Python必备开发工具
三、Python视频合集
观看零基础学习视频，看视频学习是最快捷也是最有效果的方式，跟着视频中老师的思路，从基础到深入，还是很容易入门的。
四、实战案例
光学理论是没用的，要学会跟着一起敲，要动手实操，才能将自己的所学运用到实际当中去，这时候可以搞点实战案例来学习。
五、Python练习题
检查学习结果。
六、面试资料
我们学习Python必然是为了找到高薪的工作，下面这些面试题是来自阿里、腾讯、字节等一线互联网大厂最新的面试资料，并且有阿里大佬给出了权威的解答，刷完这一套面试资料相信大家都能找到满意的工作。
最后祝大家天天进步！！
上面这份完整版的Python全套学习资料已经上传至CSDN官方，朋友如果需要可以直接微信扫描下方CSDN官方认证二维码免费领取【保证100%免费】。
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
