front.f90

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

        IF( lsavc( 3 ) .LT. 0 ) WRITE( 12, 3800 ) -lsavc( 3 )
        IF( lsavc( 3 ) .EQ. 1 ) WRITE( 12, 4000 )
        IF( lsavc( 3 ) .EQ. 2 ) WRITE( 12, 4200 )
        RETURN
2000    FORMAT( //2x, '致命错误: 内存缓冲区不足以存交换数据.' )
2200    FORMAT( //2x, '致命错误: 内外存缓冲区不足以存所有数据.', I10 )
2400    FORMAT( //20x, '一致矩阵的存储信息' // )
2600    FORMAT( 2x, '刚度阵存于内存, 首址', I10 )
2800    FORMAT( 2x, '刚度阵存于外存' )
3000    FORMAT( 2X, '集中刚度阵, 存于内存' )
3200    FORMAT( 2X, '质量阵存于内存, 首址', I10 )
3400    FORMAT( 2X, '质量阵存于外存' )
3600    FORMAT( 2X, '集中质量阵, 存于内存' )
3800    FORMAT( 2X, '阻尼阵存于内存, 首址', I10 )
4000    FORMAT( 2X, '阻尼阵存于外存' )
4200    FORMAT( 2x, '集中阻尼阵, 存于内存' )
        END

        SUBROUTINE SetRec
!........
!       模块功能
!           确定记录长度和记录对应的起终节点.
!........
!       lbcks 回代信息            lrcrd 记录信息
!       mdofn 最大自由度数              lrcrd( 1, ircrd ) 起始节点号
!       ndata 一致矩阵元素数            lrcrd( 2, ircrd ) 终止节点号
!       nexch 数据交换区尺寸            lrcrd( 3, ircrd ) 记录长度
!       nrcrd 一致矩阵记录数      npoin 节点总数
!........
        USE CtrlData
        USE FrontData
        IF( nerrc .GT. 0 ) RETURN

        IF( nexch .EQ. 0 ) THEN
!........如果为内存存储.
         nexch = ndata
         lrcrd( 1, 1 ) = 1
         lrcrd( 2, 1 ) = npoin
         lrcrd( 3, 1 ) = ndata
        ELSE
!........重新标定记录信息.
         kdata = 0
         krcrd = nrcrd
         DO ipoin = npoin, 1,-1
          npntw = lbcks( 1, ipoin )
          IF( npntw .GT. 0 ) THEN
           ndofw = npntw * mdofn
           idata = ndofw * mdofn - ( mdofn + 1 ) * mdofn / 2
           IF( kdata + idata .GT. nexch ) THEN
            lrcrd( 3, krcrd ) = kdata
            lrcrd( 1, krcrd ) = ipoin + 1
            krcrd = krcrd - 1
            kdata = idata
           ELSE
            kdata = kdata + idata
           END IF
          END IF
         END DO
         lrcrd( 1,     1 ) = 1
         lrcrd( 3,     1 ) = kdata
         lrcrd( 2, nrcrd ) = npoin
         DO ircrd = 1, nrcrd - 1
          krcrd = ircrd + 1
          lrcrd( 2, ircrd ) = lrcrd( 1, krcrd ) - 1
         END DO
        END IF
        WRITE( 12, 2000 )
        WRITE( 12, 2200 ) nrcrd
        WRITE( 12, 2400 )
        DO ircrd = 1, nrcrd
         WRITE( 12, 2600 )ircrd, lrcrd( 1, ircrd ),                            &
              lrcrd( 2, ircrd ), lrcrd( 3, ircrd )
        END DO
2000    FORMAT( // 20x, '外存记录信息' // )
2200    FORMAT( 2x, '一致阵外存记录数', I10 )
2400    FORMAT( 2x, '记录号  起始节点  终止节点  有效长度' )
2600    FORMAT( 2X, I6, 3( 2X, I8 ) )
        RETURN
        END

        SUBROUTINE DeComp( iswth, index )
!........
!       模块功能
!           一致矩阵的分解.
!........
!       buffs 缓冲区数组         diags 一致矩阵对角元
!       fbcks 工作数组           stifw 波阵区
!       lbcks 回代信息           nrcrd 一致矩阵记录总数
!       lnods 单元定义           lopts 优化单元顺序
!       lpnte 绕单元节点数       lpntw 波阵节点总体编号
!       lsavc 一致矩阵存储信息   lrcrd 记录信息
!       iswth 一致矩阵开关       mdofn 最大节点自由度
!       mdofw 最大波阵自由度     mnode 最大单元节点数
!       mpntw 最大波阵节点       mstfw 最大波阵容量
!       nbufs 缓冲区容量         ndofs 自由度总数
!       nelem 单元总数           nexch 数据交换区容量
!       npoin 节点总数
!........
        USE CtrlData
        USE MeshData
        USE GlobData
        USE SolvData
        USE FrontData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
!.......判断是否有必要分解.
        isavc = lsavc( iswth )
        IF( isavc .EQ. 0 ) RETURN
        IF( isavc .EQ. 2 ) RETURN
!.......参数初始化.
        knume = 0
        npntw = 0
        kcoun = 0
        DO ipntw = 1, mpntw
         lpntw( ipntw ) = 0
        END DO
        DO istfw = 1, mstfw
         stifw( istfw ) = 0.0D0
        END DO
!.......寻找最大对角元素
        dmaxv = 0.0D0
        IF( iswth .EQ. 1 ) THEN
          DO idofs = 1, ndofs
            dtemp = DABS( diagk( idofs ) )
            IF( dtemp .GT. dmaxv ) dmaxv = dtemp
          END DO
        ELSE IF( iswth .EQ. 2 ) THEN
          DO idofs = 1, ndofs
            dtemp = DABS( diagm( idofs ) )
            IF( dtemp .GT. dmaxv ) dmaxv = dtemp
          END DO
        ELSE IF( iswth .EQ. 3 ) THEN
          DO idofs = 1, ndofs
            dtemp = DABS( diagd( idofs ) )
            IF( dtemp .GT. dmaxv ) dmaxv = dtemp
          END DO
        END IF
!.......确定阀值，当主消元素小于该值时即认为该自由度为刚体位移，必须约束。
        IF( index .EQ. 0 ) index = -16
        IF( index .GT. 0 ) index = -index
        dswth = dmaxv * 10 ** index
!.......按记录顺序消元.
        DO ircrd = 1, nrcrd
         inump = lrcrd( 1, ircrd )
         jnump = lrcrd( 2, ircrd )
         IF( isavc .EQ. 1 ) THEN
          kdata = 0
          CALL OptRec( buffs, ircrd, 1, iswth, nexch )
         ELSE
          kdata =-isavc - 1
         END IF
         DO knump = inump, jnump
          kpoin = lbcks( 2, knump )
          kpntw = lbcks( 3, knump )
          kpnte = lpnte(    kpoin )
          kcoun = kcoun + 1
          IF( necho .NE. 0 ) CALL ShowProcess( kcoun, npoin )
!.........对绕主消节点的单元循环.
          DO ipnte = 1, kpnte
           knume = knume + 1
           kelem = lopts( knume )
           DO inode = 1, mnode
            ipoin = IABS( lnods( inode, kelem ) )
            IF( ipoin .NE. 0 ) THEN
!............当前节点为非零节点,查是否在波阵中.
             label = 0
             jpntw = 1
             ipoin = llnks( ipoin )
             DO WHILE( jpntw .LE. mpntw .AND. label .EQ. 0 )
              jpoin = lpntw( jpntw )
              IF( jpoin .EQ. ipoin ) THEN
               label = 1
              ELSE
               jpntw = jpntw + 1
              END IF
             END DO
!............如果当前节点不在波阵中,添加该节点.
             IF( label .EQ. 0 ) THEN
              jpntw = 1
              DO WHILE( label .EQ. 0 )
               IF( lpntw( jpntw ) .EQ. 0 ) THEN
                lpntw( jpntw ) = ipoin
                label = 1
               ELSE
                jpntw = jpntw + 1
               END IF
              END DO
              npntw = npntw + 1
             END IF
            END IF
           END DO
          END DO
          IF( npntw .GT. 0 ) THEN
           ndofw = npntw * mdofn
!..........组装一致矩阵对角元素.
           DO idofn = 1, mdofn
            idofw = ( kpntw - 1 ) * mdofn + idofn
            idofs = ( kpoin - 1 ) * mdofn + idofn
            istfw = idofw * ( idofw + 1 ) / 2
            IF( iswth .EQ. 1 ) THEN
              stifw( istfw ) = stifw( istfw ) + diagk( idofs )
              diagk( idofs ) = stifw( istfw )
            ELSE IF( iswth .EQ. 2 ) THEN
              stifw( istfw ) = stifw( istfw ) + diagm( idofs )
              diagm( idofs ) = stifw( istfw )
            ELSE IF( iswth .EQ. 3 ) THEN
              stifw( istfw ) = stifw( istfw ) + diagd( idofs )
              diagd( idofs ) = stifw( istfw )
            END IF
           END DO
!..........组装一致矩阵元素.
           locat = kdata
           kdofw = ( kpntw - 1 ) * mdofn
           DO idofn = 1, mdofn
            idofw = kdofw + idofn
            DO jpntw = 1, mpntw
             IF( lpntw( jpntw ) .NE. 0 ) THEN
              DO jdofn = 1, mdofn
               jdofw = ( jpntw - 1 ) * mdofn + jdofn
               IF( jpntw .LT. kpntw ) THEN
                locat = locat + 1
                istfw = idofw * ( idofw - 1 ) / 2 + jdofw
                stifw( istfw ) = stifw( istfw ) + buffs( locat )
               ELSE IF( jdofw .GT. idofw ) THEN
                locat = locat + 1
                istfw = jdofw * ( jdofw - 1 ) / 2 + idofw
                stifw( istfw ) = stifw( istfw ) + buffs( locat )
               END IF
              END DO
             END IF
            END DO
           END DO
!..........分解一致矩阵.
           DO idofn = 1, mdofn
            idofw = ( kpntw - 1 ) * mdofn + idofn
            kdofs = ( kpoin - 1 ) * mdofn + idofn
            istfw = ( idofw + 1 ) * idofw / 2
            pdiag = stifw( istfw )
            IF( DABS( pdiag ) .LT. dswth ) pdiag = dmaxv
            IF( lfixs( idofn, kpoin ) .NE. 0 ) pdiag = 1.0D0
            IF( iswth .EQ. 1 ) diagk( kdofs ) = pdiag
            IF( iswth .EQ. 2 ) diagm( kdofs ) = pdiag
            IF( iswth .EQ. 3 ) diagd( kdofs ) = pdiag
            stifw( istfw ) = 0.0D0
!...........当前行元素除以对角元素后存储.
            locat = 0
            DO jpntw = 1, mpntw
             jpoin = lpntw( jpntw )
             IF( jpoin .NE. 0 ) THEN
              DO jdofn = 1, mdofn
               jdofw = ( jpntw - 1 ) * mdofn + jdofn
               IF( jpntw .LT. kpntw ) THEN
                locat = locat + 1
                istfw = idofw * ( idofw - 1 ) / 2 + jdofw
                fbcks( locat ) = stifw( istfw ) / pdiag
                stifw( istfw ) = 0.0D0
               ELSE IF( jdofw .GT. idofw ) THEN
                locat = locat + 1
                istfw = jdofw * ( jdofw - 1 ) / 2 + idofw
                fbcks( locat ) = stifw( istfw ) / pdiag
                stifw( istfw ) = 0.0D0
               END IF
              END DO
             END IF
            END DO
            IF( iswth .EQ. 1 ) pdiag = diagk( kdofs )
            IF( iswth .EQ. 2 ) pdiag = diagm( kdofs )
            IF( iswth .EQ. 3 ) pdiag = diagd( kdofs )
            IF( DABS( pdiag ) .LT. 1.0D-16 ) THEN
!............若对角元为零,则回代信息置零.
             DO idofw = 1, ndofw - idofn
              fbcks( idofw ) = 0.0D0
             END DO
            ELSE IF( lfixs( idofn, kpoin ) .EQ. 0 ) THEN
!............波阵消元处理.
             idofw = 0
             DO ipntw = 1, mpntw
              ipoin = lpntw( ipntw )
              IF( ipoin .NE. 0 ) THEN
               DO kdofn = 1, mdofn
                jdofw = 0
                iloca = ( ipntw - 1 ) * mdofn + kdofn
                IF( ipntw .NE. kpntw .OR.                                      &
                   kdofn .GT. idofn ) THEN
                 idofw = idofw + 1
                 ifixs = lfixs( kdofn, ipoin )
                 IF( ifixs .EQ. 0 ) THEN
                  DO jpntw = 1, ipntw
                   jpoin = lpntw( jpntw )
                   IF( jpoin .NE. 0 ) THEN
                    DO jdofn = 1, mdofn
                     IF( jpntw .NE. kpntw .OR.                                 &
                        jdofn .GT. idofn ) THEN
                      jloca = ( jpntw - 1 ) * mdofn +                          &
                                jdofn
                      jdofw = jdofw + 1
                      jfixs = lfixs( jdofn, jpoin )
                      IF( jfixs .EQ. 0 ) THEN
                       IF( jloca .LE. iloca ) THEN
                        kstfw = ( iloca - 1 ) *                                &
                                  iloca / 2   + jloca
                        stifw( kstfw ) = stifw( kstfw ) -                      &
                        fbcks( idofw ) * fbcks( jdofw ) *                      &
                        pdiag
                       END IF
                      END IF
                     END IF
                    END DO
                   END IF
                  END DO
                 END IF
                END IF
               END DO
              END IF
             END DO
            END IF
            DO idofw = 1, ndofw - idofn
             locat = kdata + idofw
             buffs( locat ) = fbcks( idofw )
            END DO
            kdata = kdata + ndofw - idofn
           END DO
!..........波阵前移.
           npntw = npntw - 1
           lpntw( kpntw ) = 0
          END IF
         END DO
         IF( isavc .EQ. 1 )                                                    &
         CALL OptRec( buffs, ircrd, 2, iswth, nexch )
        END DO
        RETURN
        END

        SUBROUTINE DealFix( iswth )
        USE CtrlData
        USE MeshData
        USE GlobData
        USE SolvData
        USE FrontData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
!.......判断是否有必要分解.
        isavc = lsavc( iswth )
        IF( isavc .EQ. 0 ) RETURN
!.......
        npntw = 0
        knume = 0
        CALL Support
        DO idofs = 1, ndofs
         statf( idofs ) = 0.0D0
        END DO
!.......对角元约束处理
        DO ipoin = 1, npoin
         DO idofn = 1, mdofn
          ifixs = lfixs( idofn, ipoin )
          idofs = ( ipoin - 1 ) * mdofn + idofn
          IF( iswth .EQ. 1 ) THEN
           IF( ifixs .NE. 0 ) diagk( idofs ) = 1.0D0
           IF( nuntc .EQ. 0 ) THEN
            vunts( idofs ) = 1.0D0
           ELSE
            vuntc = DABS( diagk( idofs ) )
            IF( vuntc .LT. 1.0D-16 ) vuntc = 1.0D0
            diagk( idofs ) = diagk( idofs ) / vuntc
            vunts( idofs ) = 1.0D0 / DSQRT( vuntc )
           END IF
          ELSE IF( iswth .EQ. 2 ) THEN
           IF( ifixs .NE. 0 ) diagm( idofs ) = 1.0D0
           IF( nuntc .EQ. 0 ) THEN
            vunts( idofs ) = 1.0D0
           ELSE
            vuntc = DABS( diagm( idofs ) )
            IF( vuntc .LT. 1.0D-16 ) vuntc = 1.0D0
            diagm( idofs ) = diagm( idofs ) / vuntc
            vunts( idofs ) = 1.0D0 / DSQRT( vuntc )
           END IF
          ELSE
           IF( ifixs .NE. 0 ) diagd( idofs ) = 1.0D0
           vuntc = DABS( diagd( idofs ) )
           IF( vuntc .LT. 1.0D-16 ) vuntc = 1.0D0
           diagd( idofs ) = diagd( idofs ) / vuntc
           vunts( idofs ) = 1.0D0 / DSQRT( vuntc )
          END IF
         END DO
        END DO
!.......按记录处理非对角元约束
        IF( isavc .EQ. 2 ) RETURN
        DO ircrd = 1, nrcrd
         inump = lrcrd( 1, ircrd )
         jnump = lrcrd( 2, ircrd )
         IF( isavc .EQ. 1 ) THEN
          kdata = 0
          CALL OptRec( buffs, ircrd, 1, iswth, nexch )
         ELSE
          kdata =-isavc - 1
         END IF
         DO knump = inump, jnump
          kpoin = lbcks( 2, knump )
          kpntw = lbcks( 3, knump )
          kpnte = lpnte(    kpoin )
!.........对绕主消节点的单元循环.
          DO ipnte = 1, kpnte
           knume = knume + 1
           kelem = lopts( knume )
           DO inode = 1, mnode
            label = 0
            ipoin = IABS( lnods( inode, kelem ) )