shape.f90

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

                            0.5D0 * shape( idime,    21 )
          END DO
         END DO
        END IF
        DO idime = 1, 4
         shape( idime,  1 ) = shape( idime,  1 ) - 0.5D0 * (                   &
         shape( idime,  9 ) + shape( idime, 16 ) + shape( idime, 17 ) )
         shape( idime,  2 ) = shape( idime,  2 ) - 0.5D0 * (                   &
         shape( idime,  9 ) + shape( idime, 18 ) + shape( idime, 13 ) )
         shape( idime,  3 ) = shape( idime,  3 ) - 0.5D0 * (                   &
         shape( idime, 13 ) + shape( idime, 19 ) + shape( idime, 10 ) )
         shape( idime,  4 ) = shape( idime,  4 ) - 0.5D0 * (                   &
         shape( idime, 10 ) + shape( idime, 16 ) + shape( idime, 20 ) )
         shape( idime,  5 ) = shape( idime,  5 ) - 0.5D0 * (                   &
         shape( idime, 12 ) + shape( idime, 15 ) + shape( idime, 17 ) )
         shape( idime,  6 ) = shape( idime,  6 ) - 0.5D0 * (                   &
         shape( idime, 12 ) + shape( idime, 14 ) + shape( idime, 18 ) )
         shape( idime,  7 ) = shape( idime,  7 ) - 0.5D0 * (                   &
         shape( idime, 11 ) + shape( idime, 14 ) + shape( idime, 19 ) )
         shape( idime,  8 ) = shape( idime,  8 ) - 0.5D0 * (                   &
         shape( idime, 11 ) + shape( idime, 15 ) + shape( idime, 20 ) )
        END DO
        DO idime = 1, ndime
         DO jdime = 1, 3
          ftemp = 0.0D0
          DO inode = 1, nnode
           ftemp = ftemp + coren( idime, inode ) *                             &
                           shape( jdime, inode )
          END DO
          xjaxm( jdime, idime ) = ftemp
         END DO
        END DO
        xjaxi( 1, 1 ) = xjaxm( 2, 2 ) * xjaxm( 3, 3 ) -                        &
                        xjaxm( 2, 3 ) * xjaxm( 3, 2 )
        xjaxi( 2, 1 ) = xjaxm( 2, 3 ) * xjaxm( 3, 1 ) -                        &
                        xjaxm( 2, 1 ) * xjaxm( 3, 3 )
        xjaxi( 3, 1 ) = xjaxm( 2, 1 ) * xjaxm( 3, 2 ) -                        &
                        xjaxm( 2, 2 ) * xjaxm( 3, 1 )
        xjaco = 0.0D0
        DO idime = 1, 3
         xjaco = xjaco + xjaxm( 1, idime ) * xjaxi( idime, 1 )
        END DO
        IF( iflag .EQ. 1 ) RETURN
        xjaxi( 1, 2 ) = xjaxm( 3, 2 ) * xjaxm( 1, 3 ) -                        &
                        xjaxm( 3, 3 ) * xjaxm( 1, 2 )
        xjaxi( 2, 2 ) = xjaxm( 3, 3 ) * xjaxm( 1, 1 ) -                        &
                        xjaxm( 3, 1 ) * xjaxm( 1, 3 )
        xjaxi( 3, 2 ) = xjaxm( 3, 1 ) * xjaxm( 1, 2 ) -                        &
                        xjaxm( 3, 2 ) * xjaxm( 1, 1 )
        xjaxi( 1, 3 ) = xjaxm( 1, 2 ) * xjaxm( 2, 3 ) -                        &
                        xjaxm( 1, 3 ) * xjaxm( 2, 2 )
        xjaxi( 2, 3 ) = xjaxm( 1, 3 ) * xjaxm( 2, 1 ) -                        &
                        xjaxm( 1, 1 ) * xjaxm( 2, 3 )
        xjaxi( 3, 3 ) = xjaxm( 1, 1 ) * xjaxm( 2, 2 ) -                        &
                        xjaxm( 1, 2 ) * xjaxm( 2, 1 )
        IF( xjaco .LE. 1.0D-8 ) THEN
          WRITE( 12, 2000 ) ielem
          nerrc = 158436
          RETURN
        END IF
        DO idime = 1, 3
         DO jdime = 1, 3
          xjaxi( idime, jdime ) = xjaxi( idime, jdime ) / xjaco
         END DO
        END DO
!.....
!       form global derivatives
!.....
        DO inode = 1, nnode
         xtemp = 0.0D0
         ytemp = 0.0D0
         ztemp = 0.0D0
         DO idime = 1, 3
          xtemp = xtemp + xjaxi( 1, idime ) * shape( idime, inode )
          ytemp = ytemp + xjaxi( 2, idime ) * shape( idime, inode )
          ztemp = ztemp + xjaxi( 3, idime ) * shape( idime, inode )
         END DO
         shape( 1, inode ) = xtemp
         shape( 2, inode ) = ytemp
         shape( 3, inode ) = ztemp
        END DO
2000    FORMAT( 2x, '第', I5, '单元的雅可比行列式值为零或负值！' )
        RETURN
        END

        SUBROUTINE TShap2D( coren, shape, lnode, xloca, yloca,                 &
                            xjaco, ielem, mdime, nerrc, nnode )
!........
!       模块功能
!           计算三节点平面单元的的形函数及其导数
!........
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION shape( 3, 6 ), coren( mdime, nnode )
        IF( nnode .GT. 3 ) THEN
          WRITE( 12, 2000 )
          nerrc = 312534
          RETURN
        END IF
        IF( mdime .LT. 2 ) THEN
          WRITE( 12, 2200 )
          nerrc = 35933
          RETURN
        END IF
        xjaco = coren(1,2) * coren(2,3) - coren(1,3) * coren(2,2) +            &
                coren(1,3) * coren(2,1) - coren(1,1) * coren(2,3) +            &
                coren(1,1) * coren(2,2) - coren(1,2) * coren(2,1)
        IF( xjaco .LT. 1.0D-8 ) THEN
          WRITE( 12, 2400 ) ielem
          nerrc = 35945
          RETURN
        END IF
        shape( 1, 1 ) = ( coren( 2, 2 ) - coren( 2, 3 ) ) / xjaco
        shape( 1, 2 ) = ( coren( 2, 3 ) - coren( 2, 1 ) ) / xjaco
        shape( 1, 3 ) = ( coren( 2, 1 ) - coren( 2, 2 ) ) / xjaco
        shape( 2, 1 ) = ( coren( 1, 3 ) - coren( 1, 2 ) ) / xjaco
        shape( 2, 2 ) = ( coren( 1, 1 ) - coren( 1, 3 ) ) / xjaco
        shape( 2, 3 ) = ( coren( 1, 2 ) - coren( 1, 1 ) ) / xjaco
        shape( 3, 1 ) = xloca
        shape( 3, 2 ) = yloca
        shape( 3, 3 ) = 1.0D0 - xloca - yloca
2000    FORMAT( 2x, '错误：TShap2D只适用于三节点单元。' )
2200    FORMAT( 2x, '错误：TShap2D只适用于二维和三维问题。' )
2400    FORMAT( 2x, '错误：三角形单元的面积非正。' )
        RETURN
        END

        SUBROUTINE TShap3D( coren, shape, lnode, xloca, yloca, zloca,          &
                            xjaco, ielem, mdime, nerrc, nnode )
!.......
!       模块功能
!           计算四节点四面体形函数及其导数
!.......
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION shape( 4, 10 ), coren( mdime, nnode )
        DIMENSION vwork( 3,  4 )
        IF( nnode .NE. 4 ) THEN
          WRITE( 12, 2000 )
          nerrc = 312533
          RETURN
        END IF
        IF( mdime .NE. 3 ) THEN
          WRITE( 12, 2200 )
          nerrc = 359333
          RETURN
        END IF

        DO idime = 1, 3
          vwork( idime, 1 ) = 0.0D0
          vwork( idime, 2 ) = coren( idime, 2 ) - coren( idime, 1 )
          vwork( idime, 3 ) = coren( idime, 3 ) - coren( idime, 1 )
          vwork( idime, 4 ) = coren( idime, 4 ) - coren( idime, 1 )
        END DO

        xjaco = vwork( 1, 2 ) * vwork( 2, 3 ) * vwork( 3, 4 ) +                &
                vwork( 1, 3 ) * vwork( 2, 4 ) * vwork( 3, 2 ) +                &
                vwork( 1, 4 ) * vwork( 2, 2 ) * vwork( 3, 3 ) -                &
                vwork( 1, 4 ) * vwork( 2, 3 ) * vwork( 3, 2 ) -                &
                vwork( 1, 3 ) * vwork( 2, 2 ) * vwork( 3, 4 ) -                &
                vwork( 1, 2 ) * vwork( 2, 4 ) * vwork( 3, 3 )

        IF( xjaco .LT. 1.0D-8 ) THEN
          WRITE( 12, 2400 ) ielem
          nerrc = 35945
          RETURN
        END IF

        shape( 1, 1 ) = coren(2, 2) * ( coren(3, 4) - coren(3, 3) ) +          &
                        coren(2, 3) * ( coren(3, 2) - coren(3, 4) ) +          &
                        coren(2, 4) * ( coren(3, 3) - coren(3, 2) )
        shape( 2, 1 ) = coren(3, 2) * ( coren(1, 4) - coren(1, 3) ) +          &
                        coren(3, 3) * ( coren(1, 2) - coren(1, 4) ) +          &
                        coren(3, 4) * ( coren(1, 3) - coren(1, 2) )
        shape( 3, 1 ) = coren(1, 2) * ( coren(2, 4) - coren(2, 3) ) +          &
                        coren(1, 3) * ( coren(2, 2) - coren(2, 4) ) +          &
                        coren(1, 4) * ( coren(2, 3) - coren(2, 2) )
        shape( 4, 1 ) = xloca

        shape( 1, 2 ) = coren(2, 3) * ( coren(3, 4) - coren(3, 1) ) +          &
                        coren(2, 4) * ( coren(3, 1) - coren(3, 3) ) +          &
                        coren(2, 1) * ( coren(3, 3) - coren(3, 4) )
        shape( 2, 2 ) = coren(3, 3) * ( coren(1, 4) - coren(1, 1) ) +          &
                        coren(3, 4) * ( coren(1, 1) - coren(1, 3) ) +          &
                        coren(3, 1) * ( coren(1, 3) - coren(1, 4) )
        shape( 3, 2 ) = coren(1, 3) * ( coren(2, 4) - coren(2, 1) ) +          &
                        coren(1, 4) * ( coren(2, 1) - coren(2, 3) ) +          &
                        coren(1, 1) * ( coren(2, 3) - coren(2, 4) )
        shape( 4, 2 ) = yloca


        shape( 1, 3 ) = coren(2, 4) * ( coren(3, 2) - coren(3, 1) ) +          &
                        coren(2, 1) * ( coren(3, 4) - coren(3, 2) ) +          &
                        coren(2, 2) * ( coren(3, 1) - coren(3, 4) )
        shape( 2, 3 ) = coren(3, 4) * ( coren(1, 2) - coren(1, 1) ) +          &
                        coren(3, 1) * ( coren(1, 4) - coren(1, 2) ) +          &
                        coren(3, 2) * ( coren(1, 1) - coren(1, 4) )
        shape( 3, 3 ) = coren(1, 4) * ( coren(2, 2) - coren(2, 1) ) +          &
                        coren(1, 1) * ( coren(2, 4) - coren(2, 2) ) +          &
                        coren(1, 2) * ( coren(2, 1) - coren(2, 4) )
        shape( 4, 3 ) = zloca


        shape( 1, 4 ) = coren(2, 1) * ( coren(3, 2) - coren(3, 3) ) +          &
                        coren(2, 2) * ( coren(3, 3) - coren(3, 1) ) +          &
                        coren(2, 3) * ( coren(3, 1) - coren(3, 2) )
        shape( 2, 4 ) = coren(3, 1) * ( coren(1, 2) - coren(1, 3) ) +          &
                        coren(3, 2) * ( coren(1, 3) - coren(1, 1) ) +          &
                        coren(3, 3) * ( coren(1, 1) - coren(1, 2) )
        shape( 3, 4 ) = coren(1, 1) * ( coren(2, 2) - coren(2, 3) ) +          &
                        coren(1, 2) * ( coren(2, 3) - coren(2, 1) ) +          &
                        coren(1, 3) * ( coren(2, 1) - coren(2, 2) )
        shape( 4, 4 ) = 1.0D0 - xloca - yloca - zloca

        DO idime = 1, 3
          DO inode = 1, 4
            shape( idime, inode ) = shape( idime, inode ) / xjaco
          END DO
        END DO

2000    FORMAT( 2x, '错误：TShap3D只适用于四节点四面体单元。' )
2200    FORMAT( 2x, '错误：TShap3D只适用于三维问题。' )
2400    FORMAT( 2x, '错误：三角形单元的面积非正。' )
        END

        SUBROUTINE TCShap3D( coren, shape, lnode, xloca, yloca,                &
                             zloca, xjaco, ielem, ndime, nerrc,                &
                             nnode, iflag )
!........
!       模块功能
!           计算六节点三棱柱单元的的形函数及其导数
!........
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION xjaxm( 3, 3 ), lnode( nnode ), shape( 4, 15 )
        DIMENSION xjaxi( 3, 3 ), coren( ndime, nnode )
        IF( nnode .GT. 6 ) THEN
          WRITE( 12, 2000 )
          nerrc = 312534
          RETURN
        END IF
        zlocd = 1.0D0 - zloca
        zlocu = 1.0D0 + zloca
        ztemp = 1.0D0 - xloca - yloca
        shape( 1, 1 ) = zlocd / 2.0D0
        shape( 2, 1 ) = 0.0D0
        shape( 3, 1 ) =-xloca / 2.0D0
        shape( 4, 1 ) = xloca * zlocd / 2.0D0
        shape( 1, 2 ) = 0.0D0
        shape( 2, 2 ) = zlocd / 2.0D0
        shape( 3, 2 ) =-yloca / 2.0D0
        shape( 4, 2 ) = yloca * zlocd / 2.0D0
        shape( 1, 3 ) =-1.0D0 * zlocd / 2.0D0
        shape( 2, 3 ) =-1.0D0 * zlocd / 2.0D0
        shape( 3, 3 ) =-ztemp / 2.0D0
        shape( 4, 3 ) = ztemp * zlocd / 2.0D0
        shape( 1, 4 ) = zlocu / 2.0D0
        shape( 2, 4 ) = 0.0D0
        shape( 3, 4 ) = xloca / 2.0D0
        shape( 4, 4 ) = xloca * zlocu / 2.0D0
        shape( 1, 5 ) = 0.0D0
        shape( 2, 5 ) = zlocu / 2.0D0
        shape( 3, 5 ) = yloca / 2.0D0
        shape( 4, 5 ) = yloca * zlocu / 2.0D0
        shape( 1, 6 ) =-zlocu / 2.0D0
        shape( 2, 6 ) =-zlocu / 2.0D0
        shape( 3, 6 ) = ztemp / 2.0D0
        shape( 4, 6 ) = ztemp * zlocu / 2.0D0
        IF( iflag .EQ. 1 ) RETURN
!.......计算Jacobi矩阵
        DO idime = 1, ndime
          DO jdime = 1, 3
            ftemp = 0.0D0
            DO inode = 1, nnode
              ftemp = ftemp + coren( idime, inode ) *                          &
                              shape( jdime, inode )
            END DO
            xjaxm( jdime, idime ) = ftemp
          END DO
        END DO
!.......计算Jacobi矩阵第一行元素的代数余子式
        xjaxi( 1, 1 ) = xjaxm( 2, 2 ) * xjaxm( 3, 3 ) -                        &
                        xjaxm( 2, 3 ) * xjaxm( 3, 2 )
        xjaxi( 2, 1 ) = xjaxm( 2, 3 ) * xjaxm( 3, 1 ) -                        &
                        xjaxm( 2, 1 ) * xjaxm( 3, 3 )
        xjaxi( 3, 1 ) = xjaxm( 2, 1 ) * xjaxm( 3, 2 ) -                        &