八节点四边形等参单元有限元程序

C2 P114-116
C
C CALCUATES THE STIFFNESS MATRIX OF ELEMENT
C
C-----------------------------------------------------------------C
COMMON/CONTR/NPOIN,NELEM,NNODE,NDOFN,NDIME,NGAUS,NPROP,NMATS
KGASP=0
C
C ENTER NUMERICAL INTEGRATION
C
DO 50 IGAUS=1,NGAUS
DO 50 JGAUS=1,NGAUS
KGASP=KGASP+1
EXISP=POSGP(IGAUS)
ETASP=POSGP(JGAUS)
C
CALL SFR2(EXISP,ETASP)
CALL JACOB2(IELEM,DJACB,KGASP)
CALL LOADPS
CALL ASSEMB
C
CALL GAUSS(ASTIF,ASLOD,NTOTV)
DO 175 I=1,NTOTV
175 ASDIS(I)=ASLOD(I)
WRITE(6,180)
180 FORMAT(/1X,'NODE',5X,'X-DISP',5X,'Y-DISP')
WRITE(6,185) (I,ASDIS(2*I-1),ASDIS(2*I),I=1,NPOIN)
COMMON/WORKS/ELCOD(2,8),SHAPE(8),DERIV(2,8),CARTD(2,8)
COMMON/WORKS/POSGP(3),WEIGP(3),GPCOD(2,9),BMATX(3,16),DMATX(3,3)
COMMON/WORKS/SMATX(3,16,9),DBMAT(3,16)
COMMON/CONTR/NVFIX,NEVAB,NSTRE,NTYPE,NTOTV
COMMON/LGDAT/COORD(40,2),PROPS(9,4),PRESC(40,2),ASDIS(80)
COMMON/LGDAT/ELOAD(9,16),NOFIX(40),IFPRE(40,2),LNODS(9,8),MATNO(9)
1 'NMATS='I3,1X,'NGAUS='I2,1X,'NTYPE=',I2)
110 FORMAT(/1X,'ELEMENT',3X,'PROPERTY',6X,'NODE NUMBER')
115 FORMAT(1X,I5,I9,6X,8I4)
120 FORMAT(/24H NODAL POINT COORDINATES)
DVOLU=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS)
IF (THICK.NE.0.0) DVOLU=DVOLU*THICK
CALL BMATPS
CALL DBE
C
3
C CALCULATE THE ELEMENT STIFFNESSES
C
DO 30 IEVAB=1,NEVAB
DO 30 JEVAB=IEVAB,NEVAB
155 FORMAT(1X,I4,5X,2I1,2F10.5)
160 FORMAT(/1X,'MATERAL PROPERTIES')
165 FORMAT(/1X,'NUMBER',7X,'PROPERTIES')
170 FORMAT(1X,I4,4X,4E14.4)
C
CALL GAUSSQ
CALL STIFPS
40 SMATX(ISTRE,IEVAB,KGASP)=DBMAT(ISTRE,IEVAB)
50 CONTINUE
C
C CONSTRUCT THE LOWER TRIAGLE OF K
C
DO 60 IEVAB=1,NEVAB
DO 60 JEVAB=1,NEVAB
60 ESTIF(JEVAB,IEVAB)=ESTIF(IEVAB,JEVAB)
IF(NGAUS.GT.2) GO TO 10
POSGP(1)=-0.577350269189626
WEIGP(1)=1.0
GO TO 20
10 POSGP(1)=-0.774596669241483
POSGP(2)=0.0
WEIGP(1)=0.5555555555555556
WEIGP(2)=0.8888888888888889
COMMON/WORKS/ELCOD(2,8),SHAPE(8),DERIV(2,8),CARTD(2,8)
COMMON/WORKS/POSGP(3),WEIGP(3),GPCOD(2,9),BMATX(3,16),DMATX(3,3)
COMMON/WORKS/SMATX(3,16,9),DBMAT(3,16)
C C-----------PARAMETRIC FEM FOR PLANE PROBLEMS----------------C
COMMON/CONTR/NPOIN,NELEM,NNODE,NDOFN,NDIME,NGAUS,NPROP,NMATS COMMON/CONTR/NVFIX,NEVAB,NSTRE,NTYPE,NTOTV COMMON/LGDAT/COORD(40,2),PROPS(9,4),PRESC(40,2),ASDIS(80) COMMON/LGDAT/ELOAD(9,16),NOFIX(40),IFPRE(40,2),LNODS(9,8),MATNO(9) COMMON/WORKS/ELCOD(2,8),SHAPE(8),DERIV(2,8),CARTD(2,8) COMMON/WORKS/POSGP(3),WEIGP(3),GPCOD(2,9),BMATX(3,16),DMATX(3,3) COMMON/WORKS/SMATX(3,16,9),DBMAT(3,16) COMMON/GENEL/ASLOD(80),ASTIF(80,80) CHARACTER *12 NAME C P231-232 WRITE(*,210) 210 FORMAT(/3X,'PLEASE ENTER DATAFILE NAME:') READ(*,220) NAME 220 FORMAT(A12) REWIND 1 REWIND 3 OPEN(5,FILE=NAME,status='old') OPEN(1,FILE='STIF.dat',form='unformatted') OPEN(3,FILE='SMAT.dat',form='unformatted') OPEN(6,FILE='OUT1.txt') C NNODE=8 NDOFN=2 NEVAB=NNODE*NDOFN NDIME=2 NSTRE=3 NPROP=4 READ(5,*)NPOIN,NELEM,NVFIX,NMATS,NGAUS,NTYPE NGASP=NGAUS*NGAUS NTOTV=NPOIN*NDOFN WRITE(6,100)NPOIN,NELEM,NVFIX,NMATS,NGAUS,NTYPE WRITE(6,110) DO 10 LELEM=1,NELEM READ(5,*) IELEM,MATNO(IELEM),(LNODS(IELEM,INODE),INODE=1,NNODE) WRITE(6,115) IELEM,MATNO(IELEM),(LNODS(IELEM,INODE),INODE=1,NNODE) 10 CONTINUE WRITE(6,120) WRITE(6,125) DO 20 IPOIN=1,NPOIN 20 READ(5,*) JPOIN,(COORD(JPOIN,IDIME),IDIME=1,NDIME) c WRITE(6,135)(I,(COORD(I,IDIME),IDIME=1,NDIME),I=1,NPOIN) WRITE(6,140) WRITE(6,145) DO 80 IVFIX=1,NVFIX READ(5,*) NOFIX(IVFIX),(IFPRE(IVFIX,IDOFN),IDOFN=1,NDOFN), 1(PRESC(IVFIX,IDOFN),IDOFN=1,NDOFN) WRITE(6,155) NOFIX(IVFIX),(IFPRE(IVFIX,IDOFN),IDOFN=1,NDOFN), 1(PRESC(IVFIX,IDOFN),IDOFN=1,NDOFN) 80 CONTINUE WRITE(6,160) WRITE(6,165)
185 FORMAT(1X,I4,2E12.4,4X,I2,2E12.4)
C
CALL STREPS
STOP
END
C
SUBROUTINE GAUSSQ
C----------------------------------------------------------------C
C1 P65
C
C GIVES THE WEIGHTING FACTORS FOR GAUSS POINTS
DO 30 ISTRE=1,NSTRE
30 ESTIF(IEVAB,JEVAB)=ESTIF(IEVAB,JEVAB)+
1 BMATX(ISTRE,IEVAB)*DBMAT(ISTRE,JEVAB)*DVOLU
C
C STORE DB MATRIX
C
DO 40 ISTRE=1,NSTRE
DO 40 IEVAB=1,NEVAB
COMMON/GENEL/ASLOD(80),ASTIF(80,80)
DIMENSION ESTIF(16,16)
C
DO 70 IELEM=1,NELEM
LPROP=MATNO(IELEM)
C
C EVALUATE THE COORDINATES
C
DO 10 INODE=1,NNODE
LNODE=LNODS(IELEM,INODE)
C
WRITE(1) ESTIF
WRITE(3) SMATX,GPCOD
70 CONTINUE
RETURN
END
C
SUBROUTINE MODPS(LPROP)
C----------------------------------------------------------------C
C3 P111
20 KGAUS=NGAUS/2
DO 30 IGASH=1,KGAUS
2
JGASH=NGAUS+1-IGASH
POSGP(JGASH)=-POSGP(IGASH)
WEIGP(JGASH)=WEIGP(IGASH)
30 CONTINUE
RETURN
END
C
SUBROUTINE STIFPS
C-----------------------------------------------------------------C
1
DO 90 IMATS=1,NMATS
READ(5,*) NUMAT,(PROPS(NUMAT,IPROP),IPROP=1,NPROP)
WRITE(6,170) NUMAT,(PROPS(NUMAT,IPROP),IPROP=1,NPROP)
90 CONTINUE
100 FORMAT(1X,'NPOIN=',I3,1X,'NELEM=',I3,1X,'NVFIX=',I3,1X,
DO 10 IDIME=1,NDIME
10 ELCOD(IDIME,INODE)=COORD(LNODE,IDIME)
C
CALL MODPS(LPROP)
THICK=PROPS(LPROP,3)
C
DO 20 IEVAB=1,NEVAB
DO 20 JEVAB=1,NEVAB
20 ESTIF(IEVAB,JEVAB)=0.0
125 FORMAT(/2(7H NODE ,7X,1HX,9X,1HY,5X))
135 FORMAT(2(1X,I4,2X,2F10.3,2X))
140 FORMAT(/16HRESTRAINED NODES)
145 FORMAT(/1X,'NODE',4X,'CODE',6X,'FIXED VALUES')
C
C----------------------------------------------------------------C
COMMON/CONTR/NPOIN,NELEM,NNODE,NDOFN,NDIME,NGAUS,NPROP,NMATS
COMMON/CONTR/NVFIX,NEVAB,NSTRE,NTYPE,NTOTV
Cwenku.baidu.com