Fortran Coder

查看: 5642|回复: 8

[子程序] 求助,子程序中一个自定义函数的问题

[复制链接]

13

帖子

4

主题

0

精华

入门

F 币
74 元
贡献
43 点
发表于 2015-12-19 14:44:51 | 显示全部楼层 |阅读模式
[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

954

帖子

0

主题

0

精华

大师

F 币
184 元
贡献
75 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
发表于 2015-12-19 16:46:46 | 显示全部楼层
您有何问题?请给出文字描述。

1957

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1335 元
贡献
563 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2015-12-19 16:57:18 | 显示全部楼层
nfunc(i,j)=(j*j-j)/2+i            !*********************************问题在此行

移动到 call Openfile 前面

13

帖子

4

主题

0

精华

入门

F 币
74 元
贡献
43 点
 楼主| 发表于 2015-12-21 20:38:35 | 显示全部楼层
移了,但还是有问题
Error        1         error LNK2019: unresolved external symbol _NEDST referenced in function _FRONT        s-front.obj       
Error        2         fatal error LNK1120: 1 unresolved externals        Debug\master-plane.exe       

1957

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1335 元
贡献
563 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2015-12-22 08:46:10 | 显示全部楼层
第 107 行
200 if(nedst(kevab).gt.nfron) nfron=ndest(kevab)
是否应该是
200 if(ndest(kevab).gt.nfron) nfron=ndest(kevab)

13

帖子

4

主题

0

精华

入门

F 币
74 元
贡献
43 点
 楼主| 发表于 2015-12-22 11:35:09 | 显示全部楼层
是的,找出错误了,太感谢了!

100

帖子

0

主题

0

精华

专家

F 币
550 元
贡献
291 点

规矩勋章元老勋章

QQ
发表于 2015-12-23 23:28:17 | 显示全部楼层
implicit none 你值得拥有。

13

帖子

4

主题

0

精华

入门

F 币
74 元
贡献
43 点
 楼主| 发表于 2015-12-25 10:57:55 | 显示全部楼层
fcode 发表于 2015-12-22 08:46
第 107 行
200 if(nedst(kevab).gt.nfron) nfron=ndest(kevab)
是否应该是

大侠,编译通过了,可是,运行过程中当call Front之前,整形数组nofix(40),ifpre(40,2),lnods(25,8),matno(25)是有内容的,可是进入Front后这些数组全部回到了0值,其他非整形的数组没有问题。不知道是什么原因???(openfile中没有用implicit none)

1957

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1335 元
贡献
563 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2015-12-25 11:26:18 | 显示全部楼层
不知道,运行以后的问题属于动态问题,需要动态分析。(调试等)

而只有代码的一部分,无法进行动态分析。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-4-18 16:09

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表