input.f90

边界元程序,供力学工作者参考,希望对大家有所帮助

2600    FORMAT( ///2x, 60('=') // 20x, '节点连接信息'                          &
                 //2x, 60('=') /// )
2800    FORMAT( 16x, I5, 10x, I5 )
3000    FORMAT( //2x, '错误: 读节点连接信息错.' )
        END

        SUBROUTINE RdElm
!........
!       模块功能
!           读写单元信息.
!........
        USE CtrlData
        USE MeshData
        USE AutoMeshData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION lnodt( 10 ), ldimt( 3 ), lwork( 100 )

        IF( nerrc .GT. 0 ) RETURN
        CALL ShowMessage( '读单元定义...' )

        jelem = 0
        jgene = 0
        iflag = 1
        DO WHILE( iflag .EQ. 1 )
         READ( 11, 2000 ) ielem, imats, igene, lnodt, ldimt
         IF( ielem .EQ. 0 .OR. ielem .GT. nelem ) THEN
           iflag = 0
         ELSE
          IF( imats .LT. 1 .OR. imats .GT. nmats ) THEN
           WRITE( 12, 2200 ) imats, ielem
           nerrc = 1004
           RETURN
          END IF
          lmats( ielem ) = imats
          DO inode = 1, MIN0( 10, mnode )
           lnods( inode, ielem ) = lnodt( inode )
          END DO
          DO iloca = 11, mnode, 10
           READ( 11, 2000 ) izero, izero, izero, lnodt
           DO inode = iloca, MIN0( iloca + 9, mnode )
            lnods( inode, ielem ) = lnodt( inode - iloca + 1 )
           END DO
          END DO
          DO idime = 1, ndime
           ldims( idime, ielem ) = ldimt( idime )
          END DO
          IF( jgene .NE. 0 ) THEN
           IF( imats .NE. jmats ) THEN
            WRITE( 12, 2400 ) jelem, ielem
            nerrc = 1005
            RETURN
           END IF
           nsegm = ( ielem - jelem ) / jgene
           IF( nsegm * jgene .NE. ielem - jelem ) THEN
            WRITE( 12, 2600 ) jgene, jelem, ielem
            nerrc = 1006
            RETURN
           END IF

           DO inode = 1, mnode
            lwork( inode ) = ( lnods( inode, ielem ) -                         &
                               lnods( inode, jelem ) ) / nsegm
            IF( lwork( inode ) * nsegm .NE. lnods( inode, ielem ) -            &
                lnods( inode, jelem ) ) THEN
             WRITE( 12, 2800 ) inode, jelem, ielem
             nerrc = 1007
             RETURN
            END IF
           END DO

           DO idime = 1, ndime
            ldimt( idime ) = ( ldims( idime, ielem ) -                         &
                               ldims( idime, jelem ) ) / nsegm
            IF( ldimt( idime ) * nsegm .NE. ldims( idime, ielem ) -            &
                ldims( idime, jelem ) ) THEN
             WRITE( 12, 2900 ) idime, jelem, ielem
             nerrc = 1007
             RETURN
            END IF
           END DO

           DO kelem = jelem + jgene, ielem - jgene
            lmats( kelem ) = jmats
            DO inode = 1, mnode
             lnods( inode, kelem ) = lnods( inode, kelem - jgene ) +           &
             lwork( inode )
            END DO
            DO idime = 1, ndime
             ldims( idime, kelem ) = ldims( idime, kelem - jgene ) +           &
             ldimt( idime )
            END DO
           END DO
          END IF
          jelem = ielem
          jmats = imats
          jgene = igene
         END IF
        END DO

        DO ielem = 1, nelem
          IF( lmats( ielem ) .EQ. 0 ) WRITE( 12, 3000 ) ielem
        END DO
        IF( nprnc .NE. 0 ) THEN
         WRITE( 12, 3200 )
         DO ielem = 1, nelem
          iwrot = 0
          imats = lmats( ielem )
          DO idime = 1, ndime
           ldimt( idime ) = ldims( idime, ielem )
          END DO
          DO inode = 1, mnode
           iloca = MOD( inode, 10 )
           IF( iloca .EQ. 0 ) iloca = 10
           lnodt( iloca ) = lnods( inode, ielem )
           IF( inode .EQ. mnode ) THEN
            DO jloca = iloca + 1, 10
             lnodt( jloca ) = 0
            END DO
           END IF
           IF( iloca .EQ. 10 .OR. inode .EQ. mnode ) THEN
            iwrot = iwrot + 1
            IF( iwrot .EQ. 1 ) THEN
             WRITE( 12, 3400 ) ielem, imats, lnodt, ldimt
            ELSE
             WRITE( 12, 3600 ) lnodt
            END IF
           END IF
          END DO
         END DO
        END IF
        RETURN
2000    FORMAT( 16I5 )
2200    FORMAT( 2x, '错误: 无效的材料号', i2, ', 单元', i5 )
2400    FORMAT( 2x, '错误: 不同类型单元', I5, '和',                            &
                I5, '间不可自动生成' )
2600    FORMAT( 2x, '错误: 无效的单元步长', i5,                                &
                    ', 单元', i5, '和', i5 )
2800    FORMAT( 2x, '错误: 无效的自动步长, 节点', i5,                          &
                    ',  单元', i5, '和', i5 )
2900    FORMAT( 2x, '错误: 无效的自动步长, 维数', i5,                          &
                    ',  单元', i5, '和', i5 )
3000    FORMAT( 2x, '警告：单元', I5, ' 没有定义' )
3200    FORMAT( // 18x, ' = = = = 单 元 信 息 = = = = ' //                     &
                   2x, ' 单元 材料   节点  节点  节点  节点  节点',            &
                                  '  节点  节点  节点  节点  节点'             &
                                  '  分段  分段  分段' )
3400    FORMAT(    2x, 2I5, 1x, 13I6 )
3600    FORMAT(   13x, 10I6 )
        END

        SUBROUTINE RdFix
!........
!       模块功能
!           读写位移约束信息.
!........
        USE CtrlData
        USE MeshData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )

        IF( nerrc .GT. 0 ) RETURN
        CALL ShowMessage( '读约束信息...' )

        DO ipoin = 1, npoin
          DO idofn = 1, mdofn
            lfixs( idofn, ipoin ) = 0
            vfixs( idofn, ipoin ) = 0.0D0
          END DO
        END DO
        label = 1
        DO WHILE( label .EQ. 1 )
          READ( 11, 2000, ERR = 200 ) ipoin, idofn, itype, vcord, value
          IF( ipoin .EQ. 0 ) THEN
            label = 0
          ELSE
            IF( idofn .LE.     0 ) nerrc = 1009
            IF( ipoin .LT.-ndime ) nerrc = 1008
            IF( ipoin .GT. npoin ) nerrc = 1008
            IF( idofn .GT. mdofn ) nerrc = 1009
            IF( nerrc .EQ.  1008 ) WRITE( 12, 2200 ) ipoin
            IF( nerrc .EQ.  1009 ) WRITE( 12, 2400 ) idofn, ipoin
            IF( nerrc .GT. 0 ) RETURN
            IF( ipoin .GT. 0 ) THEN
              ipoin = llnks( ipoin )
              lfixs( idofn, ipoin ) = itype + 1
              vfixs( idofn, ipoin ) = value
            ELSE
              idime =-ipoin
              DO ipoin = 1, npoin
                vdisp = DABS( coord( idime, ipoin ) - vcord )
                IF( vdisp .LT. 1.0D-8 ) THEN
                  lfixs( idofn, llnks( ipoin ) ) = itype + 1
                  vfixs( idofn, llnks( ipoin ) ) = value
                END IF
              END DO
            END IF
          END IF
        END DO
        IF( nprnc .NE. 0 ) THEN
          WRITE( 12, 2600 )
          DO ipoin = 1, npoin
            DO idofn = 1, mdofn
              IF( lfixs( idofn, ipoin ) .NE. 0 ) THEN
                value = vfixs( idofn, ipoin )
                itype = lfixs( idofn, ipoin ) - 1
                WRITE( 12, 2800 ) ipoin, idofn, itype, value
              END IF
            END DO
          END DO
        END IF
        RETURN
200     WRITE( 12, 3000 )
        nerrc = 1010
        RETURN
2000    FORMAT( 3I5, 2F15.5 )
2200    FORMAT( 2x, '错误: 无效的节点号', i5 )
2400    FORMAT( 2x, '错误: 无效的自由度', I3, ', 节点', I5 )
2600    FORMAT( //8x, ' = = = = 位 移 约 束 信 息 = = = =' //                  &
                  5x, '   节点    自由度         给定位移 ' )
2800    FORMAT(   7x, i5, 6x, I2, 6x, I2, 8x, F15.5 )
3000    FORMAT( 2x, '错误: 读约束代码出错' )
        END

        SUBROUTINE RdMat
!........
!       模块功能
!           读写材料信息.
!........
        USE CtrlData
        USE MeshData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        ALLOCATABLE lprpt(:)

        ALLOCATE( lprpt( mdofn + 3 ), STAT = ierro )
        IF( ierro .NE. 0 ) nerrc = 4
        IF( nerrc .GT. 0 ) RETURN

        CALL ShowMessage( '读材料信息...' )
        IF( nprnc .NE. 0 ) WRITE( 12, 2000 )

        iswth = 0
        nhstr = 0
        nstrh = 0
        nmatp = mdofn + 3
        DO imats = 1, nmats
          READ( 11, 2200 ) jmats, lprpt
          IF( jmats .LT. 1 .OR. jmats .GT. nmats ) nerrc = 1011
          IF( nprnc .NE. 0 ) THEN
            WRITE( 12, 2400 ) jmats, lprpt( 1 ), lprpt( 2 )
            WRITE( 12, 2600 ) ( lprpt( imatp + 2 ), imatp = 1, mdofn )
          END IF
          IF( nerrc .GT. 0 ) DEALLOCATE( lprpt )
          IF( nerrc .GT. 0 ) RETURN
          ietyp = lprpt(     1 )
          ipara = lprpt( nmatp )
          jetyp = lprps(     1, jmats )
          jpara = lprps( nmatp, jmats )

          IF( jetyp .NE. 0 ) THEN
            IF( ietyp .NE. jetyp .OR. ipara .NE. jpara ) THEN
              WRITE( 12, 2800 ) jmats, jetyp, jpara, ietyp, ipara
              DEALLOCATE( lprpt )
              nerrc = 3401
              RETURN
            END IF
          END IF

          nchck = 0
          DO idofn = 1, mdofn
            nchck = nchck + lprpt( idofn + 2 )
          END DO
          IF( nchck .EQ. 0 ) THEN
            DO idofn = 1, mdofn
              lprpt( idofn + 2 ) = idofn
            END DO
          END IF
          DO imatp = 1, nmatp
            lprps( imatp, jmats ) = lprpt( imatp )
          END DO
          CALL ElmOpt( jmats, iswth )
          IF( nerrc .GT. 0 ) DEALLOCATE( lprpt )
          IF( nerrc .GT. 0 ) RETURN
        END DO
        DEALLOCATE( lprpt )
        RETURN
2000    FORMAT( //10x, '= = = = 材 料 信 息 = = = = ' )
2200    FORMAT( 12I5 )
2400    FORMAT( //2x, '材料类型序号', I5 / 2x, '单元类型序号', I5              &
                / 2x, '材料作用因子', I5 )
2600    FORMAT( 2x, '自由度序列表', 10I5 / 14x, 10I5 )
2800    FORMAT( //2x, '错误: 材料', I5, '应输入单元类型', I5, '-', I5,         &
                 /2x, '          ', 5x, '实输入单元类型', I5, '-', I5 )
        END

        SUBROUTINE RdFrc
!........
!       模块功能
!           读写载荷信息.
!........
        USE CtrlData
        USE LoadData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION fwork( 40 )

        IF( nerrc .GT. 0 ) RETURN
        CALL ShowMessage( '读载荷信息...' )

!.......读节点载荷.
        knodf = 0
        kedgf = 0
        kplnf = 0
        nnodg = mnodg - 1
        nnodp = mnodp - 1
        DO iload = 1, nload
          IF( nprnc .NE. 0 )                                                   &
          WRITE( 12, 2000 ) iload
          READ( 11, 2200 ) nnodf, nedgf, nplnf
          loads( 1, iload ) = nnodf
          loads( 2, iload ) = nedgf
          loads( 3, iload ) = nplnf
          IF( nnodf .GT. 0 ) THEN
            IF( nprnc .NE. 0 ) WRITE( 12, 2400 )