reslv.f90

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

        SUBROUTINE ReSolveStif
        USE CtrlData
        USE ElmtData
        USE MeshData
        USE GlobData
        USE ReSolveData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION pxgas( 6 ), wxgas( 6 )
        DIMENSION pygas( 6 ), wygas( 6 )
        DIMENSION shape( 3, 8 )

        IF( ndofr .EQ. 0 ) RETURN

        iflag = 2
        mdime = 2
        nxgas = 3
        nygas = 3
        iswth = 1

        CALL Gauss( pxgas, wxgas, nxgas, nerrc )
        CALL Gauss( pygas, wygas, nygas, nerrc )

        DO iequr = 1, nequr
          forcr( iequr ) = 0.0D0
        END DO

        DO istfr = 1, nstfr
          stifr( istfr ) = 0.0D0
        END DO

        DO ielem = 1, nelem
          nfilt = 0
          imats = lmats( ielem )
          xfilt = props( 10, imats )
          yfilt = props( 11, imats )
          CALL ElmInf( ielem, iswth )
          IF( xfilt + yfilt .LT. 1.0D-16 ) nfilt = 1
          DO ixgas = 1, nxgas
            DO iygas = 1, nygas
              xloca = pxgas( ixgas )
              yloca = pygas( iygas )
              CALL Shap2D( shape, coren, xloca, yloca, lnode, xjaco,           &
                           ndime, mdime, nnode, ielem, iflag, nerrc )
              IF( nerrc .GT. 0 ) RETURN
              dvolu = wxgas( ixgas ) * wygas( iygas ) * xjaco
              xcons = dvolu * props( 10, imats ) / 9.8D0
              ycons = dvolu * props( 11, imats ) / 9.8D0
              DO inode = 1, nnode
                DO jnode = 1, nnode
                  stife( inode, jnode ) = stife( inode, jnode ) +              &
                  shape( 1, inode ) * shape( 1, inode ) * xcons +              &
                  shape( 2, inode ) * shape( 2, inode ) * ycons
                END DO
              END DO
              IF( nfilt .EQ. 0 ) THEN
                const = 0.0D0
                DO inode = 1, nnode
                  idofu = ( lnode( inode ) - 1 ) * mdofn + 1
                  idofv = ( lnode( inode ) - 1 ) * mdofn + 2
                  const = const +                                              &
                  shape( 1, inode ) * veloc( idofu, kload ) +                  &
                  shape( 2, inode ) * veloc( idofv, kload )
                END DO
                DO inode = 1, nnode
                  force( inode ) = force(    inode ) + dvolu *                 &
                                   shape( 3, inode ) * const
                END DO
              END IF
            END DO
          END DO
          CALL ReAssemble
        END DO
        RETURN
        END

        SUBROUTINE ReAssemble
        USE CtrlData
        USE ElmtData
        USE ReSolveData
        DO inode = 1, mnode
          ipoin = lnode( inode )
          DO idofr = 1, ndofr
            iequr = lequr( idofr, ipoin )
            idofe = ( inode - 1 ) * ndofr + idofr
            IF( iequr .GT. 0 ) THEN
              DO jnode = 1, mnode
                jpoin = lnode( jnode )
                DO jdofr = 1, ndofr
                  jequr = lequr( jdofr, jpoin )
                  jdofe = ( jnode - 1 ) * ndofr + jdofr
                  IF( jequr .GT. 0 .AND. jequr .LE. iequr ) THEN
                    locat = lposr( iequr ) - iequr + jequr
                    stifr( locat ) = stifr( locat        ) +                   &
                                     stife( idofe, jdofe )
                  END IF
                END DO
              END DO
              forcr( iequr ) = forcr( iequr ) + force( idofe )
            END IF
          END DO
        END DO
        RETURN
        END

        SUBROUTINE ReDecomp
        USE ReSolveData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )

        DO iequr = 2, nequr
          jeque = iequr - 1
          jequb = lposr( jeque ) - lposr( iequr ) + iequr + 1
          IF( jequb .LE. jeque ) THEN
            DO jequr = jequb, jeque
              IF( jequr .NE. 1 ) THEN
                keque = jequr - 1
                kposr = lposr( iequr ) - iequr + jequr
                kequb = lposr( keque ) - lposr( jequr ) + jequr + 1
                IF( kequb .LT. jequb ) kequb = jequb
                IF( keque .GE. kequb ) THEN
                  DO kequr = kequb, keque
                    iposr = lposr( iequr ) - iequr + kequr
                    jposr = lposr( jequr ) - jequr + kequr
                    stifr( kposr ) = stifr( kposr ) -                          &
                    stifr( iposr ) * stifr( jposr )
                  END DO
                END IF
              END IF
            END DO

            DO jequr = jequb, jeque
              iposr = lposr( iequr )
              jposr = lposr( jequr )
              kposr = lposr( iequr ) - iequr + jequr
              stifr( kposr ) = stifr( kposr ) / stifr( jposr )
              stifr( iposr ) = stifr( iposr ) - stifr( kposr ) *               &
              stifr( kposr ) * stifr( jposr )
            END DO
          END IF
        END DO
        RETURN
        END

        SUBROUTINE ReBakSub
        USE ReSolveData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )

        DO iequr = 2, nequr
          jeque = iequr - 1
          jequb = lposr( jeque ) - lposr( iequr ) + iequr + 1
          IF( jequb .LE. jeque ) THEN
            DO jequr = jequb, jeque
              iposr = lposr( iequr ) - iequr + jequr
              forcr( iequr ) = forcr( iequr ) -                                &
              forcr( jequr ) * stifr( iposr )
            END DO
          END IF
        END DO

        DO iequr = 1, nequr
          iposr = lposr( iequr )
          IF( DABS( stifr( iposr ) ) .LT. 1.0D-16 ) THEN
            forcr( iequr ) = 0.0D0
          ELSE
            forcr( iequr ) = forcr( iequr ) / stifr( iposr )
          END IF
        END DO

        DO iequr = 2, nequr
          iposr = lposr( iequr )
          IF( DABS( stifr( iposr ) ) .GT. 1.0D-16 ) THEN
            kequr = nequr - iequr + 2
            jeque = nequr - iequr + 1
            jequb = lposr( jeque ) - lposr( kequr ) + kequr + 1
            IF( jequb .LE. jeque ) THEN
              DO jequr = jequb, jeque
                kposr = lposr( kequr ) - kequr + jequr
                forcr( jequr ) = forcr( jequr ) -                              &
                forcr( kequr ) * stifr( kposr )
              END DO
            END IF
          END IF
        END DO
        RETURN
        END

        SUBROUTINE ReExchang
        USE CtrlData
        USE GlobData
        USE ResolveData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DO ipoin = 1, npoin
          DO idofr = 1, ndofr
            idofn = ldofr( idofr )
            iequr = lequr( idofr, ipoin )
            idofs = ( ipoin - 1 ) * mdofn + idofn
            IF( iequr .GT. 0 ) THEN
              disps( idofs, kload ) = forcr( iequr )
            ELSE
              disps( idofs, kload ) = 0.0D0
            END IF
          END DO
        END DO
        RETURN
        END