flagshyp.f90

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

        eload (ig) = 0.d0
      END DO
      DO ig = 1, negdf
        pdisp (ig) = 0.d0
      END DO
!
!     Reads point loads and allocates them accordingly
!
      DO il = 1, nplds
        READ (1, *, err = 10) ip, (force (id), id = 1, ndime)
        IF ( ip > npoin ) goto 10
        DO id = 1, ndime
          ig = ldgof (id, ip)
          IF ( ig > 0 ) eload (ig) = force (id)
        END DO
      END DO
!
!     Reads in prescribed displacements. This must applied to only
!     fixed degrees of freedom
!
      DO i = 1, nprs
        READ (1, *, err = 10) ip, id, presc
        IF ( ip > npoin ) goto 10
        il = ldgof (id, ip)
        IF ( il > 0 ) goto 10
        pdisp ( - il) = presc
      END DO
!
!     Reads in pressure elements
!
      IF ( nbpel == 0 ) return
      IF ( nbpel > mbpel ) goto 10
      DO ie = 1, nbpel
        READ (1, *, err = 10) je, (lbnod (i, je), i = 1, nnodb), &
                              press (je)
      END DO
      RETURN
!
   10 WRITE (6, '(a)') ' error reading nodal forces'
      STOP 'inloads: Error reading nodal forces'
END SUBROUTINE inloads
 
!-----------------------------------------------------------------------
SUBROUTINE incontr (nincr, xlmax, dlamb, miter, cnorm, searc,    &
                    arcln, lun)
!-----------------------------------------------------------------------
!
!     Reads in parameters to control the solution process
!
!     nincr --> number of load increments
!     xlmax --> maximum load parameter
!     dlamb --> load parameter increment
!     miter --> maximum number of iterations per increment
!     cnorm --> convergence norm
!     searc --> line search parameter (=0.0 means no line search)
!     arcln --> arc-length parameter (=0.0 means no arc length)
!     lun   --> logical unit from which the above are read
!
!     Note that arc length and line search methods cannot be
!          used together
!
!-----------------------------------------------------------------------
!
  IMPLICIT DOUBLE PRECISION (a - h, o - z)
!
      IF ( lun == 5 ) write (6, 100)
        100 FORMAT(                                                  &
        &  ' Enter: No. of increments; max. load parameter;',/,      &
        &  ' load parameter increment; max. No. of iterations;',/,   &
        &  ' convergence norm; line search and arc length values: ')
      READ (lun, *, err = 10) nincr, xlmax, dlamb, miter, cnorm, &
                              searc, arcln
      IF ( searc * arcln /= 0.d0 ) goto 10
      IF ( searc == 0.d0 ) searc = 1.d5
      RETURN
!
   10 WRITE (6, '(a)') ' error reading solution control parameters'
      STOP 'incontr: Error reading solution control parameters'
END SUBROUTINE incontr
 
!-----------------------------------------------------------------------
SUBROUTINE initno (ndime, npoin, ndgof, nprof, x, x0, stifd,     &
                   stifp, resid)
!-----------------------------------------------------------------------
!
!     Initializes nodal data such as initial coordinates, stiffness
!     matrix and residual forces
!
!     ndime  -->  number of dimensions
!     npoin  -->  number of points
!     ndgof  -->  number of degrees of freedom
!     nprof  -->  number of profile entries
!     x      -->  nodal coordinates
!     x0     -->  initial nodal coordinates
!     stifd  -->  diagonal of K
!     stifp  -->  profile part of K
!     resid  -->  residual forces
!
!-----------------------------------------------------------------------
!
  IMPLICIT DOUBLE PRECISION (a - h, o - z)
  DIMENSION x (ndime, * ), x0 (ndime, * ), stifd ( * ), stifp ( * ),&
      resid ( * )
!
!     initialises the stiffness matrix to zero
!
      DO ig = 1, ndgof
        stifd (ig) = 0.d0
        resid (ig) = 0.d0
      END DO
      DO ip = 1, nprof
        stifp (ip) = 0.d0
      END DO
!
!     copies the nodal coordinates to the initial nodal coordinates
!
      DO ip = 1, npoin
        DO id = 1, ndime
          x0 (id, ip) = x (id, ip)
        END DO
      END DO
END SUBROUTINE initno
 
!-----------------------------------------------------------------------
SUBROUTINE initel (ndime, nnode, ngaus, nelem, gravt, x, eledb,  &
                   lnods, matno, matyp, props, ldgof, eload, kprof, &
                   stifd, stifp, vinc, elecd, elacd, vol0)
!-----------------------------------------------------------------------
!
!     Initializes element data such as initial stiffness matrix,
!     element loads,...
!
!     ndime  -->  number of dimensions
!     nnode  -->  number of nodes
!     ngaus  -->  number of Gauss points
!     nelem  -->  gravity vector
!     gravt  -->  gravity vector
!     x      -->  nodal coordinates
!     eledb  -->  element data base
!     lnods  -->  element connectivities
!     matno  -->  materal number
!     matyp  -->  material type
!     props  -->  material properties
!     ldgof  -->  degree of freedom number
!     eload  -->  nodal point and body forces
!     kprof  -->  profile addresses
!     stifd  -->  diagonal of K
!     stifp  -->  profile part of K
!     vinc   -->  jacobian at each gauss point
!     elecd  -->  element Cartesian derivatives at each Gauss point
!     elacd  -->  element average Cartesian derivatives
!     vol0   -->  initial element volume
!
!-----------------------------------------------------------------------
!
  IMPLICIT DOUBLE PRECISION (a - h, o - z)
  DIMENSION gravt (ndime), lnods (nnode, * ), props (8, *),         &
      ldgof (ndime, *), eload (*), kprof (*), stifd (*), stifp (*), &
      x (ndime, *), vol0 (*), eledb (ndime+1, nnode+1, ngaus),      &
      vinc (*), elacd (ndime, *), ctens (3, 3, 3, 3), matno (*),    &
      matyp (*), elecd (ndime, nnode, *)
!
!     Loops over all the elements. First finds the Cartesian
!     derivatives of the shape function for each Gauss point
!
      totv = 0.d0
      DO ie = 1, nelem
        vol0 (ie) = 0.d0
        im = matno (ie)
        mat = matyp (im)
        CALL gradsh (ie, ndime, nnode, ngaus, eledb, lnods, x, &
                     elecd, vinc)
!
!     For incompressible types of elements finds the volumetric
!     component of the tangent matrix
!
        IF (  (mat == 5) .or. (mat == 7)  ) then
          xkapp = props (4, im)
          CALL kvolume (ndime, nnode, ngaus, xkapp, vinc, elecd, &
                        elacd, lnods (1, ie), ldgof, kprof, stifd, &
                        stifp)
        END IF
!
!     Loops over the element Guass points
!
        DO ig = 1, ngaus
!
!     For plane stress materials takes the initial thickness
!     into account
!
          imm = mat / 2
          IF ( mat == 2 * imm ) vinc (ig) = vinc (ig) * props (4, im)
!
!     Finds forces due to gravity and adds them to the global
!     vector eload
!
          DO in = 1, nnode
            shape = eledb (1, in, ig) * props (1, im) * vinc (ig)
            ip = lnods (in, ie)
            DO id = 1, ndime
              if = ldgof (id, ip)
              IF ( if > 0 ) then
                eload (if) = eload (if) + gravt (id) * shape
              END IF
            END DO
          END DO
!
!     Evaluates the initial stiffness matrix and adds the volume
!     contribution
!
          IF ( mat <= 8 ) CALL cisotp (ndime, props (3, im), &
                                       props (2, im), ctens)
          CALL kconst (ndime, nnode, lnods (1, ie), ldgof, ctens, &
                       elecd (1, 1, ig), vinc (ig), kprof, stifd, stifp)
          totv = totv + vinc (ig)
          vol0 (ie) = vol0 (ie) + vinc (ig)
        END DO
      END DO
      WRITE (6, 101) totv
      101 FORMAT(' Total mesh volume is: ',g15.5)
!
END SUBROUTINE initel
 
!-----------------------------------------------------------------------
SUBROUTINE initrk (ndgof, nprof, negdf, xlamb, eload, tload,     &
                   resid, react, stifd, stifp)
!-----------------------------------------------------------------------
!
!     Initialises the residual vector and the stiffness matrices
!
!     ndgof  --> number of degrees of freedom
!     nprof  --> number of entries in the profile
!     negdf  --> negative degrees of freedom
!     xlamb  --> current load parameter value
!     resid  --> residual force
!     react  --> reactions
!     eload  --> nominal nodal loads
!     tload  --> total loads including pressure
!     stifd  --> diagonal component of the stiffness matrix
!     stifp  --> upper profile component of stiffness
!
!-----------------------------------------------------------------------
!
  IMPLICIT DOUBLE PRECISION (a - h, o - z)
  DIMENSION resid ( * ), eload ( * ), stifd ( * ), stifp ( * ),     &
      react ( * ), tload ( * )
!
!     Initialises the diagonal stiffness vector, the residual vector
!     and the total loads
!
      DO id = 1, ndgof
        stifd (id) = 0.d0
        resid (id) = xlamb * eload (id)
        tload (id) = eload (id)
      END DO
!
!     Initialises the reactions
!
      DO id = 1, negdf
        react (id) = 0.d0
      END DO
!
!     Initialises the profile of K
!
      DO id = 1, nprof
        stifp (id) = 0.d0
      END DO
END SUBROUTINE initrk
 
!-----------------------------------------------------------------------
SUBROUTINE bpress (ndime, nnodb, ngaub, nbpel, xlamb, elbdb,     &
                   lbnod, press, x, ldgof, tload, react, resid,  &
                   kprof, stifd, stifp)
!-----------------------------------------------------------------------
!
!     Evaluates the nodal forces due to boundary normal pressure and
!     assembles the corresponding component of the tangent matrix
!
!     ndime --> number of dimensions
!     nnodb --> number of nodes per boundary element
!     ngaub --> number of Gauss points per boundary element
!     nbpel --> number of boundary elements with applied pressure
!     xlamb --> force factor lambda
!     elbdb --> bondary elements information
!     lbnod --> boundary element connectivities
!     press --> applied pressure
!     x     --> nodal coordinates
!     ldgof --> degrees of freedom number of each node
!     tload --> total nominal load
!     react --> reactions
!     resid --> residual forces
!     kprof --> profile addresses
!     stifd --> diagonal stiffness
!     stifp --> profile component of tangent matrix
!
!-----------------------------------------------------------------------
!
  IMPLICIT DOUBLE PRECISION (a - h, o - z)
  DIMENSION elbdb (ndime, nnodb + 1, * ), lbnod (nnodb, * ),        &
      press ( * ), x (ndime, * ), ldgof (ndime, * ), tload ( * ),   &
      react ( * ), resid ( * ), kprof ( * ), stifd ( * ),           &
      stifp ( * ), dxis (3, 2)
  DATA ( (dxis (i, j), i = 1, 3), j = 1, 2) &
         / 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0 /
!
!     If no pressure elements are present returns
!
      IF ( nbpel == 0 ) return
!
!     otherwise loops over the elements and Gauss points and finds the
!     value of the applied pressure
!
      DO ie = 1, nbpel
        epres = press (ie)
!b      if ( epress == 0.d0 ) cycle ! to next element !b check this
        DO ig = 1, ngaub
!
!     Initialises the tangential vectors dx/dxi and dx/deta:. For 2-D
!     sets dx/deta=(0,0,-1)
!
        IF ( ndime == 2 ) dxis (3, 2) = - 1.d0
!
!     Evaluates the derivatives using equation (7.6b)
!
        DO id = 1, ndime
          DO jd = 1, ndime-1
            dxis (id, jd) = 0.d0
            DO in = 1, nnodb
              ip = lbnod (in, ie)
              dxis (id, jd) = dxis (id, jd) &
                            + x (id, ip) * elbdb (jd+1, in, ig)
            END DO
          END DO
        END DO
!
!     Finds and assembles the nodal forces
!
        CALL pforce (ndime, nnodb, epres, xlamb, elbdb (1, 1, ig),  &
                     lbnod (1, ie), ldgof, tload, resid, react, dxis)
        apres = epres * xlamb
!
!     Finds and assembles the tangent matrix component
!
        CALL kpress (ndime, nnodb, apres, elbdb (1, 1, ig),         &
                     lbnod (1, ie), ldgof, dxis, kprof, stifd, stifp)
        END DO
      END DO
END SUBROUTINE bpress
 
!-----------------------