flagshyp.f90

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

!     in the link list at position nfree
!
              nconn (2, nadd0) = nfree
              nconn (1, nfree) = jp
              nconn (2, nfree) = 0
              nconn (1, ip) = nconn (1, ip) + 1
              nfree = nfree+1
   10       END IF
          END DO
        END DO
      END DO
END SUBROUTINE nodecon
 
!-----------------------------------------------------------------------
SUBROUTINE degfrm (mdgof, npoin, ndime, nconn, ndgof, negdf,     &
                   ldgof, ndque)
!-----------------------------------------------------------------------
!
!     Determines the profile solver information, in particular degrees
!     of freedom corresponding to each node.
!     The degree of freedom numbering is performed in such as way
!     that the profile size is reasonable.
!
!     mdgof -->  maximum number of degrees of freedom allowed
!     npoin -->  number of points
!     ndime -->  number of dimensions
!     nconn -->  node to node connectivity array
!     ndgof -->  number of degrees of freedom
!     negdf -->  number of fixed degrees of freedom
!     ldgof -->  degrees of freedom of each node, array of dimensions
!                (ndime,npoin)
!     ndque -->  circular queue, storing next nodes' degrees of freedom
!                to number
!
!-----------------------------------------------------------------------
!
  IMPLICIT DOUBLE PRECISION (a - h, o - z)
  DIMENSION nconn (2, * ), ldgof (ndime, * ), ndque ( * )
!
!     Initialize the queue by puting node 1 in it.
!     numbn -->  number of nodes already dealt with
!
      ndque (1) = 1
      nconn (1, 1) = 0
      lqbeg = 1
      lqend = 1
      numbn = 0
      ndgof = 0
      negdf = 0
!
!     Loops while there are points to be numbered. If the queue is empty
!     finds the first node unnumbered.
!
      DO WHILE (numbn < npoin)
      IF ( lqend < lqbeg ) then
        ip = 2
        DO WHILE (nconn (1, ip) <= 0)
          ip = ip + 1
        END DO
!
!     Otherwise takes the first node in the queue and removes it from
!     the queue
!
      ELSE
        lqb = 1 + mod (lqbeg - 1, npoin)
        ip = ndque (lqb)
        lqbeg = lqbeg + 1
      END IF
      numbn = numbn + 1
!
!     Adds nodes connected to ip to the queue
!
      naddr = nconn (2, ip)
      DO WHILE (naddr > 0)
        jp = nconn (1, naddr)
        naddr = nconn (2, naddr)
        IF ( nconn (1, jp) > 0 ) then
          lqend = lqend+1
          ndque (lqend) = jp
          nconn (1, jp) = 0
        END IF
      END DO
!
!     Defines the degree of freedom numbers for node ip. If fixed sets
!     the corresponding position to 0.
!
        DO id = 1, ndime
          IF ( ldgof (id, ip) == 0 ) then
            ndgof = ndgof + 1
            ldgof (id, ip) = ndgof
          ELSE
            negdf = negdf + 1
            ldgof (id, ip) = - negdf
          END IF
        END DO
      END DO
!
!     Checks that the total number of degrees of freedom is OK
!
      IF ( ndgof > mdgof ) then
        WRITE (6, 100) ndgof
        100 FORMAT('degfrm: problem dimensions exceed maximum, ', &
        &          'set mdgof to:',i10)
        STOP 'degfrm: dimensions exceed mdgof'
      END IF
!
END SUBROUTINE degfrm
 
!-----------------------------------------------------------------------
SUBROUTINE profile (mprof, ndgof, nelem, nnode, ndime, lnods,    &
                    ldgof, nprof, kprof)
!-----------------------------------------------------------------------
!
!     Determines the profile solver information, in particular profile
!     heights and then profile addresses.
!
!     mprof -->  maximum number of degrees of freedom allowed
!     nelem -->  number of elements
!     ndime -->  number of dimensions
!     nnode -->  number of nodes per element
!     ndgof -->  number of degrees of freedom
!     lnods -->  element connectivities
!     ldgof -->  degrees of freedom of each node, array of dimensions
!                (ndime,npoin)
!     nprof -->  number of entries in the out of diagonal part of K
!     kprof -->  first height of each profile column, then address in
!                stifp of each column.
!
!-----------------------------------------------------------------------
!
  IMPLICIT DOUBLE PRECISION (a - h, o - z)
  DIMENSION ldgof (ndime, * ), kprof ( * ), lnods (nnode, * )
!
!     Initialise  kprof to 0
!
      DO ig = 1, ndgof
        kprof (ig) = 0
      END DO
!
!     Loops over each element and first finds the lowest degree of
!     freedom in that element (lg)
!
      DO ie = 1, nelem
        lg = 5**5
        DO in = 1, nnode
          ip = lnods (in, ie)
          DO id = 1, ndime
            ig = ldgof (id, ip)
            IF ( ig > 0 ) lg = min (lg, ig)
          END DO
        END DO
!
!     Then finds the difference between other degrees of freedom in the
!     element and lg.
!
        DO in = 1, nnode
          ip = lnods (in, ie)
          DO id = 1, ndime
            ig = ldgof (id, ip)
            IF ( ig > 0 ) kprof (ig) = max (kprof (ig), ig - lg)
          END DO
        END DO
      END DO
!
!     Transforms the height information into addresses of each column in
!     stifp, by simply accumulating the heights of each column. Also
!     finds the total number of spaces need to store stifp.
!
      kprof (1) = 0
      DO ig = 2, ndgof
        kprof (ig) = kprof (ig) + kprof (ig - 1)
      END DO
      nprof = kprof (ndgof)
!
!     Checks that there is enough storage alocatted to stifp.
!
      IF ( nprof > mprof ) then
        WRITE (6, 100) nprof
        100 FORMAT('profile: problem dimensions exceed maximum, ', &
        &          'set mprof to:',i10)
        STOP 'profile: dimensions exceed mprof'
      END IF
      RETURN
!
END SUBROUTINE profile
 
!-----------------------------------------------------------------------
SUBROUTINE matprop (ndime, nmats, matyp, props)
!-----------------------------------------------------------------------
!
!     Reads in the material type and material properties. Also
!     re-arranges the properties to save time later on
!
!     ndime --> number of dimensions
!     matyp --> material type
!     props --> vector of properties. The first is always the initial
!               density, the rest depend on the material type as:
!
!    |----------|--------|--------|--------|--------|--------|--------|
!    | material |props(1)|props(2)|props(3)|props(4)|props(5)|props(6)|
!    |==========|========|========|========|========|========|========|
!    |     1    |  rho   |   mu   | lambda |   -    |   -    |   -    |
!    |----------|--------|--------|--------|--------|--------|--------|
!    |     2    |  rho   |   -    |   -    |   -    |   -    |   -    |
!    |----------|--------|--------|--------|--------|--------|--------|
!    |     3    |  rho   |   mu   | lambda |   -    |   -    |   -    |
!    |----------|--------|--------|--------|--------|--------|--------|
!    |     4    |  rho   |   mu   | lambda | thick  | (gam)  |   -    |
!    |----------|--------|--------|--------|--------|--------|--------|
!    |     5    |  rho   |   mu   | -2mu/3 | kappa  |   -    |   -    |
!    |----------|--------|--------|--------|--------|--------|--------|
!    |     6    |  rho   |   mu   |   2mu  | thick  |   -    |   -    |
!    |----------|--------|--------|--------|--------|--------|--------|
!    |     7    |  rho   |   mu   | -2mu/3 | kappa  |   -    |   -    |
!    |----------|--------|--------|--------|--------|--------|--------|
!    |     8    |  rho   |   mu   |   2mu  | thick  |   -    |   -    |
!    |----------|--------|--------|--------|--------|--------|--------|
!
!     Material types:
!
!        1 --> plane strain (or 3d) compressible neo-Hookean
!        2 --> plane stress compressible neo-Hookean (not implemented)
!        3 --> plane strain (or 3d) hyperelastic in principal directions
!        4 --> plane stress hyperelastic in principal directions
!        5 --> plane strain (or 3d) nearly incompressible neo-Hookean
!        6 --> plane stress incompressible neo-Hookean
!        7 --> plane strain (or 3d) nearly incompressible in principal
!              directions
!        8 --> plane stress incompressible in principal directions
!
!-----------------------------------------------------------------------
!
  IMPLICIT DOUBLE PRECISION (a - h, o - z)
  DIMENSION props (8, * ), matyp ( * )
!
!     First makes sure that the number of materials coincides with the
!     number obtained from element information
!
      READ (1, *, err = 10) nmat2
      IF ( nmat2 /= nmats ) goto 10
!
!     Loops over all the material types
!
      DO jm = 1, nmats
        READ (1, *, err = 10) im, matyp (im)
!
!     Check that the material is implemented and re-arranges parameters
!
        IF ( matyp (im) == 1 ) then
          READ (1, *, err = 10) props (1, im), props (2, im), &
                                props (3, im)
        ELSEIF ( matyp (im) == 3 ) then
          READ (1, *, err = 10) props (1, im), props (2, im), &
                                props (3, im)
!
!     For plane stress material obtains the effective Lambda coefficient
!     and the parameter gamma (see equation 5.68)
!
        ELSEIF ( matyp (im) == 4 ) then
          READ (1, *, err = 10) (props (i, im), i = 1, 4)
          xmu = props (2, im)
          xla = props (3, im)
          gam = 2 * xmu / (xla + 2 * xmu)
          props (5, im) = gam
          props (3, im) = gam * xla
          IF ( props (4, im) == 0.d0 ) props (4, im) = 1.d0
!
!     Incompressible materials. sets lambda equal to minus two thirds of
!     mu for the evaluation of the initial stiffness matrix.
!
        ELSEIF ( matyp (im) == 5 ) then
          READ (1, *, err = 10) props (1, im), props (2, im), &
                                props (4, im)
          props (3, im) = - 2 * props (2, im) / 3.d0
        ELSEIF ( matyp (im) == 7 ) then
          READ (1, *, err = 10) props (1, im), props (2, im), &
                                props (4, im)
          props (3, im) = - 2 * props (2, im) / 3.d0
!
!     Plane stress incompressibility is enforced exactly.
!
        ELSEIF ( matyp (im) == 6 ) then
          READ (1, *, err = 10) props (1, im), props (2, im), &
                                props (4, im)
          props (3, im) = 2 * props (2, im)
          IF ( props (4, im) == 0.d0 ) props (4, im) = 1.d0
        ELSEIF ( matyp (im) == 8 ) then
          READ (1, *, err = 10) props (1, im), props (2, im), &
                                props (4, im)
          props (3, im) = 2 * props (2, im)
          IF ( props (4, im) == 0.d0 ) props (4, im) = 1.d0
        ELSE
          WRITE (6, 101) matyp (im)
          101 FORMAT(' Material type not implemented: ',i3)
          STOP 'matprop: Material type not implemented'
        END IF
!
!     Checks that plane stress materials are used with ndime=2
!
        imt = matyp (im) / 2
        imt = 2 * imt
        IF (  (imt == matyp (im) ) .and. (ndime /= 2)  ) then
          WRITE (6, 102)
          102 FORMAT(' Wrong material for number of dimensions')
          STOP 'matprop: Wrong material for number of dimensions'
        END IF
      END DO
      RETURN
!
   10 WRITE (6, '(a)') ' error reading material properties'
      STOP 'matprop: Error reading material properties'
!
END SUBROUTINE matprop
 
!-----------------------------------------------------------------------
SUBROUTINE inloads (ndime, npoin, ndgof, negdf, nnodb, mbpel,    &
                    ldgof, eload, pdisp, gravt, nprs, nbpel,     &
                    lbnod, press)
!-----------------------------------------------------------------------
!
!     Reads in loads and prescribed displacements
!
!     ndime --> number of dimensions
!     npoin --> number of points
!     ndgof --> number of degrees of freedom
!     negdf --> number of fixed degrees of freedom
!     nnodb --> number of nodes per boundary element
!     mbpel --> maximum number of boundary pressure elements
!     ldgof --> array containing the degrees of freedom of each node or
!               0 for fixed nodes. If the coordinate of these fixed
!               nodes is to be prescribed, the number stored will be
!               -(the prescribed condition number)
!     eload --> external load on each degree of freedom
!     pdisp --> prescribed positions of fixed nodes
!     gravt --> gravity acceleration vector
!     nprs  --> number of prescribed displacements
!     nbpel --> number of boundary pressure elements
!     lbnod --> pressure elements connectivity matrix
!     press --> applied pressures
!
!-----------------------------------------------------------------------
!
  IMPLICIT DOUBLE PRECISION (a - h, o - z)
  DIMENSION ldgof (ndime, * ), eload ( * ), pdisp ( * ), &
      gravt (ndime), force (3), press ( * ), lbnod (nnodb, * )
!
!     Reads the number of ponts loads, prescribed displacements and
!     gravity vector
!
      READ (1, *, err = 10) nplds, nprs, nbpel, gravt
!
!     Initialises the force vector and prescribed displacements
!
      DO ig = 1, ndgof