求助个fortran有限元小程序问题

zhuying727 · 发表于 2012-4-24 07:43:17

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

我做个有限元小程序，是从书上面抄下来的，fortran77的程序，可是我运行的时候就说我数组越界，我也找到了问题行，可是不知道如何修改，麻烦哪位帮我看看。代码在附件中，谢谢大家了！

!##################################################################################################
program psspap
!plane stress/strain problem analysis program
!##################################################################################################
dimension lnd(500,3),x(200),y(200),jz(50,3),pj(50,3),p(500),tl(20),ak(500,100),ake(6,6)
open (6,file='psspapin.txt')
open (8,file='psspapout.txt')
read (6,10) tl
10 format(20a4)
read(6,*) nj,ne,nz,npj,ips
write(8,10) tl
if (ips==1) write(8,20)
if (ips==2) write(8,30)
20 format(/1x,'plane stress problem')
30 format(/1x,'plane strain problem')
npj0=npj
if (npj==0) npj=1
call read_in(nj,ne,nz,nd,npj,npj0,ips,e,pr,t,v,lnd,x,y,jz,pj)
n=2*nj
do i=1,n
do j=1,nd
ak(i,j)=0.0
end do
end do
do ie=1,ne
call mke(ie,nj,ne,lnd,x,y,e,pr,t,ake)
do i=1,3
do ii=1,2
ih=2*(i-1)+ii
idh=2*(lnd(ie,i)-1)+ii
do j=1,3
do jj=1,2
l=2*(j-1)+jj
il=2*(lnd(ie,j)-1)+jj
idl=il-idh+1
if (idl>0) then
ak(idh,idl)=ak(idh,idl)+ake(ih,l)
end if
end do
end do
end do
end do
end do
call mf(nj,ne,npj,npj0,n,t,v,lnd,x,y,pj,p)
call rkr(nz,nd,n,jz,ak,p)
call slov(nj,n,nd,ak,p)
call made(ne,nj,n,e,pr,lnd,x,y,p)
close(6)
close(8)
stop
end program psspap
!##################################################################################################
subroutine read_in(nj,ne,nz,nd,npj,npj0,ips,e,pr,t,v,lnd,x,y,jz,pj)
!##################################################################################################
dimension lnd(500,3),x(nj),y(nj),jz(nz,3),pj(npj,3)
read(6,*) e,pr,t,v
read(6,*) ((lnd(i,j),j=1,3),i=1,ne)
read(6,*) (x(i),y(i),i=1,nj)
read(6,*) ((jz(i,j),j=1,3),i=1,nz)
if (npj0/=0) then
read(6,*) ((pj(i,j),j=1,3),i=1,npj)
end if
nd=0
do ie=1,ne
do i=1,3
do j=1,3
iw=iabs(lnd(ie,i)-lnd(ie,j))
if (nd<iw) then
nd=iw
end if
end do
end do
end do
nd=(nd+1)*2
if (ips/=1) then
e=e/(1.0-pr*pr)
pr=pr/(1.0-pr)
end if
end
!##################################################################################################
subroutine mke(ie,nj,ne,lnd,x,y,e,pr,t,ake)
!##################################################################################################
dimension lnd(ne,3),x(nj),y(nj),b(3,6),d(3,3),s(3,6),ake(6,6)
call ma(ie,nj,ne,lnd,x,y,ae)
call md(e,pr,d)
call mb(ie,nj,ne,lnd,x,y,ae,b)
do i=1,3
do j=1,6
s(i,j)=0.0
do k=1,3
s(i,j)=s(i,j)+d(i,k)*b(k,j)
end do
end do
end do
do i=1,6
do j=1,6
ake(i,j)=0.0
do k=1,3
ake(i,j)=ake(i,j)+s(k,i)*b(k,j)*ae*t
end do
end do
end do
end
!##################################################################################################
subroutine ma(ie,nj,ne,lnd,x,y,ae)
!##################################################################################################
dimension lnd(500,3),x(nj),y(nj)
i=lnd(ie,1)
j=lnd(ie,2)
k=lnd(ie,3)
xij=x(j)-x(i)
yij=y(j)-y(i)
xik=x(k)-x(i)
yik=y(k)-y(i)
ae=0.5*(xij*yik-xik*yij)
return
end
!##################################################################################################
subroutine md(e,pr,d)
!##################################################################################################
dimension d(3,3)
do i=1,3
do j=1,3
d(i,j)=0.0
end do
end do
d(1,1)=e/(1.0-pr*pr)
d(1,2)=e*pr/(1.0-pr*pr)
d(2,1)=d(1,2)
d(2,2)=d(1,1)
d(3,3)=0.5*e/(1.0+pr)
return
end
!##################################################################################################
subroutine mb(ie,nj,ne,lnd,x,y,ae,b)
!##################################################################################################
dimension lnd(500,3),x(nj),y(nj),b(3,6)
i=lnd(ie,1)
j=lnd(ie,2)
k=lnd(ie,3)
do ii=1,3
do jj=1,6
b(ii,jj)=0.0
end do
end do
b(1,1)=y(j)-y(k)
b(1,3)=y(k)-y(i)
b(1,5)=y(i)-y(j)
b(2,2)=x(k)-x(j)
b(2,4)=x(i)-x(k)
b(2,6)=x(j)-x(i)
b(3,1)=b(2,2)
b(3,2)=b(1,1)
b(3,3)=b(2,4)
b(3,4)=b(1,3)
b(3,5)=b(2,6)
b(3,6)=b(1,5)
do i1=1,3
do j1=1,6
b(i1,j1)=0.5/ae*b(i1,j1)
end do
end do
return
end
!##################################################################################################
subroutine mf(nj,ne,npj,npj0,n,t,v,lnd,x,y,pj,p)
!##################################################################################################
dimension lnd(500,3),x(nj),y(nj),pj(npj,3),p(n)
do i=1,n
p(i)=0.0
end do
if (npj0/=0) then
do i=1,npj
ii=pj(i,1)
p(2*ii-1)=pj(i,2)
p(2*ii)=pj(i,3)
end do
end if
if (v/=0) then
do ie=1,ne
call ma(ie,nj,ne,lnd,x,y,ae)
pe=-v*ae*t/3.0
do i=1,3
ii=lnd(ie,i)
p(2*ii)=p(2*ii)+pe
end do
end do
end if
return
end
!##################################################################################################
subroutine rkr(nz,nd,n,jz,ak,p)
!##################################################################################################
dimension p(n),jz(nz,3),ak(500,100)
do i=1,nz
ir=jz(i,1)
do j=2,3
if (jz(i,j)/=0) then
ii=2*ir+j-3
ak(ii,1)=1.0
do jj=2,nd
ak(ii,jj)=0.0
end do
if (ii>nd) jo=nd
if (ii<=nd) jo=ii
do jj=2,jo
ak(ii-jj+1,jj)=0.0
end do
p(ii)=0.0
end if
end do
end do
return
end
!##################################################################################################
subroutine slov(nj,n,nd,ak,p)
!##################################################################################################
dimension p(n),ak(500,100)
nj1=n-1
do k=1,nj1
if (n>k+nd-1) im=k+nd-1
if (n<=k+nd-1) im=n
k1=k+1
do i=k1,im
l=i-k+1
c=ak(k,l)/ak(k,1)
iw=nd-l+1
do j=1,iw
m=j+i-k
ak(i,j)=ak(i,j)-c*ak(k,m)
end do
p(i)=p(i)-c*p(k)
end do
end do
p(n)=p(n)/ak(n,1)
do i1=1,nj1
i=n-i1
if (nd>n-i+1) jo=n-i+1
if (nd>n-i+1) jo=nd
do j=2,jo
k=j+i-1
p(i)=p(i)-ak(i,j)*p(k)
end do
p(i)=p(i)/ak(i,1)
end do
write(8,50)
50 format(/5x,5('**'),'result of calculation',5('**')//1x,'nodal displacements'//3x,'node',5x,'x-disp.',8x,'y-disp.')
do i=1,nj
write(8,70) i,p(2*i-1),p(2*i)
70 format(1x,i5,2e15.6)
end do
return
end
!##################################################################################################
subroutine made(ne,nj,n,e,pr,lnd,x,y,p)
!##################################################################################################
dimension lnd(500,3),x(nj),y(nj),d(3,3),b(3,6),s(3,6),st(3),p(n),de(6)
write(8,10)
10 format(/1x,'elemente stresses'/)
call md(e,pr,d)
do ie=1,ne
call ma(ie,nj,ne,lnd,x,y,ae)
call mb(ie,nj,ne,lnd,x,y,ae,b)
do i=1,3
do j=1,6
s(i,j)=0.0
do k=1,3
s(i,j)=s(i,j)+d(i,k)*b(k,j)
end do
end do
end do
do i=1,3
do j=1,2
ih=2*(i-1)+j
iw=2*(lnd(ie,i)-1)+j
de(ih)=p(iw)
end do
end do
do i=1,3
st(i)=0.0
do j=1,6
st(i)=st(i)+s(i,j)*de(j)
end do
end do
stx=st(1)
sty=st(2)
txy=st(3)
ast=(stx+sty)/2
rst=sqrt(0.25*(stx-sty)**2+txy*txy)
stma=ast+rst
stmi=ast-rst
if (sty==stmi) ceta=0.0
if (sty/=stmi) ceta=90.0-57.29578*atan(txy/(sty-stmi))
write(8,60) ie,stx,sty,txy,stma,stmi,ceta
60 format(1x,'elemetn no.=',i5/3x,'stx=',e15.6,2x,'sty=',e15.6,2x,'txy=',e15.6/3x,'stma=',e15.6,2x,'stmi=',e15.6,2x,'ceta=',e15.6)
end do
return
end

!################################################################################################## 
program psspap 
!plane stress/strain problem analysis program 
!################################################################################################## 
dimension lnd(500,3),x(200),y(200),jz(50,3),pj(50,3),p(500),tl(20),ak(500,100),ake(6,6) 
 
open (6,file='psspapin.txt') 
open (8,file='psspapout.txt') 
read (6,10) tl 
10 format(20a4) 
read(6,*) nj,ne,nz,npj,ips 
write(8,10) tl 
 
if (ips==1) write(8,20) 
if (ips==2) write(8,30) 
 
 
20 format(/1x,'plane stress problem') 
30 format(/1x,'plane strain problem') 
 
npj0=npj 
if (npj==0) npj=1 
call read_in(nj,ne,nz,nd,npj,npj0,ips,e,pr,t,v,lnd,x,y,jz,pj) 
 
n=2*nj 
 
do i=1,n 
 do j=1,nd 
  ak(i,j)=0.0 
 end do 
end do 
 
do ie=1,ne 
 call mke(ie,nj,ne,lnd,x,y,e,pr,t,ake) 
 do i=1,3 
  do ii=1,2 
   ih=2*(i-1)+ii 
   idh=2*(lnd(ie,i)-1)+ii 
   do j=1,3 
    do jj=1,2 
         l=2*(j-1)+jj 
         il=2*(lnd(ie,j)-1)+jj 
         idl=il-idh+1 
         if (idl>0) then 
          ak(idh,idl)=ak(idh,idl)+ake(ih,l) 
         end if 
        end do 
   end do 
  end do 
 end do 
end do 
 
call mf(nj,ne,npj,npj0,n,t,v,lnd,x,y,pj,p) 
call rkr(nz,nd,n,jz,ak,p) 
call slov(nj,n,nd,ak,p) 
call made(ne,nj,n,e,pr,lnd,x,y,p) 
 
close(6) 
close(8) 
 
stop 
end program psspap 
 
!################################################################################################## 
subroutine read_in(nj,ne,nz,nd,npj,npj0,ips,e,pr,t,v,lnd,x,y,jz,pj) 
!################################################################################################## 
dimension lnd(500,3),x(nj),y(nj),jz(nz,3),pj(npj,3) 
 read(6,*) e,pr,t,v 
 read(6,*) ((lnd(i,j),j=1,3),i=1,ne) 
 read(6,*) (x(i),y(i),i=1,nj) 
 read(6,*) ((jz(i,j),j=1,3),i=1,nz) 
 
 if (npj0/=0) then 
  read(6,*) ((pj(i,j),j=1,3),i=1,npj) 
 end if 
 
 nd=0 
 
 do ie=1,ne 
  do i=1,3 
   do j=1,3 
    iw=iabs(lnd(ie,i)-lnd(ie,j)) 
        if (nd<iw) then 
         nd=iw 
        end if 
   end do 
  end do 
 end do 
 
 nd=(nd+1)*2 
 
 if (ips/=1) then 
  e=e/(1.0-pr*pr) 
  pr=pr/(1.0-pr) 
 end if 
 
end 
 
!################################################################################################## 
subroutine mke(ie,nj,ne,lnd,x,y,e,pr,t,ake) 
!################################################################################################## 
dimension lnd(ne,3),x(nj),y(nj),b(3,6),d(3,3),s(3,6),ake(6,6) 
call ma(ie,nj,ne,lnd,x,y,ae) 
call md(e,pr,d) 
call mb(ie,nj,ne,lnd,x,y,ae,b) 
 
do i=1,3 
  do j=1,6 
   s(i,j)=0.0 
   do k=1,3 
    s(i,j)=s(i,j)+d(i,k)*b(k,j) 
   end do 
  end do 
 end do 
 
 do i=1,6 
  do j=1,6 
   ake(i,j)=0.0 
   do k=1,3 
    ake(i,j)=ake(i,j)+s(k,i)*b(k,j)*ae*t 
   end do 
  end do 
 end do 
 
end 
 
!################################################################################################## 
subroutine ma(ie,nj,ne,lnd,x,y,ae) 
!################################################################################################## 
dimension lnd(500,3),x(nj),y(nj) 
 
 i=lnd(ie,1) 
 j=lnd(ie,2) 
 k=lnd(ie,3) 
 xij=x(j)-x(i) 
 yij=y(j)-y(i) 
 xik=x(k)-x(i) 
 yik=y(k)-y(i) 
 ae=0.5*(xij*yik-xik*yij) 
 
return 
end 
 
!################################################################################################## 
subroutine md(e,pr,d) 
!################################################################################################## 
dimension d(3,3) 
 
 do i=1,3 
  do j=1,3 
   d(i,j)=0.0 
  end do 
 end do 
 
 d(1,1)=e/(1.0-pr*pr) 
 d(1,2)=e*pr/(1.0-pr*pr) 
 d(2,1)=d(1,2) 
 d(2,2)=d(1,1) 
 d(3,3)=0.5*e/(1.0+pr) 
 
return 
end 
 
!################################################################################################## 
subroutine mb(ie,nj,ne,lnd,x,y,ae,b) 
!################################################################################################## 
dimension lnd(500,3),x(nj),y(nj),b(3,6) 
 
 i=lnd(ie,1) 
 j=lnd(ie,2) 
 k=lnd(ie,3) 
 
 do ii=1,3 
  do jj=1,6 
   b(ii,jj)=0.0 
  end do 
 end do 
 
 b(1,1)=y(j)-y(k) 
 b(1,3)=y(k)-y(i) 
 b(1,5)=y(i)-y(j) 
 b(2,2)=x(k)-x(j) 
 b(2,4)=x(i)-x(k) 
 b(2,6)=x(j)-x(i) 
 b(3,1)=b(2,2) 
 b(3,2)=b(1,1) 
 b(3,3)=b(2,4) 
 b(3,4)=b(1,3) 
 b(3,5)=b(2,6) 
 b(3,6)=b(1,5) 
 
 do i1=1,3 
  do j1=1,6 
   b(i1,j1)=0.5/ae*b(i1,j1) 
  end do 
 end do 
 
return 
end 
 
!################################################################################################## 
subroutine mf(nj,ne,npj,npj0,n,t,v,lnd,x,y,pj,p) 
!################################################################################################## 
dimension lnd(500,3),x(nj),y(nj),pj(npj,3),p(n) 
 
 do i=1,n 
  p(i)=0.0 
 end do 
 
 if (npj0/=0) then 
  do i=1,npj 
   ii=pj(i,1) 
   p(2*ii-1)=pj(i,2) 
   p(2*ii)=pj(i,3) 
  end do 
 end if 
 
 if (v/=0) then 
  do ie=1,ne 
   call ma(ie,nj,ne,lnd,x,y,ae) 
   pe=-v*ae*t/3.0 
   do i=1,3 
    ii=lnd(ie,i) 
        p(2*ii)=p(2*ii)+pe 
   end do 
  end do 
 end if 
 
return 
end 
 
!################################################################################################## 
subroutine rkr(nz,nd,n,jz,ak,p) 
!################################################################################################## 
dimension p(n),jz(nz,3),ak(500,100) 
 
 do i=1,nz 
  ir=jz(i,1) 
  do j=2,3 
   if (jz(i,j)/=0) then 
    ii=2*ir+j-3 
        ak(ii,1)=1.0 
        do jj=2,nd 
         ak(ii,jj)=0.0 
        end do 
        if (ii>nd) jo=nd 
        if (ii<=nd) jo=ii 
        do jj=2,jo 
         ak(ii-jj+1,jj)=0.0 
        end do 
        p(ii)=0.0 
   end if 
  end do 
 end do 
 
return 
end 
 
!################################################################################################## 
subroutine slov(nj,n,nd,ak,p) 
!################################################################################################## 
dimension p(n),ak(500,100) 
 
 nj1=n-1 
 
 do k=1,nj1 
  if (n>k+nd-1) im=k+nd-1 
  if (n<=k+nd-1) im=n 
 
  k1=k+1 
  do i=k1,im 
   l=i-k+1 
   c=ak(k,l)/ak(k,1) 
   iw=nd-l+1 
   do j=1,iw 
    m=j+i-k 
        ak(i,j)=ak(i,j)-c*ak(k,m) 
   end do 
   p(i)=p(i)-c*p(k) 
  end do 
 end do 
 
 p(n)=p(n)/ak(n,1) 
 
 do i1=1,nj1 
  i=n-i1 
  if (nd>n-i+1) jo=n-i+1 
  if (nd>n-i+1) jo=nd 
 
  do j=2,jo 
   k=j+i-1 
   p(i)=p(i)-ak(i,j)*p(k) 
  end do 
  p(i)=p(i)/ak(i,1) 
 end do 
 
 write(8,50) 
50 format(/5x,5('**'),'result of calculation',5('**')//1x,'nodal displacements'//3x,'node',5x,'x-disp.',8x,'y-disp.') 
 
 do i=1,nj 
  write(8,70) i,p(2*i-1),p(2*i) 
70 format(1x,i5,2e15.6) 
 end do 
 
return 
end 
 
!################################################################################################## 
subroutine made(ne,nj,n,e,pr,lnd,x,y,p) 
!################################################################################################## 
dimension lnd(500,3),x(nj),y(nj),d(3,3),b(3,6),s(3,6),st(3),p(n),de(6) 
 
 write(8,10) 
10 format(/1x,'elemente stresses'/) 
 call md(e,pr,d) 
  do ie=1,ne 
   call ma(ie,nj,ne,lnd,x,y,ae) 
   call mb(ie,nj,ne,lnd,x,y,ae,b) 
   do i=1,3 
    do j=1,6 
         s(i,j)=0.0 
         do k=1,3 
          s(i,j)=s(i,j)+d(i,k)*b(k,j) 
         end do 
        end do 
   end do 
 
   do i=1,3 
    do j=1,2 
         ih=2*(i-1)+j 
         iw=2*(lnd(ie,i)-1)+j 
         de(ih)=p(iw) 
        end do 
   end do 
 
   do i=1,3 
    st(i)=0.0 
    do j=1,6 
         st(i)=st(i)+s(i,j)*de(j) 
        end do 
   end do 
 
   stx=st(1) 
   sty=st(2) 
   txy=st(3) 
   ast=(stx+sty)/2 
   rst=sqrt(0.25*(stx-sty)**2+txy*txy) 
   stma=ast+rst 
   stmi=ast-rst 
 
   if (sty==stmi) ceta=0.0 
   if (sty/=stmi) ceta=90.0-57.29578*atan(txy/(sty-stmi)) 
 
   write(8,60) ie,stx,sty,txy,stma,stmi,ceta 
60 format(1x,'elemetn no.=',i5/3x,'stx=',e15.6,2x,'sty=',e15.6,2x,'txy=',e15.6/3x,'stma=',e15.6,2x,'stmi=',e15.6,2x,'ceta=',e15.6) 
 
 end do 
 
return 
end

mofangbao · 发表于 2012-4-24 08:17:29

那就看看是哪个数组越界了吧把前后的变量都输出来看看没做过有限元不清楚

言深深 · 发表于 2012-4-24 08:39:11

我表示压力比较大，朋友，你是第一次做程序？从书上抄了360行的“小程序”

如果真的如楼主所言是数组越界了，那么找到错误行看看是哪一个数组错了，在最前面定义的地方改一下就可以了。（ps：我得说你一下，你都找到出错行了，怎么不说一下是哪一行呢?）360行的程序，没有人会帮你看的

另外，有限元是数值里面的概念，建议先阅读文献，再从小程序下手

topmad · 发表于 2012-4-24 08:43:34

这个小程序有点大哇

看是哪里出了问题，标记出来看一下哇

		自动登录	找回密码
密码			立即注册

[求助] 求助个fortran有限元小程序问题

登录后查看更多精彩内容~

相关帖子