对flac-tecplot程序增加功能时，程序出错求助

cumtmbqq1

程序出现错误为 s_b函数运行时，出现了数组下标为非整形数的错误，程序的功能为导入tecplot后出现一个新的后处理项，即损伤值。。。。。损伤定义即在s_b函数中。可是一运行到函数，就出错。程序代：
;; by dynamax@simwe
;; 9/14/2004
;; Modifications:
;;     (1) change data packing mode from point to block
;;     (2) add zone-related info to be plotted
;;
;; 9/11/2004
;; Modifications:
;;     (1) add range tec_range to set plot range
;;     (2) add function plot_test to hide null zones
;;
;; 7/11/2004
;; Flac3d Mesh to Tecplot
;; Limits:
;;     (1)Zone Type ---- TETRAHEDRON or BRICK
;;     (2)Tecplot Version ---- Version 10
;;     (3)Only DISPLACEMENTS of Grid Points are outputted
;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; 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"SYY(Pa)" \n"SZZ(Pa)" \n"sunshang(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-14]=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 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
;; 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
                            buf(1) = buf(1) + string(z_syy(p_z)) + '  '
                     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

;定义损伤变量
def s_b
p_z = zone_head
loop while p_z # null
        array z_DDDD(n_zone)
         ef=0.1     !峰值应变
         eu=0.2     !极限应变
        ee=18*z_prop(p_z,'bu')*z_prop(p_z,'shear')/(2*z_prop(p_z,'shear')+6*z_prop(p_z,'bu'))      ;弹性模量
        uu=(3*z_prop(p_z,'bu')-2*z_prop(p_z,'shear'))/(2*z_prop(p_z,'shear')+6*z_prop(p_z,'bu'))     ;泊松比
        
         ; ee=z_prop(p_z,'young')
         ; uu=z_prop(p_z,'poisson')
;确定et，ec
         array eee(3)
         eee(1)=(1/ee)*(z_sig1(p_z)-uu*(z_sig2(p_z)+z_sig3(p_z)))   ;第1主营变
         eee(2)=(1/ee)*(z_sig2(p_z)-uu*(z_sig1(p_z)+z_sig3(p_z)))     ;第2主营变
         eee(3)=(1/ee)*(z_sig3(p_z)-uu*(z_sig1(p_z)+z_sig2(p_z)))   ;第3主营变
         ezong=sqrt(eee(1)*eee(1)+eee(2)*eee(2)+eee(3)*eee(3))  ;总应变
        
         et2=0
         ec2=0
         loop i (1,3)
              if eee(i)>0 then
               et2=et2+eee(i)*eee(i)
               else
               ec2=ec2+eee(i)*eee(i)
              end_if
         end_loop
         et=sqrt(et2)
         ec=sqrt(ec2)
;确定Dt和Dc

       if ezong=0 then
           z_DDDD(p_z)=0
       endif


       if ezong # 0 then



        if et>ef then
           Dt=(eu*(et-ef))/(ezong*(eu-ef))
        else
           Dt=0
        endif

        if ec>ef then
           Dc=(eu*(ec-ef))/(ezong*(eu-ef))
        else
           Dc=0
        endif

        at=(et/ezong)^2
        ac=(ec/ezong)^2
        z_DDDD(p_z)=at*Dt+ac*Dc
       endif
        p_z = z_next(p_z)

endloop
end
s_b

;; Write 损伤变量
def write_zone_info1
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
                 buf(1) = buf(1) + string(z_DDDD(p_z)) + '  '
                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     
     
            ;sunshang_bianliang
            ;write_zone_info1
     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

kays

z_DDDD(p_z)=0
p_z  是pointer 不能作为整数使用 使用其id即可

cumtmbqq1

kays 发表于 2011-10-12 16:09
z_DDDD(p_z)=0
p_z  是pointer 不能作为整数使用 使用其id即可

多谢。。。。。。。。。。。。，那怎么处理呢

kays

z_DDDD(z_id(p_z))=0

cumtmbqq1

kays 发表于 2011-10-13 08:21
z_DDDD(z_id(p_z))=0

谢谢，程序已经调通了。。。兄弟，你知道fish中如何表示应变吗？

cumtdabing

好好的学习下！