基于python的掺杂介质六面体nastran网格生成脚本

1.效果图

1.1 网格

$Mes{h}_1$	$Mes{h}_2$

1.2 应力

（使用COMSOL求解）

$misses$应力	切应力

3.程序

import numpy as np
import scipy.interpolate as si
import matplotlib.pyplot as plt

N_origin= 10 # 初始矩阵尺寸
N_inter= 50 # 插值之后矩阵尺寸
t_0= np.linspace(0,1,N_origin)
t_1= np.linspace(0,1,N_inter)
originMat= np.random.random(( N_origin, N_origin, N_origin))
interMat= np.zeros(( N_inter, N_inter, N_inter))
# mesh_01, mesh_02, mesh_03 = np.meshgrid( t_0, t_0, t_0)

# originMat[originMat> 0.3]= 1.0
# originMat[originMat<= 0.3]= 0.0
# plt.imshow(originMat[0,:,:], cmap= "gray")
# plt.colorbar()
# plt.show()

mesh_1 = np.meshgrid( t_1, t_1, t_1)
xi_list= np.zeros((N_inter*N_inter*N_inter, 3))
# [[x1, y1, z1], [x2, y2, z2],...,[xn, yn, zn]]
for i in range(N_inter):
    for j in range(N_inter):
        for k in range(N_inter):
            for m in range(3):
                xi_list[((i*N_inter)+j)*N_inter+ k, m]= mesh_1[m][i, j, k]
interList= si.interpn(points= ( t_0, t_0, t_0),values= originMat,xi= xi_list, method= "linear")
for i in range(N_inter):
    for j in range(N_inter):
        for k in range(N_inter):
            interMat[i, j, k]= interList[((i*N_inter)+j)*N_inter+ k]

# 二值化
interMat[interMat> 0.4]= 1.0
interMat[interMat<= 0.4]= 0.0

# plt.imshow(interMat[0,:,:], cmap= "gray")
# plt.colorbar()
# plt.show()
# interMat= "[读取到的矩阵]"
# N_inter = "[读取到的矩阵尺寸]"


# 开始写入网格文件
file= open("mesh.nas", mode= "w")
file.write("$ Generated by COMSOL 5.6.0.401\n$ Jul 15 2021, 15:23\nBEGIN BULK\n")
file.write("$ Grid data section\n")
for i in range(N_inter +1):
    for j in range(N_inter +1):
        for k in range(N_inter +1):
            file.write("GRID{:12d}        ".format(((i*(N_inter+1))+j)*(N_inter+1)+ k+1))
            for xi in [k/N_inter*1.0, j/N_inter*1.0, i/N_inter*1.0]:
                if abs(xi- int(xi))< 1e-6:
                    file.write("{:8.5f}".format(xi))
                else:
                    file.write("{:8.6f}".format(xi))
            file.write("\n")
file.write("$ Element data section\n")
"""
CHEXA          1       1       1       2       7       3       4       6+CONT
+CONT         15       9
"""
for i in range(N_inter):
    for j in range(N_inter):
        for k in range(N_inter):
            num_ele= (i*N_inter+j)*N_inter+ k+1
            node_1= num_ele+ (i*N_inter+ j)+(i*(N_inter+1))
            node_2= node_1+ 1
            node_3= node_1+ N_inter+1+1
            node_4= node_3- 1
            node_5= node_1 + (N_inter+1)*(N_inter+1)
            node_6= node_5 + 1
            node_7= node_5+ N_inter+1+1
            node_8= node_7- 1
            file.write("CHEXA{:11d}{:8.0f}{:8d}{:8d}{:8d}{:8d}{:8d}{:8d}+CONT\n+CONT{:11d}{:8d}\n".format(
                num_ele,interMat[i,j,k]+1,node_1, node_2, node_3, node_4, node_5, node_6, node_7, node_8 ))
file.write("$ Property data section\n")
file.write("PSOLID         1\n")
file.write("PSOLID         2\n")
file.write("ENDDATA\n")
file.close()