shape.f90

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

        SUBROUTINE Shap1D( shape, coren, xloca, ndime, nnode, lnode,           &
                           xjaco, iflag )
!.......
!       P R O G R A M
!             To calculate shape function and its derivations of one
!             local dimension.
!             if iflag == 1 derivations to local coordinate system.
!             if iflag == 2 derivations to global coordinate system.
!........
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION coren( ndime, nnode ), lnode( nnode ), shape( 2, 3 )

        xjaco = 0.0D0
        shape( 1, 1 ) = -0.5D0
        shape( 1, 2 ) =  0.5D0
        shape( 2, 1 ) = ( 1.0D0 - xloca ) / 2.0D0
        shape( 2, 2 ) = ( 1.0D0 + xloca ) / 2.0D0
        IF( nnode .GE. 3 .AND. lnode( 3 ) .NE. 0 ) THEN
          shape( 1, 3 ) = -2.0D0 * xloca
          shape( 2, 3 ) =  1.0D0 - xloca * xloca
          shape( 1, 1 ) = shape( 1, 1 ) - 0.5D0 * shape( 1, 3 )
          shape( 1, 2 ) = shape( 1, 2 ) - 0.5D0 * shape( 1, 3 )
          shape( 2, 1 ) = shape( 2, 1 ) - 0.5D0 * shape( 2, 3 )
          shape( 2, 2 ) = shape( 2, 2 ) - 0.5D0 * shape( 2, 3 )
        END IF
        IF( iflag .EQ. 1 ) RETURN
        DO idime = 1, ndime
          fwork = 0.0D0
          DO inode = 1, nnode
            fwork = fwork + shape( 1, inode ) * coren( idime, inode )
          END DO
          xjaco = xjaco + fwork * fwork
        END DO
        xjaco = DSQRT( xjaco )
        DO inode = 1, nnode
          shape( 1, inode ) = shape( 1, inode ) / xjaco
        END DO
        RETURN
        END

        SUBROUTINE Shap2D( shape, coren, xloca, yloca, lnode, xjaco,           &
                           mdime, ndime, nnode, ielem, iflag, nerrc )
!.....
!       P R O B L E M
!             To calculate shape functions and its' derivations.
!             if iflag == 1 derivations to local  coordinate sys.
!             if iflag == 2 derivations to grobal coordinate sys.
!.....
        IMPLICIT  DOUBLE PRECISION( a-h, o-z )
        DIMENSION             fwork( 2, 3 )
        DIMENSION xloct( 4 ), xjaci( 2, 2 )
        DIMENSION yloct( 4 ), shape( 3, 8 )
        DIMENSION lnode( nnode ), coren( mdime, nnode )

        DATA xloct /-0.5D0,  0.5D0, 0.5D0, -0.5D0/
        DATA yloct /-0.5D0, -0.5D0, 0.5D0,  0.5D0/

        DO inode = 1, 4
          shape( 1, inode ) = ( 0.5D0 + yloct( inode ) * yloca ) *             &
                                        xloct( inode )
          shape( 2, inode ) = ( 0.5D0 + xloct( inode ) * xloca ) *             &
                                        yloct( inode )
          shape( 3, inode ) = ( 0.5D0 + xloct( inode ) * xloca ) *             &
                              ( 0.5D0 + yloct( inode ) * yloca )
        END DO
        DO inode = 5, 8
          DO idime = 1, 3
            shape( idime, inode ) = 0.0D0
          END DO
        END DO

        xloc2 = ( 1.0D0 - xloca * xloca ) / 2.0D0
        yloc2 = ( 1.0D0 - yloca * yloca ) / 2.0D0

        IF( nnode .GE. 5 .AND. lnode( 5 ) .NE. 0 ) THEN
          shape( 2, 5 ) = -xloc2
          shape( 1, 5 ) = -xloca * ( 1.0D0 - yloca )
          shape( 3, 5 ) =  xloc2 * ( 1.0D0 - yloca )
        END IF

        IF( nnode .GE. 6 .AND. lnode( 6 ) .NE. 0 ) THEN
          shape( 1, 6 ) =  yloc2
          shape( 2, 6 ) = -yloca * ( 1.0D0 + xloca )
          shape( 3, 6 ) =  yloc2 * ( 1.0D0 + xloca )
        END IF

        IF( nnode .GE. 7 .AND. lnode( 7 ) .NE. 0 ) THEN
          shape( 2, 7 ) =  xloc2
          shape( 1, 7 ) = -xloca * ( 1.0D0 + yloca )
          shape( 3, 7 ) =  xloc2 * ( 1.0D0 + yloca )
        END IF

        IF( nnode .GE. 8 .AND. lnode( 8 ) .NE. 0 ) THEN
          shape( 1, 8 ) = -yloc2
          shape( 2, 8 ) = -yloca * ( 1.0D0 - xloca )
          shape( 3, 8 ) =  yloc2 * ( 1.0D0 - xloca )
        END IF

        DO inode = 1, 4
          jnode = inode + 3
          knode = inode + 4
          IF( inode .EQ. 1 ) jnode = 8
          DO idime = 1, 3
            shape( idime, inode ) = shape( idime, inode ) - (                  &
            shape( idime, jnode ) + shape( idime, knode ) ) / 2.0D0
          END DO
        END DO

        IF( iflag .EQ. 1 ) RETURN

        DO idime = 1, 2
          DO jdime = 1, ndime
            fwork( idime, jdime ) = 0.0D0
            DO inode = 1, nnode
              fwork( idime, jdime ) = fwork( idime, jdime ) +                  &
              shape( idime, inode ) * coren( jdime, inode )
            END DO
          END DO
        END DO

        IF( ndime .EQ. 2 ) THEN
          xjaco = fwork( 1, 1 ) * fwork( 2, 2 ) -                              &
                  fwork( 1, 2 ) * fwork( 2, 1 )
          IF( xjaco .LT. 1.0D-8 ) THEN
            WRITE( 12, 2000 ) ielem
            nerrc = 158436
            RETURN
          END IF
          xjaci( 1, 1 ) =  fwork( 2, 2 ) / xjaco
          xjaci( 2, 2 ) =  fwork( 1, 1 ) / xjaco
          xjaci( 1, 2 ) = -fwork( 1, 2 ) / xjaco
          xjaci( 2, 1 ) = -fwork( 2, 1 ) / xjaco
          DO inode = 1, nnode
            dndgx = shape( 1, inode ) * xjaci( 1, 1 ) +                        &
                    shape( 2, inode ) * xjaci( 1, 2 )
            dndgy = shape( 1, inode ) * xjaci( 2, 1 ) +                        &
                    shape( 2, inode ) * xjaci( 2, 2 )
            shape( 1, inode ) = dndgx
            shape( 2, inode ) = dndgy
          END DO
        ELSE
          dadlx = fwork( 1, 2 ) * fwork( 2, 3 ) -                              &
                  fwork( 1, 3 ) * fwork( 2, 2 )
          dadly = fwork( 1, 3 ) * fwork( 2, 1 ) -                              &
                  fwork( 1, 1 ) * fwork( 2, 3 )
          dadlz = fwork( 1, 1 ) * fwork( 2, 2 ) -                              &
                  fwork( 1, 2 ) * fwork( 2, 1 )
          xjaco = DSQRT( dadlx * dadlx + dadly * dadly +                       &
                         dadlz * dadlz )
        END IF
2000    FORMAT( 2x, '第', I5, '单元的雅可比行列式值为零或负值！' )
        RETURN
        END

        SUBROUTINE Shap3D( xloca, yloca, zloca, coren, shape, xjaco,           &
                           ndime, nnode, lnode, ielem, iflag, nerrc )
!......
!       P R O G R A M
!             To set shape function of three dimensional elements with
!             8 to 21 nodes.
!......
        IMPLICIT DOUBLE PRECISION ( a-h, o-z )
        DIMENSION xjaxi( 3,  3 ), xloct( 8 )
        DIMENSION xjaxm( 3,  3 ), yloct( 8 ), lnode( nnode )
        DIMENSION shape( 4, 21 ), zloct( 8 ), coren( ndime, nnode )
        data xloct/-1., 1., 1.,-1.,-1., 1., 1.,-1./
        data yloct/-1.,-1., 1., 1.,-1.,-1., 1., 1./
        data zloct/-1.,-1.,-1.,-1., 1., 1., 1., 1./
!......
!       for 8-node linear shape functions
!......
        DO inode = 1, 21
         DO idime = 1, 4
          shape( idime, inode ) = 0.0D0
         END DO
        END DO
        IF( nnode .lt. 8 ) RETURN
        DO inode = 1, 8
         xtemp = 1.0D0 + xloct( inode ) * xloca
         ytemp = 1.0D0 + yloct( inode ) * yloca
         ztemp = 1.0D0 + zloct( inode ) * zloca
         shape( 4, inode ) = 0.125D0 * xtemp * ytemp * ztemp
         shape( 1, inode ) = 0.125D0 * ytemp * ztemp * xloct( inode )
         shape( 2, inode ) = 0.125D0 * ztemp * xtemp * yloct( inode )
         shape( 3, inode ) = 0.125D0 * xtemp * ytemp * zloct( inode )
        END DO
!.....
!       midside nodes ( serendipity )
!.....
        xlocm =   1.0D0 - xloca
        ylocm =   1.0D0 - yloca
        zlocm =   1.0D0 - zloca
        xlocp =   1.0D0 + xloca
        ylocp =   1.0D0 + yloca
        zlocp =   1.0D0 + zloca
        xtemp = ( 1.0D0 - xloca * xloca ) / 4.
        ytemp = ( 1.0D0 - yloca * yloca ) / 4.
        ztemp = ( 1.0D0 - zloca * zloca ) / 4.
        IF( nnode .GE. 9 .AND. lnode( 9 ) .NE. 0 ) THEN
         shape( 4, 9 ) = ylocm * zlocm * xtemp
         shape( 1, 9 ) =-0.5D0 * xloca * ylocm * zlocm
         shape( 2, 9 ) =-xtemp * zlocm
         shape( 3, 9 ) =-xtemp * ylocm
        END IF
        IF( nnode .GE. 10 .AND. lnode( 10 ) .NE. 0 ) THEN
         shape( 4, 10 ) = ylocp * zlocm * xtemp
         shape( 1, 10 ) =-0.5D0 * xloca * ylocp * zlocm
         shape( 2, 10 ) = xtemp * zlocm
         shape( 3, 10 ) =-xtemp * ylocp
        END IF
        IF( nnode .GE. 11 .AND. lnode( 11 ) .NE. 0 ) THEN
         shape( 4, 11 ) = ylocp * zlocp * xtemp
         shape( 1, 11 ) =-0.5D0 * xloca * ylocp * zlocp
         shape( 2, 11 ) = xtemp * zlocp
         shape( 3, 11 ) = xtemp * ylocp
        END IF
        IF( nnode .GE. 12 .AND. lnode( 12 ) .NE. 0 ) THEN
         shape( 4, 12 ) = ylocm * zlocp * xtemp
         shape( 1, 12 ) =-0.5D0 * xloca * ylocm * zlocp
         shape( 2, 12 ) =-xtemp * zlocp
         shape( 3, 12 ) = xtemp * ylocm
        END IF
        IF( nnode .GE. 13 .AND. lnode( 13 ) .NE. 0 ) THEN
         shape( 4, 13 ) = zlocm * xlocp * ytemp
         shape( 1, 13 ) = ytemp * zlocm
         shape( 2, 13 ) =-0.5D0 * yloca * zlocm *xlocp
         shape( 3, 13 ) =-ytemp * xlocp
        END IF
        IF( nnode .GE. 14 .AND. lnode( 14 ) .NE. 0 ) THEN
         shape( 4, 14 ) = xlocp * zlocp * ytemp
         shape( 1, 14 ) = ytemp * zlocp
         shape( 2, 14 ) =-0.5D0 * yloca * xlocp * zlocp
         shape( 3, 14 ) = ytemp * xlocp
        END IF
        IF( nnode .GE. 15 .AND. lnode( 15 ) .NE. 0 ) THEN
         shape( 4, 15 ) = xlocm * zlocp * ytemp
         shape( 1, 15 ) =-zlocp * ytemp
         shape( 2, 15 ) =-0.5D0 * yloca * xlocm * zlocp
         shape( 3, 15 ) = ytemp * xlocm
        END IF
        IF( nnode .GE. 16 .AND. lnode( 16 ) .NE. 0 ) THEN
         shape( 4, 16 ) = xlocm * zlocm * ytemp
         shape( 1, 16 ) =-zlocm * ytemp
         shape( 2, 16 ) =-0.5D0 * yloca * xlocm * zlocm
         shape( 3, 16 ) =-xlocm * ytemp
        END IF
        IF( nnode .GE. 17 .AND. lnode( 17 ) .NE. 0 ) THEN
         shape( 4, 17 ) = ylocm * xlocm * ztemp
         shape( 1, 17 ) =-ylocm * ztemp
         shape( 2, 17 ) =-xlocm * ztemp
         shape( 3, 17 ) =-0.5D0 * zloca * xlocm * ylocm
        END IF
        IF( nnode .GE. 18 .AND. lnode( 18 ) .NE. 0 ) THEN
         shape( 4, 18 ) = xlocp * ylocm * ztemp
         shape( 1, 18 ) = ylocm * ztemp
         shape( 2, 18 ) =-xlocp * ztemp
         shape( 3, 18 ) =-0.5D0 * zloca * xlocp * ylocm
        END IF
        IF( nnode .GE. 19 .AND. lnode( 19 ) .NE. 0 ) THEN
         shape( 4, 19 ) = xlocp * ylocp * ztemp
         shape( 1, 19 ) = ylocp * ztemp
         shape( 2, 19 ) = xlocp * ztemp
         shape( 3, 19 ) =-0.5D0 * zloca * xlocp * ylocp
        END IF
        IF( nnode .GE. 20 .AND. lnode( 20 ) .NE. 0 ) THEN
         shape( 4, 20 ) = xlocm * ylocp * ztemp
         shape( 1, 20 ) =-ylocp * ztemp
         shape( 2, 20 ) = xlocm * ztemp
         shape( 3, 20 ) =-0.5D0 * zloca * xlocm * ylocp
        END IF
        IF( nnode .GE. 21 .AND. lnode( 21 ) .NE. 0 ) THEN
         shape( 4, 21 ) = 64. * xtemp * ytemp * ztemp
         shape( 1, 21 ) =-32. * xloca * ytemp * ztemp
         shape( 2, 21 ) =-32. * yloca * xtemp * ztemp
         shape( 3, 21 ) =-32. * xtemp * ytemp * zloca
         DO idime = 1, 4
          DO inode = 1, 8
           shape( idime, inode ) = shape( idime, inode ) -                     &
                          0.25D0 * shape( idime,    21 )
          END DO
          DO inode = 9, 20
           IF( nnode .GE. 9 .AND. lnode( inode ) .NE. 0 )                      &
            shape( idime, inode ) = shape( idime, inode ) -                    &