load.f90

非线性有限元分析程序

        SUBROUTINE Load( forcs )
!.......
!       模块功能
!           计算集中力, 线分布载荷, 面分布载荷的等效节点力, 并组装到
!       整体力矢量中.
!.......
        USE CtrlData
        USE LoadData
        USE MeshData
        USE FactorData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION forcs( ndofs )
        CALL Load01( forcs )
        CALL Load02( forcs )
        CALL Load03( forcs )
        RETURN
        END

        SUBROUTINE Load01( forcs )
!.....
!       模块功能
!           组装集中力到整体力矢量中.
!.....
        USE CtrlData
        USE LoadData
        USE MeshData
        USE FactorData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION forcs( ndofs )

        ktype = -99999
        nnodf = loads( 1, kload )
        IF( nnodf .LE. 0 ) RETURN
        IF( nerrc .GT. 0 ) RETURN

        knodf = 0
        DO iload = 1, kload - 1
          knodf = knodf + loads( 1, iload )
        END DO
        DO inodf = knodf + 1, knodf + nnodf
          itype =        lnodf( 2, inodf )
          ipoin = llnks( lnodf( 1, inodf ) )
          IF( itype .NE. ktype ) THEN
            ktype = itype
            CALL Factor( ktype )
            IF( nerrc .GT. 0 ) RETURN
          END IF
          DO idofn = 1, mdofn
            idofs = ( ipoin - 1 ) * mdofn + idofn
            forcs( idofs ) = forcs( idofs ) +                                  &
            vnodf( idofn, inodf ) * fctor
          END DO
        END DO
        RETURN
        END

        SUBROUTINE Load02( forcs )
!.......
!       模块功能
!           计算线分布载荷的等效节点力, 并组装到整体力矢量中.
!.......
        USE CtrlData
        USE LoadData
        USE MeshData
        USE FactorData
        IMPLICIT  DOUBLE PRECISION( a-h, o-z )
        DIMENSION distg( 6, 3 ), lnodg( 3 )
        DIMENSION shape( 2, 3 ), wgaus( 6 )
        DIMENSION forcg( 6, 9 ), pgaus( 6 ), forcs( ndofs )
        ALLOCATABLE coreg(:,:)

        nedgf = loads( 2, kload )
        IF( nedgf .LE. 0 ) RETURN
        IF( nerrc .GT. 0 ) RETURN
        ALLOCATE( coreg( ndime, mnodg ), STAT = ierro )
        IF( ierro .GT. 0 ) nerrc = 13245
        IF( nerrc .GT. 0 ) RETURN

        iflag = 2
        ngaus = 2
        ktype =-99999
        nnodg = mnodg - 1
        IF( nnodg .GT. 2 ) ngaus = 3
        CALL Gauss( pgaus, wgaus, ngaus, nerrc )
        IF( nerrc .GT. 0 ) RETURN

        kedgf = 0
        DO iload = 1, kload - 1
          kedgf = kedgf + loads( 2, iload )
        END DO
        DO iedgf = kedgf + 1, kedgf + nedgf
!.........设定当前线分布载荷信息
          DO inodg = 1, nnodg
            ipoin = ledgf( inodg, iedgf )
            lnodg( inodg ) = ipoin
            IF( ipoin .GT. 0 ) THEN
              DO idime = 1, ndime
                coreg( idime, inodg ) = coord( idime, ipoin )
              END DO
              DO idofn = 1, mdofn
                locat = ( inodg - 1 ) * mdofn + idofn
                distg( idofn, inodg ) = vedgf( locat, iedgf )
              END DO
            END IF
          END DO
          itype = ledgf( mnodg, iedgf )
          IF( itype .NE. ktype ) THEN
            ktype = itype
            CALL Factor( ktype )
            IF( nerrc .GT. 0 ) RETURN
          END IF
          DO idofn = 1, mdofn
            DO inodg = 1, nnodg
              distg( idofn, inodg ) = distg( idofn, inodg ) * fctor
            END DO
          END DO
!.........初时化
          DO idofn = 1, mdofn
            DO inodg = 1, nnodg
              forcg( idofn, inodg ) = 0.0d0
            END DO
          END DO
!.........计算线分布载荷的等效节点力
          DO igaus = 1, ngaus
            xloca = pgaus( igaus )
            CALL shap1d( shape, coreg, xloca, ndime, nnodg, lnodg,             &
                         xjaco, iflag )
            dvolu = wgaus( igaus ) * xjaco
            DO idofn = 1, mdofn
              fgaus = 0.0d0
              DO inodg = 1, nnodg
                fgaus = fgaus + shape(2, inodg) * distg(idofn, inodg)
              END DO
              DO inodg = 1, nnodg
                forcg( idofn, inodg ) = forcg( idofn, inodg ) +                &
                        dvolu * fgaus * shape(     2, inodg )
              END DO
            END DO
          END DO
!.........组装到整体力矢量中.
          DO inodg = 1, nnodg
            ipoin = lnodg( inodg )
            IF( ipoin .GT. 0 ) THEN
              ipoin = llnks( ipoin )
              DO idofn = 1, mdofn
                idofs = ( ipoin - 1 ) * mdofn + idofn
                forcs( idofs ) = forcs( idofs ) + forcg( idofn, inodg )
              END DO
            END IF
          END DO
        END DO
        DEALLOCATE( coreg )
        RETURN
        END

        SUBROUTINE Load03( forcs )
!.......
!       模块功能
!           计算面分布载荷的等效节点力, 并组装到整体力矢量中.
!.......
        USE CtrlData
        USE LoadData
        USE MeshData
        USE FactorData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION shape( 3, 8 ), wgaus( 6 )
        DIMENSION distp( 6, 9 ), pgaus( 6 )
        DIMENSION forcp( 6, 9 ), lnodp( 9 ), forcs( ndofs )
        ALLOCATABLE corep(:,:)

        nplnf = loads( 3, kload )
        IF( nplnf .LE. 0 ) RETURN
        IF( nerrc .GT. 0 ) RETURN
        ALLOCATE( corep( ndime, mnodp ), STAT = ierro )
        IF( ierro .GT. 0 ) nerrc = 43554
        IF( nerrc .GT. 0 ) RETURN

        iflag = 2
        ngaus = 2
        ktype =-99999
        nnodp = mnodp - 1
        IF( nnodp .GT. 4 ) ngaus = 3
        CALL gauss( pgaus, wgaus, ngaus, nerrc )
        IF( nerrc .GT. 0 ) RETURN

        kplnf = 0
        DO iload = 1, kload - 1
          kplnf = kplnf + loads( 3, iload )
        END DO
        DO iplnf = kplnf + 1, kplnf + nplnf
!.........设定当前面分布载荷信息
          DO inodp = 1, nnodp
            ipoin = lplnf( inodp, iplnf )
            lnodp( inodp ) = ipoin
            IF( ipoin .GT. 0 ) THEN
              DO idime = 1, ndime
                corep( idime, inodp ) = coord( idime, ipoin )
              END DO
              DO idofn = 1, mdofn
                locat = ( inodp - 1 ) * mdofn + idofn
                distp( idofn, inodp ) = vplnf( locat, iplnf )
              END DO
            END IF
          END DO
          itype = lplnf( mnodp, iplnf )
          IF( itype .NE. ktype ) THEN
            ktype = itype
            CALL Factor( ktype )
            IF( nerrc .GT. 0 ) RETURN
          END IF
          DO idofn = 1, mdofn
            DO inodp = 1, nnodp
              distp( idofn, inodp ) = distp( idofn, inodp ) * fctor
            END DO
          END DO
!.........初时化
          DO idofn = 1, mdofn
            DO inodp = 1, nnodp
              forcp( idofn, inodp ) = 0.0d0
            END DO
          END DO
!.........计算面分布载荷的等效节点力
          DO ixgas = 1, ngaus
            DO iygas = 1, ngaus
              xloca = pgaus( ixgas )
              yloca = pgaus( iygas )
              call shap2d( shape, corep, xloca, yloca, lnodp, xjaco,           &
                           ndime, ndime, nnodp, ielem, iflag, nerrc )
              IF( nerrc .GT. 0 ) RETURN
              dvolu = wgaus( ixgas ) * wgaus( iygas ) * xjaco
              DO idofn = 1, mdofn
                fgaus = 0.0d0
                DO inodp = 1, nnodp
                  fgaus = fgaus + shape(3, inodp) * distp(idofn, inodp)
                END DO
                DO inodp = 1, nnodp
                  forcp( idofn, inodp ) = forcp( idofn, inodp ) +              &
                  shape(     3, inodp ) * fgaus * dvolu
                END DO
              END DO
            END DO
          END DO
!.........组装到整体力矢量中.
          DO inodp = 1, nnodp
            ipoin = lnodp( inodp )
            IF( ipoin .GT. 0 ) THEN
              ipoin = llnks( ipoin )
              DO idofn = 1, mdofn
                idofs = ( ipoin - 1 ) * mdofn + idofn
                forcs( idofs ) = forcs( idofs ) + forcp( idofn, inodp )
              END DO
            END IF
          END DO
        END DO
        DEALLOCATE( corep )
        RETURN
        END