!----------------------------------------------------------------------- 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, akin@rice.edu. ! 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, mprof) !b last entry from jhaemer@stanford.com to fix "out of bounds" !b matyp (mmats), ndque (2*mdgof), nconn (2, mpoin) !b 6/17/05 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: ') READ (5, '(a)') name IF ( .not. rest ) then OPEN (2, file = name, status = 'unknown', form = 'formatted', & err = 20) ELSE OPEN (2, file = name, status = 'old', form = 'formatted', & err = 20, action = 'write', position = 'append') !b access = 'append') END IF ! ! Opens the restart file ! 30 WRITE (6, 104) 104 FORMAT(' Enter the restart file name: ') READ (5, '(a)') name OPEN (3, file = name, status = 'unknown', form = 'unformatted', & err = 30) ! END SUBROUTINE welcome !----------------------------------------------------------------------- SUBROUTINE elinfo (mnode, mgaus, ndime, nnode, ngaus, nnodb, & ngaub, eledb, elbdb, eltyp) !----------------------------------------------------------------------- ! ! Obtains element information ! ! mnode --> maximum number of nodes per element ! mgaus --> maximum number of gauss points per element ! nnode --> number of nodes per element ! ngaus --> number of gauss points per element ! ndime --> number of dimensions for the given element ! nnodb --> number of nodes per boundary element ! ngaub --> number of gauss points per boundary element ! eledb --> element data matrix of dimensions ! (ndime+1,nnode+1,ngaus): ! ! | N1, N2,.., Nnnode, wgt | ! | | ! | dN1 dN2 dNnnode | ! | ---,---,...,-------, xi | ! | dxi dxi dxi | X number of Gauss points ! | | ! | dN1 dN2 dNnnode | ! | ---,---,...,-------, eta | ! | det det det | ! ! elbdb --> boundary elements data, format as above ! eltyp --> element type ! !----------------------------------------------------------------------- ! IMPLICIT none INTEGER, INTENT (IN ) :: mnode, mgaus INTEGER, INTENT (OUT) :: ndime, nnode, ngaus, nnodb, ngaub DOUBLE PRECISION, INTENT (OUT) :: eledb (4, mnode+1, mgaus), & elbdb (3, mnode+1, mgaus) CHARACTER (LEN=10), INTENT (OUT) :: eltyp ! ! First reads the element type the finds element information ! READ (1, '(a10)', err = 10) eltyp ! ! Three noded linear triangle ! IF ( eltyp (1:5) == 'tria3' ) then ndime = 2 ; nnode = 3 ; ngaus = 1 ; nnodb = 2 ; ngaub = 1 CALL tria3db (eledb) ; CALL lin2db (elbdb) ! ! Six noded quadratic triangle ! ELSEIF ( eltyp (1:5) == 'tria6' ) then ndime = 2 ; nnode = 6 ; ngaus = 3 ; nnodb = 3 ; ngaub = 2 CALL tria6db (eledb) ; CALL qua3db (elbdb) ! ! Four noded bi-linear quadrilateral ! ELSEIF ( eltyp (1:5) == 'quad4' ) then ndime = 2 ; nnode = 4 ; ngaus = 4 ; nnodb = 2 ; ngaub = 1 CALL quad4db (eledb) ; CALL lin2db (elbdb) ! ! Four noded linear tetrahedron ! ELSEIF ( eltyp (1:5) == 'tetr4' ) then ndime = 3 ; nnode = 4 ; ngaus = 1 ; nnodb = 3 ; ngaub = 1 CALL tetr4db (eledb) ; CALL tria3db (elbdb) ! ! Ten noded quadratic tetrahedron ! ELSEIF ( eltyp (1:6) == 'tetr10' ) then ndime = 3 ; nnode = 10 ; ngaus = 5 ; nnodb = 6 ; ngaub = 3 CALL tetr10db (eledb) ; CALL tria6db (elbdb) ! ! Tri-linear brick element ! ELSEIF ( eltyp (1:5) == 'hexa8' ) then ndime = 3 ; nnode = 8 ; ngaus = 8 ; nnodb = 4 ; ngaub = 4 CALL hexa8db (eledb) ; CALL quad4db (elbdb) ! ! No more elements implemented yet ! ELSE WRITE (6, '(a)') ' Unknown element type' STOP 'elinfo: Unknown element type' END IF ! ! Checks that the maximum dimensions are OK ! IF ( ngaus > mgaus ) then WRITE (6, 100) ngaus 100 FORMAT(' problem dimensions exceed maximum, ', & & 'set mgaus to:',i10) STOP 'elinfo: dimensions exceed mgaus' ELSEIF ( nnode > mnode ) then WRITE (6, 101) nnode 101 FORMAT(' problem dimensions exceed maximum, ', & & 'set mnode to:',i10) STOP 'elinfo: dimensions exceed mnode' END IF RETURN ! 10 WRITE (6, '(a)') ' Error reading element type' STOP 'elinfo: Error reading element type' END SUBROUTINE elinfo !----------------------------------------------------------------------- SUBROUTINE lin2db (eledb) !----------------------------------------------------------------------- ! ! Two noded linear element ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION eledb (2, 3, 1) ! eledb (1, 1, 1) = 0.5d0 eledb (2, 1, 1) = - 0.5d0 eledb (1, 2, 1) = 0.5d0 eledb (2, 2, 1) = 0.5d0 eledb (1, 3, 1) = 2.d0 eledb (2, 3, 1) = 0.d0 END SUBROUTINE lin2db !----------------------------------------------------------------------- SUBROUTINE qua3db (eledb) !----------------------------------------------------------------------- ! ! Three noded quadratic element ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION eledb (2, 4, 2), gauss (2) DATA gauss / - 0.577350269189626d0, 0.577350269189626d0 / ! shap1 (xi) = xi * (xi - 1) / 2.d0 ! obsolete under f95 shap2 (xi) = 1 - xi * xi shap3 (xi) = xi * (xi + 1) / 2.d0 dshp1 (xi) = xi - 0.5d0 dshp2 (xi) = - 2 * xi dshp3 (xi) = xi + 0.5d0 ! DO ig = 1, 2 xi = gauss (ig) eledb (1, 1, ig) = shap1 (xi) eledb (2, 1, ig) = dshp1 (xi) eledb (1, 2, ig) = shap2 (xi) eledb (2, 2, ig) = dshp2 (xi) eledb (1, 3, ig) = shap3 (xi) eledb (2, 3, ig) = dshp3 (xi) eledb (1, 4, ig) = 1.d0 eledb (2, 4, ig) = xi END DO END SUBROUTINE qua3db !----------------------------------------------------------------------- SUBROUTINE tria3db (eledb) !----------------------------------------------------------------------- ! ! Three noded linear triangle ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION eledb (3, 4, 1) ! ! Gauss point and weight ! eledb (1, 4, 1) = 1.d0 / 2.d0 eledb (2, 4, 1) = 1.d0 / 3.d0 eledb (3, 4, 1) = 1.d0 / 3.d0 ! ! shape function and derivatives at gauss point ! eledb (1, 1, 1) = 1.d0 / 3.d0 eledb (2, 1, 1) = - 1.d0 eledb (3, 1, 1) = - 1.d0 eledb (1, 2, 1) = 1.d0 / 3.d0 eledb (2, 2, 1) = 1.d0 eledb (3, 2, 1) = 0.d0 eledb (1, 3, 1) = 1.d0 / 3.d0 eledb (2, 3, 1) = 0.d0 eledb (3, 3, 1) = 1.d0 END SUBROUTINE tria3db !----------------------------------------------------------------------- SUBROUTINE tria6db (eledb) !----------------------------------------------------------------------- ! ! Six noded linear triangle with midside integration ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION eledb (3, 7, 3), gauss (2, 3) DATA (gauss (i, 1), i = 1, 2) / 0.5d0, 0.0d0 / DATA (gauss (i, 2), i = 1, 2) / 0.5d0, 0.5d0 / DATA (gauss (i, 3), i = 1, 2) / 0.0d0, 0.5d0 / ! ! Shape functions and derivatives ! obsolete under f95 ! shap1 (xi, et) = (xi + et - 1) * (2 * xi + 2 * et - 1) shap2 (xi, et) = 4 * xi * (1 - xi - et) shap3 (xi, et) = xi * (2 * xi - 1) shap4 (xi, et) = 4 * xi * et shap5 (xi, et) = et * (2 * et - 1) shap6 (xi, et) = 4 * et * (1 - xi - et) dsh1x (xi, et) = 4 * xi + 4 * et - 3 dsh2x (xi, et) = 4 - 8 * xi - 4 * et dsh3x (xi, et) = 4 * xi - 1 dsh4x (xi, et) = 4 * et dsh5x (xi, et) = 0 dsh6x (xi, et) = - 4 * et dsh1e (xi, et) = 4 * xi + 4 * et - 3 dsh2e (xi, et) = - 4 * xi dsh3e (xi, et) = 0 dsh4e (xi, et) = 4 * xi dsh5e (xi, et) = 4 * et - 1 dsh6e (xi, et) = 4 - 4 * xi - 8 * et ! ! Constructs the data base. First Gauss point position and weight ! DO ig = 1, 3 eledb (1, 7, ig) = 1.d0 / 6.d0 xi = gauss (1, ig) et = gauss (2, ig) eledb (2, 7, ig) = xi eledb (3, 7, ig) = et ! ! Shape functions and derivatives ! eledb (1, 1, ig) = shap1 (xi, et) eledb (2, 1, ig) = dsh1x (xi, et) eledb (3, 1, ig) = dsh1e (xi, et) eledb (1, 2, ig) = shap2 (xi, et) eledb (2, 2, ig) = dsh2x (xi, et) eledb (3, 2, ig) = dsh2e (xi, et) eledb (1, 3, ig) = shap3 (xi, et) eledb (2, 3, ig) = dsh3x (xi, et) eledb (3, 3, ig) = dsh3e (xi, et) eledb (1, 4, ig) = shap4 (xi, et) eledb (2, 4, ig) = dsh4x (xi, et) eledb (3, 4, ig) = dsh4e (xi, et) eledb (1, 5, ig) = shap5 (xi, et) eledb (2, 5, ig) = dsh5x (xi, et) eledb (3, 5, ig) = dsh5e (xi, et) eledb (1, 6, ig) = shap6 (xi, et) eledb (2, 6, ig) = dsh6x (xi, et) eledb (3, 6, ig) = dsh6e (xi, et) END DO END SUBROUTINE tria6db !----------------------------------------------------------------------- SUBROUTINE tetr4db (eledb) !----------------------------------------------------------------------- ! ! Four noded linear tetrahedron ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION eledb (4, 5, 1) ! ! Gauss point and weight ! eledb (1, 5, 1) = 1.d0 / 6.d0 eledb (2, 5, 1) = 1.d0 / 4.d0 eledb (3, 5, 1) = 1.d0 / 4.d0 eledb (4, 5, 1) = 1.d0 / 4.d0 ! ! shape function and derivatives at gauss point ! eledb (1, 1, 1) = 1.d0 / 4.d0 eledb (2, 1, 1) = - 1.d0 eledb (3, 1, 1) = - 1.d0 eledb (4, 1, 1) = - 1.d0 eledb (1, 2, 1) = 1.d0 / 4.d0 eledb (2, 2, 1) = 1.d0 eledb (3, 2, 1) = 0.d0 eledb (4, 2, 1) = 0.d0 eledb (1, 3, 1) = 1.d0 / 4.d0 eledb (2, 3, 1) = 0.d0 eledb (3, 3, 1) = 1.d0 eledb (4, 3, 1) = 0.d0 eledb (1, 4, 1) = 1.d0 / 4.d0 eledb (2, 4, 1) = 0.d0 eledb (3, 4, 1) = 0.d0 eledb (4, 4, 1) = 1.d0 END SUBROUTINE tetr4db !----------------------------------------------------------------------- SUBROUTINE tetr10db (eledb) !----------------------------------------------------------------------- ! ! Ten noded quadratic tetrahedron ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) PARAMETER (sixth = 0.166666666666667d0) DIMENSION eledb (4, 11, 5), gauss (3, 5) ! ! Gauss point positions ! DATA (gauss (i, 1), i = 1, 3) / 0.25d0, 0.25d0, 0.25d0 / DATA (gauss (i, 2), i = 1, 3) / sixth, sixth, sixth / DATA (gauss (i, 3), i = 1, 3) / 0.5d0, sixth, sixth / DATA (gauss (i, 4), i = 1, 3) / sixth, 0.5d0, sixth / DATA (gauss (i, 5), i = 1, 3) / sixth, sixth, 0.5d0 / ! ! Shape functions and derivatives ! obsolete under f95 ! shap1 (xi, et, ze) = (1 - 2 * xi - 2 * et - 2 * ze) & * (1 - xi - et - ze) shap2 (xi, et, ze) = xi * (2 * xi - 1) shap3 (xi, et, ze) = et * (2 * et - 1) shap4 (xi, et, ze) = ze * (2 * ze - 1) shap5 (xi, et, ze) = 4 * xi * (1 - xi - et - ze) shap6 (xi, et, ze) = 4 * xi * et shap7 (xi, et, ze) = 4 * et * (1 - xi - et - ze) shap8 (xi, et, ze) = 4 * ze * (1 - xi - et - ze) shap9 (xi, et, ze) = 4 * xi * ze shap0 (xi, et, ze) = 4 * ze * et ! dsh1x (xi, et, ze) = 4 * (xi + et + ze) - 3 dsh2x (xi, et, ze) = 4 * xi - 1 dsh3x (xi, et, ze) = 0 dsh4x (xi, et, ze) = 0 dsh5x (xi, et, ze) = 4 * (1 - 2 * xi - et - ze) dsh6x (xi, et, ze) = 4 * et dsh7x (xi, et, ze) = - 4 * et dsh8x (xi, et, ze) = - 4 * ze dsh9x (xi, et, ze) = 4 * ze 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 ! 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 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 !----------------------------------------------------------------------- SUBROUTINE pforce (ndime, nnodb, press, xlamb, elbdb, lbnod, & ldgof, tload, resid, react, dxis) !----------------------------------------------------------------------- ! ! Evaluates the nodal forces due to boundary normal pressure and ! assembles it into the residual and total nominal load vector ! ! ndime --> number of dimensions ! nnodb --> number of nodes per boundary element ! press --> applied pressure ! xlamb --> force factor lambda ! elbdb --> bondary elements information ! lbnod --> boundary element connectivities ! ldgof --> degrees of freedom number of each node ! tload --> total nominal load ! react --> reactions ! resid --> residual forces ! react --> reactions ! dxis --> derivatives of geometry with respect isoparametric ! coordinates ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION elbdb (ndime, nnodb + 1), lbnod (nnodb), dxis (3, 2), & xnorm (3), ldgof (ndime, * ), tload ( * ), react ( * ), & resid ( * ) ! ! First finds the normal vector ! xnorm (1) = dxis (2, 1) * dxis (3, 2) - dxis (2, 2) * dxis (3, 1) xnorm (2) = dxis (1, 2) * dxis (3, 1) - dxis (1, 1) * dxis (3, 2) xnorm (3) = dxis (1, 1) * dxis (2, 2) - dxis (1, 2) * dxis (2, 1) ! ! Obtains and assembles the nodal forces ! DO in = 1, nnodb ip = lbnod (in) vl = press * elbdb (1, in) * elbdb (1, nnodb + 1) DO id = 1, ndime if = ldgof (id, ip) IF ( if > 0 ) then tload (if) = tload (if) + xnorm (id) * vl resid (if) = resid (if) + xnorm (id) * vl * xlamb ELSE react ( - if) = react ( - if) - xnorm (id) * vl * xlamb END IF END DO END DO END SUBROUTINE pforce !----------------------------------------------------------------------- SUBROUTINE kpress (ndime, nnodb, press, elbdb, lbnod, ldgof, dxis,& kprof, stifd, stifp) !----------------------------------------------------------------------- ! ! Obtains and assembles the pressure component of the stiffness matrix ! ! ndime --> number of dimensions ! nnodb --> number of nodes ! press --> pressure ! elbdb --> element information ! lbnod --> element connectivities ! ldgof --> nodal degrees of freedom ! dxis --> tangent vectors ! kprof --> profile addresses ! stifd --> diagonal stiffness matrix ! stifp --> profile component of K ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION lbnod ( * ), ldgof (ndime, * ), dxis (3, 2), elbdb ( & ndime, * ), kprof ( * ), stifd ( * ), stifp ( * ), eps (3, 3, 3) DATA ( (eps (i, j, 1), i = 1, 3), j = 1, 3) & / 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, - 1.d0, 0.d0, 1.d0, 0.d0 / DATA ( (eps (i, j, 2), i = 1, 3), j = 1, 3) & / 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, - 1.d0, 0.d0, 0.d0 / DATA ( (eps (i, j, 3), i = 1, 3), j = 1, 3) & / 0.d0, - 1.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0 / ! ! Loops over nodes ! val = - press * elbdb (1, nnodb + 1) / 2.d0 DO in = 1, nnodb ip = lbnod (in) DO jn = in, nnodb jp = lbnod (jn) val2 = val * (elbdb (1, in) * elbdb (2, jn) & - elbdb (1, jn) * elbdb (2, in) ) val3 = 0.d0 IF ( ndime == 3 ) val3 = val * (elbdb (1, in) * elbdb (3, jn)& - elbdb (1, jn) * elbdb (3, in) ) ! ! Loops over dimensions ! DO 10 id = 1, ndime if = ldgof (id, ip) IF ( if <= 0 ) goto 10 DO 20 jd = 1, ndime jf = ldgof (jd, jp) IF ( jf <= 0 ) goto 20 IF ( (in == jn ) .and. (jd < id ) ) goto 20 ! ! Uses equation (7.49) to evaluate the constitutive component ! of K in an indicial form ! sum = 0.d0 DO kd = 1, 3 sum = sum + eps (id, jd, kd) * (val2 * dxis (kd, 2) & - val3 * dxis (kd, 1) ) END DO ! ! Assembles sum in the right place. First case: diagonal term. ! IF ( jf == if ) then stifd (if) = stifd (if) + sum ELSE ! ! Then the off-diagonal terms into the profile. ! jc = max (jf, if) ir = min (jf, if) lp = kprof (jc) - jc + ir + 1 stifp (lp) = stifp (lp) + sum END IF 20 END DO 10 END DO END DO END DO END SUBROUTINE kpress !----------------------------------------------------------------------- SUBROUTINE elemtk (ndime, nnode, ngaus, nstrs, nelem, x, x0, eledb, & lnods, matno, matyp, props, ldgof, stres, resid, & kprof, stifd, stifp, react, vinc, elecd, elacd, vol0) !----------------------------------------------------------------------- ! ! Obtains the elemental stresses, assembles the internal equivalent ! forces and the tangent matrix ! ! ndime --> number of dimensions ! nnode --> number of nodes ! ngaus --> number of Gauss points ! nstrs --> number of stresses per Gauss point ! nelem --> gravity vector ! x --> nodal coordinates ! x0 --> initial nodal coordinates ! eledb --> element data base ! lnods --> element connectivities ! matno --> material number ! matyp --> material type ! props --> material properties ! ldgof --> degree of freedom number ! resid --> residual force ! stres --> elemetal stresses ! kprof --> profile addresses ! stifd --> diagonal of K ! stifp --> profile part of K ! react --> accumulated reactions ! vinc --> jacobian at each Gauss point ! elecd --> Cartesian derivatives at each Gauss point ! elacd --> average Cartesian derivatives ! vol0 --> intial element volume ! ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION lnods (nnode, * ), props (8, * ), ldgof (ndime, * ), & kprof ( * ), stifd ( * ), stifp ( * ), x (ndime, * ), x0 (ndime, & * ), matno ( * ), eledb (ndime+1, nnode+1, ngaus), elacd (ndime, & * ), resid ( * ), ctens (3, 3, 3, 3), sigma (3, 3), btens (3, 3),& matyp ( * ), react ( * ), stres (nstrs, ngaus, * ), stret (3), & princ (3, 3), sprin (3), vol0 ( * ), elecd (ndime, nnode, * ), & vinc ( * ) DATA ( (btens (i, j), i = 1, 3), j = 1, 3) & / 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0 / ! ! Loops over all the elements. ! First finds the gradient of the shape functions and Jacobian at ! each Gauss point ! DO ie = 1, nelem 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, theta (the volume ratio) ! and the element pressure. First for nearly incompressible ! neo-Hookean materials. ! IF ( mat == 5 ) then CALL getheta (ngaus, vol0 (ie), vinc, theta) ! Eq (6.50a) xkapp = props (4, im) press = xkapp * (theta - 1.d0) ! E2 (6.52) xkapp = xkapp * theta ! Eq (6.61) CALL kvolume (ndime, nnode, ngaus, xkapp, vinc, elecd, & elacd, lnods (1, ie), ldgof, kprof, & stifd, stifp) ! ! Also for nearly incompressible material in principal directions ! ELSEIF ( mat == 7 ) then CALL getheta (ngaus, vol0 (ie), vinc, theta) ! Eq (6.50a) xkapp = props (4, im) press = xkapp * log (theta) / theta ! Eq (5.99) xkapp = (xkapp / theta) - press CALL kvolume (ndime, nnode, ngaus, xkapp, vinc, elecd, & elacd, lnods (1, ie), ldgof, kprof, & stifd, stifp) END IF ! ! Loops over the Gauss points and first obtains the left ! Cauchy-Green tensor b ! DO ig = 1, ngaus CALL leftcg (ndime, nnode, lnods (1, ie), x0, & elecd (1, 1, ig), detf, btens) ! ! Material 1 compressible neo-Hookean ! IF ( mat == 1 ) then xlamb = props (3, im) / detf xmu = props (2, im) / detf xme = xmu - xlamb * log (detf) ! Eq (5.38) CALL stress1 (ndime, xmu, xme, btens, sigma) ! Eq (5.29) CALL cisotp (ndime, xlamb, xme, ctens) ! Eq (5.37) ! ! Material 3 hyperelasticity in principal directions ! ELSEIF ( mat == 3 ) then xlam = props (3, im) / detf xmu = props (2, im) / detf CALL jacobi (btens, stret, princ) ! for principal values CALL stress3 (ndime, xmu, xlam, detf, stret, princ, & sigma, sprin) ! Eq (5.77 & 5.90) CALL cprinc (ndime, xlam, xmu, stret, princ, sprin, & ctens) ! Eq (5.86, 5.87, 5.91) ! ! Material 4 plane stress hyperelastic in principal directions ! ELSEIF ( mat == 4 ) then det3d = detf**props (5, im) ! Eq (5.113) xlam = props (3, im) / det3d xmu = props (2, im) / det3d CALL jacobi (btens, stret, princ) ! for principal values CALL stress3 (ndime, xmu, xlam, detf, stret, princ, & sigma, sprin) ! Eq (5.114) CALL cprinc (ndime, xlam, xmu, stret, princ, sprin, ctens) vinc (ig) = vinc (ig) * props (4, im) * det3d / detf stres (4, ig, ie) = props (4, im) * det3d / detf ! ! Material 5 nearly incompressible neo-Hookean ! ELSEIF ( mat == 5 ) then xmu = props (2, im) CALL stress5 (ndime, xmu, detf, btens, sigma) ! Eq (5.51) CALL addpres (ndime, press, sigma) ! Eq (5.51) CALL cdevia (ndime, xmu, detf, btens, ctens) ! Eq (5.55a) CALL cvolum (ndime, press, ctens) ! Eq (5.55b) ! ! Material 6 plane stress incompressible neo-Hookean ! ELSEIF ( mat == 6 ) then ! (See exercise 5.1 and p. 203) xmu = props (2, im) xme = xmu / (detf * detf) xla = 2.d0 * xme ! Eq (5.37) CALL stress6 (ndime, xmu, detf, btens, sigma) CALL cisotp (ndime, xla, xme, ctens) vinc (ig) = vinc (ig) * props (4, im) / detf stres (4, ig, ie) = props (4, im) / detf ! ! Material 7 nearly incompressible in principal directions ! ELSEIF ( mat == 7 ) then xmu = props (2, im) / detf xlam = - 2.d0 * xmu / 3.d0 CALL jacobi (btens, stret, princ) ! for principal values CALL stress3 (ndime, xmu, xlam, detf, stret, princ, & sigma, sprin) ! Eq (5,90, 5.105) CALL addpres (ndime, press, sigma) CALL cprinc (ndime, xlam, xmu, stret, princ, sprin, & ctens) ! Eq (5.55a) CALL cvolum (ndime, press, ctens) ! Eq (5.55a) ! ! Material 8 plane stress incompressible in principal directions ! ELSEIF ( mat == 8 ) then xmu = props (2, im) xlam = 2.d0 * xmu CALL jacobi (btens, stret, princ) ! for principal values CALL stress3 (ndime, xmu, xlam, detf, stret, princ, & sigma, sprin) CALL cprinc (ndime, xlam, xmu, stret, princ, sprin, ctens) vinc (ig) = vinc (ig) * props (4, im) / detf stres (4, ig, ie) = props (4, im) / detf END IF ! ! Evaluates the internal equivalent forces and tangent matrix ! CALL internal (ndime, nnode, lnods (1, ie), ldgof, & sigma, elecd (1, 1, ig), vinc (ig), & resid, react) ! Eq (7.15b) CALL kconst (ndime, nnode, lnods (1, ie), ldgof, ctens, & elecd (1, 1, ig), vinc (ig), kprof, stifd, & stifp) ! Eq (7.35) CALL ksigma (ndime, nnode, lnods (1, ie), ldgof, sigma, & elecd (1, 1, ig), vinc (ig), kprof, stifd, & stifp) ! Eq (7.45c) ! ! Stores the stresses in stres ! ks = 1 DO id = 1, ndime DO jd = id, ndime stres (ks, ig, ie) = sigma (id, jd) ks = ks + 1 END DO END DO END DO ! Gauss loop END DO ! elements loop END SUBROUTINE elemtk !----------------------------------------------------------------------- SUBROUTINE gradsh (ie, ndime, nnode, ngaus, eledb, lnods, x, & elecd, vinc) !----------------------------------------------------------------------- ! ! Evaluates the Cartesian derivatives of the shape functions ! and the Jacobean at all Gauss points ! ! ie --> element number ! ndime --> number of dimensions ! nnode --> number of nodes ! ngaus --> number of Gauss points ! eledb --> element data base ! lnods --> element connectivities ! x --> nodal coordinates ! elecd --> Cartesian derivatives ! ! | dN1/dx, dN2/dx,... | ! | dN1/dy, dN2/dy,... | x ngaus ! | dN1/dz, dN2/dz,... | ! ! vinc --> jacobian at each Gauss point !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION eledb (ndime+1, nnode+1, * ), lnods (nnode, * ), & vinc ( * ), x (ndime, * ), xj (3, 3), xji (3, 3), elecd (ndime, & nnode, * ) DATA ( (xj (i, j), i = 1, 3), j = 1, 3) & / 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0 / ! ! First loops over all the gauss points ! DO ig = 1, ngaus weight = eledb (1, nnode+1, ig) ! ! Initialises the jacobian matrix. ! Puts a 1 in position 3,3 of xj so no special case is made for 2D. ! IF ( ndime == 2 ) xj (3, 3) = 1.d0 ! ! creates the Jacobian matrix using equation (7.6b) ! DO id = 1, ndime DO jd = 1, ndime xj (id, jd) = 0.d0 DO in = 1, nnode ip = lnods (in, ie) xj (id, jd) = xj (id, jd) & + x (id, ip)*eledb (jd+1, in, ig) END DO END DO END DO ! ! Determines the determinant of xj and checks for zero values ! detjb = xj (1, 1) * xj (2, 2) * xj (3, 3) & + xj (1, 2) * xj (2, 3) * xj (3, 1) & + xj (1, 3) * xj (2, 1) * xj (3, 2) & - xj (1, 3) * xj (2, 2) * xj (3, 1) & - xj (1, 2) * xj (2, 1) * xj (3, 3) & - xj (1, 1) * xj (2, 3) * xj (3, 2) vinc (ig) = detjb * weight IF ( detjb <= 0.d0 ) then WRITE (6, 100) ie 100 FORMAT(' Negative Jacobian for element: ',i6) STOP 'gradsh: Negative Jacobian in an element' END IF ! ! Determines the inverse of xj ! xji (1, 1) = (xj (2, 2) * xj (3, 3) & - xj (2, 3) * xj (3, 2) ) / detjb xji (1, 2) = (xj (1, 3) * xj (3, 2) & - xj (1, 2) * xj (3, 3) ) / detjb xji (1, 3) = (xj (1, 2) * xj (2, 3) & - xj (1, 3) * xj (2, 2) ) / detjb xji (2, 1) = (xj (2, 3) * xj (3, 1) & - xj (2, 1) * xj (3, 3) ) / detjb xji (2, 2) = (xj (1, 1) * xj (3, 3) & - xj (1, 3) * xj (3, 1) ) / detjb xji (2, 3) = (xj (1, 3) * xj (2, 1) & - xj (1, 1) * xj (2, 3) ) / detjb xji (3, 1) = (xj (2, 1) * xj (3, 2) & - xj (2, 2) * xj (3, 1) ) / detjb xji (3, 2) = (xj (1, 2) * xj (3, 1) & - xj (1, 1) * xj (3, 2) ) / detjb xji (3, 3) = (xj (1, 1) * xj (2, 2) & - xj (1, 2) * xj (2, 1) ) / detjb ! ! Finally obtains the Cartesian gradients using equation (7.6a) ! DO in = 1, nnode DO id = 1, ndime elecd (id, in, ig) = 0.d0 DO jd = 1, ndime elecd (id, in, ig) = elecd (id, in, ig) & + xji (jd, id) * eledb (jd+1, in, ig) END DO END DO END DO END DO ! for Gauss loop ! END SUBROUTINE gradsh !----------------------------------------------------------------------- SUBROUTINE leftcg (ndime, nnode, lnods, x0, gradn, detf, btens) !----------------------------------------------------------------------- ! ! Evaluates the left Cauchy-Green tensor given the initial ! nodal positions and the current shape function derivatives ! ! ndime --> number of dimensions ! nnode --> number of nodes ! lnods --> element connectivities ! x0 --> initial nodal coordinates ! gradn --> Cartesian gradient of shape functions: ! detf --> determinant of F (i.e. J) ! btens --> left Cauchy-Green tensor ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION lnods (nnode), gradn (ndime, nnode), x0 (ndime, * ), & btens (3, 3), finvr (3, 3), ftens (3, 3) DATA ( (finvr (i, j), i = 1, 3), j = 1, 3) & / 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0 / ! ! Initializes by putting a 1. in position 3,3 to cope with 2-D ! IF ( ndime == 2 ) finvr (3, 3) = 1.d0 ! ! Creates F inverse using equation (3.11) ! DO id = 1, ndime DO jd = 1, ndime finvr (id, jd) = 0.d0 DO in = 1, nnode ip = lnods (in) finvr (id, jd) = finvr (id, jd) & + x0 (id, ip) * gradn (jd, in) END DO END DO END DO ! ! Obtains the determinant ! detf = finvr (1, 1) * finvr (2, 2) * finvr (3, 3) & + finvr (1, 2) * finvr (2, 3) * finvr (3, 1) & + finvr (1, 3) * finvr (2, 1) * finvr (3, 2) & - finvr (1, 3) * finvr (2, 2) * finvr (3, 1) & - finvr (1, 2) * finvr (2, 1) * finvr (3, 3) & - finvr (1, 1) * finvr (2, 3) * finvr (3, 2) IF ( detf == 0.d0 ) detf = 1.d0 detf = 1.d0 / detf ! ! Determines F ! ftens (1, 1) = (finvr (2, 2) * finvr (3, 3) & - finvr (2, 3) * finvr (3, 2) ) * detf ftens (1, 2) = (finvr (1, 3) * finvr (3, 2) & - finvr (1, 2) * finvr (3, 3) ) * detf ftens (1, 3) = (finvr (1, 2) * finvr (2, 3) & - finvr (1, 3) * finvr (2, 2) ) * detf ftens (2, 1) = (finvr (2, 3) * finvr (3, 1) & - finvr (2, 1) * finvr (3, 3) ) * detf ftens (2, 2) = (finvr (1, 1) * finvr (3, 3) & - finvr (1, 3) * finvr (3, 1) ) * detf ftens (2, 3) = (finvr (1, 3) * finvr (2, 1) & - finvr (1, 1) * finvr (2, 3) ) * detf ftens (3, 1) = (finvr (2, 1) * finvr (3, 2) & - finvr (2, 2) * finvr (3, 1) ) * detf ftens (3, 2) = (finvr (1, 2) * finvr (3, 1) & - finvr (1, 1) * finvr (3, 2) ) * detf ftens (3, 3) = (finvr (1, 1) * finvr (2, 2) & - finvr (1, 2) * finvr (2, 1) ) * detf ! ! Finally obtains the tensor b using equation (7.9c) ! IF ( ndime == 2 ) btens (3, 3) = 1.d0 DO id = 1, ndime DO jd = 1, ndime btens (id, jd) = 0.d0 DO kd = 1, ndime btens (id, jd) = btens (id, jd) & + ftens (id, kd) * ftens (jd, kd) END DO END DO END DO END SUBROUTINE leftcg !----------------------------------------------------------------------- SUBROUTINE jacobi (btens, stret, princ) !----------------------------------------------------------------------- ! ! Evaluates the stretches and principal directions given the b matrix ! using the Jacobi iteration. Adapted from numerical recpies ! ! btens --> left Cauchy-Green tensor ! stret --> vector containing the stretches ! princ --> matrix containing the three principal column vectors ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) INTEGER, PARAMETER :: msweep = 50 DIMENSION btens (3, 3), stret (3), princ (3, 3) ! ! ! Initialise princ to the identity ! DO i = 1, 3 DO j = 1, 3 princ (i, j) = 0.d0 END DO princ (i, i) = 1.d0 stret (i) = btens (i, i) END DO ! ! Starts sweeping. ! DO is = 1, msweep sum = 0.d0 DO ip = 1, 2 DO iq = ip + 1, 3 sum = sum + abs (btens (ip, iq) ) END DO END DO ! ! If the sum of off-diagonal terms is zero evaluates the ! stretches and returns ! IF ( sum == 0.d0 ) then DO i = 1, 3 stret (i) = sqrt (stret (i) ) END DO RETURN END IF ! ! Performs the sweep in three rotations. One per off diagonal term ! DO ip = 1, 2 DO iq = ip + 1, 3 od = 100.d0 * abs (btens (ip, iq) ) IF ( (od+abs (stret (ip) ) /= abs (stret (ip) )) & .and. (od+abs (stret (iq) ) /= abs (stret (iq) ))) then hd = stret (iq) - stret (ip) ! ! Evaluates the rotation angle ! IF ( abs (hd) + od == abs (hd) ) then t = btens (ip, iq) / hd ELSE theta = 0.5d0 * hd / btens (ip, iq) t = 1.d0 / (abs (theta) + sqrt (1.d0 + theta**2) ) IF ( theta < 0.d0 ) t = - t END IF ! ! Re-evaluates the diagonal terms ! c = 1.d0 / sqrt (1.d0 + t**2) s = t * c tau = s / (1.d0 + c) h = t * btens (ip, iq) stret (ip) = stret (ip) - h stret (iq) = stret (iq) + h ! ! Re-evaluates the remaining off-diagonal terms ! ir = 6 - ip - iq g = btens (min (ir, ip), max (ir, ip) ) h = btens (min (ir, iq), max (ir, iq) ) btens (min (ir, ip), max (ir, ip) ) = g & - s * (h + g * tau) btens (min (ir, iq), max (ir, iq) ) = h & + s * (g - h * tau) ! ! Rotates the eigenvectors ! DO ir = 1, 3 g = princ (ir, ip) h = princ (ir, iq) princ (ir, ip) = g - s * (h + g * tau) princ (ir, iq) = h + s * (g - h * tau) END DO END IF btens (ip, iq) = 0.d0 END DO END DO END DO ! over sweeps ! ! If convergence is not achieved stops ! WRITE (6, 100) 100 FORMAT(' Jacobi iteration unable to converge') STOP ' Jacobi iteration unable to converge' END SUBROUTINE jacobi !----------------------------------------------------------------------- SUBROUTINE stress1 (ndime, xmu, xme, btens, sigma) !----------------------------------------------------------------------- ! ! Determines the Cauchy stress tensor for compressible neo-Hookean ! materials ! ! ndime --> number of dimensions ! xmu --> mu coefficient ! xme --> effective mu coefficient ! btens --> right Cauchy-Green tensor ! sigma --> Cauchy stress tensor ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION sigma (3, 3), btens (3, 3) ! ! Obtains the Cauchy stress using equation (5.29) ! DO id = 1, ndime DO jd = 1, ndime sigma (id, jd) = xmu * btens (id, jd) END DO sigma (id, id) = sigma (id, id) - xme END DO END SUBROUTINE stress1 !----------------------------------------------------------------------- SUBROUTINE stress3 (ndime, xmu, xlamb, detf, stret, princ, sigma,& sprin) !----------------------------------------------------------------------- ! ! Determines the Cauchy stress tensor for materials defined in ! principal directions ! ! ndime --> number of dimensions ! xmu --> mu coefficient ! xlamb --> lambda coefficient ! detf --> determinant of F, i.e. J ! stret --> vector containing the stretches ! princ --> matrix containing the the principal directions ! sigma --> Cauchy stress tensor ! sprin --> principal stresses ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION sigma (3, 3), princ (3, 3), stret (3), sprin (3) ! ! Obtains the principal stress using equation (5.90) ! DO id = 1, ndime sprin (id) = xlamb * log (detf) + 2.d0 * xmu * log (stret (id)) END DO ! ! Obtains the Cartesian Cauchy stress tensor using Eq (4.8), (5.77) ! DO jd = 1, ndime DO kd = 1, ndime sum = 0.d0 DO id = 1, ndime sum = sum + sprin (id) * princ (jd, id) * princ (kd, id) END DO sigma (jd, kd) = sum END DO END DO END SUBROUTINE stress3 !----------------------------------------------------------------------- SUBROUTINE stress5 (ndime, xmu, detf, btens, sigma) !----------------------------------------------------------------------- ! ! Determines the Cauchy stress tensor for nearly incompressible ! neo-Hookean materials ! ! ! ndime --> number of dimensions ! xmu --> mu coefficient ! detf --> determinant of F ! btens --> right Cauchy-Green tensor ! sigma --> Cauchy stress tensor ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION sigma (3, 3), btens (3, 3) ! ! Obtains the Cauchy stress using equation (5.51) ! First finds the trace ! trace = 0.d0 DO id = 1, 3 trace = trace + btens (id, id) END DO ! xm = xmu / (detf** (5.d0 / 3.d0) ) DO id = 1, ndime DO jd = 1, ndime sigma (id, jd) = xm * btens (id, jd) END DO sigma (id, id) = sigma (id, id) - xm * trace / 3.d0 END DO END SUBROUTINE stress5 !----------------------------------------------------------------------- SUBROUTINE stress6 (ndime, xmu, detf, btens, sigma) !----------------------------------------------------------------------- ! ! Determines the Cauchy stress tensor for plane stress incompressible ! neo-Hookean materials ! ! ! ndime --> number of dimensions ! xmu --> mu coefficient ! detf --> determinant of F ! btens --> right Cauchy-Green tensor ! sigma --> Cauchy stress tensor ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION sigma (3, 3), btens (3, 3) ! ! Obtains the Cauchy stress using equation on p.203 ! DO id = 1, ndime DO jd = 1, ndime sigma (id, jd) = xmu * btens (id, jd) END DO sigma (id, id) = sigma (id, id) - xmu / (detf**2) END DO END SUBROUTINE stress6 !----------------------------------------------------------------------- SUBROUTINE addpres (ndime, press, sigma) !----------------------------------------------------------------------- ! ! Adds the pressure the the deviatoric stresses ! ! ndime --> number of dimensions ! press --> pressure ! sigma --> Cauchy stress tensor ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION sigma (3, 3) ! ! Obtains the Cauchy stress from deviatoric and pressure ! components using equation (5.51 & 4.49) ! DO id = 1, ndime sigma (id, id) = sigma (id, id) + press END DO END SUBROUTINE addpres !----------------------------------------------------------------------- SUBROUTINE getheta (ngaus, vol0, vinc, theta) !----------------------------------------------------------------------- ! ! Obtains the volume ratio between the current and initial volumes ! ! ngaus --> number of Gauss points ! vol0 --> initial volume ! vinc --> determinant of Jacobian at each Gauss point ! theta --> volume ratio ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION vinc ( ngaus ) ! ! First finds the total volume ! evol = 0.d0 DO ig = 1, ngaus evol = evol + vinc (ig) END DO ! ! Then finds the volume ratio ! theta = evol / vol0 END SUBROUTINE getheta !----------------------------------------------------------------------- SUBROUTINE cisotp (ndime, xlamb, xmu, ctens) !----------------------------------------------------------------------- ! ! Determines the spatially isotropic tensor c(i,j,k,l). Valid ! for compressible neo-Hookean materials and all isotropic ! materials at the initial configuation. ! ! ndime --> number of dimensions ! xlamb --> effective lambda coefficient ! xmu --> effective mu coefficient ! ctens --> fourth order tensor c(i,j,k,l) ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION ctens (3, 3, 3, 3), delta (3, 3) DATA ( (delta (i, j), i = 1, 3), j = 1, 3) & / 1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0 / ! ! Evalutes c by using equation (5.37) ! DO i = 1, ndime DO j = 1, ndime DO k = 1, ndime DO l = 1, ndime ctens (i, j, k, l) = xlamb * delta (i, j) * delta (k, l) & + xmu * delta (i, k) * delta (j, l) & + xmu * delta (i, l) * delta (j, k) END DO END DO END DO END DO END SUBROUTINE cisotp !----------------------------------------------------------------------- SUBROUTINE cdevia (ndime, xmu, detf, btens, ctens) !----------------------------------------------------------------------- ! ! Determines the elasticity tensor for nearly incompressible ! neo-Hokean materials. ! ! ! ndime --> number of dimensions ! xmu --> effective mu coefficient ! detf --> determinant of F ! btens --> left Cauchy-Green tensor ! ctens --> fourth order tensor c(i,j,k,l) ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION ctens (3, 3, 3, 3), delta (3, 3), btens (3, 3) DATA ( (delta (i, j), i = 1, 3), j = 1, 3) & / 1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0 / ! ! First finds the trace of b and the effective mu ! xm = xmu / (3 * detf** (5.d0 / 3.d0) ) trace = 0.d0 DO id = 1, 3 trace = trace + btens (id, id) END DO ! ! Evalutes c by using equation (5.55a) ! DO i = 1, ndime DO j = 1, ndime DO k = 1, ndime DO l = 1, ndime ctens (i, j, k, l) = (2 * xm * trace / 3.d0) & * delta (i, j) * delta (k, l) & + xm * trace * delta (i, k) * delta (j, l) & + xm * trace * delta (i, l) * delta (j, k) & - 2 * xm * btens (i, j) * delta (k, l) & - 2 * xm * delta (i, j) * btens (k, l) END DO END DO END DO END DO END SUBROUTINE cdevia !----------------------------------------------------------------------- SUBROUTINE cvolum (ndime, press, ctens) !----------------------------------------------------------------------- ! ! Adds the volumetric component to the elasicity tensor for ! nearly incompressible materials ! ! ! ndime --> number of dimensions ! press --> pressure ! ctens --> fourth order tensor c(i,j,k,l) ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION ctens (3, 3, 3, 3), delta (3, 3) DATA ( (delta (i, j), i = 1, 3), j = 1, 3) & / 1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0 / ! ! Finds cvol using equation (5.55b) ! DO i = 1, ndime DO j = 1, ndime DO k = 1, ndime DO l = 1, ndime ctens (i, j, k, l) = ctens (i, j, k, l) & + press * delta (i, j) * delta (k, l) & - press * delta (i, k) * delta (j, l) & - press * delta (i, l) * delta (j, k) END DO END DO END DO END DO END SUBROUTINE cvolum !----------------------------------------------------------------------- SUBROUTINE cprinc (ndime, xlamb, xmu, stret, princ, sprin, ctens) !----------------------------------------------------------------------- ! ! Determines the constitutitve tangent modulus for materials ! defined in principal directions. ! ! ndime --> number of dimensions ! xlamb --> effective lambda coefficient ! xmu --> effective mu coefficient ! stret --> stretches ! princ --> principal directions ! sprin --> principal stresses ! ctens --> fourth order tensor c(i,j,k,l) ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION ctens (3, 3, 3, 3), stret (3), princ (3, 3), sprin (3) ! ! Evalutes c by using equation (5.56), (5.57) and (5.61) ! DO k = 1, ndime DO l = 1, ndime DO m = 1, ndime DO n = 1, ndime ! ! First the two first terms in equation (5.56) ! sum = 0.d0 DO i = 1, ndime DO j = 1, ndime sum = sum + xlamb * princ (k, i) * princ (l, i) & * princ (m, j) * princ (n, j) END DO ! j sum = sum + 2 * (xmu - sprin (i) ) * princ (k, i) & * princ (l, i) * princ (m, i) * princ (n, i) END DO ! i ! ! Then the last term. For the case where the stretches are equal ! equation (5.57) combined with equation (5.61) is used. ! DO i = 1, ndime DO j = 1, ndime !b IF ( j == i ) cycle ! to next j IF ( j == i ) goto 10 sti = stret (i) **2 stj = stret (j) **2 IF ( abs (sti - stj) >= 1.d-5 ) then val = (sprin (i) * stj - sprin (j) * sti) & / (sti - stj) ELSE val = xmu - sprin (j) END IF sum = sum + val * princ (k, i) * princ (l, j) & * (princ (m, i) * princ (n, j) & + princ (m, j) * princ (n, i) ) !b END DO ! j 10 END DO END DO ctens (k, l, m, n) = sum END DO END DO END DO END DO END SUBROUTINE cprinc !----------------------------------------------------------------------- SUBROUTINE internal (ndime, nnode, lnods, ldgof, sigma, gradn, & vinc, resid, react) !----------------------------------------------------------------------- ! ! Obtains and assembles the internal equivalent forces ! ! ndime --> number of dimensions ! nnode --> number of nodes ! lnods --> element connectivities ! ldgof --> nodal degrees of freedom ! sigma --> Cauchy stresses ! gradn --> shape function gradient ! vinc --> elemental volume ! resid --> residual forces ! react --> reactions ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION lnods (nnode), ldgof (ndime, * ), sigma (3, 3), react ( & * ), gradn (ndime, nnode), resid ( * ) ! ! Loops over nodes and degrees of freedom ! DO in = 1, nnode ip = lnods (in) DO id = 1, ndime if = ldgof (id, ip) ! ! Evaluates the component id at node ip of the internal force using ! equation (7.15b) ! sum = 0.d0 DO jd = 1, ndime sum = sum + sigma (id, jd) * gradn (jd, in) * vinc END DO ! ! If the corresponding degree of freedom is active, assembles ! the force. Otherwise acumulates the value squared as a reaction ! IF ( if > 0 ) then resid (if) = resid (if) - sum ELSE react ( - if) = react ( - if) + sum END IF END DO END DO END SUBROUTINE internal !----------------------------------------------------------------------- SUBROUTINE kconst (ndime, nnode, lnods, ldgof, ctens, gradn, vinc,& kprof, stifd, stifp) !----------------------------------------------------------------------- ! ! Obtains and assembles the constitutive component of the ! stiffness matrix ! ! ndime --> number of dimensions ! nnode --> number of nodes ! lnods --> element connectivities ! ldgof --> nodal degrees of freedom ! ctens --> fourth order constitutive tensor ! gradn --> shape function gradient ! vinc --> elemental volume ! kprof --> profile addresses ! stifd --> diagonal stiffness matrix ! stifp --> profile component of K ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION lnods (nnode), ldgof (ndime, * ), ctens (3, 3, 3, 3), & gradn (ndime, nnode), kprof ( * ), stifd ( * ), stifp ( * ) ! ! Loops over nodes and dimensions (rows) ! DO in = 1, nnode ip = lnods (in) DO id = 1, ndime if = ldgof (id, ip) !b IF ( if <= 0 ) cycle ! to next id IF ( if <= 0 ) goto 10 ! ! Loops again over nodes and dimensions (columns) ! DO jn = in, nnode jp = lnods (jn) DO jd = 1, ndime !b IF ( (in == jn ) .and. (jd < id ) ) cycle ! next jd IF ( (in == jn ) .and. (jd < id ) ) goto 20 jf = ldgof (jd, jp) !b IF ( jf <= 0 ) cycle ! next jd IF ( jf <= 0 ) goto 20 ! ! Uses equation (7.42) to evaluate the constitutive component ! of K in an indicial form ! sum = 0.d0 DO kd = 1, ndime DO ld = 1, ndime sum = sum + gradn (kd, in) * ctens (id, kd, jd, ld) & * gradn (ld, jn) * vinc END DO END DO ! ! Assembles sum in the right place. First case: diagonal term. ! IF ( jf == if ) then stifd (if) = stifd (if) + sum ELSE ! ! Then the off-diagonal terms into the profile. ! jc = max (jf, if) ir = min (jf, if) lp = kprof (jc) - jc + ir + 1 stifp (lp) = stifp (lp) + sum END IF 20 END DO END DO 10 END DO END DO END SUBROUTINE kconst !----------------------------------------------------------------------- SUBROUTINE ksigma (ndime, nnode, lnods, ldgof, sigma, gradn, vinc,& kprof, stifd, stifp) !----------------------------------------------------------------------- ! ! Obtains and assembles the initial component of the stiffness matrix ! ! ndime --> number of dimensions ! nnode --> number of nodes ! lnods --> element connectivities ! ldgof --> nodal degrees of freedom ! sigma --> Cauchy stresses ! gradn --> shape function gradient ! vinc --> elemental volume ! kprof --> profile addresses ! stifd --> diagonal stiffness matrix ! stifp --> profile component of K ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION lnods (nnode), ldgof (ndime, * ), sigma (3, 3), gradn ( & ndime, nnode), kprof ( * ), stifd ( * ), stifp ( * ) ! ! Loops over nodes ! DO in = 1, nnode ip = lnods (in) DO jn = in, nnode jp = lnods (jn) ! ! Uses equation (7.45) to evaluate the inital stress component ! of K in an indicial form ! sum = 0.d0 DO kd = 1, ndime DO ld = 1, ndime sum = sum + gradn (kd, in) * sigma (kd, ld) & * gradn (ld, jn) * vinc END DO END DO ! ! Assembles sum in the right places ! DO id = 1, ndime if = ldgof (id, ip) jf = ldgof (id, jp) !b IF ( ( if <= 0 ) .or. (jf <= 0 ) ) cycle ! next id IF ( ( if <= 0 ) .or. (jf <= 0 ) ) goto 10 ! ! First case: diagonal term. ! IF ( jf == if ) then stifd (if) = stifd (if) + sum ELSE ! ! Then the off-diagonal terms into the profile. ! jc = max (jf, if) ir = min (jf, if) lp = kprof (jc) - jc + ir + 1 stifp (lp) = stifp (lp) + sum END IF !b END DO 10 END DO END DO END DO END SUBROUTINE ksigma !----------------------------------------------------------------------- SUBROUTINE kvolume (ndime, nnode, ngaus, xkapp, vinc, elecd, & elacd, lnods, ldgof, kprof, stifd, stifp) !----------------------------------------------------------------------- ! ! Obtains and assembles the volumetric component of the stiffness ! matrix needed for nearly incompressible materials. ! ! ndime --> number of dimensions ! nnode --> number of nodes ! ngaus --> number of gauss poins ! xkapp --> effective kappa value ! vinc --> jacobian determinant per gauss point ! elecd --> Cartesian derivatives per Gauss point ! elacd --> average Cartesian derivatives ! lnods --> element connectivities ! ldgof --> nodal degrees of freedom ! kprof --> profile addresses ! stifd --> diagonal stiffness matrix ! stifp --> profile component of K ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION lnods (nnode), ldgof (ndime, * ), vinc ( * ), elacd ( & ndime, nnode), elecd (ndime, nnode, *), kprof ( *), stifd ( *),& stifp ( * ) ! ! First initilises the average cartesian derivaives ! DO in = 1, nnode DO id = 1, ndime elacd (id, in) = 0.d0 END DO END DO ! ! Loops over Gauss points and finds the current volume and average ! Cartesian derivatives ! evol = 0.d0 DO ig = 1, ngaus evol = evol + vinc (ig) DO in = 1, nnode DO id = 1, ndime elacd (id, in) = elacd (id, in) + elecd (id, in, ig) & * vinc (ig) END DO END DO END DO ! ! Constructs the volumetric tangent matrix. first loops over rows. ! DO in = 1, nnode ip = lnods (in) DO 10 id = 1, ndime !b drop 10 if = ldgof (id, ip) IF ( if <= 0 ) goto 10 !b IF ( if <= 0 ) cycle ! to next id ! ! Loops over columns ! DO jn = in, nnode jp = lnods (jn) DO 20 jd = 1, ndime !b drop 20 IF ( (in == jn) .and. (jd < id) ) goto 20 !b IF ( (in == jn) .and. (jd < id) ) cycle ! to next jd jf = ldgof (jd, jp) IF ( jf <= 0 ) goto 20 !b IF ( jf <= 0 ) cycle ! to next jd ! ! Finds the value to be assembled using equation (7.60). ! First case: diagonal term. ! sum = xkapp * elacd (id, in) * elacd (jd, jn) / evol IF ( jf == if ) then stifd (if) = stifd (if) + sum ELSE ! ! Then the off-diagonal terms into the profile. ! jc = max (jf, if) ir = min (jf, if) lp = kprof (jc) - jc + ir + 1 stifp (lp) = stifp (lp) + sum END IF 20 END DO !b drop 20 END DO 10 END DO !b drop 10 END DO END SUBROUTINE kvolume !----------------------------------------------------------------------- SUBROUTINE datri (al, au, ad, jp, neq, flg, jfile) !----------------------------------------------------------------------- ! !.... triangular decomposition of a matrix stored in profile form ! !.... input parameters ! al(jp(neq)) - lower triangular part of matrix ! au(jp(neq)) - upper part of triangular matrix ! ad(neq) - diagonals of triangular matrix ! jp(neq) - pointers to bottom of colums of al and au arrays ! neq - number of equations to be solved ! flg - if true equations are unsymmetric ! if false equations are symmetric and calling ! address of al may be same as that for au ! (i.e., au and al share same memory) ! jfile - unit number for printed output of warning ! messages. !.... output parameters ! al(jp(neq)) - lower triangular factor of matrix ! au(jp(neq)) - upper triangular factor of matrix ! ad(neq) - inverse of diagonal matrix in triangular factor ! !----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (a - h, o - z) LOGICAL flg DIMENSION al ( * ), au ( * ), ad ( * ), jp ( * ) DATA tol / 1.d-7 / ! !.... n.b. tol should be set to approximate half-word precision. !.... loop through the columns to perform the triangular decomposition ! jd = 1 DO j = 1, neq jr = jd+1 jd = jp (j) jh = jd-jr IF ( jh > 0 ) then is = j - jh ie = j - 1 ! !.... if diagonal is 0.0 compute a norm for singularity test ! IF ( ad (j) == 0.d0 ) CALL datest (au (jr), jh, daval) DO i = is, ie jr = jr + 1 id = jp (i) ih = min0 (id-jp (i - 1), i - is + 1) IF ( ih > 0 ) then jrh = jr - ih idh = id-ih + 1 au (jr) = au (jr) - dot (au (jrh), al (idh), ih) IF ( flg ) al (jr) = al (jr) & - dot (al (jrh), au (idh), ih) END IF END DO END IF ! !.... reduce the diagonal ! IF ( jh >= 0 ) then dd = ad (j) jr = jd-jh jrh = j - jh - 1 CALL dredu (al (jr), au (jr), ad (jrh), jh + 1, flg, ad (j) ) ! !.... check for possible errors and print warnings ! IF ( abs (ad (j)) < tol * abs (dd) ) write (jfile, 2000) j 2000 FORMAT('datri: Solver warning: loss of at least 7 ', & & ' digits in reducing diagonal: ',i5) IF ( dd < 0.d0 .and. ad (j) > 0.d0 ) write (jfile, 2001) j IF ( dd > 0.d0 .and. ad (j) < 0.d0 ) write (jfile, 2001) j 2001 FORMAT('datri: Solver warning: sign of diagonal ', & & 'changed when reducing equation: ',i5) IF ( ad (j) == 0.d0 ) write (jfile, 2002) j 2002 FORMAT('datri: Solver warning: reduced diagonal ', & & 'is zero for equation:',i5) ! !.... complete rank test for a 0.0 diagonal case ! IF ( dd == 0.d0 .and. jh > 0 ) then IF ( abs (ad (j) ) < tol * daval) write (jfile, 2003) j 2003 FORMAT('datri: Solver warning: rank failure for ', & & 'zero unreduced diagonal in equation:',i5) END IF END IF ! !.... store reciprocal of diagonal ! IF ( ad (j) /= 0.d0 ) ad (j) = 1.d0 / ad (j) END DO ! END SUBROUTINE datri !----------------------------------------------------------------------- SUBROUTINE dasol (al, au, ad, b, u, jp, neq, jfile, energy) !----------------------------------------------------------------------- ! !.... solution of equations stored in profile form: a * u = b !.... coefficient matrix must be decomposed into its triangular !.... factors using datri before using dasol. ! !.... input parameters ! al(jp(neq)) - lower triangular factor of matrix ! au(jp(neq)) - upper triangular factor of matrix ! (au and al have same calling address for ! symmetric matrices) ! ad(neq) - reciprocal of diagonal of triangular factor ! b(neq) - right hand side vector in equations ! jp(neq) - pointer array to bottom of columns of al and au ! neq - number of equations to be solved ! jfile - unit number for printed warning message. ! energy - energy norm for equations (rhs * soln) ! !.... output parameter ! u(neq) - solution of equations ! !----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION al ( * ), au ( * ), ad ( * ), b ( * ), jp ( * ), & u ( * ) ! !.... find the first non zero entry in the right hand side ! DO is = 1, neq u (is) = b (is) IF ( b (is) /= 0.d0 ) goto 200 END DO WRITE (jfile, 2000) 2000 FORMAT('dasol: Solver warning: zero right-hand-side vector') RETURN ! !.... reduce the right hand side ! 200 IF ( is < neq ) then DO j = is + 1, neq jr = jp (j - 1) jh = jp (j) - jr u (j) = b (j) IF ( jh > 0 ) then u (j) = u (j) - dot (al (jr + 1), u (j - jh), jh) END IF END DO END IF ! !.... multiply by inverse of diagonal elements ! energy = 0.d0 DO j = is, neq bd = u (j) u (j) = u (j) * ad (j) energy = energy + bd * u (j) END DO ! !.... backsubstitution ! IF ( neq > 1 ) then j = neq DO WHILE (j > 1) jr = jp (j - 1) jh = jp (j) - jr IF ( jh > 0 ) then CALL colred (au (jr + 1), u (j), jh, u (j - jh) ) END IF j = j - 1 END DO END IF ! END SUBROUTINE dasol !----------------------------------------------------------------------- SUBROUTINE datest (au, jh, daval) !----------------------------------------------------------------------- ! !.... test for rank ! !.... inputs ! ! au(j) - column of unreduced elements in array ! jh - number of elements in column ! !.... outputs ! ! daval - sum of absolute values ! !----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION au (jh) ! daval = 0.d0 DO j = 1, jh daval = daval + abs (au (j) ) END DO END SUBROUTINE datest !----------------------------------------------------------------------- FUNCTION dot (a, b, n) !b Replace with DOT_PRODUCT intrinsic !----------------------------------------------------------------------- ! !.... vector dot product ! !.... input parameters ! a(n),b(n) - two vectors ! n - length of vectors !.... ouput parameter ! dot - dot product of a and b ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION a ( * ), b ( * ) dot = 0.d0 DO i = 1, n dot = dot + a (i) * b (i) END DO END FUNCTION dot ! !----------------------------------------------------------------------- SUBROUTINE dredu (al, au, ad, jh, flg, dj) !----------------------------------------------------------------------- ! !.... reduce diagonal element in triangular decomposition ! !.... input parameters ! au(jh) - column of upper triangular part of matrix ! al(jh) - row of lower triangular part of matrix ! ad(jh) - reciprocal of diagonals in triangular factors ! jh - number of terms in vectors ! flg - if .true. equations are unsymmetric ! if .false. equations are symmetric ! dj - diagonal in matrix to be factored !.... output parameter ! dj - reduced diagonal of factor ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) LOGICAL flg DIMENSION al (jh), au (jh), ad (jh) ! DO j = 1, jh ud = au (j) * ad (j) dj = dj - al (j) * ud au (j) = ud END DO ! !.... finish computation of column of al for unsymmetric matrices ! IF ( flg ) then DO j = 1, jh al (j) = al (j) * ad (j) END DO END IF END SUBROUTINE dredu !----------------------------------------------------------------------- SUBROUTINE colred (au, xj, nn, b) !----------------------------------------------------------------------- ! !.... columnwise reduction for backsubstitution ! !.... input parameters ! au(nn) - column of upper triangular factor ! xj - solution of j-th equation ! nn - number of terms in current column !.... output parameter ! b(nn) - vector being modified to solution vector ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION au ( * ), b ( * ) ! DO n = 1, nn b (n) = b (n) - au (n) * xj END DO END SUBROUTINE colred !----------------------------------------------------------------------- SUBROUTINE force (ndgof, dlamb, eload, tload, resid) !----------------------------------------------------------------------- ! ! Obtains the current value of the external forces ! ! ndgof --> number of degrees of freedom ! resid --> residual vector ! eload --> nominal point load vector ! tload --> totla load vector ! dlamb --> load factor increment ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION resid ( * ), eload ( * ), tload ( * ) ! ! Evaluates loads and new residuals. Sets tload equal to eload ! DO id = 1, ndgof tload (id) = eload (id) resid (id) = resid (id) + dlamb * eload (id) END DO END SUBROUTINE force !----------------------------------------------------------------------- SUBROUTINE prescx (ndime, npoin, ldgof, pdisp, x0, x, xlamb) !----------------------------------------------------------------------- ! ! Imposes the current presicribed dispacements ! ! npoin --> number of points ! ldgof --> nodal degrees of freedom ! pdisp --> prescribed displacements ! x0 --> initial coordinates ! x --> current coordinates ! xlamb --> current load factor ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION ldgof (ndime, npoin ), pdisp ( * ), x (ndime, npoin ), & x0 ( ndime, npoin ) ! ! Interpolates current positions given lambda ! DO ip = 1, npoin DO id = 1, ndime ig = ldgof (id, ip) IF ( ig < 0 ) x (id, ip) = x0 (id, ip) + xlamb * pdisp ( - ig) END DO END DO END SUBROUTINE prescx !----------------------------------------------------------------------- SUBROUTINE update (ndime, npoin, ldgof, x, displ, eta) !----------------------------------------------------------------------- ! ! Updates the current positions according to the displacement ! vector displ and the arc length scaling factor eta ! ! ndime --> number of dimensions ! npoin --> number of points ! ldgof --> nodal degrees of freedom ! x --> nodal coordinates ! displ --> displacement vector ! resid --> residual force ! cload --> current loads ! eta --> scaling factor ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION ldgof (ndime, npoin ), x (ndime, npoin ), displ ( * ) ! ! Updates the positions ! DO ip = 1, npoin DO id = 1, ndime if = ldgof (id, ip) IF ( if > 0 ) x (id, ip) = x (id, ip) + eta * displ (if) END DO END DO END SUBROUTINE update !----------------------------------------------------------------------- SUBROUTINE arclen (ndgof, niter, arcln, displ, dispf, xincr, & xlamb, dlamb) !----------------------------------------------------------------------- ! ! Implements the arc-length method to obtain the load factor ! ! ndgof --> number of degrees of freedom ! niter --> number of iterations ! arcln --> arc length ! displ --> Newton-Raphson displacement vector (u) ! dispf --> load component of the displacement vector ! xincr --> total displacement over the load increment ! xlamb --> load factor ! dlamb --> load factor increment ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION displ ( * ), dispf ( * ), xincr ( * ) ! ! First obtains all the dot products required ! ufuf = 0.d0 ufdx = 0.d0 ufur = 0.d0 urur = 0.d0 urdx = 0.d0 dxdx = 0.d0 DO id = 1, ndgof ufuf = ufuf + dispf (id) * dispf (id) urur = urur + displ (id) * displ (id) ufur = ufur + dispf (id) * displ (id) ! ! For the first iteration initialises the total displacement to 0.0 ! IF ( niter == 1 ) then xincr (id) = 0.d0 ELSE ufdx = ufdx + dispf (id) * xincr (id) urdx = urdx + displ (id) * xincr (id) dxdx = dxdx + xincr (id) * xincr (id) END IF END DO ! ! Obtains the coefficients for equation (7.60) ! a1 = ufuf a2 = 2 * (ufdx + ufur) a3 = 2 * urdx + urur IF ( niter == 1 ) a3 = a3 - arcln * arcln * ndgof ! ! solves equation (7.60a) ! discr = sqrt (a2 * a2 - 4 * a1 * a3) IF ( a2 < 0.d0 ) discr = - discr discr = - (a2 + discr) / 2.d0 gamm1 = discr / a1 gamm2 = a3 / discr ! ! Determines which of the two solutions is better. For the first ! iteration chooses the smallest in absolute value ! gamma = gamm1 IF ( niter == 1 ) then IF ( abs (gamm2) < abs (gamm1) ) gamma = gamm2 ELSE cos1 = urdx + gamm1 * ufdx cos2 = urdx + gamm2 * ufdx IF ( cos2 > cos1 ) gamma = gamm2 END IF ! ! Updates the displacements and the accumulated displcements ! DO id = 1, ndgof displ (id) = displ (id) + gamma * dispf (id) xincr (id) = xincr (id) + displ (id) END DO ! ! Updates the values of the increment of lambda and lambda ! dlamb = dlamb + gamma xlamb = xlamb + gamma END SUBROUTINE arclen !----------------------------------------------------------------------- SUBROUTINE checkr (incrm, niter, ndgof, negdf, xlamb, resid, & eload, react, rnorm) !----------------------------------------------------------------------- ! ! Obtains the residual norm given the residual vector, the current ! vector of loads and the accumulated reactions ! ! incrm --> load increment ! niter --> Newton-Raphson iteration ! ndgof --> number of degrees of freedom ! negdf --> number of negative degrees of freedom ! xlamb --> current load parameter value ! resid --> residual force ! eload --> nominal loads ! react --> accumulated reactions ! rnorm --> residual norm ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION resid ( * ), eload ( * ), react ( * ) ! ! Evaluates the residual norm ! rnorm = 0.d0 fnorm = 0.d0 enorm = 0.d0 DO id = 1, ndgof rnorm = rnorm + resid (id) * resid (id) fnorm = fnorm + eload (id) * eload (id) END DO fnorm = fnorm * xlamb * xlamb ! ! Finds the reaction component ! DO id = 1, negdf enorm = enorm + react (id) * react (id) END DO rnorm = sqrt (rnorm / (fnorm + enorm) ) WRITE (6, 100) incrm, niter, rnorm 100 FORMAT(' Step: ',i3,' iteration: ',i3,' residual: ',g10.3) ! END SUBROUTINE checkr !----------------------------------------------------------------------- SUBROUTINE search (ndgof, resid, displ, eta0, eta, rtu0, rtu) !----------------------------------------------------------------------- ! ! Implements the line search algorithm ! ! ndgof --> number of degrees of freedom ! resid --> residual force ! displ --> displacement ! eta0 --> previous value of the parameter eta ! eta --> displacement factor eta ! rtu0 --> initial dot product of r by u ! rtu --> current dot product of r by u ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION resid ( * ), displ ( * ) ! ! First evaluates the dot product of r by u ! rtu = 0.d0 DO id = 1, ndgof rtu = rtu + resid (id) * displ (id) END DO IF ( rtu == 0.d0 ) return rtu1 = (rtu - rtu0 * (1 - eta) ) / (eta * eta) alpha = rtu0 / rtu1 eta0 = eta ! ! Takes different action depending on the value of alpha. For ! negative alpha obtains the root between (0,1). Otherwise finds ! the position of the minimum, provided that this is smaller than 1. ! IF ( alpha < 0.d0 ) then q = (alpha - sqrt (alpha * (alpha - 4) ) ) / 2.d0 eta = alpha / q ELSEIF ( alpha < 2.d0 ) then eta = alpha / 2.d0 ELSE eta = 1.d0 rtu = 0.d0 END IF END SUBROUTINE search !----------------------------------------------------------------------- SUBROUTINE output (ndime, nnode, ngaus, nstrs, npoin, nelem, eltyp, & title, icode, nincr, xlamb, x, lnods, ldgof, matno, & stres, eload, react) !----------------------------------------------------------------------- ! ! Writes formatted output ! ! ndime --> number of dimensions ! nnode --> number of nodes per element ! ngaus --> number of Gauss points per element ! nstrs --> number of stress components ! npoin --> number of points ! nelem --> number of elements ! eltyp --> element type used ! title --> problem title ! icode --> boundary code ! nincr --> increment or step number ! xlamb --> current load parameter ! x --> current nodal coordinates ! lnods --> element connectivities ! ldgof --> degree of freeedom numbers per node ! matno --> material number ! stres --> Gauss point stresses ! eload --> nominal loads ! react --> reactions ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) CHARACTER (10) eltyp CHARACTER (80) title DIMENSION x (ndime, * ), lnods (nnode, *), stres (nstrs, ngaus, * ), & icode ( * ), matno ( * ), eload ( * ), react ( * ), force (3), & ldgof (ndime, *) ! ! Writes the step number into the title, prints the title and the ! element type ! WRITE (title (41:80), 100) nincr, xlamb 100 FORMAT(' at increment: ',i4,', load: ',g10.3) WRITE (2, '(a)') title WRITE (2, '(a)') eltyp ! ! Writes the number of nodes and nodal coordinates ! WRITE (2, '(i5)') npoin DO ip = 1, npoin DO id = 1, ndime im = ldgof (id, ip) IF ( im > 0 ) force (id) = xlamb * eload (im) IF ( im < 0 ) force (id) = react ( - im) END DO WRITE (2, '(2i4,6e12.4)') ip, icode (ip), (x (id, ip), & id = 1, ndime) , (force (id) , id = 1, ndime) END DO ! ! Writes the number of elements and nodal connectivities ! WRITE (2, '(i5)') nelem DO ie = 1, nelem WRITE (2, '(16i5)') ie, matno (ie), & (lnods (in, ie) , in = 1, nnode) END DO ! ! Writes in the stresses per Gauss point ! DO ie = 1, nelem DO ig = 1, ngaus WRITE (2, '(6g15.5)') (stres (is, ig, ie) , is = 1, nstrs) END DO END DO ! END SUBROUTINE output !----------------------------------------------------------------------- SUBROUTINE 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) !----------------------------------------------------------------------- ! ! Dumps onto file all information required to restart the program ! ! title --> program title ! eltyp --> element type ! ndime --> number of dimensions ! npoin --> number of points ! nnode --> number of nodes per element ! ngaus --> number of Gauss points per element ! nstrs --> number of stresses per Gauss point ! nelem --> number of elements ! ndgof --> number of degrees of freedom ! negdf --> number of fixed dof ! nprs --> number of prescribed coordinates ! nprof --> number of entries in the profile of K ! nmats --> number of materials ! incrm --> increment number ! xlamb --> current load factor ! nbpel --> number of boundary elements ! nnodb --> number of nodes per boundary element ! ngaub --> number of gauss points per boundary element ! matyp --> material types ! props --> material properties ! matno --> material number of each element ! lnods --> element connectivities ! x --> nodal coordinates ! x0 --> initial nodal coordinates ! kprof --> profile addresses ! stifd --> diagonal stiffness ! stifp --> profile part of stiffness ! resid --> residual force ! eload --> nodal forces ! ldgof --> degree of freedom numbers ! icode --> nodal boundary codes ! eledb --> element database ! pdisp --> prescribed coordinates ! vol0 --> initial element volumes ! elbdb --> boundary elements information ! lbnod --> boundary elements nodal connectivities ! press --> external applied pressures ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) CHARACTER (80) title CHARACTER (10) eltyp ! DIMENSION x (ndime, npoin), x0 (ndime, npoin), icode (npoin), & ldgof (ndime, npoin) DIMENSION lnods (nnode, nelem), eledb (ndime+1, nnode+1, ngaus), & matno (nelem), vol0 (nelem), elbdb (ndime, nnodb + 1, ngaub), & lbnod (nnodb, nbpel), press (nbpel) DIMENSION stifd (ndgof), stifp (nprof), kprof (ndgof), eload ( & ndgof), pdisp (negdf), resid (ndgof) DIMENSION props (8, nmats), matyp (nmats) ! ! Rewinds the file and writes scalar variables ! REWIND 3 WRITE (3, err = 10) title, eltyp, ndime, npoin, nnode, ngaus, & nstrs, nelem, ndgof, negdf, nprs, nprof, nmats, incrm, xlamb, & nbpel, nnodb, ngaub ! ! Writes arrays ! WRITE (3, err = 10) matyp, props, eledb WRITE (3, err = 10) x WRITE (3, err = 10) x0 WRITE (3, err = 10) icode WRITE (3, err = 10) ldgof WRITE (3, err = 10) lnods WRITE (3, err = 10) matno WRITE (3, err = 10) stifd WRITE (3, err = 10) stifp WRITE (3, err = 10) kprof WRITE (3, err = 10) eload WRITE (3, err = 10) resid WRITE (3, err = 10) pdisp WRITE (3, err = 10) vol0 WRITE (3, err = 10) elbdb IF ( nbpel > 0 ) then WRITE (3, err = 10) lbnod WRITE (3, err = 10) press END IF RETURN ! 10 WRITE (6, 100) 100 FORMAT(' Error writing restart file') STOP 'dump: Error writing restart file' END SUBROUTINE dump !----------------------------------------------------------------------- SUBROUTINE restar1 (title, eltyp, ndime, npoin, nnode, ngaus, nstrs, & nelem, ndgof, negdf, nprs, nprof, nmats, incrm, & xlamb, nbpel, nnodb, ngaub) !----------------------------------------------------------------------- ! ! reads from file all scalar information required to restart ! the program from the last converged position ! ! title --> program title ! eltyp --> element type ! ndime --> number of dimensions ! npoin --> number of points ! nnode --> number of nodes per element ! ngaus --> number of Gauss points per element ! nstrs --> number of stresses per Gauss point ! nelem --> number of elements ! ndgof --> number of degrees of freedom ! negdf --> number of fixed degrees of freedom ! nprs --> number of prescribed coordinates ! nprof --> number of entries in the profile of K ! nmats --> number of materials ! incrm --> increment number ! xlamb --> current load factor ! nbpel --> number of boundary elements ! nnodb --> number of nodes per boundary element ! ngaub --> number of gauss points per boundary element ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) CHARACTER (80) title CHARACTER (10) eltyp ! ! Rewinds the file and reads scalar variables ! REWIND 3 READ (3, err = 10) title, eltyp, ndime, npoin, nnode, ngaus, & nstrs, nelem, ndgof, negdf, nprs, nprof, nmats, incrm, xlamb, & nbpel, nnodb, ngaub RETURN ! 10 WRITE (6, 100) 100 FORMAT('restar1: Error reading restart file') STOP 'restar1: Error reading restart file' END SUBROUTINE restar1 !----------------------------------------------------------------------- SUBROUTINE 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) !----------------------------------------------------------------------- ! ! Reads from file all vector information required to restart ! the program from the last converged position ! ! ndime --> number of dimensions ! npoin --> number of points ! nnode --> number of nodes per element ! ngaus --> number of Gauss points per element ! nelem --> number of elements ! ndgof --> number of degrees of freedom ! negdf --> number of fixed degrees of freedom ! nprof --> number of entries in the profile of K ! nmats --> number of materials ! nbpel --> number of boundary elements ! nnodb --> number of nodes per boundary element ! ngaub --> number of gauss points per boundary element ! matyp --> material types ! props --> material properties ! matno --> material number of each element ! lnods --> element connectivities ! x --> nodal coordinates ! x0 --> initial nodal coordinates ! kprof --> profile addresses ! stifd --> diagonal stiffness ! stifp --> profile part of stiffness ! resid --> residual force ! eload --> nodal forces ! ldgof --> degree of freedom numbers ! icode --> nodal boundary codes ! eledb --> element database ! pdisp --> prescribed coordinates ! vol0 --> initial volume of each element ! elbdb --> boundary elements information ! lbnod --> boundary elements nodal connectivities ! press --> external applied pressures ! !----------------------------------------------------------------------- ! IMPLICIT DOUBLE PRECISION (a - h, o - z) DIMENSION x (ndime, npoin), x0 (ndime, npoin), icode (npoin), & ldgof (ndime, npoin) DIMENSION lnods (nnode, nelem), eledb (ndime+1, nnode+1, ngaus), & matno (nelem), vol0 (nelem), elbdb (ndime, nnodb + 1, ngaub), & lbnod (nnodb, nbpel), press (nbpel) DIMENSION stifd (ndgof), stifp (nprof), kprof (ndgof), & eload (ndgof), pdisp (negdf), resid (ndgof) DIMENSION props (8, nmats), matyp (nmats) ! ! Reads arrays ! READ (3, err = 10) matyp, props, eledb READ (3, err = 10) x READ (3, err = 10) x0 READ (3, err = 10) icode READ (3, err = 10) ldgof READ (3, err = 10) lnods READ (3, err = 10) matno READ (3, err = 10) stifd READ (3, err = 10) stifp READ (3, err = 10) kprof READ (3, err = 10) eload READ (3, err = 10) resid READ (3, err = 10) pdisp READ (3, err = 10) vol0 READ (3, err = 10) elbdb IF ( nbpel > 0 ) then READ (3, err = 10) lbnod READ (3, err = 10) press END IF RETURN ! 10 WRITE (6, 100) 100 FORMAT('restar2: Error reading restart file') STOP 'restar2: Error reading restart file' END SUBROUTINE restar2