elmt501.f90

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

        SUBROUTINE Elmt501( smaxv, ielem, imats, iswth, iuses )
!.......
!       P R O G R A M
!             To get stiffen matrix of panel elements.
!.......
        USE CtrlData
        USE ElmtData
        USE MeshData

        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION fwork( 4, 4 ), bound( 4 ), iwork( 6 )
        IF( iswth .EQ. 0 ) THEN
          IF( nprnc .NE. 0 ) WRITE( 12, 2000 )
          READ( 11, 2200 ) icode
          IF( nprnc .NE. 0 ) WRITE( 12, 2400 ) icode
          props( 2, imats ) = icode
        ELSE IF( iswth .EQ. 8 ) THEN
          RETURN
        ELSE IF( iswth .EQ. 9 ) THEN
          IF( iuses .EQ. 0 ) RETURN
          RETURN
        ELSE IF( iswth .EQ. 21 ) THEN
          WRITE( 25 ) 0
          RETURN
        END IF
!.......
!       initialize the arrays.
!.......
        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 .NE. 1 .AND. iswth .NE. 5 ) RETURN
!.....
!       get the No. of base points.
!.....
        nbasp = 0
        DO inode = 1, 3
          IF( lnode( inode ) .gt. 0 ) THEN
            nbasp = nbasp + 1
            lnode( nbasp ) = lnode( inode )
          END IF
        END DO
        DO inode = nbasp + 1, 3
          lnode( inode ) = 0
        END DO
!.......
!       code for boundary conditions.
!.......
        icode = props( 2, imats )
        IF( icode .LE. 0 ) THEN
          RETURN
        ELSE IF( icode .GT. 0 ) THEN
          DO idofn = 1, ndofn
            iwork( idofn ) = icode / 10 ** ( ndofn - idofn )
            icode = icode - iwork( idofn ) * 10 ** ( ndofn - idofn )
          END DO
        END IF
!.......
!       get boundary stiffness matrix.
!.......
        panel = 1000.D0
        DO inode = 4, nnode
          IF( lnode( inode ) .NE. 0 ) THEN
            CALL Bounc( bound, inode, nbasp, ielem )
            IF( nerrc .GT. 0 ) RETURN
            DO jnode = 1, 4
              DO knode = 1, 4
                fwork( jnode, knode ) = bound( jnode ) *                       &
                        panel * smaxv * bound( knode )
              END DO
            END DO
            DO idofn = 1, ndofn
              IF( iwork( idofn ) .GT. 0 ) THEN
                DO jnode = 1, 4
                  DO knode = 1, 4
                    jloca = jnode
                    kloca = knode
                    IF( jnode .EQ. 4 ) jloca = inode
                    IF( knode .EQ. 4 ) kloca = inode
                    jloca = ( jloca - 1 ) * ndofn + idofn
                    kloca = ( kloca - 1 ) * ndofn + idofn
                    stife( jloca, kloca ) = stife( jloca, kloca ) +            &
                                            fwork( jnode, knode )
                  END DO
                END DO
              END IF
            END DO
          END IF
        END DO
        IF( nincc .EQ. 2 ) THEN
          DO idofe = 1, ndofe
            DO jdofe = 1, ndofe
              force( idofe ) = force( idofe ) -                                &
              dispe( jdofe ) * stife( idofe, jdofe )
            END DO
          END DO
        END IF
        RETURN
2000    FORMAT( 2x, '单元类型:                罚单元' )
2200    FORMAT( I10 )
2400    FORMAT( 2x, '约束代码: ', I12  )
        END
        SUBROUTINE bounc( bound, inode, nbasp, ielem )
        USE CtrlData
        USE ElmtData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION bound( 4 ), fwork( 6 )
        DO ivare = 1, 4
         bound( ivare ) = 0.0d0
        END DO
        IF( nbasp .eq. 1 ) then
         bound( 1 ) =  1.0d0
         bound( 4 ) = -1.0d0
        ELSE IF( nbasp .eq. 2 ) then
         DO idime = 1, ndime
          fwork( idime ) =                                                     &
            ( coren( idime, inode ) - coren( idime, 1 ) ) /                    &
            ( coren( idime,     2 ) - coren( idime, 1 ) )
         END DO
         ratio = fwork( 1 )
         DO idime = 2, ndime
          IF( dabs( ratio - fwork( idime ) ) .gt. 1.d-8 ) then
           WRITE( 12, 2000 ) ielem, inode
           nerrc = 1051
           RETURN
          END IF
         END DO
         bound( 1 ) = 1.0d0 - ratio
         bound( 2 ) =         ratio
         bound( 4 ) = 1.0d0
        ELSE IF( nbasp .eq. 3 ) then
         DO idime = 1, ndime
          apara = 0.5d0 * ( coren( idime, 1 ) + coren( idime, 2 ) -            &
                                        2.0d0 * coren( idime, 3 ) )
          bpara =-0.5d0 * ( coren( idime, 2 ) - coren( 1, idime ) )
          cpara = coren( idime, 3 ) - coren( idime, inode )
          cpara = dsqrt( bpara * bpara - 4.0d0 * apara * cpara )
          root1 = 0.5d0 * ( bpara + cpara ) / apara
          root2 = 0.5d0 * ( bpara - cpara ) / apara
          IF( root1 .ge. -1.d0 .and. root1 .le. 1.0d0 )                        &
          fwork( idime ) = root1
          IF( root2 .ge. -1.d0 .and. root2 .le. 1.0d0 )                        &
          fwork( idime ) = root2
         END DO
         xloca = fwork( 1 )
         DO idime = 2, ndime
          IF( dabs( xloca - fwork( idime ) ) .gt. 1.d-8 ) then
           WRITE( 12, 2000 ) ielem, inode
           nerrc = 1051
           RETURN
          END IF
         END DO
         bound( 1 ) = -xloca * ( 1.0d0 - xloca ) / 2.0d0
         bound( 2 ) =  xloca * ( 1.0d0 + xloca ) / 2.0d0
         bound( 3 ) =  1.0d0 - xloca * xloca
         bound( 4 ) = -1.0d0
        END IF
2000    FORMAT( 2x, 'Error of element ', i5, ' node ', i3,                     &
                1x, 'not compatiable with the base nodal points ' )
        RETURN
        END