[Fortran] 纯文本查看 复制代码
subroutine Front
dimension fixed(160),equat(60),vecrv(160),gload(60),gstif(1830),estif(16,16),&
iffix(160),nacva(60),locel(16),ndest(16)
common/Contro/npoin,nelem,nnode,ndofn,ndime,nstre,ntype,ngaus,nprop,nmats,nvfix,nevab,&
icase,ncase,itemp,iprob,nprob
common/Lgdata/coord(80,2),props(10,5),presc(40,2),asdis(160),eload(25,16),strin(3,225),&
nofix(40),ifpre(40,2),lnods(25,8),matno(25)
common/Work/elcod(2,8),shape(8),deriv(2,8),dmatx(3,3),cartd(2,8),dbmat(3,16),bmatx(3,16),&
smatx(3,16,9),posgp(3),weigp(3),gpcod(2,9),neror(24)
!integer nfunc
call Openfile
nfunc(i,j)=(j*j-j)/2+i !*********************************问题在此行
mfron=60
mstif=1830
!
! interpret fixity data in vector form
!
ntotv=npoin*ndofn
do 100 itotv=1,ntotv
iffix(itotv)=0
100 fixed(itotv)=0.0
do 110 idofn=1,ndofn
ngash=nloca+idofn
iffix(ngash)=ifpre(ivfix,idofn)
110 fixed(ngash)=presc(ivfix,jdofn)
!
! change the sign of the last appearance
! of each node
!
do 140 ipoin=1,npoin
klast=0
do 130 ielem=1,nelem
do 120 inode=1,nnode
if(lnods(ielem,inode).ne.ipoin) goto 120
klast=ielem
nlast=inode
120 continue
130 continue
if(klast.ne.0) lnods(klast,nlast)=-ipoin
140 continue
!
! start by initializing everything that
! matters to zero
!
do 150 istif=1,mstif
150 gstif(istif)=0.0
do 160 ifron=1,mfron
gload(ifron)=0.0
equat(ifron)=0.0
vecrv(ifron)=0.0
160 nacva(ifron)=0
!
! and prepare for disc reading and writing
! operations
!
rewind 2
rewind 4
rewind 6
rewind 8
!
! enter main element assembly-reduction loop
!
nfron=0
kelva=0
do 380 ielem=1,nelem
kevab=0
read(2,*) estif
do 170 inode=1,nnode
do 170 idofn=1,ndofn
nposi=(inode-1)*ndofn+idofn
locno=lnods(ielem,inode)
if(locno.gt.0) locel(nposi)=(locno-1)*ndofn-idofn
if(locno.lt.0) locel(nposi)=(locno+1)*ndofn-idofn
170 continue
!
! start by looking for existing destinations
!
do 210 ievab=1,nevab
nikno=iabs(locel(ievab))
kexis=0
do 180 ifron=1,nfron
if(nikno.ne.nacva(ifron)) goto 180
kevab=kevab+1
kexis=1
ndest(kevab)=ifron
180 continue
if(kexis.ne.0) goto 210
!
! we now seek new empty places for
! destination vector
!
do 190 ifron=1,mfron
if(nacva(ifron).ne.0) goto 190
nacva(ifron)=nikno
kevab=kevab+1
ndest(kevab)=ifron
190 continue
!
! the new places may demand an increase
! in current frontwidth
!
200 if(nedst(kevab).gt.nfron) nfron=ndest(kevab)
210 continue
!
! assemble element loads
!
do 240 ievab=1,nevab
idest=ndest(ievab)
gload(idest)=gload(idest)+eload(ielem,ievab)
!
! assemble the element stifnesses
! but not in resolution
!
if(icase.gt.1) goto 230
do 220 jevab=1,ievab
jdest=ndest(jevab)
ngash=nfunc(idest,jdest)
ngish=nfunc(jdest,idest)
if(jdest.ge.idest) gstif(ngash)=gstif(ngash)+estif(ievab,jevab)
if(jdest.lt.idest) gstif(ngish)=gstif(ngish)+estif(ievab,jevab)
220 continue
230 continue
240 continue
!
! re-examine each element node, to
! enquire which can be eliminated
!
do 370 ievab=1,nevab
nikno=-locel(ievab)
if(nikno.le.0) goto 370
!
! find positions of variables ready
! for elimination
!
do 350 ifron=1,nfron
if(nacva(ifron).ne.nikno) goto 350
!
! extract the cofficients of the
! new equation for elimination
!
if(icase.gt.1) goto 260
do 250 jfron=1,mfron
if(ifron.lt.jfron) nloca=nfunc(ifron,jfron)
if(ifron.ge.jfron) nloca=nfunc(jfron,ifron)
equat(jfron)=gstif(nloca)
250 gstif(nloca)=0.0
260 continue
!
! and fxtract the corresponding right
! hand sides
!
eqrhs=gload(ifron)
gload(ifron)=0.0
kelva=kelva+1
!
! write equations to disc or to tape
!
if(icase.gt.1) goto 270
write(4,*) equat,eqrhs,ifron,nikno
goto 280
270 write(8,*) eqrhs
read(4,*) wquat,dummy,idumm,nikno
280 continue
!
! deal with pivot
!
pivot=equat(ifron)
equat(ifron)=0.0
!
! enquire whether present varizable is
! free or prescribed
!
if(iffix(nikno).eq.0) goto 300
!
! deal with a prescribed deflection
!
do 290 jfron=1,nfron
290 gload(jfron)=gload(jfron)-fixed(nikno)*equat(jfron)
goto 340
!
! eliminate a free variable-deal with
! the right hand side first
!
300 do 330 jfron=1,nfron
gload(jfron)=gload(jfron)-equat(jfron)*eqrhs/pivot
!
! now deal with the coefficients in core
!
if(icase.gt.1) goto 320
if(equat(jfron).eq.0.0) goto 330
nloca=nfunc(0,jfron)
do 310 lfron=1,jfron
ngash=lfron+nloca
310 gstif(ngash)=gstif(ngash)-equat(jfron)*equat(lfron)/pivot
320 continue
330 continue
340 equat(ifron)=pivot
!
! record the new vacant space, and reduce
! frontwidth if possible
!
nacva(ifron)=0
goto 360
!
! complete the element loop in the forward
! elimination
!
350 continue
360 if(nacva(nfron).ne.0) goto 370
nfron=nfron-1
if(nfron.gt.0) goto 360
370 continue
380 continue
!
! enter back-substitution phase, loop
! backwards through variables
!
do 410 ielva=1,kelva
!
! read a new equation
!
backspace 4
read(4,*) equat,eqrhs,ifron,nikno
backspace 4
if(icase.eq.1) goto 390
backspace 8
read(8,*) eqrhs
backspace 8
390 continue
!
! prepare to back-substitute from the
! current equation
!
pivot=equat(ifron)
if(iffix(nikno).eq.1) vecrv(ifron)=fixed(nikno)
if(iffix(nikno).eq.0) equat(ifron)=0.0
!
! back-substitute in the current equation
!
do 400 jfron=1,mfron
400 eqrhs=eqrhs-vecrv(jfron)*equat(jfron)
!
! put the final values where they belong
!
if(iffix(nikno).eq.0) vecrv(ifron)=eqrhs/pivot
if(iffix(nikno).eq.1) fixed(nikno)=-eqrhs
asdis(nikno)=vecrv(ifron)
410 continue
write(12,900)
900 format(1ho,5x,13hDisplacements)
if(ndofn.ne.2) goto 430
if(ndime.ne.1) goto 420
write(12,905)
905 format(1ho,5x,4hnode,6x,5hDisp.,7x,8hRotation)
goto 440
420 write(12,910)
910 format(1ho,5x,4hNode,5x,7hX-disp.,7x,7hY-disp.)
430 write(12,915)
915 format(1ho,5x,4hNode,5x,7hX-disp.,7x,7hY-disp.)
440 continue
do 450 ipoin=1,npoin
ngash=ipoin*ndofn
ngish=ngash-ndofn+1
450 write(12,920) ipoin,(asdis(igash),igash=ngish,ngash)
920 format(i10,3e14.6)
write(12,925)
925 format(1ho,5x,9hReactions)
if(ndofn.ne.2) goto 470
if(ndime.ne.1) goto 460
write(12,930)
930 format(1ho,5x,4hNode,6x,5hForce,8x,6hMomnet)
goto 480
460 write(12,935)
935 format(1ho,5x,4hNode,5x,7hX-force,7x,7hY-froce)
goto 480
470 write(12,940)
940 format(1ho ,5x,4hNode,6x,5hForce,6x,9hXz-moment,5x,9hYz-moment)
480 continue
do 510 ipoin=1,npoin
nloca=(ipoin-1)*ndofn
do 490 idofn=1,ndofn
ngush=nloca+idofn
if(iffix(ngush).gt.0) goto 500
490 continue
goto 510
500 ngash=nloca+ndofn
ngish=nloca+1
write(12,945) ipoin,(fixed(igash),igash=ngish,ngash)
510 continue
945 format(i10,3e14.6)
!
! post front-reset all element connection numbers to positive
! values for subsequent use in stress calculation
!
do 520 ielem=1,nelem
do 520 inode=1,nnode
520 lnods(ielem,inode)=iabs(lnods(ielem,inode))
return
end