input.f90

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

        SUBROUTINE Input
        USE CtrlData
        IF( nfmtc .EQ. 0 ) CALL FixInput
        IF( nfmtc .NE. 0 ) CALL MacroInput
        RETURN
        END

        SUBROUTINE FixInput
        CALL RdCrd
        CALL RdLnk
        CALL RdElm
        CALL RdFix
        CALL RdMat
        CALL RdFrc
        CALL RdCnt
        CALL RdSlnt
        CALL RdOptAxis
        RETURN
        END

        SUBROUTINE MacroInput
        USE CtrlData
        USE MeshData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        CHARACTER strng * 80, mwork * 10
        IF( nerrc .GT. 0 ) RETURN

        label = 1
        ifcrd = 0
        ifelm = 0
        iffix = 0
        ifmat = 0
        iffrc = 0
        ifcnt = 0
        ifsln = 0
        iflnk = 0
        ifaut = 0
        ifmid = 0

!.......节点连接信息初时化.
        DO ipoin = 1, npoin
          llnks( ipoin ) = ipoin
        END DO

        DO WHILE( label .EQ. 1 )
          READ( 11, '( a80 )', END = 400 ) strng
          CALL ReLocate( strng, nerrc )
          IF( nerrc .GT. 0 ) RETURN
          mwork = strng( 1:10 )
          IF( mwork .EQ. 'end       ' ) THEN
            RETURN
          ELSE IF( mwork .EQ. 'coordinate' ) THEN
            CALL RdCrd
            ifcrd = ifcrd + 1
          ELSE IF( mwork .EQ. 'element   ' ) THEN
            CALL RdElm
            ifelm = ifelm + 1
          ELSE IF( mwork .EQ. 'fix       ' ) THEN
            CALL RdFix
            iffix = iffix + 1
          ELSE IF( mwork .EQ. 'material  ' ) THEN
            CALL RdMat
            ifmat = ifmat + 1
          ELSE IF( mwork .EQ. 'force     ' ) THEN
            CALL RdFrc
            iffrc = iffrc + 1
          ELSE IF( mwork .EQ. 'link      ' ) THEN
            CALL RdLnk
            iflnk = iflnk + 1
          ELSE IF( mwork .EQ. 'center    ' ) THEN
            CALL RdCnt
            ifcnt = ifcnt + 1
          ELSE IF( mwork .EQ. 'slent     ' ) THEN
        CALL RdSlnt
            ifsln = ifsln + 1
          ELSE IF( mwork .EQ. 'autolink  ' ) THEN
            CALL AutoLink( strng )
            ifaut = ifaut + 1
          ELSE IF( mwork .EQ. 'midnode   ' ) THEN
            CALL MidNode
            ifmid = ifmid + 1
          END IF
          IF( nerrc .NE. 0 ) RETURN
        END DO
        RETURN
400     WRITE( 12, 2000 )
        nerrc = 3004
        RETURN
2000    FORMAT( '' )
        END

        SUBROUTINE RdCrd
!........
!       模块功能
!           读写节点坐标.
!........
        USE CtrlData
        USE MeshData
        USE AutoMeshData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION corpi( 3 ), corpj( 3 ), cstep( 3 )
        IF( nerrc .GT. 0 ) RETURN
        pai = datan( 1.0D0 ) / 45.0D0
        CALL ShowMessage( '读节点坐标...' )

        jgene = 0
        iflag = 1
        DO WHILE( iflag .EQ. 1 )
         READ( 11, 2000 ) ipoin, ictyp, igene, corpi, tempi, vdimi
         IF( ipoin .EQ. 0 .OR. ipoin .GT. npoin ) THEN
           iflag = 0
         ELSE
          llnks( ipoin ) = llnks( ipoin ) + 1
          temps( ipoin ) = tempi
          vdims( ipoin ) = vdimi
          IF( ictyp .EQ. 0 ) THEN
           DO idime = 1, ndime
            coord( idime, ipoin ) = corpi( idime )
           END DO
          ELSE IF( ictyp .EQ. 1 ) THEN
           theta = corpi( 2 ) * pai
           coord( 1, ipoin ) = corpi( 1 ) * dcos( theta )
           coord( 2, ipoin ) = corpi( 1 ) * dsin( theta )
           IF( ndime .EQ. 3 ) coord( 3, ipoin ) = corpi( 3 )
          ELSE IF( ictyp .EQ. 2 ) THEN
           alpha = corpi( 2 ) * pai
           theta = corpi( 3 ) * pai
           coord( 1, ipoin ) = corpi(1) * dsin( alpha ) * dcos( theta )
           coord( 2, ipoin ) = corpi(1) * dsin( alpha ) * dsin( theta )
           coord( 3, ipoin ) = corpi(1) * dcos( alpha )
          ELSE
           WRITE( 12, 2200 ) ictyp, ipoin
           nerrc = 1001
           RETURN
          END IF
          IF( jgene .NE. 0 ) THEN
           IF( ictyp .NE. jctyp ) THEN
            WRITE( 12, 2400 ) jpoin, ipoin
            nerrc = 1002
            RETURN
           END IF
           nsegm = ( ipoin - jpoin ) / jgene
           IF( nsegm * jgene .NE. ipoin - jpoin ) THEN
            WRITE( 12, 2600 ) jgene, jpoin, ipoin
            nerrc = 1003
            RETURN
           END IF
           IF( nsegm .NE. 0 ) THEN
            tmpst = ( tempi - tempj ) / nsegm
            dimst = ( vdimi - vdimj ) / nsegm
            DO idime = 1, ndime
             cstep( idime ) = (corpi( idime ) - corpj( idime )) / nsegm
            END DO
            DO kpoin = jpoin + jgene, ipoin - jgene, jgene
             llnks( kpoin ) = llnks( kpoin ) + 1
             tempj = tempj + tmpst
             vdimj = vdimj + dimst
             DO idime = 1, ndime
              corpj( idime ) = corpj( idime ) + cstep( idime )
             END DO
             temps( kpoin ) = tempj
             vdims( kpoin ) = vdimj
             IF( jctyp .EQ. 0 ) THEN
              DO idime = 1, ndime
               coord( idime, kpoin ) = corpj( idime )
              END DO
             ELSE IF( jctyp .EQ. 1 ) THEN
              theta = corpj( 2 ) * pai
              coord( 1, kpoin ) = corpj( 1 ) * dcos( theta )
              coord( 2, kpoin ) = corpj( 1 ) * dsin( theta )
              IF( ndime .EQ. 3 ) coord( 3, kpoin ) = corpj( 3 )
             ELSE IF( jctyp .EQ. 2 ) THEN
              alpha = corpj( 2 ) * pai
              theta = corpj( 3 ) * pai
              coord( 1, kpoin ) = corpj( 1 ) * dsin( alpha ) *                 &
                                               dcos( theta )
              coord( 2, kpoin ) = corpj( 1 ) * dsin( alpha ) *                 &
                                               dsin( theta )
              coord( 3, kpoin ) = corpj( 1 ) * dcos( alpha )
             END IF
            END DO
           END IF
          END IF
          jpoin = ipoin
          jgene = igene
          jctyp = ictyp
          tempj = tempi
          vdimj = vdimi
          corpj( 1 ) = corpi( 1 )
          corpj( 2 ) = corpi( 2 )
          corpj( 3 ) = corpi( 3 )
         END IF
        END DO

        DO ipoin = 1, npoin
         IF( llnks( ipoin ) .EQ. 0 ) THEN
           WRITE( 12, 2800 ) ipoin
         ELSE IF( llnks( ipoin ) .GT. 1 ) THEN
           WRITE( 12, 3000 ) ipoin
         END IF
         llnks( ipoin ) = 0
        END DO

        IF( nprnc .NE. 0 ) THEN
         WRITE( 12, 3200 )
          corpi( 1 ) = 0.0D0
          corpi( 2 ) = 0.0D0
          corpi( 3 ) = 0.0D0
         DO ipoin = 1, npoin
          DO idime = 1, ndime
           corpi( idime ) = coord( idime, ipoin )
          END DO
          tempi = temps( ipoin )
          vdimi = vdims( ipoin )
          WRITE( 12, 3400 ) ipoin, corpi, tempi, vdimi
         END DO
        END IF
2000    FORMAT( 3i5, 5f15.5 )
2200    FORMAT( 2x, '错误: 无效的坐标类型', i2, ', 节点', i5 )
2400    FORMAT( 2x, '错误: 自动坐标不适用于不同坐标类型, ',                    &
                    '节点', i5, '和', i5 )
2600    FORMAT( 2x, '错误: 无效的自动步长', I5, '节点', i5, '和', i5 )
2800    FORMAT( 2x, '警告：节点', I5, ' 没有定义' )
3000    FORMAT( 2x, '警告：节点', I5, ' 重复定义' )
3200    FORMAT( //10x, '= = = = 节点坐标 = = = ='                              &
                // 3x, '节点    x-坐标       y-坐标       z-坐标'              &
                       '      节点温度     单元尺度' )
3400    FORMAT(   2X, I5, 5F13.5 )
        RETURN
        END

        SUBROUTINE RdLnk
!........
!       模块功能
!           读取节点连接信息.
!........
        USE CtrlData
        USE MeshData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION lwork( 10 )
        IF( nerrc .GT. 0 ) RETURN

!.......节点连接信息初时化.
        DO ipoin = 1, npoin
          llnks( ipoin ) = ipoin
        END DO
!.......坐标相同的节点自动连接.
        IF( nlnkc .NE. 0 ) THEN
          DO ipoin = 1, npoin
            IF( ipoin .EQ. llnks( ipoin ) ) THEN
              DO jpoin = ipoin + 1, npoin
                IF( jpoin .EQ. llnks( jpoin ) ) THEN
                  displ = DABS( coord(1, ipoin) - coord(1, jpoin) )
                  IF( displ .LT. 1.5D-5 ) THEN
                    nflag = 0
                    DO idime = 2, ndime
                      displ = DABS( coord( idime, ipoin ) -                    &
                                    coord( idime, jpoin ) )
                      IF( displ .GT. 1.5D-5 ) nflag = 1
                    END DO
                    IF( nflag .EQ. 0 ) llnks( jpoin ) = ipoin
                  END IF
                END IF
              END DO
            END IF
          END DO
        END IF
!.......读入节点连接信息
        nflag = 1
        DO WHILE( nflag .EQ. 1 )
          kpoin = 0
          ipoin = npoin
          READ( 11, 2000, ERR = 100 ) lwork
          DO inode = 1, 10
            jpoin = lwork( inode )
            IF( jpoin .NE. 0 ) THEN
              IF( jpoin .GT. npoin .OR. jpoin .LT. 0 ) THEN
                WRITE( 12, 2200 ) jpoin
                nerrc = 24356
                RETURN
              ELSE
                ipoin = jpoin
                kpoin = kpoin + 1
              END IF
            END IF
          END DO
          IF( kpoin .EQ. 0 ) THEN
            nflag = 0
          ELSE IF( kpoin .EQ. 1 ) THEN
            WRITE( 12, 2400 )
            nerrc = 24356
            RETURN
          ELSE
            DO inode = 1, 10
              jpoin = lwork( inode )
              IF( jpoin .GT. 0 .AND. jpoin .LE. npoin ) THEN
                DO WHILE( ipoin .NE. llnks( ipoin ) )
                  ipoin = llnks( ipoin )
                END DO
                DO WHILE( jpoin .NE. llnks( jpoin ) )
                  jpoin = llnks( jpoin )
                END DO
                IF( ipoin .LT. jpoin ) llnks( jpoin ) = ipoin
                IF( ipoin .GT. jpoin ) llnks( ipoin ) = jpoin
              END IF
            END DO
          END IF
        END DO
        IF( nprnc .NE. 0 ) THEN
          WRITE( 12, 2600 )
          DO ipoin = 1, npoin
            WRITE( 12, 2800 ) ipoin, llnks( ipoin )
          END DO
        END IF
        RETURN
100     WRITE( 12, 3000 )
        nerrc = 59759
        RETURN
2000    FORMAT( 10I5 )
2200    FORMAT( 2x, '错误: 无效节点 ', I5 )
2400    FORMAT( 2x, '错误: 至少有两个节点才能连接在一起.' )