flagshyp.f90

来自「大应变超弹性平面分析程序」· F90 代码 · 共 1,766 行 · 第 1/5 页

      dsh0x (xi, et, ze) = 0
!
      dsh1e (xi, et, ze) = 4 * (xi + et + ze) - 3
      dsh2e (xi, et, ze) = 0
      dsh3e (xi, et, ze) = 4 * et - 1
      dsh4e (xi, et, ze) = 0
      dsh5e (xi, et, ze) = - 4 * xi
      dsh6e (xi, et, ze) = 4 * xi
      dsh7e (xi, et, ze) = 4 * (1 - xi - 2 * et - ze)
      dsh8e (xi, et, ze) = - 4 * ze
      dsh9e (xi, et, ze) = 0
      dsh0e (xi, et, ze) = 4 * ze
!
      dsh1z (xi, et, ze) = 4 * (xi + et + ze) - 3
      dsh2z (xi, et, ze) = 0
      dsh3z (xi, et, ze) = 0
      dsh4z (xi, et, ze) = 4 * ze - 1
      dsh5z (xi, et, ze) = - 4 * xi
      dsh6z (xi, et, ze) = 0
      dsh7z (xi, et, ze) = - 4 * et
      dsh8z (xi, et, ze) = 4 * (1 - xi - et - 2 * ze)
      dsh9z (xi, et, ze) = 4 * xi
      dsh0z (xi, et, ze) = 4 * et
!
!     Constructs the data base. First Gauss point position and weight
!
      DO ig = 1, 5
        IF ( ig == 1 ) then
          eledb (1, 11, ig) = - 0.1333333333333333d0
        ELSE
          eledb (1, 11, ig) = 0.075d0
        END IF
        xi = gauss (1, ig)
        et = gauss (2, ig)
        ze = gauss (3, ig)
        eledb (2, 11, ig) = xi
        eledb (3, 11, ig) = et
        eledb (4, 11, ig) = ze
!
!     Shape functions and derivatives
!
        eledb (1, 1, ig) = shap1 (xi, et, ze)
        eledb (2, 1, ig) = dsh1x (xi, et, ze)
        eledb (3, 1, ig) = dsh1e (xi, et, ze)
        eledb (4, 1, ig) = dsh1z (xi, et, ze)
        eledb (1, 2, ig) = shap2 (xi, et, ze)
        eledb (2, 2, ig) = dsh2x (xi, et, ze)
        eledb (3, 2, ig) = dsh2e (xi, et, ze)
        eledb (4, 2, ig) = dsh2z (xi, et, ze)
        eledb (1, 3, ig) = shap3 (xi, et, ze)
        eledb (2, 3, ig) = dsh3x (xi, et, ze)
        eledb (3, 3, ig) = dsh3e (xi, et, ze)
        eledb (4, 3, ig) = dsh3z (xi, et, ze)
        eledb (1, 4, ig) = shap4 (xi, et, ze)
        eledb (2, 4, ig) = dsh4x (xi, et, ze)
        eledb (3, 4, ig) = dsh4e (xi, et, ze)
        eledb (4, 4, ig) = dsh4z (xi, et, ze)
        eledb (1, 5, ig) = shap5 (xi, et, ze)
        eledb (2, 5, ig) = dsh5x (xi, et, ze)
        eledb (3, 5, ig) = dsh5e (xi, et, ze)
        eledb (4, 5, ig) = dsh5z (xi, et, ze)
        eledb (1, 6, ig) = shap6 (xi, et, ze)
        eledb (2, 6, ig) = dsh6x (xi, et, ze)
        eledb (3, 6, ig) = dsh6e (xi, et, ze)
        eledb (4, 6, ig) = dsh6z (xi, et, ze)
        eledb (1, 7, ig) = shap7 (xi, et, ze)
        eledb (2, 7, ig) = dsh7x (xi, et, ze)
        eledb (3, 7, ig) = dsh7e (xi, et, ze)
        eledb (4, 7, ig) = dsh7z (xi, et, ze)
        eledb (1, 8, ig) = shap8 (xi, et, ze)
        eledb (2, 8, ig) = dsh8x (xi, et, ze)
        eledb (3, 8, ig) = dsh8e (xi, et, ze)
        eledb (4, 8, ig) = dsh8z (xi, et, ze)
        eledb (1, 9, ig) = shap9 (xi, et, ze)
        eledb (2, 9, ig) = dsh9x (xi, et, ze)
        eledb (3, 9, ig) = dsh9e (xi, et, ze)
        eledb (4, 9, ig) = dsh9z (xi, et, ze)
        eledb (1, 10, ig) = shap0 (xi, et, ze)
        eledb (2, 10, ig) = dsh0x (xi, et, ze)
        eledb (3, 10, ig) = dsh0e (xi, et, ze)
        eledb (4, 10, ig) = dsh0z (xi, et, ze)
      END DO
END SUBROUTINE tetr10db
 
!-----------------------------------------------------------------------
SUBROUTINE quad4db (eledb)
!-----------------------------------------------------------------------
!
!     Four noded quadrilateral definintion
!
!-----------------------------------------------------------------------
!
 IMPLICIT DOUBLE PRECISION (a - h, o - z)
 PARAMETER (gauss = 0.577350269189626d0)
 DIMENSION eledb (3, 5, 4), node (2, 4)
!
!     Positions of the nodes (and Gauss points if multiplied by 0.577...
!     in the parametric plane
!
 DATA (node (i, 1), i = 1, 2) / - 1, - 1 /
 DATA (node (i, 2), i = 1, 2) /   1, - 1 /
 DATA (node (i, 3), i = 1, 2) /   1,   1 /
 DATA (node (i, 4), i = 1, 2) / - 1,   1 /
!
!     shape functions and derivatives  ! obsolete under f95
!
      shap (xi, et, in) = (1 + node (1, in) * xi) &
                        * (1 + node (2, in) * et) / 4.d0
      dshx (xi, et, in) = node (1, in) * (1 + node (2, in) * et) / 4.d0
      dshe (xi, et, in) = node (2, in) * (1 + node (1, in) * xi) / 4.d0
!
!     constructs the data base. first positions of Gauss points
!     and weights
!
      DO ig = 1, 4
        xi = float (node (1, ig) ) * gauss
        et = float (node (2, ig) ) * gauss
        eledb (1, 5, ig) = 1.d0
        eledb (2, 5, ig) = xi
        eledb (3, 5, ig) = et
!
!     Shape functions and derivatives
!
        DO in = 1, 4
          eledb (1, in, ig) = shap (xi, et, in)
          eledb (2, in, ig) = dshx (xi, et, in)
          eledb (3, in, ig) = dshe (xi, et, in)
        END DO
      END DO
END SUBROUTINE quad4db
 
!-----------------------------------------------------------------------
SUBROUTINE hexa8db (eledb)
!-----------------------------------------------------------------------
!
!     Eight noded hexahedron definintion
!
!-----------------------------------------------------------------------
!
 IMPLICIT DOUBLE PRECISION (a - h, o - z)
 PARAMETER (gauss = 0.577350269189626d0)
 DIMENSION eledb (4, 9, 8), node (3, 8)
!
!     Positions of the nodes (and Gauss points if multiplied by 0.577...
!     in the parametric plane
!
 DATA (node (i, 1), i = 1, 3) / - 1, - 1, - 1 /
 DATA (node (i, 2), i = 1, 3) /   1, - 1, - 1 /
 DATA (node (i, 3), i = 1, 3) /   1,   1, - 1 /
 DATA (node (i, 4), i = 1, 3) / - 1,   1, - 1 /
 DATA (node (i, 5), i = 1, 3) / - 1, - 1,   1 /
 DATA (node (i, 6), i = 1, 3) /   1, - 1,   1 /
 DATA (node (i, 7), i = 1, 3) /   1,   1,   1 /
 DATA (node (i, 8), i = 1, 3) / - 1,   1,   1 /
!
!     shape functions and derivatives  ! obsolete under f95
!
      shap (xi, et, ze, in) = (1 + node (1, in) * xi) &
                            * (1 + node (2, in) * et) &
                            * (1 + node (3, in) * ze) / 8.d0
      dshx (xi, et, ze, in) = node (1, in) * (1 + node (2, in) * et)  &
                            * (1 + node (3, in) * ze) / 8.d0
      dshe (xi, et, ze, in) = node (2, in) * (1 + node (1, in) * xi)  &
                            * (1 + node (3, in) * ze) / 8.d0
      dshz (xi, et, ze, in) = node (3, in) * (1 + node (2, in) * et)  &
                            * (1 + node (1, in) * xi) / 8.d0
!
!     constructs the data base. First positions of Gauss points
!     and weights
!
      DO ig = 1, 8
        eledb (1, 9, ig) = 1.d0
        xi = float (node (1, ig) ) * gauss
        et = float (node (2, ig) ) * gauss
        ze = float (node (3, ig) ) * gauss
        eledb (2, 9, ig) = xi
        eledb (3, 9, ig) = et
        eledb (4, 9, ig) = ze
 
!
!     Shape functions and derivatives
!
        DO in = 1, 8
          eledb (1, in, ig) = shap (xi, et, ze, in)
          eledb (2, in, ig) = dshx (xi, et, ze, in)
          eledb (3, in, ig) = dshe (xi, et, ze, in)
          eledb (4, in, ig) = dshz (xi, et, ze, in)
        END DO
      END DO
END SUBROUTINE hexa8db
 
!-----------------------------------------------------------------------
SUBROUTINE innodes (mpoin, ndime, npoin, x, icode, ldgof)
!-----------------------------------------------------------------------
!
!     Reads initial nodal coordinates
!
!     mpoin    -->  maximum number of nodes
!     ndime    -->  number of dimensions
!     npoin    -->  number of points
!     x        -->  nodal coordinates
!     icode    -->  boundary codes:
!
!                    0: free
!                    1: x fixed
!                    2: y fixed
!                    3: x+y fixed
!                    4: z fixed
!                    5: x+z fixed
!                    6: y+z fixed
!                    7: x+y+z fixed
!
!     ldgof    -->  boundary code (1 or 0) for each degree of freedom
!                   corresponding to each node. Later degree of freedom
!                   numbers for each node.
!
!-----------------------------------------------------------------------
!
 IMPLICIT DOUBLE PRECISION (a - h, o - z)
 !b DIMENSION x (ndime, * ), icode ( * ), ldgof (ndime, * )
 INTEGER,          INTENT (IN)  :: mpoin, ndime
 INTEGER,          INTENT (OUT) :: npoin
 INTEGER,          INTENT (OUT) :: icode (mpoin), ldgof (ndime, mpoin)
 DOUBLE PRECISION, INTENT (OUT) :: x (ndime, mpoin)
!
!     First reads in the total number of nodes
!
      READ (1, *, err = 10) npoin
      IF ( npoin > mpoin ) then
        WRITE (6, 100) npoin
        100 FORMAT('innodes: problem dimensions exceed maximum, ', &
        &          'set mpoin to:',i10)
        STOP 'innodes: dimensions exceed mpoin'
      END IF
!
!     Reads nodal coordinates and boundary codes
!
      DO jp = 1, npoin
        READ (1, *, err = 10) ip, icode (ip), &
                              (x (id, ip), id = 1, ndime)
      END DO
!
!     Initialises the degree of freedom array ldgof
!
      DO ip = 1, npoin
        ic = icode (ip)
        DO id = 1, ndime
          ic0 = ic
          ic = ic / 2
          ldgof (id, ip) = ic0 - ic * 2
        END DO
      END DO
      RETURN
!
      10 WRITE (6, '(a)') ' error reading nodal coordinates'
      STOP 'innodes: dimensions exceed mpoin'
!
END SUBROUTINE innodes
 
!-----------------------------------------------------------------------
SUBROUTINE inelems (melem, nelem, nnode, lnods, matno, nmats)
!-----------------------------------------------------------------------
!
!     Reads in the element connectivities
!
!     melem   --> maximum number of elements
!     nelem   --> number of elements
!     nnode   --> number of nodes per element
!     lnods   --> nodal conectivities of dimensions (nnode,nelem)
!     matno   --> material number of each element
!     nmats   --> number of materials
!
!-----------------------------------------------------------------------
!
 IMPLICIT DOUBLE PRECISION (a - h, o - z)
 DIMENSION lnods (nnode, * ), matno ( * )
!
!     First reads the number of elements
!
      READ (1, *, err = 10) nelem
      IF ( nelem > melem ) then
        WRITE (6, 100) nelem
        100 FORMAT('inelems: problem dimensions exceed maximum, ', &
        &          'set melem to:',i10)
        STOP 'inelems: dimensions exceed melem'
      END IF
      nmats = 0
!
!     Reads in element connectivities and material types
!
      DO je = 1, nelem
        READ (1, *, err = 10) ie, matno (ie), &
                              (lnods (in, ie), in = 1, nnode)
        nmats = max (nmats, matno (ie) )
      END DO
      RETURN
!
      10 WRITE (1, '(a)') ' Error reading element connectivities'
      STOP 'inelems: dimensions exceed melem'
END SUBROUTINE inelems
 
!-----------------------------------------------------------------------
SUBROUTINE nodecon (npoin, nelem, nnode, lnods, nconn)
!-----------------------------------------------------------------------
!
!     Determines the node to node connectivities and stores them
!     in nconn as a linked list
!
!     npoin  -->  number of mesh nodes
!     nelem  -->  number of elements
!     nnode  -->  number of nodes per element
!     lnods  -->  element nodal connectivities
!     nconn  -->  nodal connectivities as a linked list of dimensions
!                 (2,total number of node to node connexions):
!
!                 | number of nodes connected to node 1,.... |
!                 | position of next node conected to 1,.... |
!
!-----------------------------------------------------------------------
!
 IMPLICIT DOUBLE PRECISION (a - h, o - z)
 DIMENSION lnods (nnode, * ), nconn (2, * )
!
!     Initialises the matrix nconn to zero and sets the first
!     free position in nconn to npoin+1
!
      DO ip = 1, npoin
        nconn (1, ip) = 0
        nconn (2, ip) = 0
      END DO
      nfree = npoin + 1
!
!     Loop over the elements to examine all nodal connexions
!
      DO ie = 1, nelem
        DO in = 1, nnode
          ip = lnods (in, ie)
          DO jn = 1, nnode
            jp = lnods (jn, ie)
            IF ( jp /= ip ) then
!
!     Node jp is connected to node ip. Jumps along the linked list
!     until no more connections are found
!
              naddr = nconn (2, ip)
              nadd0 = ip
              DO WHILE (naddr > 0)
                kp = nconn (1, naddr)
                IF ( kp == jp) goto 10
                nadd0 = naddr
                naddr = nconn (2, nadd0)
              END DO
!
!     Puts the new node jp as a connection to ip and creaes a new cell