以前用3.0版本的时候有一个程序段,可以输出一个dat文件,导入tecplot可以看到能量,换成5.0版本之后不能用了,哪位大神能帮忙修改一下?
老程序如下:
;; Initialization
def ini_mesh2tec
IO_READ = 0
IO_WRITE = 1
IO_FISH = 0
IO_ASCII = 1
N_RECORD = 5
ZONE_NGP = z_numgp(zone_head)
array buf(1)
tec_file = 'tec10.dat'
;; Edit the tec_range to set plot range
command
ran name tec_range
endcommand
end
ini_mesh2tec
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; If plotit = 1, plot the zone
def plot_test
plotit = 1
if z_model(p_z) = 'NULL' then
plotit = 0
endif
if inrange('tec_range',p_z) = 0 then
plotit = 0
endif
end
;; Get number of zones to plot
def get_nzone
n_zone = 0
p_z = zone_head
loop while p_z # null
plot_test
if plotit = 1 then
n_zone = n_zone + 1
endif
p_z = z_next(p_z)
endloop
end
get_nzone
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Write Tecplot File Head
def write_head
buf(1) = 'TITLE = "FLAC3D to Tecplot 10"\n'
buf(1) = buf(1) + 'VARIABLES = "X(m)" \n"Y(m)" \n"Z(m)" \n'
buf(1) = buf(1) + '"DISP(m)" \n"XDISP(m)" \n"YDISP(m)" \n"ZDISP(m)" \n'
buf(1) = buf(1) + '"SIG1(Pa)" \n"SIG2(Pa)" \n"SIG3(Pa)" \n'
buf(1) = buf(1) + '"SXX(Pa)" \n"ENERGY(J)" \n"SZZ(Pa)" \n'
buf(1) = buf(1) + 'ZONE T="GLOBAL" \n'
buf(1) = buf(1) + ' N=' + string(ngp) + ','
buf(1) = buf(1) + ' E=' + string(n_zone) + ','
if ZONE_NGP = 4 then
buf(1) = buf(1) + ' ZONETYPE=FETETRAHEDRON \n'
else
buf(1) = buf(1) + ' ZONETYPE=FEBrick \n'
endif
buf(1) = buf(1) + ' DATAPACKING=BLOCK \n'
buf(1) = buf(1) + ' VARLOCATION=([8-13]=CELLCENTERED) \n'
buf(1) = buf(1) + ' DT=(SINGLE SINGLE SINGLE'
buf(1) = buf(1) + ' SINGLE SINGLE SINGLE SINGLE'
buf(1) = buf(1) + ' SINGLE SINGLE SINGLE'
buf(1) = buf(1) + ' SINGLE SINGLE SINGLE )'
status = write(buf,1)
end
;; Calculate displacement magnitude
def get_gp_disp
gp_disp = gp_xdisp(p_gp)*gp_xdisp(p_gp)
gp_disp = gp_disp + gp_ydisp(p_gp)*gp_ydisp(p_gp)
gp_disp = gp_disp + gp_zdisp(p_gp)*gp_zdisp(p_gp)
gp_disp = sqrt(gp_disp)
end
;弹性模量和泊松比需在此输入
def get_z_energy
bosong=0.28
tanmo=350000000
z_energy = z_sig1(p_z)*z_sig1(p_z)
z_energy = z_energy + z_sig2(p_z)*z_sig2(p_z)
z_energy = z_energy + z_sig3(p_z)*z_sig3(p_z)
z_energy = z_energy -2*bosong* z_sig1(p_z)*z_sig2(p_z)
z_energy = z_energy -2*bosong* z_sig1(p_z)*z_sig3(p_z)
z_energy = z_energy -2*bosong* z_sig2(p_z)*z_sig3(p_z)
z_energy = z_energy/2
z_energy = z_energy/tanmo
end
;; Write gp-related data,such as Coordinates and Displacements
def write_gp_info
p_gp = gp_head
loop while p_gp # null
buf(1) = ''
loop i(1,N_RECORD)
if p_gp # null then
caseof info_flag
case 0
buf(1) = buf(1) + string(gp_xpos(p_gp)) + ' '
case 1
buf(1) = buf(1) + string(gp_ypos(p_gp)) + ' '
case 2
buf(1) = buf(1) + string(gp_zpos(p_gp)) + ' '
case 4
get_gp_disp
buf(1) = buf(1) + string(gp_disp) + ' '
case 8
buf(1) = buf(1) + string(gp_xdisp(p_gp)) + ' '
case 16
buf(1) = buf(1) + string(gp_ydisp(p_gp)) + ' '
case 32
buf(1) = buf(1) + string(gp_zdisp(p_gp)) + ' '
endcase
p_gp = gp_next(p_gp)
endif
endloop
status = write(buf,1)
endloop
end
;; Write zone-related data,such as Stresses
def write_zone_info
p_z = zone_head
loop while p_z # null
buf(1) = ''
loop i(1,N_RECORD)
if p_z # null then
plot_test
if plotit = 1 then
caseof info_flag
case 0
buf(1) = buf(1) + string(z_sig1(p_z)) + ' '
case 1
buf(1) = buf(1) + string(z_sig2(p_z)) + ' '
case 2
buf(1) = buf(1) + string(z_sig3(p_z)) + ' '
case 4
buf(1) = buf(1) + string(z_sxx(p_z)) + ' '
case 8
get_z_energy
buf(1) = buf(1) + string(z_energy) + ' '
case 16
buf(1) = buf(1) + string(z_szz(p_z)) + ' '
endcase
endif
p_z = z_next(p_z)
endif
endloop
status = write(buf,1)
endloop
end
;; Write Zone Connectivity
def write_zone
p_z = zone_head
loop while p_z # null
buf(1) = ' '
if ZONE_NGP = 4 then
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 1))) + ' '
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 2))) + ' '
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 3))) + ' '
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 4))) + ' '
else
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 1))) + ' '
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 2))) + ' '
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 5))) + ' '
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 3))) + ' '
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 4))) + ' '
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 7))) + ' '
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 8))) + ' '
buf(1) = buf(1) + string(gp_id(z_gp(p_z, 6))) + ' '
endif
plot_test
if plotit = 1 then
status = write(buf,1)
endif
p_z = z_next(p_z)
endloop
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Main Function
def mesh2tec
status = close
status = open(tec_file,IO_WRITE,IO_ASCII)
if status = 0 then
write_head
info_flag = 0
write_gp_info
info_flag = 1
write_gp_info
info_flag = 2
write_gp_info
info_flag = 4
write_gp_info
info_flag = 8
write_gp_info
info_flag = 16
write_gp_info
info_flag = 32
write_gp_info
info_flag = 0
write_zone_info
info_flag = 1
write_zone_info
info_flag = 2
write_zone_info
info_flag = 4
write_zone_info
info_flag = 8
write_zone_info
info_flag = 16
write_zone_info
write_zone
status = close
ii = out('Successfully Write ' + string(n_zone) + ' Zone Data Into File ' + tec_file)
else
ii=out('Open File Error! Status = ' + string(status))
endif
end
mesh2tec