fortran编的有限元计算程序

tsy

自己编三节点三角形单元计算程序
供fortran初学者参考参考：）
ps：fortran还是可以用动态数组，还可以是全局的
刚在前面看到好象有人讨论说fortran不行，呵呵

才完成80%，忙其他的去了，就丢下了
需完善的地方：
1、均布荷载的移置计算可能有问题，得先推推计算公式
2、没有考虑个地震影响系数
3、增加对平面应变问题的处理
好玩编的，也没有啥意义了，呵呵
只有等啥时候心血来潮了写着玩玩：）

MODULE DYNA
CHARACTER:UES*70
!     QUES---------问题名称
INTEGER::NELEM, NNODE, NFIX, NP, NMATERIAL, MFREE, NSP
! NELEM--------单元总数
! NNODE--------节点总数
! NFIX --------约束节点总数
! NP   --------荷载作用数
! NMATERIAL----材料类型数
! MFREE--------自由度数
INTEGER, DIMENSION(, ALLOCATABLE,SAVE::NFREE, NFIXX, NFIXY
! NFREE()------自由度对应的节点信息
! NFIXX(I)-----I号节点X方向约束信息；NFIXY(I)---I号节点Y方向约束信息
REAL, DIMENSION(:), ALLOCATABLE,SAVE:X, PY, P, BL
! PX(I)--------I节点X方向荷载;  PY(I)---I节点Y方向荷载
! P()    ------自由度荷载列阵;  BL(I)---位移计算结果
REAL, DIMENSION(:,:), ALLOCATABLE,SAVE::MELEM, MNODE, MMATERIAL, D
! MELEM()   -------单元信息数组
! MELEM(I,1)-------I号单元的I节点；MELEM(I,2)---I号单元的J节点；
! MELEM(I,3)-------I号单元的M接点；MELEM(I,4)---I号单元的材料类型
! MNODE()坐标信息
! MNODE(I,1)------I号节点X坐标；MNODE(I,2)---I号节点Y坐标
! MMATERIAL()-----材料信息
! MMATERIAL(I,1)---材料弹模；MMATERIAL(I,2)---材料横截面积
! D(I,J)        ---整体刚度矩阵;
REAL, DIMENSION(:), ALLOCATABLE, SAVE::UU,VV
!UU()-------节点X向位移； VV()--------节点Y向位移
REAL, DIMENSION(:,:), ALLOCATABLE, SAVE::STRESS
REAL, DIMENSION(:,:), ALLOCATABLE, SAVE::SPXY
END MODULE DYNA

PROGRAM MAIN
USE DYNA
CHARACTER::NAME*8, FILE_IN*12, FILE_OUT*12
WRITE(*,'(8/1X,A)')
WRITE(*,100)
WRITE(*, *)
WRITE(*,200)
WRITE(*, *)
WRITE(*,300)
WRITE(*, *)
WRITE(*,400)
100 FORMAT('                     **大学*******专业       ')
200 FORMAT('                              **研  TSy                ')
300 FORMAT('                            2006-4-18                 ')
400 FORMAT('                   请输入数据文件名[<=8字符]:'$         )
READ(*,*)NAME
FILE_IN = TRIM(NAME)//".TXT"
FILE_OUT = TRIM(NAME)//".DAT"
CALL INPUT(FILE_IN)  
CALL SFREE
CALL LOAD
CALL K
CALL GAUSSJ
CALL STRE_STRA
CALL OUTPUT(FILE_OUT)
END PROGRAM MAIN

SUBROUTINE INPUT(IN)
USE DYNA
CHARACTER::IN*12
INTEGER::NTEMP
OPEN(1,FILE=IN)
!读入项目名称
READ(1,'(A)')
READ(1,*) QUES  
READ(1,'(A)')
!读入单元、节点数
READ(1,*)NELEM, NNODE
ALLOCATE(MELEM(NELEM,4), STAT = NTEMP)
IF(NTEMP /= 0) STOP "内存空间不够"
ALLOCATE(MNODE(NNODE,2), STAT = NTEMP)
IF(NTEMP /= 0) STOP "内存空间不够"
ALLOCATE(NFIXX(NNODE), NFIXY(NNODE),PX(NNODE), PY(NNODE), STAT = NTEMP)
IF(NTEMP /= 0) STOP "内存空间不够"
NFIXX = 1
NFIXY = 1      !设定初始节点均自由
PX = 0
PY = 0         !设定初始节点荷载为0
!     读入单元信息
READ(1,'(A)')
DO I=1,NELEM
READ(1,*) II, (MELEM(II,J),J=1,4)
END DO
!     读入节点坐标信息
READ(1,'(A)')
DO I=1,NNODE
READ(1,*) II, (MNODE(II,J),J=1,2)
END DO
! 填充节点约束信息
READ(1,'(A)')
READ(1,*) NFIX
READ(1,'(A)')
DO I=1,NFIX
READ(1,*) II, NFIXX(II), NFIXY(II)  
END DO
! 读入材料信息
READ(1,'(A)')
READ(1,*) NMATERIAL
ALLOCATE(MMATERIAL(NMATERIAL, 3), STAT = NTEMP)
IF(NTEMP /= 0) STOP "内存不足"
READ(1,'(A)')
DO  I=1,NMATERIAL
READ(1,*) II, (MMATERIAL(II,J), J=1, 3)
END DO
! 读入点荷载信息
READ(1,'(A)')
READ(1,*) NP
READ(1,'(A)')
DO I=1,NP
READ(1,*) II,PX(II), PY(II)
END DO
!读入面荷载信息
READ(1,'(A)')
READ(1,*) NSP
ALLOCATE(SPXY(NSP,6), STAT = NTEMP)
IF(NTEMP /= 0) STOP "内存不足"
READ(1,'(A)')
DO I=1,NSP
READ(1,*) (SPXY(I,J),J=1,6)
END DO
CLOSE(1)
RETURN
END SUBROUTINE INPUT

SUBROUTINE SFREE
USE DYNA
INTEGER::NTEMP
!求最大自由度
MFREE = 0
DO I = 1, NNODE
IF(NFIXX(I).EQ.1)THEN
  MFREE = MFREE + 1
ENDIF
IF(NFIXY(I).EQ.1)THEN
  MFREE = MFREE + 1
ENDIF
ENDDO
ALLOCATE(NFREE(MFREE), P(MFREE), STAT = NTEMP)
IF(NTEMP /= 0) STOP "内存空间不够"
!进行自由度编号
MFREE = 0
DO I = 1, NNODE
IF(NFIXX(I).EQ.1)THEN
  MFREE = MFREE + 1
  NFIXX(I) = MFREE
  NFREE(MFREE) = I
  P(MFREE) = -PX(I)
ENDIF
IF(NFIXY(I).EQ.1)THEN
  MFREE = MFREE + 1
  NFIXY(I) = MFREE
  NFREE(MFREE) = -I
  P(MFREE) = -PY(I)
ENDIF
ENDDO
RETURN
END SUBROUTINE SFREE

SUBROUTINE K
USE DYNA
INTEGER::NTEMP
INTEGER::R, S, II
REAL::E, U, AREA, BR, CR, BS, CS, KK
REAL::A(3), B(3), C(3)
REAL:K(6,6), IJM(5),KDIR(6)
!DK(6,6)---单元刚度矩阵; IJM(5)----下标运算矩阵
ALLOCATE(D(MFREE,MFREE), STAT = NTEMP)
IF(NTEMP /= 0) STOP "内存空间不够"
D = 0
DO I = 1, NELEM
IJM(1) = MELEM(I, 1)
IJM(2) = MELEM(I, 2)
IJM(3) = MELEM(I, 3)
IJM(4) = MELEM(I, 1)
IJM(5) = MELEM(I, 2)
II = MELEM(I, 4)
E = MMATERIAL(II, 1)
U = MMATERIAL(II, 2)
DO J = 1, 3
A(J)=MNODE(IJM(J+1),1)*MNODE(IJM(J+2),2)-MNODE(IJM(J+1),2)*MNODE(IJM(J+2),1)


  B(J) = MNODE(IJM(J+1),2) - MNODE(IJM(J+2),2)
  C(J) = MNODE(IJM(J+2),1) - MNODE(IJM(J+1),1)
END DO
AREA = (A(1)+A(2)+A(3))/2
KK = E/(4*(1-U**2)*AREA)
DO  R = 1, 3
  DO  S = 1, 3
   DK(2*R-1,2*S-1)=KK*(B(R)*B(S)+(1-U)*C(R)*C(S)/2)
   DK(2*R-1,2*S)=KK*(U*B(R)*C(S)+(1-U)*C(R)*B(S)/2)
   DK(2*R,2*S-1)=KK*(U*C(R)*B(S)+(1-U)*B(R)*C(S)/2)
   DK(2*R,2*S)=KK*(C(R)*C(S)+(1-U)*B(R)*B(S)/2)
  END DO
END DO
J = MELEM(I, 1)
KDIR(1) = NFIXX(J)
KDIR(2) = NFIXY(J)
J = MELEM(I, 2)
KDIR(3) = NFIXX(J)
KDIR(4) = NFIXY(J)
J = MELEM(I, 3)
KDIR(5) = NFIXX(J)
KDIR(6) = NFIXY(J)
DO IP = 1, 6
  DO IQ = 1, 6
  IF((KDIR(IP)*KDIR(IQ)).NE.0)THEN
   D(KDIR(IP),KDIR(IQ)) = D(KDIR(IP),KDIR(IQ)) + DK(IP,IQ)
  ENDIF
  END DO
END DO
END DO
END SUBROUTINE K


SUBROUTINE LOAD
USE DYNA
INTEGER::J1,J2
REAL::A(3),IJM(5)
!J1、J2面力端节点
REAL::COS, SIN, TL, L
DO I = 1, NSP
J1 = SPXY(I,1)
J2 = SPXY(I,4)
L = SQRT((MNODE(J1,1)-MNODE(J2,1))**2+(MNODE(J1,2)-MNODE(J2,2))**2)
COS = (MNODE(J2,1)-MNODE(J1,1))/L
SIN = (MNODE(J2,2)-MNODE(J1,2))/L

! PP1=MLP1(I)*L/3+MLP2(I2)*L/6+MLJ1(I)*L/3+MLJ2(I)*L/6
! PP2=MLP2(I)*L/3+MLP1(I2)*L/6+MLJ2(I)*L/3+MLJ1(I)*L/6

PP1 = SPXY(I,2)*L/3+SPXY(I,6)*L/6+SPXY(I,3)*L/3+SPXY(I,5)*L/6
PP2 = SPXY(I,3)*L/3+SPXY(I,5)*L/6+SPXY(I,2)*L/3+SPXY(I,6)*L/6
JX1 = NFIXX(J1)
JX2 = NFIXX(J2)
JY1 = NFIXY(J1)
JY2 = NFIXY(J2)
IF(JX1.NE.0)THEN
  WRITE(*,*)'JX1=',JX1,'P(JX1)=',P(JX1)
  P(JX1) = P(JX1)+COS*PP1
ENDIF
IF(JX2.NE.0)THEN
  P(JX2) = P(JX2)+COS*PP2
ENDIF
IF(JY1.NE.0)THEN
  P(JY1) = P(JY1)+SIN*PP1
ENDIF
IF(JY2.NE.0)THEN
  P(JY2) = P(JY2)+SIN*PP2
ENDIF
ENDDO
DO I = 1, NELEM
IJM(1) = MELEM(I, 1)
IJM(2) = MELEM(I, 2)
IJM(3) = MELEM(I, 3)
IJM(4) = MELEM(I, 1)
IJM(5) = MELEM(I, 2)
II = MELEM(I, 4)
DO J = 1, 3
  A(J)=MNODE(IJM(J+1),1)*MNODE(IJM(J+2),2)-MNODE(IJM(J+1),2)*MNODE(IJM(J+2),1)


END DO
AREA = (A(1)+A(2)+A(3))/2
TL = MMATERIAL(II,3)*AREA/3
II = NFIXY(IJM(1))
JJ = NFIXY(IJM(2))
MM = NFIXY(IJM(3))
IF(II.NE.0)THEN
  P(II) = P(II)+TL
ENDIF
IF(JJ.NE.0)THEN
  P(JJ) = P(JJ)+TL
ENDIF
IF(MM.NE.0)THEN
  P(MM) = P(MM)+TL
ENDIF
ENDDO
END SUBROUTINE LOAD



SUBROUTINE GAUSSJ
USE DYNA
INTEGER::NTEMP
INTEGER::I,ICOL,IROW,J,K,L,LL
REAL, DIMENSION(:,:), ALLOCATABLE::DL
REAL, DIMENSION(:), ALLOCATABLE::INDXC,INDXR,IPIV
REAL::BIG,DUM,PIVINV
ALLOCATE(INDXC(MFREE),INDXR(MFREE),IPIV(MFREE),BL(MFREE),STAT = NTEMP)
IF(NTEMP /= 0) STOP "内存空间不够"
ALLOCATE(DL(MFREE,MFREE),STAT = NTEMP)
IF(NTEMP /= 0) STOP "内存空间不够"
BL = P
DL = D
DO J=1,MFREE
   IPIV(J)=0
END DO
DO I=1,MFREE
   BIG=0
   DO J=1,MFREE
      IF(IPIV(J).NE.1)THEN
      DO  K=1,MFREE
   IF (IPIV(K).EQ.0) ThEN
      IF (ABS(DL(J,K)).GE.BIG)THEN
       BIG=ABS(DL(J,K))
                IROW=J
                ICOL=K
            ENDIF
            ELSE IF (IPIV(K).GT.1) THEN
                PAUSE 'SINGULAR MATRIX IN GAUSSJ'
            ENDIF
  END DO
     ENDIF
END DO
    IPIV(ICOL)=IPIV(ICOL)+1
    IF (IROW.NE.ICOL) THEN
       DO L=1,MFREE
          DUM=DL(IROW,L)
          DL(IROW,L)=DL(ICOL,L)
          DL(ICOL,L)=DUM
    END DO
       DUM=BL(IROW)
       BL(IROW)=BL(ICOL)
       BL(ICOL)=DUM
    ENdIF
    INDXR(I)=IROW
    INDXC(I)=ICOL
    IF (DL(ICOL,ICOL).EQ.0.) PAUSE 'SINGULAR MATRIX IN GAUSSJ'
       PIVINV=1./DL(ICOL,ICOL)
       DL(ICOL,ICOL)=1.
       DO L=1,MFREE
           DL(ICOL,L)=DL(ICOL,L)*PIVINV
       END DO
       BL(ICOL)=BL(ICOL)*PIVINV
       DO LL=1,MFREE
          IF(LL.NE.ICOL)THEN
             DUM=DL(LL,ICOL)
             DL(LL,ICOL)=0
             DO L=1,MFREE
                DL(LL,L)=DL(LL,L)-DL(ICOL,L)*DUM
             END DO
             BL(LL)=BL(LL)-BL(ICOL)*DUM
          ENdIF
    END DO
END DO
DO L=MFREE,1,-1
   IF(INDXR(L).NE.INDXC(L))THEN
      DO K=1,MFREE
         DUM=DL(K,INDXR(L))
         DL(K,INDXR(L))=DL(K,INDXC(L))
         DL(K,INDXC(L))=DUM
      END DO
   ENdIF
END DO
DEALLOCATE(DL)
DEALLOCATE(D)
DEALLOCATE(INDXC)
DEALLOCATE(INDXR)
DEALLOCATE(IPIV)
RETURN
END SUBROUTINE GAUSSJ

!计算应力、应变
SUBROUTINE STRE_STRA
USE DYNA
INTEGER::NTEMP
REAL::AREA, TEMP1, TEMP2, K1, K2
REAL::IJM(5), A(3), B(3), C(3), U1(3), V(3)
ALLOCATE(UU(NNODE),VV(NNODE),STAT = NTEMP)
IF(NTEMP /= 0) STOP "内存空间不够"
ALLOCATE(STRESS(NELEM,6),STAT = NTEMP)
IF(NTEMP /= 0) STOP "内存空间不够"
STRESS = 0
DO  I = 1, NNODE
IF(NFIXX(I).EQ.0)THEN
  UU(I) = 0.0
ELSE
  UU(I) = BL(NFIXX(I))
ENDIF
IF(NFIXY(I).EQ.0)THEN
  VV(I) = 0.0
ELSE
  VV(I) = BL(NFIXY(I))
ENDIF
END DO
DO I = 1, NELEM
IJM(1) = MELEM(I, 1)
IJM(2) = MELEM(I, 2)
IJM(3) = MELEM(I, 3)
IJM(4) = MELEM(I, 1)
IJM(5) = MELEM(I, 2)
III = MELEM(I, 4)
E = MMATERIAL(III, 1)
U = MMATERIAL(III, 2)
DO J = 1, 3
A(J)=MNODE(IJM(J+1),1)*MNODE(IJM(J+2),2)-MNODE(IJM(J+1),2)*MNODE(IJM(J+2),1)


  B(J) = MNODE(IJM(J+1),2) - MNODE(IJM(J+2),2)
  C(J) = MNODE(IJM(J+2),1) - MNODE(IJM(J+1),1)
END DO
AREA = (A(1)+A(2)+A(3))/2
K1 = E/(2*(1-U**2)*AREA)
K2 = E/(4*(1+U)*AREA)
DO II = 1, 3
  U1(II) = UU(IJM(II))
  V(II) = VV(IJM(II))
END DO
X = 0
Y = 0
XY = 0
DO II = 1, 3
  X = X + U1(II)*B(II)
  Y = Y + +C(II)*V(II)
  XY = XY + U1(II)*C(II) + V(II)*B(II)
END DO
STRESS(I,1) = K1*(X+U*Y)
STRESS(I,2) = K1*(Y+U*X)
STRESS(I,3) = K2*XY
TEMP1=STRESS(I,1)+STRESS(I,2)
TEMP2=SQRT((STRESS(I,1)-STRESS(I,2))*(STRESS(I,1)-STRESS(I,2))+4.0*STRESS(I,3

)*STRESS(I,3))
STRESS(I,4)=(TEMP1+TEMP2)/2.0
STRESS(I,5)=(TEMP1-TEMP2)/2.0
IF(ABS(STRESS(I,3)).GT.1E-4)THEN
  STRESS(I,6)=ATAN((STRESS(I,4)-STRESS(I,1))/STRESS(I,3))*57.29578
ELSE
     STRESS(I,6)=90.0
END IF
END DO
END SUBROUTINE STRE_STRA


!输出程序
SUBROUTINE OUTPUT(OUT)
USE DYNA
CHARACTER::OUT*12
OPEN(2, FILE = OUT)
WRITE(2,*)"题目"
WRITE(2,*)QUES
WRITE(2,*)"单元数       节点数"
WRITE(2,*) NELEM, NNODE
WRITE(2,*)"单元号     I节点号    J节点号     M节点号    材料类型"
DO I = 1, NELEM
WRITE(2, *)I, (MELEM(I, J),J= 1, 4)
END DO
WRITE(2,*)"节点号      X坐标      Y坐标"
DO I = 1, NNODE
WRITE(2, *)I, (MNODE(I, J),J= 1,2)
END DO
WRITE(2, *)"约束节点数"
WRITE(2, *)NFIX
WRITE(2, *)"节点号    X向     Y向(0--约束，1--自由)"
DO I = 1, NNODE
WRITE(2, *)I, NFIXX(I), NFIXY(I)
END DO
WRITE(2, *)"材料类型数"
WRITE(2, *)NMATERIAL
WRITE(2, *)"材料类型号      弹性模量    泊松比  重度"
DO I = 1, NMATERIAL
WRITE(2, *)I, (MMATERIAL(I,J), J=1, 3)
END DO
WRITE(2, *)"点荷载作用数"
WRITE(2, *)NP
WRITE(2, *)"节点号   X向荷载    Y向荷载"
DO I = 1, NNODE
WRITE(2, *)I, PX(I), PY(I)
END DO
WRITE(2, *)"面荷载作用数"
WRITE(2, *)NSP
WRITE(2, *)"节点号1   X向荷载    Y向荷载   节点号2   X向荷载    Y向荷载"
DO I = 1, NSP
WRITE(2, *)(SPXY(I,J),J=1,6)
END DO
WRITE(2,*)'最大自由度数：'
WRITE(2,*)MFREE
WRITE(2,*)'荷载向量：'
WRITE(2,*)P
!WRITE(2,*)'D矩阵'
!DO I = 1, MFREE
!  WRITE(2,*)(D(I,J),J = 1, MFREE)
!END DO
WRITE(2,*)'BL矩阵'
WRITE(2,*)BL
WRITE(2,*)'**************************************'
WRITE(2,*)'单元号    σx      σy      τxy      σ1     σ2   σ1与X轴夹角'


DO I = 1, NELEM
WRITE(2,*)I,STRESS(I,1),STRESS(I,2),STRESS(I,3),STRESS(I,4),STRESS(I,5),STRES

S(I,6)
ENDDO
END SUBROUTINE OUTPUT

iwtry

不错不错，

feiyang52013

我就是需要这样的
谢谢