⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 elmt300.f90

📁 边界元程序,供力学工作者参考,希望对大家有所帮助
💻 F90
📖 第 1 页 / 共 5 页
字号:
        DIMENSION xaxis( 3 ), dmatx( 6, 6 )
        DIMENSION yaxis( 3 ), corep( 3, 8 )
        DIMENSION zaxis( 3 ), tmatx( 6, 6 ), dwork( 6, 6 )

        CALL Axis300( xaxis, yaxis, zaxis, xloca, yloca, corep, ielem )
        CALL TMatrix300( xaxis, yaxis, zaxis, tmatx )

        DO istre = 1, 6
          DO jstre = 1, 6
            dwork( istre, jstre ) = 0.0D0
            DO kstre = 1, 6
              dwork( istre, jstre ) = dwork( istre, jstre ) +                  &
              dmatx( istre, kstre ) * tmatx( kstre, jstre )
            END DO
          END DO
        END DO

        DO istre = 1, 6
          DO jstre = 1, 6
            dmatx( istre, jstre ) = 0.0D0
            DO kstre = 1, 6
              dmatx( istre, jstre ) = dmatx( istre, jstre ) +                  &
              tmatx( kstre, istre ) * dwork( kstre, jstre )
            END DO
          END DO
        END DO

        RETURN
        END

        SUBROUTINE Axis300( xaxis, yaxis, zaxis, xloca, yloca, corep,          &
                            ielem )
        USE CtrlData
        USE ElmtData
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION xaxis( 3 ), lnodp( 8 )
        DIMENSION yaxis( 3 ), corep( 3, 8 )
        DIMENSION zaxis( 3 ), shape( 3, 8 )

        iflag = 1
        DO inode = 1, 8
          lnodp( inode ) = inode
        END DO
        CALL Shap2D( shape, corep, xloca, yloca, lnodp, xjaco,                 &
                     ndime, ndime,     8, ielem, iflag, nerrc )
        IF( nerrc .NE. 0 ) RETURN
        DO idime = 1, ndime
          xaxis( idime ) = 0.0D0
          yaxis( idime ) = 0.0D0
          DO jnode = 1, 8
            xaxis(    idime ) = xaxis( idime        ) +                        &
            shape( 1, jnode ) * corep( idime, jnode )
            yaxis(    idime ) = yaxis( idime        ) +                        &
            shape( 2, jnode ) * corep( idime, jnode )
          END DO
        END DO
        zaxis( 1 ) = xaxis( 2 ) * yaxis( 3 ) - xaxis( 3 ) * yaxis( 2 )
        zaxis( 2 ) = xaxis( 3 ) * yaxis( 1 ) - xaxis( 1 ) * yaxis( 3 )
        zaxis( 3 ) = xaxis( 1 ) * yaxis( 2 ) - xaxis( 2 ) * yaxis( 1 )
        yaxis( 1 ) = zaxis( 2 ) * xaxis( 3 ) - zaxis( 3 ) * xaxis( 2 )
        yaxis( 2 ) = zaxis( 3 ) * xaxis( 1 ) - zaxis( 1 ) * xaxis( 3 )
        yaxis( 3 ) = zaxis( 1 ) * xaxis( 2 ) - zaxis( 2 ) * xaxis( 1 )
        xleng = xaxis( 1 ) * xaxis( 1 ) + xaxis( 2 ) * xaxis( 2 ) +            &
                xaxis( 3 ) * xaxis( 3 )
        yleng = yaxis( 1 ) * yaxis( 1 ) + yaxis( 2 ) * yaxis( 2 ) +            &
                yaxis( 3 ) * yaxis( 3 )
        zleng = zaxis( 1 ) * zaxis( 1 ) + zaxis( 2 ) * zaxis( 2 ) +            &
                zaxis( 3 ) * zaxis( 3 )
        DO idime = 1, 3
          xaxis( idime ) = xaxis( idime ) / DSQRT( xleng )
          yaxis( idime ) = yaxis( idime ) / DSQRT( yleng )
          zaxis( idime ) = zaxis( idime ) / DSQRT( zleng )
        END DO
        RETURN
        END

        SUBROUTINE TMatrix300( xaxis, yaxis, zaxis, tmatx )
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION xaxis( 3 ), yaxis( 3 )
        DIMENSION zaxis( 3 ), tmatx( 6, 6 )
        tmatx( 1, 1 ) = xaxis( 1 ) * xaxis( 1 )
        tmatx( 1, 2 ) = xaxis( 2 ) * xaxis( 2 )
        tmatx( 1, 3 ) = xaxis( 3 ) * xaxis( 3 )
        tmatx( 1, 4 ) = xaxis( 2 ) * xaxis( 3 )
        tmatx( 1, 5 ) = xaxis( 3 ) * xaxis( 1 )
        tmatx( 1, 6 ) = xaxis( 1 ) * xaxis( 2 )

        tmatx( 2, 1 ) = yaxis( 1 ) * yaxis( 1 )
        tmatx( 2, 2 ) = yaxis( 2 ) * yaxis( 2 )
        tmatx( 2, 3 ) = yaxis( 3 ) * yaxis( 3 )
        tmatx( 2, 4 ) = yaxis( 2 ) * yaxis( 3 )
        tmatx( 2, 5 ) = yaxis( 3 ) * yaxis( 1 )
        tmatx( 2, 6 ) = yaxis( 1 ) * yaxis( 2 )

        tmatx( 3, 1 ) = zaxis( 1 ) * zaxis( 1 )
        tmatx( 3, 2 ) = zaxis( 2 ) * zaxis( 2 )
        tmatx( 3, 3 ) = zaxis( 3 ) * zaxis( 3 )
        tmatx( 3, 4 ) = zaxis( 2 ) * zaxis( 3 )
        tmatx( 3, 5 ) = zaxis( 3 ) * zaxis( 1 )
        tmatx( 3, 6 ) = zaxis( 1 ) * zaxis( 2 )

        tmatx( 4, 1 ) = yaxis( 1 ) * zaxis( 1 ) * 2.0D0
        tmatx( 4, 2 ) = yaxis( 2 ) * zaxis( 2 ) * 2.0D0
        tmatx( 4, 3 ) = yaxis( 3 ) * zaxis( 3 ) * 2.0D0
        tmatx( 4, 4 ) = yaxis( 2 ) * zaxis( 3 ) +                              &
                        zaxis( 2 ) * yaxis( 3 )
        tmatx( 4, 5 ) = yaxis( 3 ) * zaxis( 1 ) +                              &
                        zaxis( 3 ) * yaxis( 1 )
        tmatx( 4, 6 ) = yaxis( 1 ) * zaxis( 2 ) +                              &
                        zaxis( 1 ) * yaxis( 2 )

        tmatx( 5, 1 ) = xaxis( 1 ) * zaxis( 1 ) * 2.0D0
        tmatx( 5, 2 ) = xaxis( 2 ) * zaxis( 2 ) * 2.0D0
        tmatx( 5, 3 ) = xaxis( 3 ) * zaxis( 3 ) * 2.0D0
        tmatx( 5, 4 ) = xaxis( 2 ) * zaxis( 3 ) +                              &
                        zaxis( 2 ) * xaxis( 3 )
        tmatx( 5, 5 ) = xaxis( 3 ) * zaxis( 1 ) +                              &
                        zaxis( 3 ) * xaxis( 1 )
        tmatx( 5, 6 ) = xaxis( 1 ) * zaxis( 2 ) +                              &
                        zaxis( 1 ) * xaxis( 2 )

        tmatx( 6, 1 ) = xaxis( 1 ) * yaxis( 1 ) * 2.0D0
        tmatx( 6, 2 ) = xaxis( 2 ) * yaxis( 2 ) * 2.0D0
        tmatx( 6, 3 ) = xaxis( 3 ) * yaxis( 3 ) * 2.0D0
        tmatx( 6, 4 ) = xaxis( 2 ) * yaxis( 3 ) +                              &
                        yaxis( 2 ) * xaxis( 3 )
        tmatx( 6, 5 ) = xaxis( 3 ) * yaxis( 1 ) +                              &
                        yaxis( 3 ) * xaxis( 1 )
        tmatx( 6, 6 ) = xaxis( 1 ) * yaxis( 2 ) +                              &
                        yaxis( 1 ) * xaxis( 2 )
        RETURN
        END

        SUBROUTINE PlotMesh300( ielem, imats )
        USE CtrlData
        USE MeshData
        USE ElmtData
        USE ExtraMesh
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION lplot( 8, 6 ), lwork( 8 )

        DATA lplot/ 4, 3, 2, 1, 10, 13,  9, 16,                                &
                    5, 6, 7, 8, 12, 14, 11, 15,                                &
                    1, 2, 6, 5,  9, 18, 12, 17,                                &
                    2, 3, 7, 6, 13, 19, 14, 18,                                &
                    3, 4, 8, 7, 10, 20, 11, 19,                                &
                    4, 1, 5, 8, 16, 17, 15, 20 /

        knode = 8
        kelem = lprtx( ielem )
        IF( nnode .EQ. 8 ) knode = 4
        nbloc = props( 2, imats ) + 0.5D0
        DO ibloc = 1, nbloc
          DO iplan = 1, 6
            DO inode = 1, knode
              jnode = lplot( inode, iplan )
              IF( jnode .LE. nnode ) THEN
                lwork( inode ) = lnodx( jnode, kelem )
                IF( lwork( inode ) .GT. 0 ) THEN
                  lwork( inode ) = lwork( inode ) + npoin
                ELSE IF( lwork( inode ) .LT. 0 ) THEN
                  lwork( inode ) =-lwork( inode )
                END IF
              ELSE
                lwork( inode ) = 0
              END IF
            END DO
            WRITE( 15, 2000 ) 4, ielem, imats, knode,                          &
                  ( lwork( inode ), inode = 1, knode )
          END DO
          kelem = kelem + 1
        END DO
2000    FORMAT( I3, 1X, I5, 1X, I4, 1X, I3, 8( 1X, I5 ) )
        END

        SUBROUTINE ExtraMesh300( imats )
        USE CtrlData
        USE MeshData
        USE ElmtData
        USE ExtraMesh
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION shape( 4, 21 ), point(  3 ), lplot( 8, 6 )
        DIMENSION coreb( 3, 21 ), lnodb( 20 )
        DATA lplot/ 4, 3, 2, 1, 10, 13,  9, 16,                                &
                    5, 6, 7, 8, 12, 14, 11, 15,                                &
                    1, 2, 6, 5,  9, 18, 12, 17,                                &
                    2, 3, 7, 6, 13, 19, 14, 18,                                &
                    3, 4, 8, 7, 10, 20, 11, 19,                                &
                    4, 1, 5, 8, 16, 17, 15, 20 /

        iflag = 1
        nbloc = props( 2, imats ) + 0.5D0
        WRITE( 25 ) nbloc
        DO ibloc = 1, nbloc
          nelmx = nelmx + 1
          igeob = props( ibloc * 2 + 2, imats ) + 0.5D0
          imatb = props( ibloc * 2 + 3, imats ) + 0.5D0
          CALL Bloc300( coreb, lnodb, igeob, nnodb )
          WRITE( 25 ) lnode
          DO inode = 1, nnode
            IF( lnode( inode ) .NE. 0 ) THEN
              IF( lnodb( inode ) .NE. 0 ) THEN
                xloca = coreb( 1, inode )
                yloca = coreb( 2, inode )
                zloca = coreb( 3, inode )
              ELSE
                DO iplan = 1, 6
                  DO imidp = 5, 8
                    IF( lplot( imidp, iplan ) .EQ. inode ) THEN
                      inodb = imidp - 4
                      inoda = inodb + 1
                      IF( inoda .GT. 4 ) inoda = inoda - 4
                      inodb = lplot( inodb, iplan )
                      inoda = lplot( inoda, iplan )
                      xloca = coreb( 1, inodb ) + coreb( 1, inoda )
                      yloca = coreb( 2, inodb ) + coreb( 2, inoda )
                      zloca = coreb( 3, inodb ) + coreb( 3, inoda )
                      xloca = xloca / 2.0D0
                      yloca = yloca / 2.0D0
                      zloca = zloca / 2.0D0
                    END IF
                  END DO
                END DO
              END IF
              CALL Shap3D( xloca, yloca, zloca, coren, shape, xjaco,           &
                           ndime, nnode, lnode, ielem, iflag, nerrc )
              IF( nerrc .GT. 0 ) RETURN
              DO idime = 1, ndime
                point( idime ) = 0.0D0
                DO jnode = 1, nnode
                  point(    idime ) = point( idime        ) +                  &
                  shape( 4, jnode ) * coren( idime, jnode )
                END DO
              END DO
              nnodx = nnodx + 1
              WRITE( 25 ) ( point( idime ), idime = 1, ndime )
            END IF
          END DO
        END DO
        END

        SUBROUTINE ExtraDisp300( ielem, imats )
        USE CtrlData
        USE MeshData
        USE ElmtData
        USE GlobData
        USE ExtraMesh
        IMPLICIT DOUBLE PRECISION( a-h, o-z )
        DIMENSION shape( 4, 21 ), lplot( 8, 6 )
        DIMENSION coreb( 3, 20 ), lnodb(   20 )

        DATA lplot/ 4, 3, 2, 1, 10, 13,  9, 16,                                &
                    5, 6, 7, 8, 12, 14, 11, 15,                                &
                    1, 2, 6, 5,  9, 18, 12, 17,                                &
                    2, 3, 7, 6, 13, 19, 14, 18,                                &
                    3, 4, 8, 7, 10, 20, 11, 19,                                &
                    4, 1, 5, 8, 16, 17, 15, 20 /
        iflag = 1
        kelem = lprtx( ielem )
        nbloc = props( 2, imats ) + 0.5D0

        DO ibloc = 1, nbloc
          igeob = props( ibloc * 2 + 2, imats ) + 0.5D0
          imatb = props( ibloc * 2 + 3, imats ) + 0.5D0
          CALL Bloc300( coreb, lnodb, igeob, nnodb )
          DO inode = 1, nnode
            inodx = lnodx( inode, kelem )
            IF( inodx .GT. 0 ) THEN
              IF( lnodb( inode ) .NE. 0 ) THEN
                xloca = coreb( 1, inode )
                yloca = coreb( 2, inode )
                zloca = coreb( 3, inode )
              ELSE
                DO iplan = 1, 6
                  DO imidp = 5, 8
                    IF( lplot( imidp, iplan ) .EQ. inode ) THEN
                      inodb = imidp - 4
                      inoda = inodb + 1
                      IF( inoda .GT. 4 ) inoda = inoda - 4
                      inodb = lplot( inodb, iplan )
                      inoda = lplot( inoda, iplan )
                      xloca = coreb( 1, inodb ) + coreb( 1, inoda )
                      yloca = coreb( 2, inodb ) + coreb( 2, inoda )
                      zloca = coreb( 3, inodb ) + coreb( 3, inoda )
                      xloca = xloca / 2.0D0
                      yloca = yloca / 2.0D0
                      zloca = zloca / 2.0D0
                    END IF
                  END DO
                END DO
              END IF

              CALL Shap3D( xloca, yloca, zloca, coren, shape, xjaco,           &
                           ndime, nnode, lnode, ielem, iflag, nerrc )
              IF( nerrc .GT. 0 ) RETURN
              DO idofn = 1, mdofn
                idofx = ( inodx - 1 ) * mdofn + idofn
                IF( ntypc .EQ. 1 .OR. ntypc .EQ. 3 ) THEN
                  DO imode = 1, nmode
                    smodx( idofx, imode ) = 0.0D0
                    DO jnode = 1, nnode
                      ipoin = lnode( jnode )
                      IF( ipoin .GT. 0 ) THEN
                        idofs = ( ipoin - 1 ) * mdofn + idofn
                        smodx( idofx, imode ) = smodx( idofx, imode ) +        &
                        shape(     4, jnode ) * shmod( idofs, imode )
                      END IF
                    END DO
                  END DO
                END IF
                IF( ntypc .NE. 1 ) THEN
                  dispx( idofx, kload ) = 0.0D0
                  IF( idofn .LE. ndofn ) THEN

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -