front.f90

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

            IF( ipoin .NE. 0 ) THEN
!............当前节点为非零节点,查是否在波阵中.
             ipoin = llnks( ipoin )
             DO jpntw = 1, mpntw
              IF( label .EQ. 0 ) THEN
               jpoin = lpntw( jpntw )
               IF( ipoin .EQ. jpoin ) label = 1
              END IF
             END DO
!............如果当前节点不在波阵中,添加该节点.
             IF( label .EQ. 0 ) THEN
              jpntw = 1
              npntw = npntw + 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
             END IF
            END IF
           END DO
          END DO
          ndofw = npntw * mdofn
          IF( ndofw .GT. 0 ) THEN
           DO idofn = 1, mdofn
            locat = kdata
            idofs = ( kpoin - 1 ) * mdofn + idofn
            IF( lfixs( idofn, kpoin ) .NE. 0 ) THEN
             DO ipntw = 1, mpntw
              ipoin = lpntw( ipntw )
              IF( ipoin .NE. 0 ) THEN
               DO jdofn = 1, mdofn
                jdofs = ( ipoin - 1 ) * mdofn + jdofn
                IF( ipntw .NE. kpntw .OR. jdofn .GT. idofn ) THEN
                 locat = locat + 1
                 IF( lfixs( jdofn, ipoin ) .EQ. 0 )                            &
                     statf( jdofs ) = statf( jdofs ) -                         &
                     buffs( locat ) * sprts( idofn, kpoin )
                 IF( nuntc .NE. 0 )                                            &
                     buffs( locat ) = buffs( locat ) *                         &
                     vunts( idofs ) * vunts( jdofs )
                END IF
               END DO
              END IF
             END DO
            ELSE
             DO ipntw = 1, mpntw
              ipoin = lpntw( ipntw )
              IF( ipoin .NE. 0 ) THEN
               DO jdofn = 1, mdofn
                jdofs = ( ipoin - 1 ) * mdofn + jdofn
                IF( ipntw .NE. kpntw .OR. jdofn .GT. idofn ) THEN
                 locat = locat + 1
                 IF( lfixs( jdofn, ipoin ) .NE. 0 )                            &
                     statf( idofs ) = statf( idofs ) -                         &
                     buffs( locat ) * sprts( jdofn, ipoin )
                 IF( nuntc .NE. 0 )                                            &
                     buffs( locat ) = buffs( locat ) *                         &
                     vunts( idofs ) * vunts( jdofs )
                END IF
               END DO
              END IF
             END DO
            END IF
            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 BckSub( iswth )
!........
!        模块功能
!            右端消元处理后回代求解.
!........
        USE CtrlData
        USE MeshData
        USE SolvData
        USE GlobData
        USE FrontData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        INTEGER * 2 year, month, date
!.......设定常数.
        isavc = lsavc( iswth )
        IF( isavc .EQ. 0 ) RETURN
!.......计算残差误差
        ftolc = 0.0D0
        vsumm = 0.0D0
        DO ipoin = 1, npoin
         DO idofn = 1, mdofn
          IF( lfixs( idofn, ipoin ) .EQ. 0 ) THEN
           idofs = ( ipoin - 1 ) * mdofn + idofn
           vtemp = vslvs( idofs ) + fstrs( idofs )
           ftolc = ftolc + vtemp * vtemp
           vsumm = vsumm + vslvs( idofs ) * vslvs( idofs )
          END IF
         END DO
        END DO
        vsumm = DSQRT( vsumm )
        ftolc = DSQRT( ftolc )
        IF( vsumm .GT. 1.0D-16 ) ftolc = ftolc / vsumm
!.......动载迭加静载(体积力和约束作用).
        DO idofs = 1, ndofs
         vslvs( idofs ) = vslvs( idofs ) +                                     &
         statf( idofs ) + fstrs( idofs )
        END DO
!.......处理给定的位移信息.
        DO ipoin = 1, npoin
         DO idofn = 1, mdofn
          IF( lfixs( idofn, ipoin ) .NE. 0 ) THEN
           idofs = ( ipoin - 1 ) * mdofn + idofn
           vslvs( idofs ) = sprts( idofn, ipoin )
          END IF
         END DO
        END DO
!.......右端项处理
        IF( nuntc .NE. 0 ) THEN
         DO idofs = 1, ndofs
          vslvs( idofs ) = vslvs( idofs ) * vunts( idofs )
         END DO
        END IF
!.......右端消元.
        npntw = 0
        knume = 0
        kcoun = 0
        DO ipntw = 1, mpntw
         lpntw( ipntw ) = 0
        END DO
        DO ircrd = 1, nrcrd
         inump = lrcrd( 1, ircrd )
         jnump = lrcrd( 2, ircrd )
         IF( isavc .EQ. 1 ) kdata = 0
         IF( isavc .LT. 0 ) kdata =-isavc - 1
         IF( isavc .EQ. 1 )                                                    &
         CALL OptRec( buffs, ircrd, 1, iswth, nexch )
         DO knump = inump, jnump
          kpoin = lbcks( 2, knump )
          kpntw = lbcks( 3, knump )
          npnte = lpnte(    kpoin )
          kcoun = kcoun + 1
          IF( necho .NE. 0 ) CALL ShowProcess( kcoun, npoin )


          DO ipnte = 1, npnte
           knume = knume + 1
           kelem = lopts( knume )
           DO inode = 1, mnode
            label = 0
            ipoin = IABS( lnods( inode, kelem ) )
            IF( ipoin .NE. 0 ) THEN
             ipoin = llnks( ipoin )
             DO ipntw = 1, mpntw
              IF( ipoin .EQ. lpntw( ipntw ) ) label = 1
             END DO
             IF( label .EQ. 0 ) THEN
              ipntw = 1
              npntw = npntw + 1
              DO WHILE( label .EQ. 0 )
               IF( lpntw( ipntw ) .EQ. 0 ) THEN
                lpntw( ipntw ) = ipoin
                label = 1
               ELSE
                ipntw = ipntw + 1
               END IF
              END DO
             END IF
            END IF
           END DO
          END DO
          ndofw = npntw * mdofn
          IF( npntw .GT. 0 ) THEN
!..........如果波阵中的节点数不为零，回代
           DO idofn = 1, mdofn
            idofs = ( kpoin - 1 ) * mdofn + idofn
            IF( iswth .EQ. 1 ) pdiag = diagk( idofs )
            IF( iswth .EQ. 2 ) pdiag = diagm( idofs )
            IF( iswth .EQ. 3 ) pdiag = diagd( idofs )
            IF( DABS( pdiag ) .LT. 1.0D-16 ) THEN
!............若对角元为零,则右端项置零.
             IF( nincc .LE. 1 ) THEN
              vslvs( idofs ) = disps( idofs, kload )
             ELSE
              vslvs( idofs ) = 0.0D0
             END IF
            ELSE IF( lfixs( idofn, kpoin ) .EQ. 0 ) THEN
!............右端消元运算.
             kdofw = 0
             DO jpntw = 1, mpntw
              IF( lpntw( jpntw ) .NE. 0 ) THEN
               DO jdofn = 1, mdofn
                IF( jpntw .NE. kpntw .OR. jdofn .GT. idofn ) THEN
                 kdofw = kdofw + 1
                 jpoin = lpntw( jpntw )
                 jdofs = ( jpoin - 1 ) * mdofn + jdofn
                 IF( lfixs( jdofn, jpoin ) .EQ. 0 )                            &
                 vslvs( jdofs ) = vslvs( jdofs ) -                             &
                 vslvs( idofs ) * buffs( kdata + kdofw )
                END IF
               END DO
              END IF
             END DO
            END IF
            kdata = kdata + ndofw - idofn
           END DO
           npntw = npntw - 1
           lpntw( kpntw ) = 0
          ELSE
!..........如果波阵中的节点数为零，右端项置零.
           DO idofn = 1, mdofn
            idofs = ( kpoin - 1 ) * mdofn + idofn
            IF( nincc .LE. 1 ) THEN
              vslvs( idofs ) = disps( idofs, kload )
            ELSE
              vslvs( idofs ) = 0.0D0
            END IF
           END DO
          END IF
         END DO
        END DO
!.......右端项除以对角元.
        DO ipoin = 1, npoin
         DO idofn = 1, mdofn
          IF( lfixs( idofn, ipoin ) .EQ. 0 ) THEN
           idofs = ( ipoin - 1 ) * mdofn + idofn
           IF( iswth .EQ. 1 ) pdiag = diagk( idofs )
           IF( iswth .EQ. 2 ) pdiag = diagm( idofs )
           IF( iswth .EQ. 3 ) pdiag = diagd( idofs )
           IF( DABS( pdiag ) .GT. 1.0D-16 )                                    &
           vslvs( idofs ) = vslvs( idofs ) / pdiag
          END IF
         END DO
        END DO
!........
!       递推法形成当步波阵节点序号与节点总体编号对照表,并按回代信息
!       求解位移向量.
!........
        DO ircrd = nrcrd, 1, -1
         inump = lrcrd( 1, ircrd )
         jnump = lrcrd( 2, ircrd )
         kdata = lrcrd( 3, ircrd )
         IF( isavc .LT. 0 ) kdata = kdata - isavc - 1
         IF( isavc .EQ. 1 )                                                    &
             CALL OptRec( buffs, ircrd, 1, iswth, nexch )
         DO knump = jnump, inump, -1
!.........取当步波阵信息
          npntw = lbcks( 1, knump )
          ipoin = lbcks( 2, knump )
          kpntw = lbcks( 3, knump )
          IF( npntw .GT. 0 ) THEN
!..........恢复波前中的总体节点号.
           ndofw = npntw * mdofn
           lpntw( kpntw ) = ipoin
           DO ipntw = 1, mpntw
             jpoin = lpntw( ipntw )
             IF( jpoin .GT. 0 ) THEN
               IF( lbcks( 4, jpoin ) .GT. knump ) lpntw( ipntw ) = 0
             END IF
           END DO
!..........回代从此开始.
           DO idofn = mdofn, 1, -1
!...........从缓冲区中取回代信息,回代运算.
            kdofw = 0
            kdata = kdata - ndofw + idofn
            idofs = ( ipoin - 1 ) * mdofn + idofn
            DO ipntw = 1, mpntw
             jpoin = lpntw( ipntw )
             IF( jpoin .NE. 0 ) THEN
              idofw = ( jpoin - 1 ) * mdofn
              DO jdofn = 1, mdofn
               idofw = idofw + 1
               IF( ipntw .NE. kpntw .OR. jdofn .GT. idofn ) THEN
                kdofw = kdofw + 1
                IF( lfixs( idofn, ipoin ) .EQ. 0 .AND.                         &
                    lfixs( jdofn, jpoin ) .EQ. 0 )                             &
                vslvs( idofs ) = vslvs( idofs ) -                              &
                vslvs( idofw ) * buffs( kdata + kdofw )
               END IF
              END DO
             END IF
            END DO
           END DO
          END IF
         END DO
        END DO

        IF( nuntc .NE. 0 ) THEN
         DO idofs = 1, ndofs
          vslvs( idofs ) = vslvs( idofs ) * vunts( idofs )
         END DO
        END IF

        RETURN
        END

        SUBROUTINE DotMat( vmatx, vresv, iswth, kmatx, nmatx )
        USE CtrlData
        USE MeshData
        USE GlobData
        USE SolvData
        USE FrontData
        IMPLICIT DOUBLE PRECISION( A-H, O-Z )
        DIMENSION vmatx( ndofs, nmatx ), vresv( ndofs )
        isavc = lsavc( iswth )
        IF( isavc .EQ. 0 ) THEN
         nerrc = 5001
         RETURN
        END IF
        IF( iswth .EQ. 1 ) THEN
          DO idofs = 1, ndofs
            vresv( idofs ) = vresv( idofs ) +                                  &
            diagk( idofs ) * vmatx( idofs, kmatx )
          END DO
        ELSE IF( iswth .EQ. 2 ) THEN
          DO idofs = 1, ndofs
            vresv( idofs ) = vresv( idofs ) +                                  &
            diagm( idofs ) * vmatx( idofs, kmatx )
          END DO
        ELSE IF( iswth .EQ. 3 ) THEN
          DO idofs = 1, ndofs
            vresv( idofs ) = vresv( idofs ) +                                  &
            diagd( idofs ) * vmatx( idofs, kmatx )
          END DO
        END IF
        IF( isavc .EQ. 2 ) RETURN
        npntw = 0
        knume = 0
        DO ipntw = 1, mpntw
         lpntw( ipntw ) = 0
        END DO
        DO ircrd = 1, nrcrd
         IF( isavc .EQ. 1 ) THEN
          kdata = 0
          CALL OptRec( buffs, ircrd, 1, iswth, nexch )
         ELSE
          kdata =-isavc - 1
         END IF
         inump = lrcrd( 1, ircrd )
         jnump = lrcrd( 2, ircrd )
         DO knump = inump, jnump
          kpoin = lbcks( 2, knump )
          kpntw = lbcks( 3, knump )
          kpnte = lpnte(    kpoin )
          DO ipnte = 1, kpnte
           knume = knume + 1
           ielem = lopts( knume )
           DO inode = 1, mnode
            ipoin = IABS( lnods( inode, ielem ) )
            IF( ipoin .NE. 0 ) THEN
             label = 0
             ipoin = llnks( ipoin )
             DO ipntw = 1, mpntw
               IF( ipoin .EQ. lpntw( ipntw ) ) label = 1
             END DO
             IF( label .EQ. 0 ) THEN
              ipntw = 1
              npntw = npntw + 1
              DO WHILE( label .EQ. 0 )
               IF( lpntw( ipntw ) .EQ. 0 ) THEN
                lpntw( ipntw ) = ipoin
                label = 1
               ELSE
                ipntw = ipntw + 1
               END IF
              END DO
             END IF
            END IF
           END DO
          END DO
          IF( npntw .GT. 0 ) THEN
           DO idofn = 1, mdofn
            idofs = ( kpoin - 1 ) * mdofn + idofn
            DO ipntw = 1, mpntw
             ipoin = lpntw( ipntw )
             IF( ipoin .GT. 0 ) THEN
              DO jdofn = 1, mdofn
               IF( ipntw .NE. kpntw .OR. jdofn .GT. idofn ) THEN
                kdata = kdata + 1
                jdofs = ( ipoin - 1 ) * mdofn + jdofn
                vresv( idofs ) = vresv( idofs ) +                              &
                buffs( kdata ) * vmatx( jdofs, kmatx )
                vresv( jdofs ) = vresv( jdofs ) +                              &
                buffs( kdata ) * vmatx( idofs, kmatx )
               END IF
              END DO
             END IF
            END DO
           END DO
           npntw = npntw - 1
           lpntw( kpntw ) = 0
          END IF
         END DO
        END DO
        RETURN
        END