elmt401.f90

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

        SUBROUTINE Elmt401( ielem, imats, iswth, iuses )
!.......
!       模块功能
!           任意位置弹簧桩4-8节点单元.
!.......
        USE CtrlData
        USE MeshData
        USE ElmtData
        USE GlobData
        USE FactorData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION vnorm( 3 ), vdisp( 3 )
        DIMENSION vcoor( 3 ), vforc( 3 ), shape( 3, 8 )

        IF( iswth .EQ. 0 ) THEN
          CALL RdMat401( imats )
          RETURN
        ELSE IF( iswth .EQ. 8 ) THEN
          nnode = 8
          ndofn = mdofn
          RETURN
        ELSE IF( iswth .EQ. 9 ) THEN
          IF( iuses .EQ. 0 ) RETURN
!          knode = MIN0( nnode, 8 )
!          WRITE( 15, 2200 ) 4, ielem, imats, knode,
!     +          ( lnode( inode ), inode = 1, knode )
          RETURN
        ELSE IF( iswth .EQ. 21 ) THEN
          WRITE( 25 ) 0
          RETURN
        END IF

        epara = props( 2, imats )
        hpara = props( 3, imats )
        xloca = props( 4, imats )
        yloca = props( 5, imats )
        DO idime = 1, ndime
          vnorm( idime ) = props( idime + 5, imats )
        END DO
        DO idime = 1, ndime
          vdisp( idime ) = props( ndime + idime + 5, imats )
        END DO
        DO idofe = 1, ndofe
          DO jdofe = 1, ndofe
            stife( idofe, jdofe ) = 0.0D0
          END DO
          force( idofe ) = 0.0D0
        END DO
        IF( iuses .EQ. 0 ) RETURN

        IF( iswth .EQ. 1 .OR. iswth .EQ. 5 .OR. iswth .EQ. 10 ) THEN
!.........计算单元刚度矩阵.
          CALL Shap2D( shape, coren, xloca, yloca, lnode, xjaco,               &
                       ndime,     2, nnode, ielem,     1, nerrc )
          IF( nerrc .GT. 0 ) RETURN
          DO inode = 1, nnode
            DO jnode = 1, nnode
              DO idime = 1, ndime
                DO jdime = 1, ndime
                  const = vnorm( idime ) * vnorm( jdime )
                  iloca = ( inode -     1 ) * ndofn + idime
                  jloca = ( jnode -     1 ) * ndofn + jdime
                  const = ( epara - hpara ) * const
                  IF( idime .EQ. jdime ) const = const + hpara
                  stife( iloca, jloca ) = shape( 3, inode ) *                  &
                                          shape( 3, jnode ) * const
                END DO
              END DO
            END DO
          END DO
!.........计算内力等效节点力
          IF( nincc .EQ. 2 .OR. iswth .EQ. 10 ) THEN
            DO idofe = 1, ndofe
              DO jdofe = 1, ndofe
                force( idofe        ) = force( idofe ) -                       &
                stife( idofe, jdofe ) * dispe( jdofe )
              END DO
            END DO
          END IF
        ELSE IF( iswth .EQ. 2 ) THEN
        ELSE IF( iswth .EQ. 3 ) THEN
        ELSE IF( iswth .EQ. 4 ) THEN
        ELSE IF( iswth .EQ. 5 ) THEN
        ELSE IF( iswth .EQ. 6 ) THEN
          CALL Shap2D( shape, coren, xloca, yloca, lnode, xjaco,               &
                       ndime,     2, nnode, ielem,     1, nerrc )
          DO inode = 1, nnode
            DO idime = 1, ndime
              idofe = ( inode - 1 ) * ndofn + idime
              DO jdime = 1, ndime
                force(    idofe ) = force( idofe ) +                           &
                shape( 3, inode ) * vnorm( idime ) *                           &
                vnorm(    jdime ) * vdisp( jdime ) * (                         &
                epara  -  hpara )
              END DO
              force(    idofe ) = force( idofe ) +                             &
              shape( 3, inode ) * vdisp( idime ) * hpara
            END DO
          END DO
        ELSE IF( MOD( iswth, mswth ) .EQ. 7 ) THEN
          CALL Shap2D( shape, coren, xloca, yloca, lnode, xjaco,               &
                       ndime,     2, nnode, ielem,     1, nerrc )
          DO idime = 1, ndime
            vdisp( idime ) =-vdisp( idime )
            DO inode = 1, nnode
              idofe = ( inode - 1 ) * ndofn + idime
              vdisp( idime ) = vdisp(    idime ) +                             &
              dispe( idofe ) * shape( 3, inode )
            END DO
          END DO

          DO idime = 1, ndime
            vcoor( idime ) = 0.0D0
            DO inode = 1, nnode
              vcoor(    idime ) = vcoor( idime ) +                             &
              shape( 3, inode ) * coren( idime, inode )
            END DO
          END DO

          DO idime = 1, ndime
            vforc( idime ) = 0.0D0
            DO jdime = 1, ndime
              vforc( idime ) = vforc( idime ) +                                &
              vnorm( idime ) * vnorm( jdime ) *                                &
              vdisp( jdime ) * ( epara - hpara )
            END DO
            vforc( idime ) = vforc( idime ) +                                  &
            vdisp( idime ) * hpara
          END DO
          WRITE( 12, 2400 ) ielem, vcoor
          WRITE( 12, 2600 ) vforc
        END IF
2200    FORMAT( 20I5 )
2400    FORMAT( I5, 5x, 3F15.5 )
2600    FORMAT( 10X, 3E15.5 )
        RETURN
        END

        SUBROUTINE RdMat401( imats )
        USE CtrlData
        USE MeshData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )

        mprop = ndime * 2 + 5
        IF( nprop .LT. mprop  ) THEN
          WRITE( 12, 2000 ) mprop
          nerrc = 56450
          RETURN
        END IF

        READ( 11, 2200 ) ( props( iprop, imats ), iprop = 2, mprop )
        IF( nprnc .NE. 0 ) THEN
          WRITE( 12, 2400 )
          WRITE( 12, 2600 ) ( props( iprop, imats ), iprop = 2, 5 )
          WRITE( 12, 2800 ) ( props( iprop, imats ), iprop = 6,                &
                              ndime + 5 )
          WRITE( 12, 3000 ) ( props( iprop, imats ), iprop =                   &
                              ndime + 6, 2 * ndime + 5 )
        END IF

2000    FORMAT( 2x, '致命错误：材料参数数目至少为', I5 )
2200    FORMAT( 4E15.5 / 4E15.5 / 4E15.5 )
2400    FORMAT( 2x, '单元类型：      弹簧桩单元' )
2600    FORMAT( 2x, '轴向模量:',  E15.5 / 2x, '横向模量:',  E15.5 /            &
                2x, '支点坐标:', 2F15.5 )
2800    FORMAT( 2x, '轴线方向:', 3F15.5 )
3000    FORMAT( 2x, '初始位移:', 3E15.5 )
        RETURN
        END