. 46
( 54 .)


initel ............finds the initial tangent matrix and gravity loads
gradsh ........finds the Cartesian derivatives of the shape functions
kvolume .......finds the mean dilatation component of the stiffness
cisotp ........finds the isotropic elasticity tensor
kconst ........obtains and assembles the constitutive component of K
initrk ............initializes the tangent matrix and residual to zero
bpress ............deals with boundary pressure elements
pforce ........evaluates the forces due to normal pressure
kpress ........evaluates the external pressure component of K
elemtk ............organises the evaluation of internal forces and K
gradsh ........finds the Cartesian derivatives of the shape functions
leftcg ........evaluates the left Cauchy--Green tensor
jacobi ........finds principal directions and stretches
stress1 .......evaluates stresses for material type 1
stress3 .......evaluates stresses for material type 3
stress5 .......evaluates stresses for material type 5
stress6 .......evaluates stresses for material type 6
addpres .......adds the internal pressure to deviatoric stresses
getjbar .......evaluates the average Jacobian or volume ratio
cisotp ........finds the isotropic elasticity tensor
cdevia ........finds the deviatoric component of the elasticity tensor
cvolum ........finds the pressure component of the elasticity tensor
cprinc ........finds the elasticity tensor for materials in principal
internal ......obtains and assembles internal forces
kconst ........obtains and assembles the constitutive component of K
ksigma ........obtains and assembles the initial stress component of K
kvolume .......finds the mean dilatation component of the stiffness
datri .............performs the LDU decomposition
dot ...........dot product of two vectors
datest ........tests the rank of the matrix
dredu .........reduces the diagonal terms
dasol .............performs the backsubstitution solution process
colred ........reduces a column
dot ...........dot product of two vectors
force .............evaluates the current force increment
prescx ............enforces prescribed displacements
update ............updates the nodal coordinates
arclen ............implements the arc-length method
checkr ............evaluates the residual norm
search ............implements the line search process
output ............outputs current converged state
dump ..............dumps the current state to a binary file
restar1 ...........prepares to restart by reading scalar values
restar2 ...........prepares to restart by reading arrays


In order to give an overview of the structure of the program FLagSHyP,
Box 8.3 lists all routines contained in the program indented according to
the calling sequence. The routines in italic typeface are described in details
in the following sections. The remainder either are routines common to
standard linear elasticity ¬nite element codes or are relatively minor and
self-evident utility routines that are a direct FORTRAN translation of some
of the equations described in the text.

BOX 8.5: Solution algorithm
r INPUT geometry, material properties, and solution parameters (seg-
ment 2)
r INITIALIZE F = 0, x = X (initial geometry), R = 0 (segment 2)
r FIND initial K (segment 3)
r LOOP over load increments (segment 5)
¯ ¯
r SET » = » + ∆», F = »F, R = R ’ ∆»F (segment 5)
r UPDATE GEOMETRY (segment 5)
r FIND F, T, K, R = T ’ F (segment 5)
r END IF (segment 5)
r DO WHILE ( R / F > tolerance ) (segment 6)
r SOLVE Ku = ’R (segment 6)
r IF ARC-LENGTH FIND » (segment 6)
r UPDATE GEOMETRY x = x + u (segment 7)
r FIND F, T, K, R = T ’ »F (segment 7)
r FIND · (segment 7)
r END LOOP (segment 7)
r END DO (segment 7)
r ENDLOOP (segment 8)

8.8 MAIN ROUTINE flagshyp

The main routine closely follows the algorithm described in Box 7.1, which
is repeated in Box 8.4 for convenience. The items in parentheses refer to
segments contained in the main routine flagshyp. Each of these program
segments is presented in the following pages in a box followed by a short de-
scription. The concatenation of these boxes comprises the main routine. To
avoid duplication, comments have been removed and only actual FORTRAN
instructions are shown. A dictionary of main variables used throughout the
code is provided as an appendix to this chapter.

¬‚agshyp segment 1 “ Dimensions
program flagshyp
implicit double precision (a-h,o-z)
parameter (mpoin=100)
parameter (melem=100)
parameter (mdgof=300)
parameter (mprof=10000)
parameter (mnode=10)
parameter (mgaus=8)
parameter (mmats=10)
parameter (mbpel=100)
character*80 title
character*10 eltyp
logical rest
dimension x(3,mpoin),x0(3,mpoin),icode(mpoin),
& ldgof(3,mpoin)
dimension lnods(mnode,melem),eledb(4,mnode+1,mgaus),
& stres(6,mgaus,melem),matno(melem),
& elecd(4,mnode,mgaus),vinc(mgaus),vol0(melem),
& elacd(4,mnode),elbdb(3,mnode+1,mgaus),
& lbnod(mnode,mbpel),press(mbpel)
dimension stifd(mdgof),stifp(mprof),kprof(mdgof),
& eload(mdgof),
& pdisp(mdgof),resid(mdgof),displ(mdgof),
& xincr(mdgof),
& react(mdgof),tload(mdgof)
dimension props(8,mmats),gravt(3),matyp(mmats)
8.8 MAIN ROUTINE flagshyp

The parameters listed above control the maximum problem size. Large
examples may require some of these dimensions to be reset. It can be seen
that these parameters control the size of various vectors and matrices whose
purpose will be described when required.

¬‚agshyp segment 2 “ Input

call welcome(title,rest)
if(.not.rest) then
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)
if(ndime.eq.3) nstrs=6
call nodecon(npoin,nelem,nnode,lnods,stifp)
call degfrm (mdgof,npoin,ndime,stifp,ndgof,negdf,
& ldgof,stifd)
call profile(mprof,ndgof,nelem,nnode,ndime,lnods,
& ldgof,nprof,kprof)
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)

The program can start either from an input data ¬le or, when convenient,
using data from a restart ¬le written during a previous incomplete analy-
sis. The mode of operation is controlled by the user interactively from the
screen and is recorded in the logical variable rest read by the welcome rou-
tine. In the following it is assumed that the data are read for the ¬rst time.
Subroutine elinfo reads the element type eltyp and establishes the arrays
eledb and elbdb discussed in Section 8.4. Subroutine innodes reads nodal
input data items 3 and 4, and subroutine inelems reads element input data
items 5 and 6. Routine nodecon determines the node to node connectivities
required by the Cuthill“McKee algorithm implemented in degfrm, which
at the same time allocates degree of freedom numbers to nodes and stores
them in ldgof. In profile the pro¬le pointers are obtained and stored

in the vector array kprof. Subroutine matprop reads in material property
data items 7 and 8. Routine inloads inputs loading and prescribed dis-
placement items 9, 10, 11, and 12, and incontr reads the solution control
item 13.

¬‚agshyp segment 3 “ Initialization and dump

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)
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)

Routine initno sets the residuals resid, the tangent matrix compo-
nents stifd and stifp to zero, and de¬nes the initial geometry x0 as the
current geometry x. Routine initel evaluates the equivalent nodal forces
due to gravity and the tangent sti¬ness matrix at the unstressed con¬gura-
tion (this coincides with the small strain linear elastic sti¬ness matrix). It
also evaluates the total mesh volume and prints this on the screen. All the
relevant information is now dumped to a restart ¬le.

¬‚agshyp segment 4 “ Restart

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)
8.8 MAIN ROUTINE flagshyp

call incontr(nincr,xlmax,dlamb,miter,cnorm,searc,
& arcln,5)

If the program is restarting from a previous incomplete solution, a restart
¬le is read and the solution control parameters are reset by reading them
interactively from the screen. This enables the user to reset the convergence
tolerance, the maximum number of iterations per step, etc.

¬‚agshyp segment 5 “ Increment loop

do while((xlamb.lt.xlmax).and.(incrm.lt.nincr))
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(nprs.gt.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)

The do while controls the load or prescribed displacement incremen-
tation. This remains active while the load scaling parameter xlamb is less
than the maximum xlmax and the increment number is smaller than the
total number of increments nincr.
The imposition of an increment of point or gravity loads immediately
creates a residual force. This is added to any small residual carried over from


. 46
( 54 .)