flagshyp.f90

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

!-----------------------------------------------------------------------
PROGRAM flagshyp ! (flagshyp.f is f90 syntax)
!-----------------------------------------------------------------------
!
! Finite element LArGe Strain HYperelasticity Program
!
! Written by: J.Bonet
!             Civil Engineering Department
!             University of Wales, Swansea
!
! This program has been written in order to demonstrate the
! concepts and theory explained in the book "Introduction to
! Nonlinear Continuum Mechanics for Finite Element Analysis", 1997, by 
! J. Bonet & R.D. Wood, Cambridge University Press, ISBN 0-521-57272-X.
!
! Copyright (c) J.Bonet & R.D. Wood
!
! Swansea, 1994,95
!-----------------------------------------------------------------------
!
! Translated from f77 syntax to f90 syntax by J. E. Akin,
! Rice University, June 2003, [e-mail]. 
! See "Object Oriented Programming via F90", by Ed Akin,
!     Cambridge University Press, 2002, ISBN 0-521-52408-3.
! 
! flagshyp.f   is f90 syntax, without dynamic memory management.  It has
!              been tested on f90 and f95 compilers.
! flagshyp.f90 will be f90 style, with dynamic memory management.
!
!-----------------------------------------------------------------------
!
! THIS PROGRAM IS LICENSED FREE OF CHARGE. THERE IS NO WARRANTY
! FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.
! EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS
! AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT
! WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING,
! BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
! AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
! TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.
! SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL
! NECESSARY SERVICING, REPAIR OR CORRECTION.
!
!-----------------------------------------------------------------------
!
! Version 1.1 6/9/1995
! ======================
!
! Materials implemented:
!
!    1 --> plane strain (or 3d) compressible neo-Hookean
!    2 --> plane stress compressible neo-Hookean (*)
!    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
!
!
! Elements implemented:
!
! Linear     3-noded triangle (tria3)
! Quadratic  6-noded triangle (tria6)
! Bi-linear  4-noded quadrilateral (quad4)
! Linear     4-noded tetrahedron  (tetr4)
! Quadratic  10-noded tetrahedron (tetr10)
! Tri-linear 8-noded hexahedron  (hexa8)
!
!-----------------------------------------------------------------------
!
! The following parameters are used to adjust the maximum size
! of the problem that the program can handle:
!
! mpoin --> maximum number of nodes
! melem --> maximum number of elements
! mdgof --> maxiumu number of degrees of freedom
! mprof --> maximum number of off-diagonal terms in tangent matrix
! mnode --> maximum number of nodes per element
! mgaus --> maximum number of Gauss points per element
! mmats --> maximum number of materials
! mbpel --> maximum number of pressure elements
! msearch > maximum number of line search iterations
!
!-----------------------------------------------------------------------
!
 IMPLICIT none  ! Always best
 INTEGER, PARAMETER :: mpoin = 100, melem = 100, mdgof = 300 
 INTEGER, PARAMETER :: mprof = 10000, mnode = 10, mgaus = 8
 INTEGER, PARAMETER :: mmats = 10, mbpel = 100, msearch = 5
 
 CHARACTER(80) title
 CHARACTER(10) eltyp
 LOGICAL rest
 
!     Dimensions nodal arrays, elemental arrays; degree of freedom
!     arrays and material properties vector
 
 INTEGER :: icode (mpoin), ldgof (3, mpoin), lnods (mnode, melem),  &
   matno (melem), lbnod (mnode, mbpel), kprof (mdgof),              &
   matyp (mmats), ndque (2*mdgof), nconn (2, mpoin)
 INTEGER :: ndime, nnode, ngaus, nnodb, ngaub, npoin, nelem, nmats, &
   nstrs, ndgof, negdf, nprof, nprs, nbpel, nincr, miter, incrm,    &
   niter, nsear

 DOUBLE PRECISION :: x (3, mpoin), x0 (3, mpoin),              &
   eledb (4, mnode+1, mgaus), stres (6, mgaus, melem),         &
   elecd (4, mnode, mgaus), vinc (mgaus), vol0 (melem),        &
   elacd (4, mnode), elbdb (3, mnode+1, mgaus), press (mbpel), &
   stifd (mdgof), stifp (mprof), eload (mdgof),                &
   pdisp (mdgof), resid (mdgof), displ (mdgof), xincr (mdgof), &
   react (mdgof), tload (mdgof), props (8, mmats), gravt (3)
 DOUBLE PRECISION :: xlmax, dlamb, cnorm, searc, arcln, xlamb, &
   rnorm, rtu0, r, eta0, eta, rtu  
 
!     Welcomes the user and determines whether the problem is being
!     restarted or a data file is to be read.
 
   CALL welcome (title, rest)
   IF ( .not. rest ) then
 
!     Reads in the initial nodal positions and boundary codes;
!     determines the element type to use and reads in the element
!     connectivities
 
   CALL elinfo (mnode, mgaus, ndime, nnode, ngaus, nnodb, ngaub,  &
                eledb, elbdb, eltyp)
   CALL innodes (mpoin, ndime, npoin, x, icode, ldgof)
   CALL inelems (melem, nelem, nnode, lnods, matno, nmats)
   nstrs = 4
   IF ( ndime == 3 ) nstrs = 6
 
!     Obtains degree of freedom numbers and profile information
!     in an optimum manner, by first finding the node to node
!     connections, then numbering the degrees of feedom following the
!     node to node connections and finally finds the profile addresses.
 
     !f95 CALL nodecon (npoin, nelem, nnode, lnods, stifp)
     !f95 CALL degfrm (mdgof, npoin, ndime, stifp, ndgof, negdf, &
     !f95              ldgof, stifd)
     CALL nodecon (npoin, nelem, nnode, lnods, nconn)
     CALL degfrm (mdgof, npoin, ndime, nconn, ndgof, negdf, ldgof,  &
                  ndque)
     CALL profile (mprof, ndgof, nelem, nnode, ndime, lnods, ldgof, &
                   nprof, kprof)
 
!     Reads in the external loads and prescribed displacements; the
!     material parameters and the iteration control information.
 
     CALL matprop (ndime, nmats, matyp, props)
     CALL inloads (ndime, npoin, ndgof, negdf, nnodb, mbpel, ldgof, &
                   eload, pdisp, gravt, nprs, nbpel, lbnod, press)
     CALL incontr (nincr, xlmax, dlamb, miter, cnorm, searc, &
                   arcln, 1)
 
!     Initializes elemental values, obtains element forces and
!     the initial stiffness matrix
 
     xlamb = 0.d0
     incrm = 0
     CALL initno (ndime, npoin, ndgof, nprof, x, x0, stifd, stifp,  &
                  resid)
     CALL initel (ndime, nnode, ngaus, nelem, gravt, x, eledb, &
                  lnods, matno, matyp, props, ldgof, eload, kprof, &
                  stifd, stifp, vinc, elecd, elacd, vol0)
 
!     Initialises load incrent variables and dumps everything to file
 
      CALL dump (title, eltyp, ndime, npoin, nnode, ngaus, nstrs,    &
                 nelem, ndgof, negdf, nprs, nprof, nmats, incrm,  &
                 xlamb, nbpel, nnodb, ngaub, matyp, props, matno, &
                 lnods, x, x0, kprof, stifd, stifp, resid, eload, &
                 ldgof, icode, eledb, pdisp, vol0, elbdb, lbnod, press)
    ELSE
 
!     Re-start the analysis from previous values
 
      CALL restar1 (title, eltyp, ndime, npoin, nnode, ngaus, nstrs, &
                    nelem, ndgof, negdf, nprs, nprof, nmats, incrm,  &
                    xlamb, nbpel, nnodb, ngaub)
      CALL restar2 (ndime, npoin, nnode, ngaus, nelem, ndgof, negdf, &
                    nprof, nmats, nnodb, ngaub, nbpel, matyp, props, &
                    matno, lnods, x, x0, kprof, stifd, stifp, resid, &
                    eload, ldgof, icode, eledb, pdisp, vol0, elbdb,  &
                    lbnod, press)
      CALL incontr (nincr, xlmax, dlamb, miter, cnorm, searc, arcln, 5)
    END IF
!
!     Starts the load increment loop
!
    DO WHILE ( (xlamb < xlmax) .and. (incrm < nincr) )
      incrm = incrm + 1
      xlamb = xlamb + dlamb
      CALL force (ndgof, dlamb, eload, tload, resid)
      CALL bpress (ndime, nnodb, ngaub, nbpel, dlamb, elbdb, lbnod,    &
                   press, x, ldgof, tload, react, resid, kprof, stifd, &
                   stifp)
!
!     If required imposes prescribed displacements, in which case the
!     stiffness matrix and the residuals need to be re-evaluated
!
      IF ( nprs > 0 ) then
        CALL prescx (ndime, npoin, ldgof, pdisp, x0, x, xlamb)
        CALL initrk (ndgof, nprof, negdf, xlamb, eload, tload, resid,  &
                     react, stifd, stifp)
        CALL bpress (ndime, nnodb, ngaub, nbpel, xlamb, elbdb, lbnod,  &
                     press, x, ldgof, tload, react, resid, kprof, &
                     stifd, stifp)
        CALL elemtk (ndime, nnode, ngaus, nstrs, nelem, x, x0, eledb,  &
                     lnods, matno, matyp, props, ldgof, stres, resid,  &
                     kprof, stifd, stifp, react, vinc, elecd, elacd, &
                     vol0)
      END IF
!
!     Starts the Newton-Raphson iteration
!
      niter = 0
      rnorm = 2 * cnorm
      DO WHILE ( (rnorm > cnorm) .and. (niter < miter) )
        niter = niter + 1
!
!     Calls the profile solver routines to obtain the displacement.
!     Also obtains the product r.u used for line searches.
!
        CALL datri (stifp, stifp, stifd, kprof, ndgof, .false., 6)
        CALL dasol (stifp, stifp, stifd, resid, displ, kprof, &
                    ndgof, 6, rtu0)                        ! Eq (7.80c)
!
!     If the arc length method is to be used, obtains the force
!     component of the displacement
!
        IF ( arcln /= 0.d0 ) then
          CALL dasol (stifp, stifp, stifd, tload, resid, kprof, &
                      ndgof, 6, r)
          CALL arclen (ndgof, niter, arcln, displ, resid, xincr, &
                       xlamb, dlamb)
        END IF
!
!     Starts the line search iteration. The total number of line search
!     iterations is limited to msearch
!
        eta0 = 0.d0 ; eta = 1.d0
        nsear = 0   ; rtu = rtu0 * searc * 2

        DO WHILE ( (abs (rtu)  > abs (rtu0 * searc) ) &
                   .and. (nsear < msearch) )
          nsear = nsear + 1
!
!     Updates the geometry and obtains the new residual forces and if
!     required implements the line search method
!
          CALL update (ndime, npoin, ldgof, x, displ, eta - eta0)
          CALL initrk (ndgof, nprof, negdf, xlamb, eload, tload, &
                       resid, react, stifd, stifp)
          CALL bpress (ndime, nnodb, ngaub, nbpel, xlamb, elbdb, &
                       lbnod, press, x, ldgof, tload, react, resid, &
                       kprof, stifd, stifp)
          CALL elemtk (ndime, nnode, ngaus, nstrs, nelem, x, x0, &
                       eledb, lnods, matno, matyp, props, ldgof, &
                       stres, resid, kprof, stifd, stifp, react, &
                       vinc, elecd, elacd, vol0)
          CALL search (ndgof, resid, displ, eta0, eta, rtu0, rtu)
        END DO
!
!     Checks for equilibrium convergence
!
          CALL checkr (incrm, niter, ndgof, negdf, xlamb, resid, &
                       tload, react, rnorm)
        END DO
!
!     If convergence was not achieved restarts from previous results
!
        IF ( niter >= miter ) then
          WRITE (6, 100)
          100 FORMAT(' Solution not converged, ', &
                   & 'restarting from previous step')
          CALL restar1 (title, eltyp, ndime, npoin, nnode, ngaus, &
                        nstrs, nelem, ndgof, negdf, nprs, nprof,  &
                        nmats, incrm, xlamb, nbpel, nnodb, ngaub)
          CALL restar2 (ndime, npoin, nnode, ngaus, nelem, ndgof, &
                        negdf, nprof, nmats, nnodb, ngaub, nbpel, &
                        matyp, props, matno, lnods, x, x0, kprof, &
                        stifd, stifp, resid, eload, ldgof, icode, &
                        eledb, pdisp, vol0, elbdb, lbnod, press)
          CALL incontr (nincr, xlmax, dlamb, miter, cnorm, searc, &
                        arcln, 5)
!
!     Otherwise writes results and dumps information to enable a
!     future restart if required
!
        ELSE
          CALL output (ndime, nnode, ngaus, nstrs, npoin, nelem, &
                       eltyp, title, icode, incrm, xlamb, x, lnods, &
                       ldgof, matno, stres, tload, react)
          CALL dump (title, eltyp, ndime, npoin, nnode, ngaus, &
                     nstrs, nelem, ndgof, negdf, nprs, nprof, nmats, &
                     incrm, xlamb, nbpel, nnodb, ngaub, matyp, props, &
                     matno, lnods, x, x0, kprof, stifd, stifp, resid, &
                     eload, ldgof, icode, eledb, pdisp, vol0, elbdb,  &
                     lbnod, press)
        END IF
      END DO
      STOP 'Normal end of PROGRAM flagshyp'
!
      END PROGRAM flagshyp
!
!-----------------------------------------------------------------------
SUBROUTINE welcome (title, rest)
!-----------------------------------------------------------------------
!
!     Opens input, output and restart files
!
!     title  --> example title
!     rest   --> logical variable: .true. is problem is restarted,
!                .false. if problem started from scratch.
!
!
!-----------------------------------------------------------------------
!
  IMPLICIT DOUBLE PRECISION (a - h, o - z)
  LOGICAL rest
  CHARACTER (20) name
  CHARACTER (80) title
  CHARACTER (1) ans
  DATA ans / ' ' /
!
!     Writes initial tile
!
      WRITE (6, 100)
      100 FORMAT(//////,                                            &
      &      10x,'          P R O G R A M    F L a g S H y P',///,  &
      &      10x,'Finite element LArGe Strain HYperelasticity ',    &
      &      'Program',////)
      WRITE (6, 101)
      101 FORMAT(' Is the problem starting from scratch (y/n) ?: ')
      READ (5, '(a1)') ans
      rest = .true.
!
!     Reads in the input file name
!
      IF ( ans == 'y' ) then
   10   WRITE (6, 102)
        102 FORMAT(' Enter the data file name   : ')
        READ (5, '(a)') name
        OPEN (1, file = name, status = 'old', form = 'formatted', &
              err = 10)
        READ (1, '(a)', err = 10) title
        rest = .false.
      END IF
!
!     Opens the output file
!
   20 WRITE (6, 103)
      103 FORMAT(' Enter the results file name: ')