主题：小女急需求助，麻烦各位帮忙看看是什么问题
lynn0976
[专家分：0] 发布于 2012-08-31 12:09:00
程序算出的雅克比为0，但其实我人工计算发现根本不为0，查找错误时发现程序中的子程序SFR2和JACOB2在进行到第四次循环时，有些变量值在输出文件中显示为空白，不能输出值，导致这些量与其他量相乘后，雅克比结果为0.
这是怎么回事，以前没遇到过这种情况。部分主程序和子程序贴在下面：
沙发
lynn0976 [专家分：0] 发布于 2012-08-31 12:10:00
DO 170 IPOIN=1,NPOIN
170 TEMPE(IPOIN)=0.0
WRITE(6,917)
917 FORMAT(1H0,5X,29HPRESCRIBED NODAL TEMPERATURES)
180 READ (5,*) NODPT,TEMPE(NODPT)
WRITE(6,916) NODPT,TEMPE(NODPT)
916 FORMAT(I5,F10.3)
IF(NODPT.LT.NPOIN) GO TO 180
KGAST=0
! MDIME=3
!
!*** LOOP OVER EACH ELEMENT
!
DO 280 IELEM=1,NELEM !开始对单元进行循环
LPROP=MATNO(IELEM)
DO 200 INODE=1,NNODE
LNODE=IABS(LNODS(IELEM,INODE))
!
!*** IDENTIFY THE COORDINATES AND TEMPERATURE OF EACH ELEMENT NODE POINT
!
DO 190 IDIME=1,2
190 ELCOD(IDIME,INODE)=COORD(LNODE,IDIME)
200 FLCOD(INODE)=TEMPE(LNODE)
!
!*** SET UP MATERIAL PROPERTIES
!
CALL MODPS (DMATX,LPROP,MMATS,NTYPE,PROPS) !计算弹性矩阵D
YOUNG=PROPS(LPROP,1)
POISS=PROPS(LPROP,2)
THICK=PROPS(LPROP,3)
ALPHA=PROPS(LPROP,8)
!
!*** ENTER LOOPS FOR AREA NUMERICAL INTEGRATION
!
!
KGASP=0
DO 270 IGAUS=1,NGAUS
EXISP=POSGP(IGAUS)
DO 270 JGAUS=1,NGAUS !进入在单元面积上高斯积分的循环
KGAST=KGAST+1 !对KGAST增量，以便给出此时高斯点的值
KGASP=KGASP+1 !
ETASP=POSGP(JGAUS) !形成高斯取样点的坐标ξ、 η
!
!*** EVALUATE THE SHAPE FUNCTIONS AND TDIPERATURE AT THE SAMPLING POINTS
! ,ELEPIENTAL VOLUME AND CARTESIAN DERIVATIVES
!
[color=FFFF00]CALL SFR2 [/color](DERIV,ETASP,EXISP,NNODE,SHAPE) !计算形函数及其导数
[color=FFFF00]CALL JACOB2 [/color](CARTD,DERIV,DJACB,ELCOD,GPCOD,IELEM,KGASP,NNODE,SHAPE) !计算雅可比矩阵及其行列式和逆、形函数导数
THERM=0.0
DO 210 INODE=1,NNODE
210 THERM=THERM+FLCOD(INODE)*SHAPE(INODE) !计算高斯点的温度值THERM
DVOLU=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS)
IF(NTYPE.EQ.1) DVOLU=DVOLU*THICK
IF(NTYPE.EQ.3) DVOLU=DVOLU*TWOPI*GPCOD(1,KGASP)
!
!*** EVALUATE THE INITIAL THERMAL STRAINS
!计算温度引起的高斯点的初始应变
EIGEN=THERM*ALPHA
IF(NTYPE.EQ.2) GO To 220
STRAN(1)=-EIGEN
STRAN(2)=-EIGEN
STRAN(3)=0.0
GO To 230
220 STRAN(1)=-(1.0+POISS)*EIGEN
STRAN(2)=-(1.0+POISS)*EIGEN
STRAN(3)=0.0
!
!*** AND THE CORRESPONDING INITJAL STRESSES
!
230 Do 250 ISTRE=1,NSTRE
STRES(ISTRE)=0.0
Do 240 JSTRE=1,NSTRE
240 STRES(ISTRE)=STRES(ISTRE)+DMATX(ISTRE,JSTRE)*STRAN(JSTRE) !计算相对应的初始应力
STRSG(ISTRE,KGAST)=STRES(ISTRE)
250 STRIN(ISTRE,KGAST)=STRES(ISTRE)
IF(NTYPE.EQ.2) STRIN(4,KGAST)=-YOUNG*EIGEN
IF(NTYPE.EQ.1) STRIN(4,KGAST)=0.0
!
!*** CALCULATE THE EQUIVALENT NODAL FORCES AND ASSOCIATE WITH THE
! ELEMENT NODES
!计算等效节点力并存储在RLOAD数组中的正确位置（p式7.21）
EXTRA=0.0
Do 260 INODE=1,NNODE
IF(NTYPE.EQ.3) EXTRA=DVOLU*SHAPE(INODE)*STRES(4)/GPCOD(1,KGASP)
NGASH=(INODE-1)*NDOFN+1
MGASH=(INODE-1)*NDOFN+2
RLOAD(IELEM,NGASH)=RLOAD(IELEM,NGASH)+EXTRA&
-(CARTD(1,INODE)*STRES(1)+CARTD(2,INODE)*STRES(3))*DVOLU
260 RLOAD(IELEM,MGASH)=RLOAD(IELEM,MGASH)&
-(CARTD(1,INODE)*STRES(3)+CARTD(2,INODE)*STRES(2))*DVOLU
270 CONTINUE !一个单元内的高斯点循环结束
280 CONTINUE !所有单元循环结束
800 CONTINUE !温度载荷部分结束
WRITE(6,907)
907 FORMAT(1H ,5X,36H TOTAL NODAL FORCES FOR EACH ELEMENT)
DO 290 IELEM=1,NELEM
290 WRITE(6,905) IELEM,(RLOAD(IELEM,IEVAB),IEVAB=1,NEVAB)
905 FORMAT(1X,I4,5X,8E12.4/(10X,8E12.4))
RETURN
END
板凳
lynn0976 [专家分：0] 发布于 2012-08-31 12:11:00
[color=FFFF00]SUBROUTINE SFR2(DERIV,ETASP,EXISP,NNODE,SHAPE[/color])
!*************************************************************************
!
!**** THIS SUBROUTINE EVALUATES SHAPE FUNCTIONS AND THEIR DERIVATIVES
! FOR LINEAR,QUADRATIC LAGRAGIAN AND SERENDIPITY
! ISOPARAMETRIC 2-D ELEMENTS
!
!**************************************************************************
DIMENSION DERIV(2,9),SHAPE(9)
S=EXISP
T=ETASP
WRITE(6,*) ETASP,EXISP
ST=S*T
S2=S*2.0
T2=T*2.0
SS=S*S
TT=T*T
ST=S*T
SST=S*S*T
STT=S*T*T
ST2=S*T*2.0
!
!*** SHAPE FUNCTIONS FOR 8 NODED ELEMENT
!
SHAPE(1)=(-1.0+ST+SS+TT-SST-STT)/4.0
SHAPE(2)=(1.0-T-SS+SST)/2.0
SHAPE(3)=(-1.0-ST+SS+TT-SST+STT)/4.0
SHAPE(4)=(1.0+S-TT-STT)/2.0
SHAPE(5)=(-1.0+ST+SS+TT+SST+STT)/4.0
SHAPE(6)=(1.0+T-SS-SST)/2.0
SHAPE(7)=(-1.0-ST+SS+TT+SST-STT)/4.0
SHAPE(8)=(1.0-S-TT+STT)/2.0
!
!*** SHAPE FUNCTION DERIVATIVES
!
DERIV(1,1)=(T+S2-ST2-TT)/4.0
DERIV(1,2)=-S+ST
DERIV(1,3)=(-T+S2-ST2+TT)/4.0
DERIV(1,4)=(1.0-TT)/2.0
DERIV(1,5)=(T+S2+ST2+TT)/4.0
DERIV(1,6)=-S-ST
DERIV(1,7)=(-T+S2+ST2-TT)/4.0
DERIV(1,8)=(-1.0+TT)/2.0
DERIV(2,1)=(S+T2-SS-ST2)/4.0
DERIV(2,2)=(-1.0+SS)/2.0
DERIV(2,3)=(-S+T2-SS+ST2)/4.0
DERIV(2,4)=-T-ST
DERIV(2,5)=(S+T2+SS+ST2)/4.0
DERIV(2,6)=(1.0-SS)/2.0
DERIV(2,7)=(-S+T2+SS-ST2)/4.0
DERIV(2,8)=-T+ST
! WRITE(6,*) ((DERIV(IDIME,INODE),IDIME=1,2),INODE=1,NNODE) [color=FFFF00] 循环到第四次时此处输出DERIV是一行空白[/color]
RETURN
END
!
[color=FFFF00] [color=000000]SUBROUTINE JACOB2(CARTD,DERIV,DJACB,ELCOD,GPCOD,IELEM,KGASP,&
NNODE,SHAPE）
!***********************************************************************
!
!**** THIS SUBROUTINE EVALUATES THE JACOBIAN MATRIX AND CARTESIAN
! SHAPE FUNCTION DERIVATIVES
!
!***********************************************************************
DIMENSION CARTD(2,9),DERIV(2,9),ELCOD(2,9),GPCOD(2,9),SHAPE(9),&
XJACI(2,2),XJACM(2,2)
!
!*** CALCULATE COORDINATES OF SAMPLING POINT
!
! WRITE(6,*) "ELCOD"
! WRITE(6,*) ((ELCOD(IDIME,INODE),IDIME=1,2),INODE=1,NNODE) [color=FFFF00]循环到第四次时此处输出ELCOD也是一行空白 [/color]
DO 2 IDIME=1,2
GPCOD(IDIME,KGASP)=0.0
DO 2 INODE=1,NNODE
GPCOD(IDIME,KGASP)=GPCOD(IDIME,KGASP)+ELCOD(IDIME,INODE)&
*SHAPE(INODE)
2 CONTINUE
!
!*** CREATE JACOBIAN MATRIX XJACM
!
DO 4 IDIME=1,2
DO 4 JDIME=1,2
XJACM(IDIME,JDIME)=0.0
DO 4 INODE=1,NNODE
XJACM(IDIME,JDIME)=XJACM(IDIME,JDIME)+DERIV(IDIME,INODE)*&
ELCOD(JDIME,INODE)
4 CONTINUE
!
!*** CALCULATE DETERMINANT AND INVERSE OF JACOBIAN MATRIX
!
DJACB=ABS(XJACM(1,1)*XJACM(2,2)-XJACM(1,2)*XJACM(2,1))
! WRITE(6,*) DJACB,(GPCOD(IDIME,KGASP),IDIME=1,2),&
! SHAPE(1),SHAPE(2),SHAPE(3),SHAPE(4),SHAPE(5),&
! SHAPE(6),SHAPE(7),SHAPE(8),SHAPE(9),XJACM(1,1),&
! XJACM(1,2),XJACM(2,1),XJACM(2,2)
! WRITE(6,*) ((DERIV(IDIME,INODE),IDIME=1,2),INODE=1,NNODE)
IF(DJACB.GT.0.0) GO TO 8
6 WRITE(6,600) IELEM
STOP
8 CONTINUE
XJACI(1,1)=XJACM(2,2)/DJACB
XJACI(2,2)=XJACM(1,1)/DJACB
XJACI(1,2)=-XJACM(1,2)/DJACB
XJACI(2,1)=-XJACM(2,1)/DJACB
!
!*** CALCULATE CARTESIAN DERIVATIVES
!
DO 10 IDIME=1,2
DO 10 INODE=1,NNODE
CARTD(IDIME,INODE)=0.0
DO 10 JDIME=1,2
CARTD(IDIME,INODE)=CARTD(IDIME,INODE)+XJACI(IDIME,JDIME)*&
DERIV(JDIME,INODE)
10 CONTINUE
600 FORMAT(//,36H PROGRAM HALTED IN SUBROUTINE JACOB2,/,11X,&
22H ZERO OR NEGATIVE AREA,/,10X,16H ELEMENT NUMBER ,I5)
RETURN
END
您所在位置：社区首页 — Fortran讨论区 — 小女急需求助，麻烦各位帮忙看看是什么问题
idealguy [专家分：110] 发布于 2012-09-02 20:49:00