无网格伽辽金方法计算源程序（fortran实现探讨）

红黑22

1. 问题描述
如图 B1 .1 所示,一个 2 mm×2 mm 的矩形板,一端固定,一端承受均匀载荷的拉伸。在板内排布 9 个结点,采用不规则排布,其中 5 号结点的坐标是( 0 .5, 0 .5 )。材料参数: 弹性模量 E =
100 GPa,泊松比μ= 0 .3,平面应力问题。采用 3×3 的背景网格,4×4 高斯积分。选择式(4 .11)的幂权函数,规定 c=2/2mm,因为结点布局很不规则,取 a = 10c。采用拉格朗日乘子强制满足位移边界条件。

2. 参量说明（略）
3..源程序及说明
   $DEBUG
C  主程序
    COMMON/ CONTRO/NPOIN,NCELL,NDOFN,NDIME,NSTRE ,
     &NTYPE,NGAUS,NNODE ,
     &NPROP, NVFIX, NBNOD, NEVAB, NNMAX , NECELX, NECELY,
     &NGASP, NLAGR, NINFL ,
     &NBASE , NBASB, NDIFL
      DIMENSION BCOOD(16,2) ,COORD(9,2) , COX( 9,2) , IFPRE(3 ,2) ,
     &LNODS(9 ) , LBNOD(9 ,4 ) , NOFIX(3) , PROPS(3 )
      DIMENSION BMATX(3 ,18) ,DERIV(2 ,9) , DMATX( 3,3 ) ,
     &DBMAT (3 ,18 ) , DEWEIG(2 ,9) , POSGP(4 ) , ESTIF(18, 18) ,
     &SMATX(3 ,18,16, 9) ,SH AB( 3,9 ) , SH APE( 9) ,WEIGP(4) ,
     &WEIGHT(9 )
      DIMENSION AA( 3, 3) ,DETA(3 ,3) ,DERIVA(3 ,3)
      DIMENSION PT (3 ) , DEP T( 3)
      DIMENSION DERIVB(3 ,9) , DERIV1 (9) , DERIV2 (9 ) , DERIV3 (9 ) ,
      &B2(3 ,9 )
      DIMENSION DSH APE( 2, 6) , GIK(6 ,6 )
      DIMENSION DIP(24) , ESTIFQ(24, 24) , ELOADQ( 24)
      OPEN( 1, FILE =′ INIT .TXT′)
      OPEN( 2, FILE =′ INPUT .TXT′)
      OPEN( 3, FILE =′ STIFPS .TXT′ )
      OPEN( 4, FILE =′ GNIK .TXT′ )
      OPEN( 5, FILE =′ LOADPS .TXT′)
      OPEN( 6, FILE =′ EFGDIS .TXT′ )
      OPEN( 7, FILE =′ ESTIFQ .TXT′)
      WRITE( * ,5)
5     FORMAT (1X,′ READING THE MESSAGE FROM INIT .TXT . . . ′ / )
      READ(1 , * )NPOIN, NECELX, NECELY, NDOFN , NDIME, NSTRE,
      &NTYPE , NPROP , NVFIX, NGAUS, NNODE , NINFL , NBASE ,
      &NBASB
C * * * * * * * * * * * * * * *
      NCELL = NECELX * NECELY
      NBNOD = ( NECELX + 1 ) * ( NECELY + 1)
      NEVAB = NDIME * NPOIN
      NLAGR = NDIME * NVFIX
      NNMAX = NEVAB+ NLAGR
      NDIFL = NDIME * NINFL
C * * * * * * * * * * * * * * *
      WRITE( * ,7) NPOIN , NECELX , NECELY, NDOFN , NDIME ,
      &NSTRE , NTYPE , NCELL , NPROP , NVFIX, NBASE , NBASB,
      &NBNOD, NEVAB, NNMAX, NGAUS, NNODE , NLAGR
7     FORMAT (1X ,′ NPOIN =′ , I3 ,1X ,′ NECELX =′ , I3, 1X,′ NECELY
      & =′ , I3, 1X, &′ NDOFN =′ , I3 ,1X ,′ NDIME =′ , I3 ,1X ,′ NSTRE
      & =′ , I3 / 1X,′ NTYPE =′ , I3 ,1X ,′ NCELL =′ , I3 ,1X,′ NPROP =′,
      &I3 ,1X ,′ NVFIX =′ , I3 ,1X,′ NBASE =′ , I3,1X ,′ NBASB =′ , I3,
      &1X,/ 1X,′ NBNOD =′ , I3,1X ,′ NEVAB =′ , I3, 1X,′ NNMAX =′
      &I3 ,1X ,′ NGAUS =′ , I3 ,1X,′ NNODE =′ , I3 ,1X,′ NLAGR =′ , I3 /)
       WRITE( * ,10 )
10    FORMAT (1X,′ THE DATA IS INPUTTING FROM INPUT .TXT . . . ′/)
      CALL INPUT(BCOOD,COORD ,DIMAX , SURPC, SURPG, IFPRE ,
      &LBNOD , NOFIX, PROPS, QLOAD)
      WRITE( * ,20 )
20    FORMAT(1X,′ CALCULATING THE STIFFNESS OF THIS
     &PROBLEM . . .′ / )
      CALL STIFPS(BCOOD ,COORD ,COX, DIMAX, LBNOD, PROPS,
      &ESTIFQ , SURPC , SH APE , ESTIF , LNODS, DERIV ,DMATX,
      &DBMAT ,BMATX, SMATX, POSGP, WEIGP ,WEIGHT , DEWEIG ,
      &SH AB, DERIVB,DERIV1 ,DERIV2 ,DERIV3 ,B2 )
C * * * * * * * * * * * * * * *
      WRITE( * ,30 )
30    FORMAT (1X,′ CALCULATING THE LAGRANGE MULTIPLIERS . . .′/)
      CALL GNIK(COORD ,COX ,DSH APE ,DIMAX, ESTIFQ , GIK ,
      &LNODS, LBNOD, SURPG, SH AB, SHAPE , POSGP ,WEIGP ,B2,
      &WEIGHT , NOFIX , IFPRE)
C * * * * * * * * * * * * * * *
      WRITE(7 ,90)
90    FORMAT(1X,′ ESTIFQ =′
      DO I = 1 , NNMAX
      WRITE(7 ,80) I , ESTIFQ( I , I )
      END DO
80    FORMAT(1X, I4,1X , E10.4)
C * * * * * * * * * * * * * * *
      WRITE( * ,40 )
40    FORMAT(1X,′ CALCULATING THE LOADS OF BOUNDARY
      &NODES . . .′ / )
      CALL LOADPS(COORD ,COX, DIMAX, LNODS , LBNOD , SHAPE ,
      &POSGP ,WEIGP , WEIGHT , SH AB, QLOAD ,SURPG ,B2, NOFIX ,
      &ELOADQ)
C * * * * * * * * * * * * * * *
      WRITE( * ,50)
40    FORMAT(1X,′ CALCULATING THE LOADS OF BOUNDARY
      &NODES . . .′ / )
      CALL LOADPS(COORD ,COX, DIMAX, LNODS , LBNOD , SHAPE ,
      &POSGP ,WEIGP , WEIGHT , SH AB, QLOAD ,SURPG ,B2, NOFIX ,
      &ELOADQ)
C * * * * * * * * * * * * * * *
      WRITE( * ,50 )
50    FORMAT(1X,′ CALCULATING THE DISPLACEMENTS . . .′ / )
      CALL GAUDIP( ESTIFQ , ELOADQ , DIP)
C * * * * * * * * * * * * * * *
      WRITE( * ,60 )
60    FORMAT(1X,′ OK , THE RESULTS HAVE BEEN WRITTEN IN
      &FILE EFGDIS .TXT .′ )
      END
C 数据输入子程序
  SUBROUTINE INPUT(BCOOD,COORD, DIMAX, SURPC, SURPG,
&IFPRE , LBNOD, NOFIX, PROPS, QLOAD)
  COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE
&NBASB, NDIFL
DIMENSION BCOOD(NBNOD, NDIME) ,COORD(NPOIN, NDIME) ,
&IFPRE( NVFIX, NDOFN) , LBNOD( NCELL , NNODE) , NOFIX
&(NVFIX) , PROPS( NPROP)
C INPUT THE YOUNG′ S MODULUS , POISSION′ S RATIO , THICK
C &OF MODEL
   DO 15 I = 1 , NPROP
   READ(2 , * ) PROPS( I )
15 CONTINUE
C INPUT THE FIX NUMBER OF DISPLACEMENT BOUNDARY
  DO 20 I = 1 , NVFIX
  READ(2 , * ) NOFIX( I )
20 CONTINUE
C INPUT THE FIX DIMENSION OF DISPLACEMENT BOUNDARY
  DO 25 I = 1 , NVFIX
  READ(2 , * ) IFPRE( I ,1) , IFPRE( I , 2)
25 CONTINUE
C INPUT THE COORDINATE OF NODES
  DO 30 I = 1 , NPOIN
READ(2 , * )COORD( I ,1 ) , COORD( I ,2)
30 CONTINUE
C * * * * * * * * * * * * * * *
C DO 60 I = 1, NPOIN
C WRITE(3 ,103) I , COORD( I ,1 ) , I ,COORD( I ,2)
C103 FORMAT(1X,′ COORD(′ , I3 ,′ ,1) =′ , F8 . 2, 2X,′ COORD(′ , I3,′ ,2
C & =′ , F8 . 2 ,2X/ )
C60 CONTINUE
C * * * * * * * * * * * * * * *
C INPUT THE NODE NUMBER OF EVERY CELL
  DO 35 I = 1 , NCELL
  READ(2 , * ) LBNOD( I ,1) , LBNOD( I ,2) , LBNOD( I ,3) , LBNOD( I,4)
35 CONTINUE
C  INPUT THE COORDINATE OF POINTS OF BAKGROUND CELLS
   DO 40 I = 1 , NBNOD
   READ(2 , * )BCOOD( I ,1) ,BCOOD( I ,2)
40 CONTINUE
C * * * * * * * * * * * * * * *
C DO 70 I = 1, NBNOD
C WRITE(3 ,104) I , BCOOD( I ,1) , I ,BCOOD( I , 2)
C104 FORMAT(1X,′ BCOOD(′ , I3 ,′ ,1) =′ , F8 . 2,2X,′ BCOOD(′ , I3 ,′ ,2)
C & =′ , F8 . 2 ,2X/ )
C70 CONTINUE
C * * * * * * * * * * * * * * *
   READ(2 , * )SURPC ,SURPG , DIMAX
   READ(2 , * ) QLOAD
   RETURN
   END
C 高斯点坐标及加权因子计算子程序
  SUBROUTINE GAUSSQ(POSGP ,WEIGP)
  COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
  &NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NVEAB,
  &NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
  &NBASB, NDIFL
  DIMENSION POSGP(NGAUS) ,WEIGP(NGAUS)
  IF(NGAUS .EQ . 2 ) THEN
  POSGP( 1) = - 0 . 577350269189626
  WEIGP(1 ) = 1 . 0
  POSGP( 2) = - POSGP( 1)
  WEIGP(2 ) = WEIGP(1 )
ELSE IF(NGAUS .EQ .3) THEN
POSGP( 1) = - 0 . 774596669241483
POSGP( 2) = 0 . 0
WEIGP(1 ) = 0 . 5555555555555556
WEIGP(2 ) = 0 . 8888888888888889
POSGP( 3) = - POSGP( 1)
WEIGP(3 ) = WEIGP(1 )
ELSE IF(NGAUS .EQ .4) THEN
POSGP( 1) = - 0 . 861136311594053
POSGP( 2) = - 0 . 339981043584856
WEIGP(1 ) = 0 . 347854845137454
WEIGP(2 ) = 0 . 652145154862546
POSGP( 3) = - POSGP( 2)
WEIGP(3 ) = WEIGP(2 )
POSGP( 4) = - POSGP( 1)
WEIGP(4 ) = WEIGP(1 )
END IF
RETURN
END
C 系统刚度计算及组装子程序
     SUBROUTINE STIFPS(BCOOD, COORD,COX ,DIMAX , LBNOD ,
     &PROPS, ESTIFQ,SURPC,SHAPE , ESTIF , LNODS ,DERIV ,DMATX ,
     &DBMAT ,BMATX, SMATX, POSGP, WEIGP ,WEIGHT , DEWEIG ,
     &SHAB, DERIVB,DERIV1 ,DERIV2 ,DERIV3 ,B2 )
      COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
     &NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
     &NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
     &NBASB, NDIFL
     DIMENSION BCOOD(NBNOD, NDIME) ,COORD(NPOIN, NDIME) ,
     &PROPS(NPROP) , LNODS( NINFL) , LBNOD(NCELL , NNODE)
     DIMENSION SHAPE(NINFL) ,DERIV(NDIME , NINFL) , DMAT
     &(NSTRE, NSTRE) , DBMAT( NSTRE , NDIFL) , BMATX(NSTRE
     &NDIFL) , DERIVB( NBASE , NINFL) , SMATX(NSTRE , NDIFL ,
     &NGASP, NCELL) , POSGP(NGAUS) ,WEIGP(NGAUS) , WEIGH
     &(NINFL) ,DEWEIG(NDIME , NINFL) , SHAB(NBASE , NINFL) ,
     &B2(NBASE , NINFL)
     DIMENSION DERIV1(NINFL) , DERIV2(NINFL) , DERIV3( NINFL)
     DIMENSION ESTIFQ( NNMAX, NNMAX)
DIMENSION ESTIF( NDIFL , NDIFL)
DIMENSION COX( NINFL , NDIME)
DO 70 ICELL = 1 , NCELL
CALL MODPS(DMATX, PROPS)
THICK = PROPS( 3)
KGASP = 0
NGAUS = 4
CALL GAUSSQ(POSGP ,WEIGP)
DO 50 IGAUS = 1 , NGAUS
DO 50 JGAUS = 1, NGAUS
DO 20 IEVAB = 1, NDIFL
DO 20 JEVAB = 1 , NDIFL
ESTIF( IEVAB, JEVAB) = 0 .0
20 CONTINUE
KGASP = KGASP + 1
EXISP = POSGP( IGAUS)
ETASP = POSGP( JGAUS)
CALL SFR(BCOOD,COORD,COX,DIMAX, EXISP, ETASP , ICELL ,
&LBNOD , LNODS, NOCE1, SURPC,DERIV , DERIVB,DERIV1 ,
&DERIV2 ,DERIV3, SH APE ,WEIGHT , DEWEIG ,SH AB, B2 )
CALL BMATPS(BMATX, DERIV , NOCE1)
CALL DBE(BMATX, DBMAT ,DMATX, NOCE1)
CALL GAUSSQ(POSGP ,WEIGP)
DVOLU = (SURPC/ 2 . 0) * (SURPC/ 2 . 0) * WEIGP( IGAUS) *
&WEIGP( JGAUS)
WRITE(3 , * )WEIGP( IGAUS) ,WEIGP( JGAUS) , DVOLU
IF( THICK .NE . 0 . 0) DVOLU = DVOLU * THICK
DO 30 IEVAB = 1, NDIME * NOCE1
DO 30 JEVAB = 1 , NDIME * NOCE1
DO 30 ISTRE = 1 , NSTRE
ESTIF( IEVAB, JEVAB) = ESTIF( IEVAB, JEVAB) +
&BMATX( ISTRE , IEVAB) * DBMAT ( ISTRE , JEVAB) * DVOLU
30 CONTINUE
DO 80 I = 1 , NOCE1
DO 80 J = 1 , NOCE1
ESTIFQ(2 * LNODS( I ) - 1,2 * LNODS( J) - 1) = ESTIFQ(2 * LNODS
&( I ) - 1 ,2 * LNODS( J ) - 1) + ESTIF( 2 * I - 1 ,2 * J - 1)
ESTIFQ(2 * LNODS( I ) - 1 ,2 * LNODS( J) ) = ESTIFQ(2 * LNODS
&( I ) - 1 ,2 * LNODS( J ) ) + ESTIF(2 * I - 1 ,2 * J)
ESTIFQ(2 * LNODS( I ) ,2 * LNODS( J ) - 1) = ESTIFQ(2 * LNODS
&( I ) ,2 * LNODS( J) - 1) + ESTIF(2 * I ,2 * J - 1)
ESTIFQ(2 * LNODS( I ) ,2 * LNODS( J ) ) = ESTIFQ(2 * LNODS
&( I ) ,2 * LNODS( J) ) + ESTIF(2 * I ,2 * J )
80 CONTINUE
C DO 40 ISTRE = 1, NSTRE
C DO 40 IEVAB = 1 , NDIME * NOCE1
C40 SMATX ( ISTRE , IEVAB , KGASP , ICELL) = DBMAT ( ISTRE, IEVAB)
50 CONTINUE
OPEN( 27 , FILE =′ ESF .TXT′ )
DO I = 1 , NOCE1
WRITE(27, 111 ) I , ESTIF( I , I )
END DO
111 FORMAT (1X, I4 ,1X, E10.5)
70 CONTINUE
RETURN
END
C 弹性矩阵计算子程序
SUBROUTINE MODPS(DMATX, PROPS)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION PROPS(NPROP)
DIMENSION DMATX(NSTRE , NSTRE)
Y = PROPS( 1)
P = PROPS(2)
DO 10 I = 1 , NSTRE
DO 10 J = 1 , NSTRE
DMATX( I , J ) = 0 . 0
10 CONTINUE
IF(NTYPE .EQ . 1) THEN
CONST = Y/ ( 1 . 0 - P * P)
DMATX( 1, 1) = CONST
DMATX( 2, 2) = CONST
DMATX( 1, 2) = CONST * P
DMATX( 2, 1) = CONST * P
DMATX( 3, 3) = ( 1 . 0 - P) * CONST * 0 .5
20 ELSE IF( NTYPE .EQ .2) THEN
CONST = Y * (1 . 0 - P)/ ( ( 1 . 0 + P) * (1 . 0 - 2 . 0 * P) )
DMATX( 1, 1) = CONST
DMATX( 2, 2) = CONST
DMATX( 1, 2) = CONST * P/ ( 1 .0 - P)
DMATX( 2, 1) = CONST * P/ ( 1 .0 - P)
DMATX( 3, 3) = CONST * ( 1 . 0 - 2 .0 * P)/ ( 2 .0 * (1 . 0 - P) )
END IF
RETURN
END
C 应变矩阵计算子程序
SUBROUTINE BMAT PS(BMATX , DERIV, NOCE1 )
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION BMATX(NSTRE , NDIFL) , DERIV(NDIME , NINFL)
NGASH = 0
DO 10 I = 1 , NOCE1
MGASH = NGASH + 1
NGASH = MGASH + 1
BMATX( 1, MGASH) = DERIV( 1, I )
BMATX( 1, NGASH) = 0 . 0
BMATX( 2, MGASH) = 0 . 0
BMATX( 2, NGASH) = DERIV(2, I )
BMATX( 3, MGASH) = DERIV( 2, I )
BMATX( 3, NGASH) = DERIV(1, I )
10 CONTINUE
RETURN
END
C 矩阵 DB相乘子程序
SUBROUTINE DBE(BMATX , DBMAT , DMATX , NOCE1 )
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION DMATX (NSTRE , NSTRE) ,DBMAT (NSTRE , NDIFL) ,
&BMATX(NSTRE , NDIFL)
DO 10 ISTRE = 1 , NSTRE
DO 10 IEVAB = 1, NDIME * NOCE1
DBMAT( ISTRE , IEVAB) = 0 . 0
DO 10 JSTRE = 1 , NSTRE
DBMAT( ISTRE , IEVAB) = DBMAT( ISTRE , IEVAB) +
&DMATX( ISTRE , JSTRE) * BMATX( JSTRE , IEVAB)
10 CONTINUE
RETURN
END
C 形函数及其导数计算子程序
SUBROUTINE SFR(BCOOD ,COORD ,COX, DIMAX, S , T , ICELL
&LBNOD , LNODS, NOCE1, SURPC,DERIV , DERIVB,DERIV1 ,
&DERIV2 ,DERIV3, SH APE ,WEIGHT , DEWEIG ,SH AB, B2 )
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION PT(NBASE) ,DEPT(NBASE) ,DERIVA(NBASE, NBASE) ,
&DERIVB(NBASE , NINFL)
DIMENSION DERIV1(NINFL) , DERIV2(NINFL) , DERIV3( NINFL)
DIMENSION COORD(NPOIN, NDIME) ,BCOOD(NBNOD, NDIME) ,
&LNODS(NINFL) , LBNOD(NCELL , NNODE)
DIMENSION SHAPE(NINFL) ,DERIV(NDIME , NINFL) ,
&WEIGHT(NINFL) ,DETA(NBASE , NBASE) , DEWEIG( NDIME ,
&NINFL) , SHAB(NBASE , NINFL) ,B2(NBASE , NINFL)
DIMENSION COX( NINFL , NDIME)
DO 10 I = 1 , NPOIN
SHAPE( I ) = 0 . 0
DERIV1 ( I ) = 0 . 0
DERIV2 ( I ) = 0 . 0
DERIV3 ( I ) = 0 . 0
10 CONTINUE
DO 35 I = 1 , NDIME
DO 35 J = 1 , NPOIN
DERIV( I , J) = 0 . 0
35 CONTINUE
DO 15 I = 1 , NPOIN
WEIGHT( I ) = 0 . 0
DEWEIG( 1, I ) = 0 . 0
DEWEIG( 2, I ) = 0 . 0
15 CONTINUE
SURPD = DIMAX * SURPC
CALL CHAGAU(S, T , EGACX, EGACY , EX, EY, ICELL , SURPC ,
&BCOOD, LBNOD)
CALL SEARCH(COORD, EX , EY , SURPD , NOCE1 , LNODS)
C-------------------------------------------------------------
2 FORMAT (1X,′ NOCE1 =′ ,1X, I4)
WRITE(3 ,2 )NOCE1
C-------------------------------------------------------------
CALL CHACOORD(COORD,BCOOD,COX, NOCE1, ICELL, LNODS,
&LBNOD)
P T(1 ) = 1 . 0
P T(2 ) = EGACX
P T(3 ) = EGACY
WRITE(3 ,400) ICELL
400 FORMAT (1X,′ ICELL =′ , I3 )
C------------------------------------------------------------
DO 20 INODE = 1 , NOCE1
DEGAC = SQRT( ( EGACX - COX( INODE , 1) ) * * 2 +
&( EGACY - COX( INODE ,2) ) * * 2)
WEIGHT( INODE) = WEIGHTT (DEGAC, SURPD, SURPC)
DEWEIG( 1, INODE) = DEWEI (DEGAC, SURPD, SURPC,
&EGACX - COX( INODE ,1) )
DEWEIG( 2, INODE) = DEWEI (DEGAC, SURPD, SURPC,
&EGACY - COX( INODE ,2 ) )
20 CONTINUE
C---------------OUT PUT WEIGHT FUNCTION-----------
DO I = 1 , NOCE1
WRITE(3 ,101) I , LNODS( I ) ,WEIGHT( I )
END DO
101 FORMAT (1X, I4 ,1X, I4, 1X,′ WEIGHT =′ , F8 . 5)
C WRITE(3 ,199) ( (DEWEIG( I , J ) , J = 1, NOCE1) , I = 1, NDIME)
C199 FORMAT(1X,′ DEWEIG =′ ,100( 1X, F8 . 5) )
C---------------------------------------------------------
CALL DETA1(COX, LNODS, NOCE1 ,WEIGHT , DETA)
CALL SHAB1 (COX , LNODS, NOCE1, WEIGHT , SHAB)
CALL PDB( PT ,DETA , SHAB,B2, SHAPE , NOCE1)
C------------------------------------------------------
DO I = 1 , NOCE1
WRITE(3 ,1001 ) I , LNODS( I ) , SH APE( I )
END DO
1001 FORMAT( 1X, I4, 1X, I4,1X ,′ SHAPE( I ) =′ , F8 .4)
C---------------------------------------------------
DO 60 I = 1 , NDIME
DO 25 J = 1 ,3
DEP T( J ) = 0 . 0
25 CONTINUE
DEP T( I + 1) = 1 . 0
CALL DERIVA1( I ,COX,DERIVA , DETA , DEWEIG, NOCE1, LNODS)
CALL DERIVB1( I , COX, DERIVB, DEWEIG , LNODS, NOC
CALL PDB(DEPT ,DETA , SHAB,B2, DERIV1, NOCE1)
CALL PDB( PT ,DERIVA , SHAB,B2, DERIV2, NOCE1)
CALL PDB( PT ,DETA , DERIVB,B2 ,DERIV3 , NOCE1 )
DO 80 K = 1 , NOCE1
DERIV( I , K) = DERIV1(K) + DERIV2 (K) + DERIV3 (K)
80 CONTINUE
60 CONTINUE
c-------------------------------------------------
DO J = 1, NOCE1
WRITE(3 ,1003 ) J ,DERIV(1 , J)
END DO
1003 FORMAT( 1X, I4, 1X,′ DERIV =′ ,100F8 . 4)
WRITE(3 ,1004 )
1004 FORMAT( 1X,/ )
RETURN
END
C 将模型坐标转换为积分背景网格的局部坐标子程序
SUBROUTINE CH ACOORD(COORD ,BCOOD, COX, NOCE1 ,
&ICELL , LNODS , LBNOD)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COORD( NPOIN , NDIME) , LNODS(NINFL) ,
&COX( NINFL , NDIME)
DIMENSION BCOOD( NBNOD, NDIME) , LBNOD(NCELL ,4 )
DO I = 1 , NINFL
DO J = 1, NDIME
COX( I , J) = 0
ENDDO
ENDDO
DO I = 1 , NOCE1
IORNO = ABS( LBNOD( ICELL , 1) )
IORNI = ABS( LBNOD( ICELL ,3) )
CENX = (BCOOD( IORNI , 1) + BCOOD( IORNO ,1 ) )/ 2
CENY = (BCOOD( IORNI ,2 ) + BCOOD( IORNO, 2) )/ 2
COX( I ,1) = COORD( LNODS( I ) ,1) - CENX
COX( I ,2) = COORD( LNODS( I ) ,2) - CENY
ENDDO
RETURN
END
C 将标准单元的高斯点坐标转换为模型坐标的子程序
SUBROUTINE CH AGAU(S , T , EGACX , EGACY , EX , EY , ICELL
&SURPC,BCOOD , LBNOD)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE
&NBASB, NDIFL
DIMENSION BCOOD(NBNOD, NDIME) , LBNOD(NCELL , NNODE)
EGACX = ( (S / 2 . 0 ) * SURPC)
EGACY = ( ( T/ 2 . 0) * SURPC)
IORNO = ABS( LBNOD( ICELL , 2) )
EX = ( ( ( 1 + S)/ 2 . 0 ) * SURPC) + BCOOD( IORNO ,1 )
EY = ( ( (1 + T)/ 2 . 0) * SURPC) + BCOOD( IORNO ,2)
WRITE(3 ,101) EGACX, EGACY
101 FORMAT (1X,′ EGACX =′ , F10 . 4, 2X,′ EGACY =′ , E10 . 4 , E10 . 4)
RETURN
END
C 高斯点形成的权函数影响域内的结点搜索子程序
SUBROUTINE SEARCH(COORD, EX, EY , DM,M, LNODS)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COORD( NPOIN , NDIME) , LNODS(NINFL)
M = 0
DO 20 I = 1 , NPOIN
LNODS( I ) = 0
20 CONTINUE
DO 10 IPOIN = 1 , NPOIN
DEGACA = SQRT ( ( EX - COORD( IPOIN , 1) ) * * 2 +
&( EY - COORD( IPOIN ,2) ) * * 2)
IF(DEGACA .LE .DM) THEN
M = M + 1
LNODS( M) = IPOIN
END IF
10 CONTINUE
WRITE(3 ,12) ( LNODS( I ) , I = 1, NPOIN)
12 FORMAT(1X,′ LNODS =′ ,100I3)
RETURN
END
C 权函数计算子函数
REAL FUNCTION WEIGHTT(D ,DM,C)
WEIGH U = EXP( - 1 * (D/ C) * * 2) - EXP( - 1 * (DM/ C) * * 2)
WEIGHD = 1 - EXP( - 1 * (DM/ C) * * 2 )
WEIGHTT = WEIGHU/ WEIGHD
END
C 权函数导数计算子函数
REAL FUNCTION DEWEI (D ,DM, C, X)
DEWEIU = - 2 * X * EXP( - 1 * (D/ C) * * 2)
DEWEID = (C * * 2 ) * (1 - EXP( - 1 * (DM/ C) * * 2) )
DEWEI = DEWEIU/ DEWEID
END
C 权函数导数计算子函数
REAL FUNCTION DEWEI (D ,DM, C, X)
DEWEIU = - 2 * X * EXP( - 1 * (D/ C) * * 2)
DEWEID = (C * * 2 ) * (1 - EXP( - 1 * (DM/ C) * * 2) )
DEWEI = DEWEIU/ DEWEID
END
C 计算形函数中矩阵 A的逆矩阵子程序
SUBROUTINE DETA1 (COX , LNODS , NOCE1, WEIGHT ,DETA)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COX( NINFL , NDIME) , LNODS(NINFL)
DIMENSION WEIGHT(NINFL) , DETA(NBASE , NBASE)
DIMENSION AA( NBASE , NBASE)
DO 5 I = 1, NBASE
DO 5 J = 1, NBASE
AA( I , J) = 0.0
5 CONTINUE
DO 10 I = 1 , NOCE1
AA(1 ,1) = AA(1 ,1 ) + WEIGHT( I )
AA(1 ,2) = AA(1 ,2 ) + COX( I , 1) * WEIGHT( I )
AA(1 ,3) = AA(1 ,3 ) + COX( I , 2) * WEIGHT( I )
AA(2 ,1) = AA(2 ,1 ) + COX( I , 1) * WEIGHT( I )
AA(2 ,2) = AA(2 ,2 ) + COX( I , 1) * COX( I ,1) * WEIGHT ( I )
AA(2 ,3) = AA(2 ,3 ) + COX( I , 1) * COX( I ,2) * WEIGHT ( I )
AA(3 ,1) = AA(3 ,1 ) + COX( I , 2) * WEIGHT( I )
AA(3 ,2) = AA(3 ,2 ) + COX( I , 1) * COX( I ,2) * WEIGHT ( I )
AA(3 ,3) = AA(3 ,3 ) + COX( I , 2) * COX( I ,2) * WEIGHT ( I )
10 CONTINUE
CALL AMATX(AA , DETA)
END
SUBROUTINE AMATX( AA, DETA)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION DETA(NBASE , NBASE)
DIMENSION AA( NBASE , NBASE) , AMAX(NBASE , NBASE)
AMAX( 1,1 ) = AA(2 ,2) * AA(3 ,3) - AA( 2,3 ) * AA(3, 2)
AMAX( 1,2 ) = - 1 * ( AA( 2,1 ) * AA(3, 3) - AA(2 ,3) * AA( 3,1 ) )
AMAX( 1,3 ) = AA(2 ,1) * AA(3 ,2) - AA( 2,2 ) * AA(3, 1)
AMAX( 2,1 ) = - 1 * ( AA( 1,2 ) * AA(3, 3) - AA(1 ,3) * AA( 3,2 ) )
AMAX( 2,2 ) = AA(1 ,1) * AA(3 ,3) - AA( 1,3 ) * AA(3, 1)
AMAX( 2,3 ) = - 1 * ( AA( 1,1 ) * AA(3, 2) - AA(1 ,2) * AA( 3,1 ) )
AMAX( 3,1 ) = AA(1 ,2) * AA(2 ,3) - AA( 1,3 ) * AA(2, 2)
AMAX( 3,2 ) = - 1 * ( AA( 1,1 ) * AA(2, 3) - AA(1 ,3) * AA( 2,1 ) )
AMAX( 3,3 ) = AA(1 ,1) * AA(2 ,2) - AA( 1,2 ) * AA(2, 1)
DDA1 = AA( 1,1 ) * AA(2 ,2) * AA(3 ,3)
DDA2 = AA( 1,3 ) * AA(2 ,1) * AA(3 ,2)
DDA3 = AA( 3,1 ) * AA(1 ,2) * AA(2 ,3)
DDA4 = AA( 1,3 ) * AA(2 ,2) * AA(3 ,1)
DDA5 = AA( 1,1 ) * AA(3 ,2) * AA(2 ,3)
DDA6 = AA( 1,2 ) * AA(2 ,1) * AA(3 ,3)
DDA = DDA1 + DDA2 + DDA3 - DDA4 - DDA5 - DDA6
DO 10 I = 1 , NBASE
DO 10 J = 1 , NBASE
DETA( I , J ) = AMAX( J , I )/ DDA
10 CONTINUE
END
C 计算形函数中矩阵 C的子程序
SUBROUTINE SHAB1(COX, LNODS, NOCE1 ,WEIGHT , SHAB)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COX( NINFL , NDIME) , LNODS(NINFL)
DIMENSION WEIGHT(NINFL) , SHAB( NBASE , NINFL)
DO 10 I = 1 , NOCE1
SHAB(1 , I ) = WEIGHT ( I )
SHAB(2 , I ) = COX( I ,1 ) * WEIGHT( I )
SHAB(3 , I ) = COX( I ,2 ) * WEIGHT( I )
10 CONTINUE
END
C 计算形函数中矩阵 A的逆矩阵导数的子程序
SUBROUTINE DERIVA1 ( I I ,COX , DERIVA , DETA , DEWEIG ,
&NOCE1 , LNODS)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COX( NINFL , NDIME) , LNODS(NINFL)
DIMENSION DEWEIG( NDIME , NINFL) , DETA( NBASE , NBASE)
DIMENSION DERIVA (NBASE , NBASE) ,DEAA( NBASE , NBASE) ,
&A2(NBASE , NBASE)
DO 30 I = 1 , NBASE
DO 30 J = 1 , NBASE
DERIVA( I , J ) = 0.0
DEAA( I , J) = 0.0
30 CONTINUE
DO 80 I = 1 , NOCE1
DEAA(1 ,1) = DEAA( 1,1 ) + DEWEIG( II , I )
DEAA(1 ,2) = DEAA( 1,2 ) + COX( I , 1) * DEWEIG( I I , I )
DEAA(1 ,3) = DEAA( 1,3 ) + COX( I , 2) * DEWEIG( I I , I )
DEAA(2 ,1) = DEAA( 2,1 ) + COX( I , 1) * DEWEIG( I I , I )
DEAA(2,2) = DEAA(2,2) + COX( I ,1) * COX( I ,1 ) * DEWEIG( I I, I )
DEAA(2,3) = DEAA(2,3) + COX( I ,1) * COX( I ,2 ) * DEWEIG( I I, I )
DEAA(3 ,1) = DEAA( 3,1 ) + COX( I , 2) * DEWEIG( I I , I )
DEAA(3,2) = DEAA(3,2) + COX( I ,1) * COX( I ,2 ) * DEWEIG( I I, I )
DEAA(3,3) = DEAA(3,3) + COX( I ,2) * COX( I ,2 ) * DEWEIG( I I, I )
80 CONTINUE
DO 5 I = 1, NBASE
DO 5 J = 1, NBASE
A2( I , J ) = 0 . 0
5 CONTINUE
DO 10 I = 1 , NBASE
DO 10 J = 1 , NBASE
DO 10 K = 1 , NBASE
A2( I , J ) = A2( I , J ) + DETA( I , K) * DEAA(K , J)
10 CONTINUE
DO 20 I = 1 , NBASE
DO 20 J = 1 , NBASE
DO 20 K = 1 , NBASE
DERIVA( I , J ) = DERIVA( I , J) - A2 ( I , K) * DETA(K , J)
20 CONTINUE
END
C 计算形函数中矩阵 C导数的子程序
SUBROUTINE DERIVB1 ( II , COX , DERIVB, DEWEIG , LNODS,
&NOCE1 )
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COX( NINFL , NDIME) , LNODS(NINFL)
DIMENSION DEWEIG( NDIME , NINFL)
DIMENSION DERIVB(NBASE , NINFL)
DO 30 I = 1 , NBASE
DO 30 J = 1 , NPOIN
DERIVB( I , J ) = 0.0
30 CONTINUE
DO 10 I = 1 , NOCE1
DERIVB( 1, I ) = DEWEIG( II , I )
DERIVB( 2, I ) = COX( I ,1) * DEWEIG( II , I )
DERIVB( 3, I ) = COX( I ,2) * DEWEIG( II , I )
10 CONTINUE
END
C 计算形函数导数各项的子程序
SUBROUTINE PDB(P , D,B, B2 , S, NOCE1)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION P(NBASE) ,D(NBASE , NBASE) ,B(NBASE , NPOIN) ,
&B2(NBASE , NINFL) , S(NPOIN)
DO 5 I = 1, NBASE
DO 5 J = 1, NPOIN
B2( I , J) = 0 . 0
5 CONTINUE
DO 12 I = 1 , NPOIN
S( I ) = 0
12 CONTINUE
DO 10 I = 1 , NBASE
DO 10 J = 1 , NOCE1
DO 10 K = 1 , NBASE
B2( I , J) = B2(I,J) + D( I , K) * B(K , J )
10 CONTINUE
DO 20 I = 1 , NOCE1
DO 20 J = 1 , NBASE
S( I ) = S( I ) + P( J ) * B2( J, I )
20 CONTINUE
END
C 计算拉格朗日乘子法中矩阵 G 及组装子程序
SUBROUTINE GNIK(COORD, COX, DSHAPE , DIMAX , ESTIFQ,
&GIK , LNODS , LBNOD , SURPG , SHAB, SHAPE , POSGP , WEIGP,
&B2, WEIGHT , NOFIX, IFPRE)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COORD( NPOIN , NDIME) , LNODS(NINFL) ,
&LBNOD(NCELL , NNODE)
DIMENSION NOFIX( NVFIX) , IFPRE(NVFIX, NDOFN)
DIMENSION SHAPE(NINFL) , POSGP(NGAUS) ,WEIGP(NGAUS) ,
&WEIGHT(NINFL) , SH AB( NBASB, NINFL) , B2 (NBASB, NVFIX)
DIMENSION DSH APE( NDOFN , NLAGR) ,DINPOL( 2,4 )
DIMENSION ESTIFQ( NNMAX, NNMAX)
DIMENSION GIK( NLAGR, NLAGR)
DIMENSION COX( NINFL , NDIME)
DO 70 ICELL = 1 , NVFIX - 1
DO 20 I = 1 , NLAGR
DO 20 J = 1 , NLAGR
GIK( I , J) = 0.0
20 CONTINUE
NGAUS = 4
CALL GAUSSQ(POSGP ,WEIGP)
DO 50 IGAUS = 1 , NGAUS
EXISP = - 1 . 0
ETASP = POSGP( IGAUS)
CALL SFR1 (COORD,COX, DIMAX, SURPG, NOCE1, EXISP, ETASP
&ICELL , NOFIX, LNODS, LBNOD, WEIGHT , SHAB,SH APE ,B2)
CALL DSH APE1 (DSHAPE , SHAPE , IFPRE)
CALL DINPOL1(COORD, DINPOL , IFPRE , ETASP , ICELL, LBNOD,
&SURPG)
CALL GAUSSQ(POSGP ,WEIGP)
DVOLU = WEIGP( IGAUS) * (SURPG/ 2 . 0)
DO 30 IEVAB = 1, NDIME * NOCE1
DO 30 JEVAB = 1 ,4
DO 30 ISTRE = 1 ,2
GIK( IEVAB, JEVAB + 2 * ( ICELL - 1) ) = GIK( IEVAB, JEVAB +
&2 * ( ICELL - 1) ) - DSH APE( ISTRE , IEVAB) *
&DINPOL( IST RE , JEVAB) * DVOLU
30 CONTINUE
50 CONTINUE
C * * * * * * * * * * * * * * *
WRITE(4 ,1 ) ( ( GIK( I , J) , J = 1 , NLAGR) , I = 1 , NLAGR)
1 FORMAT (1X,′ GIK =′ ,10E10 . 3)
WRITE(4 , * )
C * * * * * * * * * * * * * * *
DO 80 I = 1 , NOCE1
DO 80 J = 1 , NVFIX
ESTIFQ(2 * LNODS( I ) - 1 ,2 * J - 1 + NEVAB) = ESTIFQ(2 * LNODS
&( I ) - 1 ,2 * J - 1 + NEVAB) + GIK( 2 * I - 1 ,2 * J - 1)
ESTIFQ(2 * LNODS( I ) - 1 ,2 * J + NEVAB) = ESTIFQ( 2 * LNODS
&( I ) - 1 ,2 * J + NEVAB) + GIK(2 * I - 1 ,2 * J)
ESTIFQ(2 * LNODS( I ) ,2 * J - 1 + NEVAB) = ESTIFQ( 2 * LNODS
&( I ) ,2 * J - 1 + NEVAB) + GIK(2 * I , 2 * J - 1)
ESTIFQ(2 * LNODS( I ) ,2 * J + NEVAB) = ESTIFQ(2 * LNODS
&( I ) ,2 * J + NEVAB) + GIK(2 * I ,2 * J )
80 CONTINUE
70 CONTINUE
DO 90 I = 1 , NEVAB
DO 90 J = NEVAB+ 1 , NNMAX
ESTIFQ( J, I ) = ESTIFQ( I , J)
90 CONTINUE
RETURN
END
C 位移边界结点形函数计算子程序
SUBROUTINE SFR1(COORD, COX, DIMAX, SURPG, NOCE1 , S,
&T , ICELL , NOFIX, LNODS, LBNOD,WEIGHT , SHAB,SHAPE,B2)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COORD( NPOIN , NDIME) , NOFIX(NVFIX) ,
&LNODS(NINFL) , LBNOD(NCELL , NNODE)
DIMENSION SHAPE(NINFL) ,WEIGHT (NINFL) ,
&DETA(NBASB, NBASB) , SHAB(NBASB, NINFL)
DIMENSION PT (NBASB) ,B2(NBASB, NVFIX)
DIMENSION COX( NINFL , NDIME)
DO 10 I = 1 , NPOIN
SHAPE( I ) = 0 . 0
10 CONTINUE
DO 15 I = 1 , NPOIN
WEIGHT( I ) = 0 . 0
15 CONTINUE
SURPD = DIMAX * SURPG
CALL CHAGAU1(COORD, S , T , EGACX , EGACY , EX , EY ,
&LBNOD , SURPG, ICELL)
CALL SEARCH1(COORD, EX, EY, SURPD, NOCE1 , LNODS, NOFIX)
c-----------------------------------------------------
C2 FORMAT( 1X,′ NOCE1 =′ ,1X , I3
C WRITE(4 ,2) NOCE1
c-----------------------------------------------------
CALL CHACOORD1(COORD,COX, NOCE1, ICELL , LNODS, LBNOD)
PT(1 ) = 1 . 0
PT(2 ) = EGACY
C WRITE(4 ,400) ICELL
C400 FORMAT(1X,′ ICELL =′ , I3)
c-----------------------------------------------------
DO 20 INODE = 1 , NOCE1
DEGAC = SQRT( ( EGACX - COX( INODE , 1) ) * * 2 +
&( EGACY - COX( INODE ,2) ) * * 2)
WEIGHT( INODE) = WEIGHTT (DEGAC, SURPD, SURPG)
20 CONTINUE
c-------------------------------------------------------
C WRITE(4,100) (WEIGHT ( I ) , I = 1, NOCE1)
C100 FORMAT( 1X,′ WEIGHT =′ ,100F8.5)
CALL DETA2(COX,WEIGHT , LNODS , NOCE1, DETA)
CALL SH AB2 (COX , LNODS, NOCE1, WEIGHT , SHAB)
CALL PDB1 (PT , SHAPE ,B2 ,DETA , SHAB)
c-----------------------------------------------------
C WRITE(4,1001) (SHAPE( I ) , I = 1, NOCE1 )
C1001 FORMAT(1X,′ SHAPE( I ) =′ ,100F8 . 4)
c------------------------------------------------------
RETURN
END
C 将位移边界上的模型坐标转换为边界积分背景网格的局部坐标子程序
SUBROUTINE CHACOORD1(COORD,COX, NOCE1 , ICELL , LNODS,
&LBNOD)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COORD( NPOIN , NDIME) , LNODS(NINFL) ,
&COX( NINFL , NDIME)
DIMENSION LBNOD(NCELL ,4)
DO I = 1 , NINFL
DO J = 1, NDIME
COX( I , J) = 0 .
ENDDO
ENDDO
DO I = 1 , NOCE1
IORNO = ABS( LBNOD( ICELL , 1) )
IORNI = ABS( LBNOD( ICELL ,2) )
CENX = (COORD( IORNI ,1) + COORD( IORNO ,1 ) )/ 2 .
CENY = (COORD( IORNI ,2 ) + COORD( IORNO ,2) )/ 2 .
COX( I ,1) = COORD( LNODS( I ) ,1) - CENX
COX( I ,2) = COORD( LNODS( I ) ,2) - CENY
ENDDO
RETURN
END
C 将位移边界上标准单元高斯点坐标转换为模型坐标子程序
SUBROUTINE CH AGAU1 (COORD , S, T , EGACX, EGACY, EX,
&EY , LBNOD, SURPG, ICELL)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COORD(NPOIN, NDIME) , LBNOD(NCELL , NNODE)
EGACX = ( (S / 2 . 0 ) * SURPG)
EGACY = ( ( T/ 2 . 0) * SURPG)
IORNO = ABS( LBNOD( ICELL , 2) )
EX = ( ( ( 1 + S)/ 2 . 0 ) * SURPG) + COORD( IORNO ,1 )
EY = ( ( (1 + T)/ 2 . 0) * SURPG) + COORD( IORNO ,2)
C WRITE(4 ,101) EGACX , EGACY
C101 FORMAT(1X,′ EGACX =′ , F10 . 4 ,2X ,′ EGACY =′ , F10 . 4 )
RETURN
END
C 位移边界上高斯点形成的权函数影响域内的结点搜索子程序
SUBROUTINE SEARCH1(COORD, EX, EY,DM,M, LNODS, NOFIX)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COORD( NPOIN , NDIME) , NOFIX(NVFIX) ,
&LNODS(NINFL)
M = 0
DO 20 I = 1 , NPOIN
LNODS( I ) = 0
20 CONTINUE
DO 10 IPOIN = 1 , NVFIX
DEGACA = SQRT ( ( EX - COORD( NOFIX( IPOIN) ,1 ) ) * * 2 +
&( EY - COORD( NOFIX( IPOIN) ,2) ) * * 2)
IF(DEGACA .LE .DM) THEN
M = M + 1
LNODS( M) = NOFIX( IPOIN)
END IF
10 CONTINUE
C WRITE(4 ,12 ) ( LNODS( I ) , I = 1 , NVFIX)
C12 FORMAT (1X,′ LNODS =′ ,1X ,20I3 )
RETURN
END
C 位移边界结点生成的矩阵 A 的逆矩阵计算子程序
SUBROUTINE DETA2 (COX ,WEIGHT , LNODS, NOCE1 ,DETA)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COX( NINFL , NDIME) , LNODS(NINFL)
DIMENSION WEIGHT(NINFL) , DETA(NBASB, NBASB)
DIMENSION AA( NBASB, NBASB)
DO 10 I = 1 , NBASB
DO 10 J = 1 , NBASB
AA( I , J) = 0 . 0
10 CONTINUE
DO 20 I = 1 , NOCE1
AA(1 ,1) = AA(1 ,1 ) + WEIGHT( I )
AA(1 ,2) = AA(1 ,2 ) + COX( I , 2) * WEIGHT( I )
AA(2 ,1) = AA(2 ,1 ) + COX( I , 2) * WEIGHT( I )
AA(2 ,2) = AA(2 ,2 ) + COX( I , 2) * COX( I ,2) * WEIGHT ( I )
20 CONTINUE
CALL AMATX1( AA, DETA)
END
SUBROUTINE AMATX1 (AA ,DETA)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION AA( NBASB, NBASB) , AMAX(NBASB, NBASB) ,
&DETA(NBASB, NBASB)
DO 5 I = 1, NBASB
DO 5 J = 1, NBASB
DETA( I , J ) = 0 . 0
5 CONTINUE
AMAX( 1,1 ) = AA(2 ,2)
AMAX( 1,2 ) = - 1 * AA(2 ,1)
AMAX( 2,1 ) = - 1 * AA(1 ,2)
AMAX( 2,2 ) = AA(1 ,1)
DDA1 = AA( 1,1 ) * AA(2 ,2)
DDA2 = AA( 1,2 ) * AA(2 ,1)
DDA = DDA1 - DDA2
DO 10 I = 1 , NBASB
DO 10 J = 1 , NBASB
DETA( I , J ) = AMAX( J , I )/ DDA
10 CONTINUE
END
C 位移边界结点形函数中矩阵 C的计算子程序
SUBROUTINE SHAB2(COX, LNODS, NOCE1 ,WEIGHT , SHAB)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COX( NINFL , NDIME) , LNODS(NINFL)
DIMENSION WEIGHT(NINFL) , SHAB( NBASB, NINFL)
DO 5 I = 1, NBASB
DO 5 J = 1, NPOIN
SHAB( I , J) = 0 . 0
5 CONTINUE
DO 10 I = 1 , NOCE1
SHAB(1 , I ) = WEIGHT ( I )
SHAB(2 , I ) = COX( I ,2 ) * WEIGHT( I )
10 CONTINUE
END
C 位移边界形函数计算子程序
SUBROUTINE PDB1( P, S ,B2, DETA, SH AB)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION P(NBASB) ,B2(NBASB, NVFIX) , S(NVFIX) ,
&DETA(NBASB, NBASB) , SHAB(NBASB, NINFL)
DO 5 I = 1, NBASB
DO 5 J = 1, NVFIX
B2( I , J) = 0 . 0
5 CONTINUE
DO 30 I = 1 , NVFIX
S( I ) = 0 . 0
30 CONTINUE
DO 10 I = 1 , NBASB
DO 10 J = 1 , NVFIX
DO 10 K = 1 , NBASB
B2( I , J) = B2( I , J) + DETA( I , K) * SHAB(K , J)
10 CONTINUE
DO 20 I = 1 , NVFIX
DO 20 J = 1 , NBASB
S( I ) = S( I ) + P( J ) * B2( J, I )
20 CONTINUE
END
C 位移边界上结点形函数的组装子程序
SUBROUTINE DSHAPE1(DSHAPE , SHAPE , IFPRE)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION SHAPE(NINFL) , IFPRE(NVFIX, NDOFN)
DIMENSION DSH APE( NDOFN , NLAGR)
DO 10 I = 1 , NDOFN
DO 10 J = 1 , NLAGR
DSHAPE( I , J) = 0 . 0
10 CONTINUE
DO 20 I = 1 , NVFIX
DSHAPE(1 ,2 * I - 1) = IFPRE( I ,1) * SHAPE( I )
DSHAPE(2 ,2 * I ) = IFPRE( I ,2 ) * SHAPE( I )
20 CONTINUE
END
C 位移边界上的拉格朗日插值函数计算子程序
SUBROUTINE DINPOL1(COORD, DINPOL , IFPRE , T , ICELL ,
&LBNOD , SURPG)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COORD( NPOIN , NDIME) , IFPRE( NVFIX, NDOFN) ,
&LBNOD(NCELL , NNODE)
DIMENSION AINPOL(2 ) , DINPOL(2 ,4)
CALL CHAGAU1(COORD, - 1 . 0 , T , EGACX , EGACY , EX , EY ,
&LBNOD , SURPG, ICELL)
CALL INPOL1 (AINPOL , COORD, EY, ICELL , LBNOD)
DO 10 I = 1 ,2
DO 10 J = 1 ,4
DINPOL( I , J) = 0 . 0
10 CONTINUE
DO 20 I = 1 ,2
DINPOL( 1, 2 * I - 1 ) = IFPRE( I + ( ICELL - 1) ,1) * AINPOL( I )
DINPOL( 2, 2 * I ) = IFPRE( I + ( ICELL - 1 ) , 2) * AINPOL( I )
20 CONTINUE
END
SUBROUTINE INPOL1( AINPOL ,COORD , T , ICELL , LBNOD)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COORD( NPOIN , NDIME) , LBNOD(NCELL ,4)
DIMENSION AINPOL(2 )
AINPOL( 1) = ( T - COORD( LBNOD( ICELL , 2) ,2) )/
&(COORD( LBNOD( ICELL ,1) ,2) - COORD(LBNOD( ICELL ,2) ,2) )
AINPOL( 2) = ( T - COORD( LBNOD( ICELL , 1) ,2) )/
&(COORD( LBNOD( ICELL ,2) ,2) - COORD(LBNOD( ICELL ,1) ,2) )
10 CONTINUE
C WRITE(4 ,100) ( AINPOL( I ) , I = 1, 2)
C100 FORMAT(1X,′ AINPOL =′ ,2F10 .4)
END
C 结点载荷计算子程序
SUBROUTINE LOADPS(COORD,COX , DIMAX , LNODS, LBNOD,
&SH APE , POSGP , WEIGP , WEIGHT , SHAB, QLOAD , SURPG ,
&B2, NOFIX , ELOADQ)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION COORD( NPOIN , NDIME) , LNODS(NINFL) ,
&LBNOD(NCELL , NNODE)
DIMENSION SHAPE(NINFL) , POSGP(NGAUS) ,WEIGP(NGAUS) ,
&WEIGHT(NINFL) , SH AB( NBASB, NINFL) , NOFIX(NVFIX)
DIMENSION B2( NBASB, NVFIX)
DIMENSION ELOADQ(NNMAX)
DIMENSION SHP(NVFIX)
DIMENSION COX( NINFL , NDIME)
DO I = 1 , NVFIX
SHP( I ) = 0 . 0
END DO
DO 60 I = 1 , NNMAX
ELOADQ( I ) = 0 . 0
60 CONTINUE
DO 70 ICELL = 1 , NVFIX - 1
NGAUS = 4
CALL GAUSSQ(POSGP ,WEIGP)
DO 50 IGAUS = 1 , NGAUS
EXISP = - 1 . 0
ETASP = POSGP( IGAUS)
CALL SFR1(COORD,COX ,DIMAX ,SURPG , NOCE1, EXISP ,
&ETASP , ICELL , NOFIX, LNODS, LBNOD ,WEIGHT , SHAB,
&SH APE ,B2)
DO I = 1 , NOCE1
SHP( LNODS( I ) ) = SH APE( I )
END DO
CALL GAUSSQ(POSGP ,WEIGP)
DVOLU = WEIGP( IGAUS) * (SURPG/ 2 . 0)
DO 30 IEVAB = 1, NVFIX
ELOADQ( 2 * ( IEVAB + 6 ) - 1) = ELOADQ(2 * ( IEVAB + 6) - 1 ) +
&SHP( IEVAB) * DVOLU * QLOAD
30 CONTINUE
50 CONTINUE
70 CONTINUE
WRITE(5 ,1 ) ( ELOADQ( I ) , I = 1 , NNMAX)
1 FORMAT (1X ,′ ELOADQ =′ , E10 . 4 )
RETURN
END
C 高斯选主元消去法求位移计算子程序
SUBROUTINE GAUDIP( ESTIFQ , ELOADQ ,DIP)
COMMON/ CONTRO/ NPOIN , NCELL , NDOFN , NDIME , NSTRE ,
&NTYPE , NGAUS, NNODE , NPROP , NVFIX, NBNOD, NEVAB,
&NNMAX , NECELX , NECELY, NGASP , NLAGR, NINFL , NBASE ,
&NBASB, NDIFL
DIMENSION ESTIFQ( NNMAX, NNMAX) , ELOADQ( NNMAX) ,
&DIP(NNMAX)
DO 10 N = 1 , NNMAX - 1
DO 20 I = N + 1, NNMAX
IF(ABS( ESTIFQ( I , N) ) .GE .ABS( ESTIFQ( N, N) ) ) THEN
TT = ELOADQ( N)
ELOADQ( N) = ELOADQ( I )
ELOADQ( I ) = TT
DO 30 J = 1 , NNMAX
TT = ESTIFQ(N , J)
ESTIFQ(N , J) = ESTIFQ( I , J)
ESTIFQ( I , J) = TT
30 CONTINUE
END IF
20 CONTINUE
DO 50 I = N + 1, NNMAX
DO 40 J = N + 1, NNMAX
ESTIFQ( I , J) = ESTIFQ( I , J) - ( ESTIFQ( I , N)/ ESTIFQ(N , N) ) *
&ESTIFQ( N, J )
40 CONTINUE
ELOADQ( I ) = ELOADQ( I ) - (ESTIFQ( I , N)/ ESTIFQ( N , N) ) *
&ELOADQ(N)
ESTIFQ( I , N) = 0 .0
50 CONTINUE
10 CONTINUE
DO 60 I = NNMAX,1 , - 1
A = 0 . 0
DO 70 J = I + 1 , NNMAX
A = A + ESTIFQ( I , J) * DIP( J)
70 CONTINUE
C WRITE(6 , * )A , ESTIFQ( I , I ) , ELOADQ( I )
DIP( I ) = ( ELOADQ( I ) - A)/ ESTIFQ( I , I )
60 CONTINUE
WRITE(6 ,100)
100 FORMAT(/ 1X ,′- DISPLACEMENTS - ′ )
DO 80 I = 1 , NNMAX/ 2
WRITE(6 ,200) I , DIP(2 * I - 1) ,DIP(2 * I )
200 FORMAT( 5X, I3, 3X, E12 . 6, 2X, E12 . 6)
80 CONTINUE
END

weikuangming

顶 11111!!!!!!!!!

王雨森

拷回去看看，顶起来！

jinsaixiang

从哪抄袭的？