线弹性小变形空间板壳静力有限元计算程序.docx
《线弹性小变形空间板壳静力有限元计算程序.docx》由会员分享,可在线阅读,更多相关《线弹性小变形空间板壳静力有限元计算程序.docx(56页珍藏版)》请在冰豆网上搜索。
线弹性小变形空间板壳静力有限元计算程序
元计算有限元自动生成系统所开发源代码系列
线弹性小变形空间板壳静力有限元计算程序
1.简介
元计算(www.ectec.asia)公司所开发的并行有限元程序自动生成系统(pFEPG)可根据用户需要开发出各种有限元计算程序源代码。
该源代码系列即为pFEPG所开发出来的求解各学科典型问题的有限元计算程序。
该组程序为线弹性小变形空间板壳静力有限元计算程序。
2.starta.for,对位移场的数据进行初始化;
implicitreal*8(a-h,o-z)
character*12fname,filename(20)
common/aa/ia(250000000)
common/bb/ib(125000000)
c....opendisp0filetogetthenumbersofnodesanddegreeoffreedom
c....knode....numberofnodes,kdgof....numberofd.o.f.
open(1,file='',form='unformatted')
read
(1)knode,kdgof
close
(1)
kvar=knode*kdgof
write(*,*)'knode,kdgof,kvar='
write(*,'(1x,4i7)')knode,kdgof,kvar
kvar1=kvar+1
kcoor=3
kelem=31250000
knb1=kdgof*knode*1
if(knb1/2*2.lt.knb1)knb1=knb1+1
kna4=kcoor*knode*2
kna1=kdgof*knode*2
kna2=kdgof*knode*2
kna3=kdgof*knode*2
kna5=knode*1
if(kna5/2*2.lt.kna5)kna5=kna5+1
knb4=kelem*1
if(knb4/2*2.lt.knb4)knb4=knb4+1
knb2=kvar1*1
if(knb2/2*2.lt.knb2)knb2=knb2+1
knb3=kvar1*1
if(knb3/2*2.lt.knb3)knb3=knb3+1
kna0=1
kna1=kna1+kna0
kna2=kna2+kna1
kna3=kna3+kna2
kna4=kna4+kna3
kna5=kna5+kna4
if(kna5-1.gt.250000000)then
write(*,*)'exceedmemoryofarrayia'
write(*,*)'memoryofia=250000000'
write(*,*)'memoryneeded=',kna5,'inprgramstart'
stop55555
endif
knb0=1
knb1=knb1+knb0
knb2=knb2+knb1
knb3=knb3+knb2
knb4=knb4+knb3
if(knb4-1.gt.125000000)then
write(*,*)'exceedmemoryofarrayib'
write(*,*)'memoryofib=125000000'
write(*,*)'memoryneeded=',knb4,'inprgramstart'
stop55555
endif
callstart(knode,kdgof,kcoor,kvar,
*kelem,maxt,kvar1,ia(kna0),ia(kna1),ia(kna2),
*ia(kna3),ia(kna4),ib(knb0),ib(knb1),ib(knb2),
*ib(knb3),
*filename)
end
subroutinestart(knode,kdgof,kcoor,kvar,
*kelem,maxt,kvar1,u0,u1,u2,
*coor,inodvar,nodvar,numcol,lm,node,
*filename)
implicitreal*8(a-h,o-z)
character*12filename(20)
DIMENSIONNODVAR(KDGOF,KNODE),COOR(KCOOR,KNODE),R(3),
*U0(KDGOF,KNODE),U1(KDGOF,KNODE),U2(KDGOF,KNODE),
*INODVAR(KNODE),node(kelem)
DIMENSIONNUMCOL(KVAR1),LM(KVAR1)
CHARACTER*1MATERIAL
logicalfilflg
C.................................................................
C.....KDGOFNUMBEROFD.O.F
C.....KNODENUMBEROFNODES
C.....INODVARIDDATA
C.....NODVARDENOTETHEEQUATIONNUMBERCORRESPONDINGTHED.O.F
C.....U0U1U2INITIALVALUE
C.....COORCOORDINATES
C.....NODEELEMENTNODALCONNECTION
C.................................................................
6FORMAT(1X,15I4)
7FORMAT(1X,8F9.3)
C.......OPENIDfile
OPEN(1,FILE='',FORM='UNFORMATTED',STATUS='OLD')
READ
(1)NUMNOD,NODDOF,((NODVAR(I,J),I=1,NODDOF),J=1,NUMNOD)
CLOSE
(1)
callchms(kdgof,knode,NODVAR)
cWRITE(*,*)'NUMNOD=',NUMNOD,'NODDOF=',NODDOF
cWRITE(*,*)'ID='
cWRITE(*,6)((NODVAR(I,J),I=1,NODDOF),J=1,NUMNOD)
C.....GETTHENATURALNODALORDER
DO12N=1,KNODE
INODVAR(N)=N
12CONTINUE
C.....OPENORDER.NODFILEANDREADTHENODALORDERIFTHEFILEEXIST
inquire(file='ORDER.NOD',exist=filflg)
if(filflg)then
OPEN(1,FILE='ORDER.NOD',FORM='UNFORMATTED',STATUS='OLD')
READ
(1)(INODVAR(I),I=1,NUMNOD)
CLOSE
(1)
WRITE(*,*)'NODORDER='
WRITE(*,6)(INODVAR(I),I=1,NUMNOD)
endif
C.....GETNVBYID
NEQ=0
DO20JNOD=1,NUMNOD
J=INODVAR(JNOD)
DO18I=1,NODDOF
IF(NODVAR(I,J).NE.1)GOTO18
NEQ=NEQ+1
NODVAR(I,J)=NEQ
18CONTINUE
20CONTINUE
DO30JNOD=1,NUMNOD
J=INODVAR(JNOD)
DO28I=1,NODDOF
IF(NODVAR(I,J).GE.-1)GOTO28
N=-NODVAR(I,J)-1
NODVAR(I,J)=NODVAR(I,N)
28CONTINUE
30CONTINUE
C.....OPENANDWRITETHENVFILE
OPEN(8,STATUS='unknown',FILE='',FORM='UNFORMATTED')
WRITE(8)((NODVAR(I,J),I=1,NODDOF),J=1,NUMNOD)
CLOSE(8)
cWRITE(*,*)'NUMNOD=',NUMNOD,'NODDOF=',NODDOF
cWRITE(*,6)((NODVAR(I,J),I=1,NODDOF),J=1,NUMNOD)
C....WRITETHEBOUNDAYCONDITIONFILEBFDACCORDINGTOTHEDISP0FILE
C....OPENDISP0FILE
OPEN(1,FILE='',FORM='UNFORMATTED',STATUS='OLD')
READ
(1)NUMNOD,NODDOF,((U0(I,J),I=1,NODDOF),J=1,NUMNOD)
CLOSE
(1)
C....OPENBFDFILE
OPEN(1,FILE='',FORM='UNFORMATTED',STATUS='unknown')
WRITE
(1)((U0(I,J),I=1,NODDOF),J=1,NUMNOD)
CLOSE
(1)
C......GETTHEINITIALTIMEFROMTIME0FILE
C.......OPENTIME0File
OPEN(1,FILE='',FORM='FORMATTED')
READ(1,*)T0,TMAX,DT
TIME=T0
IT=0
WRITE(*,*)'TMAX,DT,TIME,IT=',TMAX,DT,TIME,IT
CLOSE
(1)
C.......OPENTIMEFile
OPEN(1,FILE='',FORM='UNFORMATTED',STATUS='unknown')
WRITE
(1)TMAX,DT,TIME,IT
CLOSE
(1)
C.......OPENCOORfile
OPEN(1,FILE='',FORM='UNFORMATTED',STATUS='OLD')
READ
(1)NUMNOD,NCOOR,((COOR(I,J),I=1,NCOOR),J=1,NUMNOD)
CLOSE
(1)
cWRITE(*,*)'COOR='
cWRITE(*,7)((COOR(I,J),I=1,NCOOR),J=1,NUMNOD)
C......GETTHEINITIALVALUEFROMTHEDATAFILESBYPREPROCESSOR
inquire(file='disp1',exist=filflg)
if(filflg)then
open(16,file='disp1',form='unformatted',status='old')
read(16)numnod,noddof,((U0(J,N),J=1,NODDOF),N=1,NUMNOD)
close(16)
endif
inquire(file='disp2',exist=filflg)
if(filflg)then
open(16,file='disp2',form='unformatted',status='old')
read(16)numnod,noddof,((U1(J,N),J=1,NODDOF),N=1,NUMNOD)
close(16)
endif
inquire(file='disp3',exist=filflg)
if(filflg)then
open(16,file='disp3',form='unformatted',status='old')
read(16)numnod,noddof,((U2(J,N),J=1,NODDOF),N=1,NUMNOD)
close(16)
endif
cWRITE(*,*)'U0='
cWRITE(*,'(6F13.3)')((U0(J,N),J=1,NODDOF),N=1,NUMNOD)
CWRITE(*,*)'U1='
CWRITE(*,'(6F13.3)')((U1(J,N),J=1,NODDOF),N=1,NUMNOD)
C......COMPUTETHEINITIALVALUEBYBOUND.FOR
zo=0.0d0
cDO321N=1,NUMNOD
cDO100J=1,NCOOR
c100R(J)=COOR(J,N)
cDO200J=1,NODDOF
cU0(J,N)=BOUND(R,zo,J)
cU1(J,N)=BOUND1(R,zo,J)
cU2(J,N)=BOUND2(R,zo,J)
c200CONTINUE
c321CONTINUE
C.......OPENANDWRITETHEINITIALVALUEFILEUNOD
OPEN(1,FILE='',FORM='UNFORMATTED',STATUS='unknown')
WRITE
(1)((U0(I,J),J=1,NUMNOD),I=1,NODDOF),
*((U1(I,J),J=1,NUMNOD),I=1,NODDOF),
*((U2(I,J),J=1,NUMNOD),I=1,NODDOF),
*((U0(I,J),J=1,NUMNOD),I=1,NODDOF)
CLOSE
(1)
c....openIOfile
open(21,file='',form='formatted',status='old')
read(21,'(1a)')material
read(21,*)numtyp
close(21)
DOI=1,NEQ
NUMCOL(i)=1
ENDDO
C.......OPENELEM0file
OPEN(3,FILE='',FORM='UNFORMATTED',STATUS='OLD')
NUMEL=0
KELEM=0
KEMATE=0
DO2000ITYP=1,NUMTYP
C.......INPUTENODE
READ(3)NUM,NNODE,
*((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM)
ccWRITE(*,*)'NUM=',NUM,'NNODE=',NNODE
ccWRITE(*,*)'NODE='
ccWRITE(*,6)((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM)
IF(KELEM.LT.NUM*NNODE)KELEM=NUM*NNODE
NNE=NNODE
IF(MATERIAL.EQ.'Y'.OR.MATERIAL.EQ.'y')THEN
READ(3)MMATE,NMATE
IF(KEMATE.LT.MMATE*NMATE)KEMATE=MMATE*NMATE
NNE=NNE-1
ENDIF
WRITE(*,*)'MMATE=',MMATE,'NMATE=',NMATE
ccWRITE(*,*)'NUM=',NUM,'NNODE=',NNODE
ccWRITE(*,*)'NODE='
ccWRITE(*,6)((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM)
DO1000NE=1,NUM
L=0
DO700INOD=1,NNE
NODI=NODE((NE-1)*NNODE+INOD)
DO600IDGF=1,KDGOF
INV=NODVAR(IDGF,NODI)
IF(INV.LE.0)GOTO600
L=L+1
LM(L)=INV
600CONTINUE
700CONTINUE
NUMEL=NUMEL+1
CWRITE(*,*)'L,LM=',L
CWRITE(*,'(1X,15I5)')(LM(I),I=1,L)
if(l.gt.0)callACLH(NEQ,NUMCOL,l,lm)
1000continue
2000CONTINUE
cCLOSE
(1)
CLOSE(3)
callBCLH(NEQ,NUMCOL)
MAXA=NUMCOL(NEQ)
C.......OPENSYSFile
OPEN(2,FILE='',FORM='UNFORMATTED',STATUS='unknown')
WRITE
(2)NUMEL,NEQ,NUMTYP,MAXA,KELEM,KEMATE
CLOSE
(2)
OPEN(2,FILE='',FORM='UNFORMATTED',STATUS='unknown')
write
(2)(NUMCOL(I),I=1,NEQ)
CLOSE
(2)
cwrite(*,*)'NEQ,NUMCOL=',NEQ
cwrite(*,6)(NUMCOL(i),i=1,NEQ)
END
subroutinechms(kdgof,knode,id)
dimensionid(kdgof,knode),ms(1000),is(1000)
do1000k=1,kdgof
m=0
do800n=1,knode
if(id(k,n).le.-1)id(k,n)=-1
if(id(k,n).le.1)goto800
j=id(k,n)
j0=0
if(m.gt.0)then
doi=1,m
if(j.eq.ms(i))j0=is(i)
enddo
endif
if(j0.eq.0)then
m=m+1
ms(m)=j
is(m)=n
id(k,n)=1
else
id(k,n)=-j0-1
endif
800continue
1000continue
return
end
SUBROUTINEACLH(NEQ,NUMCOL,ND,LM)
implicitreal*8(a-h,o-z)
DIMENSIONLM(ND),NUMCOL(NEQ)
LS=LM
(1)+1
DO100I=1,ND
110IF(LM(I)-LS)120,100,100
120LS=LM(I)
100CONTINUE
DO200I=1,ND
II=LM(I)
ME=II-LS
IF(ME.GT.NUMCOL(II))NUMCOL(II)=ME
200CONTINUE
RETURN
END
SUBROUTINEBCLH(NEQ,NUMCOL)
implicitreal*8(a-h,o-z)
DIMENSIONNUMCOL(NEQ)
CNUMCOL
(1)=1
DO490I=2,NEQ
490NUMCOL(I)=NUMCOL(I)+NUMCOL(I-1)+1
RETURN
END
3.eshell3da.for,Galerkin法求解位移场的主程序
implicitreal*8(a-h,o-z)
character*12fname,filename(20)
common/aa/ia(250000000)
common/bb/ib(125000000)
common/cc/ic(62500000)
open(1,file='',form='unformatted',status='old')
read
(1)knode,kdgof
close
(1)
MAXT=250000000/2/2
C.......OPENSYSFile
OPEN(2,FILE='',FORM='UNFORMATTED',STATUS='OLD')
read
(2)NUMEL,NEQ,NUMTYP,MAXA,KELEM,KEMATE
CLOSE
(2)
IF(MAXA.GT.MAXT)THEN
WRITE(*,*)'MATRIXAEXCEEDCOREMEMERY....',MAXA
WRITE(*,*)'REQUIREDCOREMEMERY...........',MAXT
STOP0000
ENDIF
KVAR=KNODE*KDGOF
KCOOR=3
CKELEM=31250000
WRITE(*,*)'KNODE,KDGOF,KVAR,KCOOR,KELEM='
WRITE(*,'(1X,6I7)')KNODE,KDGOF,KVAR,KCOOR,KELEM
kna1=kdgof*knode*1
if(kna1/2*2.lt.kna1)kna1=kna1+1
knc1=kdgof*knode*2
knc2=kcoor*knode*2
knc7=kdgof*knode*2
knc3=neq*2
knb1=maxa*2
knb2=maxa*2
kna2=neq*1
if(kna2/2*2.lt.kna2)kna2=kna2+1
knc6=kemate*2
kna3=kelem*1
if(kna3/2*2.lt.kna3)kna3=kna3+1
knc8=100000*2
knc5=neq*2
knc4=kdgof*knode*2
kna0=1
kna1=kna1+kna0
kna2=kna2+kna1
kna3=kna3+kna2
if(kna3-1.gt.125000000)then
write(*,*)'exceedmemoryofarrayib'
write(*,*)'memoryofib=125000000'
write(*,*)'memoryneeded=',kna3,'inprgrameshell3da'
stop55555
endif
knb0=1
knb1=knb1+knb0
knb2=knb2+knb1
if(knb2-1.gt.250000000)then
write(*,*)'exc