C :
C :    ******************************************************
C :    *****            COPYRIGHT  (C)  1996            *****
C :    *****    BY STRUCTURAL DYNAMICS RESEARCH CORP.   *****
C :    *****             ALL RIGHTS RESERVED            *****
C :    ******************************************************
C : 
C :
C : NBBPLSS.PRG
C :     Plasticity program file
C :
C :---------------------------------------------------------------------
C :
C : matrix and variable names:
C :
C : ALPHA    - back stress or kinematic hardening tensor (at element
C :            integration points)
C :
C : CONFAC   - element condensation factors
C :
C : DALPHA   - kinematic hardening tensor increment (at integration 
C :            points)
C : DELDI    - independent dof partition of the incremental displacement
C : DELDK    - kinematic partition of the incremental displacement
C : DELDSP   - incremental displacements in nodal space
C : DELEN    - strain energy increment
C : DEPSCR   - creep strain increments (at integration points)
C : DEPSPL   - plastic strains increments (at integration points)
C : DKAPPA   - isotropic hardening variable increment (for plasticity)
C : DISCRC   - displacements at contact reference configuration
C : DISP     - displacements in nodal space
C : DISPRV   - nodal displacements from previous solution point
C : DLAMDA   - increment in contact tractions for current load case
C : DLAMDO   - last iteration increment in contact tractions for current
C :            load case
C :
C : EBCDIS   - contact boolean for all COMPLETE contact elements
C : EBOCFU   - contact boolean for COMPLETE-ACTIVE contact elements
C : EBOLFU   - boolean from element space to full nodal space
C : ECSTAT   - contact element status flags
C : ECTOFF   - contact element offset list
C : EDSPRC   - element condensed displacements, previous solution point
C : EKOFF    - element stiffness offset table
C : ELAMDA   - contact tractions
C : ELCFOF   - element condensed freedom offset table
C : ELCFOR   - contact force in contact element space
C : ELCNST   - element normal contact stiffness
C : ELCTST   - element tangential contact stiffness
C : ELDLDR   - element incremental displacements on retained dof
C : ELDIPR   - element displacements at previous solution point
C : ELDISP   - element displacements
C : ELDLDC   - element incremental displacements on condensed dof
C : ELDFST   - element stress stiffness matrix
C : ELDSPC   - element condensed displacements
C : ELEDIS   - element displacements at previous solution point for 
C :            beams
C : ELIDIS   - element displacements at current iteration for beams
C : ELLD     - element internal forces at current iteration
C : ELLDPR   - element internal forces at previous solution point 
C : ELLHST   - element (material) stiffness matrix
C : ELNTMP   - element temperature matrix
C : ELNTRF   - element temperatures for warp (not implemented)
C : ELTPOF   - element temperature offset table
C : ELSTRN   - element strain for convergence checking
C : ELSTRR   - element stress resultants for stress stiffness
C : ELSTRS   - element stress for stress stiff and convergence checking
C : EOFF     - element offset table
C : EPSCRE   - creep strains (at element integration points)
C : EPSPLA   - plastic strains (at element integration points)
C :
C : FDHAT    - load vector, dependent partition
C : FEQC     - element internal forces at condensed dof
C : FEQTOT   - equivalent end forces for beams
C : FEXT     - external (applied) forces in nodal space
C : FEXTC    - element external forces at condensed dof
C : FEXTSV   - saved external forces (prior to adding contact forces)
C : FI       - load vector, independent partition
C : FK       - load vector, kinemetic partition
C : FINT     - internal force in nodal space
C : FP       - reaction force on prscribed dof
C : FPEXSV   - saved external forces on prescribed dof (prior to adding
C :            contact forces)
C : FPEXT    - external plus contact forces on prescribed dof
C : FPINT    - internal force at prescribed dof
C : FRESCK   - force residual for force convergence checking
C : FRESID   - residual force in nodal space, FEXT - FINT
C : FRESIC   - residual force on condensed dof, element space
C : FS       - reaction force on suppressed dof
C : FSEXSV   - saved external forces on suppressed dof (prior to adding
C :            contact forces)
C : FSEXT    - external plus contact forces on suppressed dof
C : FSINT    - internal force at suppressed dof
C :
C : KAPPA    - element isotropic hardening variable for plasticity at 
C :            integration points
C : KIK      - off diagonal partition of stiffness matrix associated
C :            with independent and kinematic dof
C : KIP      - off diagonal partition of stiffness matrix associated
C :            with independent and prescribed dof
C : KIS      - off diagonal partition of stiffness matrix associated
C :            with independent and suppressed dof
C : KKK      - diagonal partition of stiffness matrix associated with
C :            kinematic dof
C : KPK      - off diagonal partition of stiffness matrix associated
C :            with prescribed and kinematic dof
C : KPS      - off diagonal partition of stiffness matrix associated
C :            with prescribed and suppressed dof
C : KSK      - off diagonal partition of stiffness matrix associated
C :            with suppressed and kinematic dof
C : KPP      - diagonal partition of stiffness matrix associated with
C :            prescribed dof
C :
C : NBOLPR   - boolean maps from full nodal space to prescribed dof
C : NBOLSU   - boolean maps from full nodal space to suppressed dof
C : NDCFOR   - contact force in nodal space
C : NMULT1   - nodal multiplicity table, number of elements connected to
C :            each node
C : NMULT2   - nodal multiplicity table, packed integer of element label
C :            and type
C : NOOF     - nodal dof offset table
C :
C : OALPHA   - previous solution point back stress or kinematic 
C :            hardening tensor prior to incremental update
C : OEPSPL   - previous solution point plastic strains prior to 
C :            incremental updates
C : OLDSTN   - previous iteration strain for convergence checking
C : OLDSTS   - previous iteration stress for convergence checking
C :
C : PI       - independent multiplier partitions of the rigid body
C :            modes, orthogonalized with respect to the mass matrix
C : PIDUAL   - dual basis vectors, independent partition
C : PK       - kinematic partition of rigid body vectors
C : PKDUAL   - dual basis vectors, kinematic partition
C :
C : SCVOFF   - scalar value offset array (Gauss point)
C :
C : TEPSPL   - last solution point plus incremental plastic strain 
C : TNVOFF   - element tensor variable offset table (for material 
C :            history)
C : TOTEN    - total strain energy increment for this load step
C :
C : Y_CHST   - global flag to reform stiffness, 0 -> do not form,
C :            1 -> reform stiffness
C : Y_CIEQ   - equilibrium iteration counter
C : Y_CONFLG - contact flag, 0 -> not contact solve, 1 -> contact solve
C : Y_CVG    - convergence flag, 0 -> not converged, 1 -> converged,
C :            -1 -> diverging
C : Y_CURK   - current stiffness parameter for the step 
C : Y_CURL   - time increment for load control or load increment for arc
C :            length control of the current step
C : Y_CURSTF - current iteration stiffness parameter
C : Y_DELTT  - current time increment
C : Y_DIS    - =1, input parameter for displacement convergence check
C : Y_DMIN   - minimum value for denominator in nonlinear convergence 
C :            check set in nbbnls.prg
C : Y_DSPL   - flag to check displacement convergence, 0 -> off, 1 -> on
C : Y_ENG    - =5, input parameter for energy convergence check
C : Y_ENRG   - flag to check energy convergence, 0 -> off, 1 -> on
C : Y_FOR    - =2, input parameter for force convergence check
C : Y_FRAT   - Contact force convergence ratio
C : Y_FRMD   - "Stress Stiffening" flag from the "Stiffness Control"
C :            form, 0 -> do not calculate, 1 -> calculate
C : Y_GLBC   - global convergence counter, incremented by 1 for each 
C :            convergence criteria satisfied 
C : Y_GLBM   - number of convergence checks required
C : Y_GNL    - "Geometric Nonlinear" flag, 0 -> off, 1 -> on
C : Y_IAUSTP - adaptive search flag, 0 -> off, 1 -> on
C : Y_IDIVRG - divergence flag
C : Y_LATLIM - limit point flag, 0 -> normal load control, 1 -> below 
C :            limit point arc length on, 2 -> limit point exceeded arc
C :            lenghth on, -1 -> limit point reached arc length off 
C : Y_LUPK   - flag to reform stiffness based on stiffness change,
C :            0 -> off, 1 -> on
C : Y_MAXI   - user input maximum number of contact force iterations
C : Y_MIEQ   - "Maximum Number of Iterations" from "Add/Modify Time 
C :            Interval" form
C : Y_MODTYP - model type;
C :            0 -> 3D
C :            1 -> 2D plane stress/strain
C :            2 -> 2D axisymmetric
C : Y_NCTE   - number of potential contact elements
C : Y_NGAPS  - number of gap elements
C : Y_NITRN  - number of equilibrium iteration from program file 
C :            nbbgnls or nbbplss 
C : Y_NOCONV - not converged flag
C : Y_NSTN   - flag to check strain convergence, 0 -> off, 1 -> on
C : Y_NSTS   - flag to check stress convergence, 0 -> off, 1 -> on
C : Y_PLAS   - plasticity flag, 0 -> off, 1 -> on
C : Y_PMOD   - verification flag, 0 -> verification run, 1 -> not a
C :            verification run, -1 -> ????
C : Y_PRESTF - previous iteration stiffness parameter
C : Y_RATD   - displacement convergence ratio
C : Y_RATE   - energy convergence ratio
C : Y_RATF   - force convergence ratio
C : Y_RATS   - strain convergence ratio
C : Y_RATT   - stress convergence ratio
C : Y_RSCH   - flag to check force convergence, 0 -> off, 1 -> on
C : Y_STN    - =4, input parameter for strain convergence check 
C : Y_STS    - =3, input parameter for stress convergence check
C : Y_ERR1     = 0 no surface to surface contact elements on special 
C :                elements (axi, plane_stress, etc.)
C :           != 0 gap elements on special elements and surface contact
C :
C : Z_IRGD   - = 0 - no rigid body dof present,
C :           != 0 - rigid body dof present 
C :
C :---------------------------------------------------------------------
C :
K :   #ON_ERROR GOTO PERROR
C :
C :---------------------------------------------------------------------
C :
C : get the fraction of initially open active contacts 
C :
K :   /GLOBAL_PARAMETER GET "Z_MPNL"
K :   #Y_NOFAC = Z_MPNL/100.0
C :
C :---------------------------------------------------------------------
C :
C : get the plasticity flag 
C :
K :   /GLOBAL_PARAMETER GET "Y_PLAS"
C :
C :---------------------------------------------------------------------
C :
C :     Initialize equilibrium iteration counter, Y_CIEQ.
C :
K :   #Y_CIEQ = 0
K :   #Y_ONE = 1
C :
C :*********************************************************************
C :---------------------------------------------------------------------
C :
C :     Initialize the contact iteration counter, Y_CITI.
C :
K :   #Y_CITI = 0
K :   #Y_CONTOL = 0.20
C :
C :     If contact analysis, initialize gap traction increment,
C :     DLAMDA, previous solution point nodal displacements, DISPRV,
C :     and previous solution point condensed displacements, EDSPRC.
C :     Also get maximum number of contact force iterations, Z_MAXI.
C :
K :   #IF ( Y_CONFLG EQ 0 ) THEN GOTO ITERLP
C :
K :   /HMPACK MATRIX KILL  DLAMDA ELDDEL
K :   /HMPACK MATRIX INQUIRE  DISP "Z_EXIST"
K :   #IF ( Z_EXIST EQ 1 ) THEN /HMPACK MATRIX COPY  DISP DISPRV
K :   /HMPACK MATRIX INQUIRE  ELDSPC "Z_EXIST"
K :   #IF ( Z_EXIST EQ 1 ) THEN /HMPACK MATRIX COPY  ELDSPC EDSPRC
K :   /GLOBAL_PARAMETER GET "Z_MAXI"
K :   #IF ( Z_MAXI LT 1 ) THEN #Z_MAXI = 1
K :   #Y_MAXI = Z_MAXI
C :
C :---------------------------------------------------------------------
C :*********************************************************************
C :
C :     Begin equilibrium iteration / contact status loop.
C :
K : #ITERLP:
C :
C :     Initialize global convergence counter, Y_GLBC.
C :
K :   #Y_GLBC = 0
C :
C :     Increment iteration counter, Y_CIEQ, store in global common.
C :
K :   #Y_CIEQ = Y_CIEQ+1
K :   /GLOBAL_PARAMETER PUT  "Y_CIEQ"
C :
C :     If the maximum number of iterations, Y_MIEQ, has not been 
C :     exceeded, proceed with the next iteration.
C :
K :   #IF ( Y_CIEQ LE Y_MIEQ ) THEN GOTO NSTOP1
C :
C :     Maximum number of iterations, Y_MIEQ, exceeded, set the not
C :     converged flag, Y_NOCONV.
C :
K :   #Y_NOCONV = 1
C :
C :     If the adaptive search flag, Y_IAUSTP, is not set, issue an
C :     error message before exiting the iteration loop.
C :
K :   #IF ( Y_IAUSTP EQ 1 ) THEN GOTO QUIT
K :   #Y_ERNO = 17611
K :   #Y_NMAR = 1
K :   #Y_SVER = 3
K :   /LOG ERROR  "Y_ERNO" "Y_SVER" "Y_NMAR" "Y_MIEQ"
K :   #GOTO QUIT
C :
K : #NSTOP1:
C :
K :   #X_CIEQ = "Equilibrium Iteration Number:  "+Y_CIEQ
K :   /LOG MESSAGE " "
K :   /LOG MESSAGE  X_CIEQ
K :   /LOG MESSAGE " "
K :   /LOG MESSAGE  "----- Begin Iteration -----"
C :
C :     Initialize the number of contact status changes.
C :
K :   #Y_NCSC = 0
C :
C :*********************************************************************
C :---------------------------------------------------------------------
C :
C :     Check for the existence of contact and gap elements
C :
K :   #IF ( Y_CONFLG EQ 0 ) THEN GOTO GAP2
C :
C :     If this is a geometric linear solve, we only need to form 
C :     contact elements once. Check to see if this has been done yet
C :     by checking for the existence of the contact offset matrix,
C :     ECTOFF.
C :
K :   /HMPACK MATRIX INQUIRE  ECTOFF "Z_EXIST"
K :   #IF ( Z_EXIST EQ 1 AND Y_GNL EQ 0 ) THEN GOTO GAP1
C : 
C :     Form/update contact elements with new geometry, return Y_NCTE.
C :
K :   /FORM CONTACT ELEMENTS  "Y_NCTE" "Y_ERR1" NOOF DISP
K :   #IF ( Y_ERR1 NE 0 ) THEN GOTO QUIT
K :   #NCCR = "Number of potential contact elements :  "+Y_NCTE
K :   /LOG MESSAGE  NCCR
C :
C :     Update displacements for the contact definition configuration,
C :     DISCRC. This is required only for a geometric nonlinear solve.
C :
K :   /HMPACK MATRIX KILL  DISCRC
K :   #IF ( Y_GNL NE 0 ) THEN /HMPACK MATRIX COPY  DISP DISCRC
C :
C :     Add gap elements to the list of contact elements and return the
C :     number of gap elements, Y_NGAPS.
C :
K :   /FORM CONTACT GAP_ELEMENTS  "Y_NGAPS"
C : 
C :     Form contact element offset list, ECTOFF.
C :
K :   /FORM ELEMENT OFFSET CONTACT ECTOFF "Y_MODTYP"
K :   /LOG MESSAGE "Contact Element Offsets Formed"
C :
K :   #IF ( Y_GNL EQ 0 OR Y_PLAS EQ 1 ) THEN GOTO GAP1X
K :   #IF ( Y_MODTYP EQ 0) THEN GOTO GAP1X
C :
C :     If ECDISP exists, expand it to full size from current EBCDIS
C :     and call it ECDISF in case the EBCDIS boolean below changes
C :     from the previously 2D purely geometric nonlinear pass.
C :
K :   /HMPACK MATRIX INQUIRE ECDISP "Z_EXIST"
K :   #IF ( Z_EXIST EQ 0 ) THEN GOTO GAP1X
K :   /HMPACK MATRIX KILL  ECDISF
K :   /HMPACK BOOLEAN_RHS PARTITION  EBCDIS ECDISP ECDISF
K :   /HMPACK MATRIX KILL  ECDISP
C :
K : #GAP1X:
C :
C :     Form the boolean for all COMPLETE contact elements, EBCDIS.
C :
K :   /HMPACK MATRIX KILL  EBCDIS ECDUM
K :   /FORM ELEMENT BOOLEAN CONTACT_ALL  ECTOFF FREE NOOF EBOCIN EBCDIS ECDUM
K :   /HMPACK MATRIX KILL EBOCIN
C :
K :   #IF ( Y_GNL EQ 0 OR Y_PLAS EQ 1 ) THEN GOTO GAP1Y
K :   #IF ( Y_MODTYP EQ 0) THEN GOTO GAP1Y
C :
C :     If ECDISF exists, boolean it to create ECDISP of the correct size
C :     using the new EBCDIS for the current 2D purely geometric nonlinear pass.
C :
K :   /HMPACK MATRIX INQUIRE ECDISF "Z_EXIST"
K :   #IF ( Z_EXIST EQ 0 ) THEN GOTO GAP1Y
K :   /HMPACK MATRIX KILL  ECDISP
K :   /HMPACK BOOLEAN_RHS EXPANSION  EBCDIS ECDISF ECDISP
K :   /HMPACK MATRIX KILL  ECDISF
C :
K : #GAP1Y:
C :
C :     Form contact element stiffness matrices, normal, ELCNST, and
C :     tangential, ELCTST.
C :
K :   /FORM ELEMENT LHS_MATRICES CONTACT_STIFFNESS  ECTOFF EBCDIS FREE,
K :                                                 ELCNST ELCTST EKCOFF
K :   /LOG MESSAGE "Contact Stiffness Matrix Formation Complete"
C :
K : #GAP1:
C :
K :   #IF ( Y_GNL EQ 0 ) THEN GOTO GAP1A
C :
C :     Boolean contact element displacement increments, ECDISO, from
C :     previous solution point nodal displacement to last contact 
C :     element update, DISPRV - DISCRC
C :
K :   #IF ( Y_PLAS EQ 1 OR Y_MODTYP EQ 0 ) THEN /HMPACK MATRIX KILL  ECDISP
K :   /HMPACK MATRIX KILL  ECIDIS ECDISO DISDUM
C :
K :   /HMPACK MATRIX INQUIRE  DISPRV "Z_EXIST"
K :   #IF ( Z_EXIST EQ 0 ) THEN GOTO GAP1B
K :   /HMPACK MATRIX COPY  DISPRV DISDUM
K :   /HMPACK MATRIX INQUIRE  DISCRC "Z_EXIST"
K :   #IF (Z_EXIST EQ 1)THEN /HMPACK MATHEMATICS SUBTRACT  DISCRC DISDUM
K :   /HMPACK BOOLEAN_RHS EXPANSION  EBCDIS DISDUM ECDISO
K :   /HMPACK MATRIX KILL  DISDUM
C :
K :   #GOTO GAP1B
C :
K : #GAP1A:
C :
C :     Geometric linear contact solution:
C :     Boolean contact element displacements, ECDISP, from nodal
C :     displacements, DISP, using the full boolean of all COMPLETE
C :     contact elements. Also form the previous solution point
C :     contact element displacements, ECDISO.
C :
K :   /HMPACK MATRIX KILL  ECIDIS ECDISP
K :   /HMPACK MATRIX INQUIRE  DISP "Z_EXIST"
K :   #IF ( Z_EXIST EQ 1 ) THEN /HMPACK BOOLEAN_RHS EXPANSION  EBCDIS,
K :                              DISP ECDISP
K :   #IF ( Z_EXIST EQ 1 ) THEN /HMPACK MATRIX COPY  ECDISP ECIDIS
C :
K :   /HMPACK MATRIX KILL  ECDISO
K :   /HMPACK MATRIX INQUIRE  DISPRV "Z_EXIST"
K :   #IF ( Z_EXIST EQ 1 ) THEN /HMPACK BOOLEAN_RHS EXPANSION  EBCDIS,
K :                              DISPRV ECDISO
C :
K : #GAP1B:
C :
C :     Form / update contact element status flag, ECSTAT.
C :
K :   /LOG MESSAGE  "Form Contact Status"
C :
K :   /FORM CONTACT STATUS  ECTOFF ECSTAT ELAMDA DLAMDA ECDISP "Y_CIEQ",
K :                         "Y_NCS0" "Y_NCS1" "Y_NCS2" "Y_NCS3" "Y_NCSC",
K :                         "Y_NOFAC" "Y_SLP" "Y_MODTYP"
C :
K :   #NCSC = "Number of contact status changes:  "+Y_NCSC
K :   /LOG MESSAGE NCSC
K :   #NCS0 = "Number of inactive contacts:       "+Y_NCS0
K :   /LOG MESSAGE  NCS0
K :   #NCS1 = "Number of active open contacts:    "+Y_NCS1
K :   /LOG MESSAGE  NCS1
K :   #NCS2 = "Number of sticking contacts:       "+Y_NCS2
K :   /LOG MESSAGE  NCS2
K :   #NCS3 = "Number of sliding contacts:        "+Y_NCS3
K :   /LOG MESSAGE  NCS3
C :
C :     Note that for a contact solution to be converged, the nonlinear 
C :     convergence criteria must be satisfied and the the number of 
C :     status changes, Y_NCSC, must be zero.
C :
C :     Combine normal, ELCNST, and tangential, ELCTST, contact 
C :     stiffnesses to get total contact stiffness, ELCSTF.
C :
K :   #Y_ZERO = 0
K :   /FORM CONTACT COMBINED_STIFFNESS ECTOFF ECSTAT ELCNST ELCTST,
K :                                    ELCSTF "Y_ZERO" EKCOFF "Y_CIEQ"
C :
C :     Form active contact element boolean, EBOCFU. 
C :     ECSTAT controls the Boolean formation
C :
K :   /FORM ELEMENT BOOLEAN CONTACT_ALL ECTOFF FREE NOOF EBOCIN EBOCFU ECSTAT
K :   /HMPACK MATRIX KILL EBOCIN
C :
K : #GAP2:
C :
C :---------------------------------------------------------------------
C :*********************************************************************
C :
C :     Compute the element (material) stiffness matrix if required.
C :                  
K :   /LOG MESSAGE  "Begin Stiffness Matrix Formation"
C :
C :     For the first iteration of the solution point the tangent
C :     material stiffness is computed from the material state at the
C :     end of the last solution point using previous solution point
C :     plastic strain, OEPSPL, previous solution point back stress, 
C :     OALPHA, last computed plastic strain increment, DEPSPL, and
C :     last computed back stress increment, DALPHA. Also use previous
C :     material state if it is geometric nonlinear and it is the
C :     second iteration, since new plastic strain and back stress
C :     increments have not yet been computed.
C :
K :   #IF ( Y_CIEQ GT 1 ) THEN GOTO NOTFST
K :   /FO EL LHS_MATRICES STIFF  EOFF   ELLHST ELNTMP ELTPOF ELNTRF,
K :                              ELDISP TNVOFF OEPSPL EPSCRE OALPHA, 
K :                              DEPSPL DEPSCR DALPHA ELCFOF CONFAC,
K :                              ELDSPC "Y_LUPK" EKOFF "Y_NCSC"
K :   /HMPACK MATRIX KILL DALPHA DEPSPL DKAPPA
K :   #GOTO ENDSTF
C :
K : #NOTFST:
C :
C :     For subsequent iterations, use beginning of load step values 
C :     of plastic strain, EPSPLA, and kinematic hardening variable, 
C :     ALPHA, in the stiffness calculation.
C :
K :   /FO EL LHS_MATRICES STIFF  EOFF ELLHST ELNTMP ELTPOF ELNTRF,
K :                              ELDISP TNVOFF EPSPLA EPSCRE ALPHA,
K :                              DEPSPL DEPSCR DALPHA ELCFOF CONFAC,
K :                              ELDSPC "Y_LUPK" EKOFF "Y_NCSC"
C :
K : #ENDSTF:
C :
C :----------------------------------------------------------------------
C : 
C :   Form externally applied loads, FEXT, and partition
C :
K :   /HMPACK PROGRAM_FILE RUN
K :                                "nlhlod.prg"    
C :
C :     Save the external force vector, FEXT, in FEXTSV. Save prescribed
C :     and suppressed partitions of the external force vector, FPEXT 
C :     and FSEXT in FPEXSV and FSEXSV
C :
K :   /HMPACK MATRIX COPY  FEXT FEXTSV
K :   /HMPACK MATRIX COPY  FSEXT FSEXSV
K :   /HMPACK MATRIX COPY  FPEXT FPEXSV
C :
C :---------------------------------------------------------------------
C :
C :     If this is the first equilibrium iteration, do not copy element 
C :     stress, ELSTRS, and element strain, ELSTRN, into old element 
C :     stress, OLDSTS, and old element strain, OLDSTN. These are used 
C :     for convergence checking.
C :    
K :   #IF ( Y_CIEQ EQ 1 ) THEN GOTO NCNVG
K :   /HMPACK MATRIX NEGATE_COPY  ELSTRS OLDSTS
K :   /HMPACK MATRIX NEGATE_COPY  ELSTRN OLDSTN
C :
K : #NCNVG:
C :
C :---------------------------------------------------------------------
C : 
K :   /LOG MESSAGE  "Begin Internal Force Formation"
C :       
C :     Form / update the following matrices: ELSTRS, ELSTRN, ELSTRR, 
C :     ELLD, FEQC, ELIDIS.
C :       
K :   #IF ( Y_CIEQ EQ 1 ) THEN GOTO NOADD
C :
C :     Add plastic strain increment to plastic strain for use in the
C :     element internal force calculation.
C :  
K :   /HMPACK MATRIX COPY  EPSPLA TEPSPL
K :   /HMPACK MATHEMATICS ADD  DEPSPL TEPSPL
K :   /FO EL F IN  EOFF   ELCFOF TNVOFF ELTPOF ELNTMP,
K :                TEPSPL EPSCRE ELDISP ELSTRS ELSTRN,
K :                ELSTRR ELLD   ELDSPC CONFAC FEQC,
K :                FEQTOT ELDIPR ELLDPR ELEDIS ELIDIS
K :   /HMPACK MATRIX KILL  TEPSPL
K :   #GOTO ENDINF
C :
K : #NOADD:
C :
K :   /FO EL F IN  EOFF   ELCFOF TNVOFF ELTPOF ELNTMP,
K :                EPSPLA EPSCRE ELDISP ELSTRS ELSTRN,
K :                ELSTRR ELLD   ELDSPC CONFAC FEQC,
K :                FEQTOT ELDIPR ELLDPR ELEDIS ELIDIS
C :
K : #ENDINF:
C :       
K :   /LOG MESSAGE  "Internal Force Formation Complete"
C :
C :---------------------------------------------------------------------
C : 
C :     Get the global change stiffness flag, Y_CHST, to see if the
C :     stiffness is to be reformed this iteration
C :
K :   /GLOBAL_PARAMETER GET  "Y_CHST"
K :   #IF ( Y_CHST EQ 0 ) THEN GOTO NODCMP
C :
C :     Check the stress stiffening flag, Y_FRMD, equilibrium iteration 
C :     number, Y_CIEQ and contact force iteration number, Y_CITI.
C :
K :   #IF ( Y_FRMD NE 1 ) THEN GOTO NOSSTF
C :
K :   #IF ( Y_CIEQ EQ 1 AND Y_CITI LT 2 ) THEN GOTO SSTF1
C :
K :   /LOG MESSAGE  "Begin Stress Stiffness Matrix Formation"
C :
K :   /FO EL LS DIFF_STIFFNESS  EOFF   ELDISP FEQTOT ELDFST ELTPOF,
K :                             ELNTMP TNVOFF ELSTRS ELSTRR ELCFOF,
K :                             ELDSPC EKOFF
C :
K :   #GOTO SSTF2
C :
K : #SSTF1:
C :
K :   /HMPACK MATRIX INQUIRE  TLSTRS "Z_EXIST"
K :   #IF (Z_EXIST EQ 0) THEN GOTO SSTF2
C :
K :   /LOG MESSAGE  "Begin Stress Stiffness Matrix Formation"
C :
K :   /FO EL LS DIFF_STIFFNESS  EOFF   ELDISP FEQTOT ELDFST ELTPOF,
K :                             ELNTMP TNVOFF TLSTRS TLSTRR ELCFOF,
K :                             ELDSPC EKOFF
C :
K : #SSTF2:
C :
K :   /HMPACK MATRIX COPY  ELSTRS TLSTRS
K :   /HMPACK MATRIX COPY  ELSTRR TLSTRR
K :   /HMPACK MATRIX KILL  ELSTRR
C : 
C :     Sum the material and stress stiffness matrices (actually subtract
C :     the negative of the stress stiffness from the material stiffness)
C :
K :   /HMPACK MATHEMATICS SUBTRACT  ELDFST ELLHST
K :   /HMPACK MATRIX KILL  ELDFST
C :
K : #NOSSTF:
C :
C :     Delete the stress resultants, ELSTRR, which are no longer needed.
C :
K :   /HMPACK MATRIX KILL  ELSTRR
C :
C :---------------------------------------------------------------------
C : 
K :   /LOG MESSAGE  "Begin Stiffness Matrix Partitioning"
C :
C :     Assemble stiffness partitions from element stiffness.
C :
K :   /HMPACK MATRIX KILL  KKK KIP KIK KPK
K :   /HMPACK MATRIX KILL  KPP KPS KSK KIS
K :   /FORM SYSTEM_MATRIX STATIC_STIFFNESS  EKOFF NOOF FREE NMULT1 NMULT2,
K :                                         ELLHST UI KIP KIS KIK KPP KPS,
K :                                         KPK KSK KKK
C :
K :   /HMPACK MATRIX KILL  ELLHST
C :
C :*********************************************************************
C :---------------------------------------------------------------------
C :
C :     Assemble the contact stiffness if this is a contact analysis.
C :
K :   #IF ( Y_CONFLG EQ 0 ) THEN GOTO GAP4
C :
C :     Form contact stiffness partitions.
C :
K :   /FORM SYSTEM_MATRIX CONTACT_STIFFNESS ECTOFF EKCOFF NOOF FREE ECSTAT,
K :                                         ELCSTF UI KIPC KISC KPPC KPSC
K :   /HMPACK MATRIX KILL  ELCSTF
K :   /LOG MESSAGE  "Contact Stiffness Partitions Formed"
C :
K : #GAP4:
C :
C :---------------------------------------------------------------------
C :*********************************************************************
C :
K :   /LOG MESSAGE  "Stiffness Matrix Formation Complete"
C :
C :----------------------------------------------------------------------
C :
C :     Setup data for sparse solver.
C :
K :   /FORM SPARSE_STATIC  UI
C : 
C :     Decompose the stiffness matrix.
C : 
K :   /SELECT DECOMPOSITION AUTO_RESTRAINT
K :   /SELECT DECOMPOSITION STATISTICS
K :   /HMPACK MATRIX KILL  UKIK UICI UG UICIUG
K :   /HMPACK MATRIX KILL  UGV KPLUS
C :
K :   /HMPACK PROGRAM_FILE RUN
K :                                "utility.prg"
C :
K :   /GLOBAL_PARAMETER GET  "Y_PMOD"
C :
C :     Check for a negative stiffness, Y_PMOD = -1
C :
K :   #IF ( Y_PMOD NE -1 ) THEN GOTO SKNGK
K :   #Y_PMOD = 1
K :   /GLOBAL_PARAMETER PUT  "Y_PMOD"
C :
C :     Set the flag for limit point exceeded, Y_LATLIM = -1, and get
C :     out of the iteration loop (this program file)
C :
K :   #Y_LATLIM = -1
K :   #GOTO QUIT
C :
K : #SKNGK:
C :
K : #NODCMP:
C :                                                             
C :*********************************************************************
C :---------------------------------------------------------------------
C :
C :     Initialize contact inner (force) loop counter, Y_CITI
C :
K :   #Y_CITI = 0
C :
K :   #IF ( Y_CONFLG EQ 0 ) THEN GOTO GAP3
C :
C :     Begin inner contact loop, initialize traction force increment 
C :     matrix, DLAMDA, and inner loop convergence flag, Y_CNVI.
C :
K :   #Y_CNVI = 0
C :
K : #DOGAPI:
C :
K :   #Y_CITI = Y_CITI+1
K :   #IF ( Y_CITI GT 1 ) THEN GOTO GAP3
K :   #X_CITI = "Begin Contact Force Calculation"
K :   /LOG MESSAGE  X_CITI
C :
K : #GAP3:
C :
C :---------------------------------------------------------------------
C :*********************************************************************
C :
C :     If this is a not a contact solve, or if it is the first 
C :     iteration of a contact solve, the internal forces have already 
C :     been calculated.
C :
K :   #IF ( Y_CONFLG EQ 0 ) THEN GOTO ENDINF2
K :   #IF ( Y_CITI EQ 1 ) THEN GOTO ENDINF2
C :
C :     Form / update the following matrices: ELSTRS, ELSTRN, ELSTRR, 
C :     ELLD, FEQC, ELIDIS.
C :       
C :     Add plastic strain increment to plastic strain for use in the
C :     element internal force calculation.
C :  
K :   /HMPACK MATRIX COPY  EPSPLA TEPSPL
K :   /HMPACK MATHEMATICS ADD  DEPSPL TEPSPL
K :   /FO EL F IN  EOFF   ELCFOF TNVOFF ELTPOF ELNTMP,
K :                TEPSPL EPSCRE ELDISP ELSTRS ELSTRN,
K :                ELSTRR ELLD   ELDSPC CONFAC FEQC,
K :                FEQTOT ELDIPR ELLDPR ELEDIS ELIDIS
K :   /HMPACK MATRIX KILL  TEPSPL
C :
K : #ENDINF2:
C :
C :     Form internal force partitions for independent, FINT, 
C :     prescribed, FPINT, and suppressed, FSINT, dof.
C :
K :   /HMPACK MATRIX KILL  FINT FPINT FSINT
K :   /HMPACK BOOLEAN_RHS PARTITION  EBOLFU ELLD FINT
K :   /HMPACK BOOLEAN_RHS PARTITION  NBOLPR FINT FPINT
K :   /HMPACK BOOLEAN_RHS PARTITION  NBOLSU FINT FSINT
C :
C :*********************************************************************
C :---------------------------------------------------------------------
C :
K :   #IF ( Y_CONFLG EQ 0 ) THEN GOTO GAP5
C :
C :     Form contact Lagrange multipliers and contact element forces.
C :     If converged, also update the total tractions (ELAMDA).
C :
C :     Set Y_CONV = 0, (to not update ELAMDA for the next time point)
C :     this update will be done after convergence.
C :
K :   #Y_CONV = 0
C :
K :   /FORM ELEMENT STATIC CONTACT_FORCE  ECTOFF ECSTAT ECDISP ECDISO, 
K :                            ELAMDA DLAMDA ELCFOR "Y_CONV" "Y_CIEQ",
K :                            "Y_MODTYP" ECIDIS
C :
C :     Boolean contact forces from element space to nodal space.
C :
K :   /HMPACK MATRIX KILL  NDCFOR
K :   /HMPACK BOOLEAN_RHS PARTITION  EBOCFU ELCFOR NDCFOR
K :   /HMPACK MATRIX KILL  ELCFOR FEXT FSEXT FPEXT
K :   /HMPACK MATRIX COPY  FEXTSV FEXT
K :   /HMPACK MATRIX COPY  FSEXSV FSEXT
K :   /HMPACK MATRIX COPY  FPEXSV FPEXT
C :
C :     Partition contact force to prescribed and suppressed dof.
C :
K :   /HMPACK BOOLEAN_RHS PARTITION  NBOLPR NDCFOR FPEXT
K :   /HMPACK BOOLEAN_RHS PARTITION  NBOLSU NDCFOR FSEXT
C :
C :     Add contact forces into applied forces.
C :
K :   /HMPACK MATHEMATICS ADD  NDCFOR FEXT
C :
K : #GAP5:
C :
C :---------------------------------------------------------------------
C :*********************************************************************
C : 
C :     Compute residual force and residual force partitions. Create
C :     FI, FK, FDHAT, FRESID, FRESIC.
C :                 
K :   /HMPACK PROGRAM_FILE RUN 
K :                                "nlslod.prg" 
C :
C :----------------------------------------------------------------------
C :
C :     Form rigid body basis.
C :
K :   #Y_FI = "FI"
K :   #Y_FK = "FK"
K :   #IF ( Z_IRGD EQ 0 ) THEN GOTO NRGD
K :   /HMPACK PROGRAM_FILE RUN 
K :                             "rigids.prg"
C :
K :   #IF ( Y_CONFLG NE 0 ) THEN GOTO NOMESS1
K :   /LOG MESSAGE  "Rigid Body Modes Constructed"
C :
K : #NOMESS1:
C :
C :----------------------------------------------------------------------
C : 
C :     Form inertia loads.
C :                            
K :   /HMPACK PROGRAM_FILE RUN
K :                                "inertia.prg"
C :
K :   #IF ( Y_CONFLG NE 0 ) THEN GOTO NRGD
K :   /LOG MESSAGE  "Inertia Loads Constructed"
C :
K : #NRGD:
C :
C :----------------------------------------------------------------------
C : 
C :     Compute displacement increments.
C :                                  
K :   #Y_DI = "DELDI"
K :   #Y_DK = "DELDK"
C :
K :   /HMPACK PROGRAM_FILE RUN
K :                                "power.prg"
C :                                                 
C :----------------------------------------------------------------------
C :
C :     Filter out rigid body modes.
C :
K :   #IF ( Z_IRGD EQ 0 ) THEN GOTO NRGD3
K :   /HMPACK PROGRAM_FILE RUN
K :                                "filter.prg"
C :
K :   #IF ( Y_CONFLG NE 0 ) THEN GOTO NOMESS2
K :   /LOG MESSAGE  "Displacements Filtered for Rigid Body Components"
C :
K : #NOMESS2:
C :
K :   /HMPACK MATRIX KILL  PI PK
K :   /HMPACK MATRIX KILL  PIDUAL PKDUAL
C :
K : #NRGD3:
C :
C :----------------------------------------------------------------------
C :                              
C :     Form nodal dispalcement increment, DELDSP, element displacement
C :     increment, ELDLDR, and element displacement, ELDISP. Update
C :     displacement partitions DELDI and DELDK, nodal displacements,
C :     DISP, and reaction force partitions FS and FP.
C :
K :   /HMPACK PROGRAM_FILE RUN
K :                                "nlinc.prg"
C :
C :----------------------------------------------------------------------
C : 
C :     Compute incremental condensed dof, ELDLDC and update total
C :     displacements on condensed dof, ELDSPC.
C :                                                      
K :   /HMPACK PROGRAM_FILE RUN
K :                                "nlcinc.prg"
C :
C :----------------------------------------------------------------------
C :
K :   #IF ( Y_CONFLG EQ 0 ) THEN GOTO SKIPIT1
K :   #Y_RAT = 1.0
K :   #IF ( Y_GNL EQ 0 ) THEN GOTO SKIPIT1
C :
C :     For geometric nonlinear analysis with contact, form the Euclidean
C :     norm of the iterative displacement and the Euclidean norm of the
C :     incremental displacement. Find the square root of the ratio (for 
C :     deciding if contact force iterations should be done).
C :
K :   /HMPACK MATRIX INQUIRE  ELDDEL "Z_EXIST"
K :   #IF ( Z_EXIST EQ 0 ) THEN /HMPACK MATRIX COPY  ELDLDR ELDDEL
K :   #IF ( Z_EXIST NE 0 ) THEN /HMPACK MATHEMATICS ADD  ELDLDR ELDDEL
K :   /HMPACK MATRIX COPY  ELDLDR TEMP1
K :   /HMPACK MATHEMATICS TRANSPOSE_MULTIPLY  ELDLDR TEMP1 TEMP2
K :   /HMPACK MATRIX EXTRACT EXTRACT_VALUE  TEMP2 "Y_ONE" "Y_ONE" "Y_NUM"
K :   /HMPACK MATRIX KILL  TEMP1 TEMP2
C :
K :   /HMPACK MATRIX COPY  ELDDEL TEMP1
K :   /HMPACK MATHEMATICS TRANSPOSE_MULTIPLY  ELDDEL TEMP1 TEMP2
K :   /HMPACK MATRIX EXTRACT EXTRACT_VALUE  TEMP2 "Y_ONE" "Y_ONE" "Y_DEN"
K :   /HMPACK MATRIX KILL  TEMP1 TEMP2
K :   #Y_NUM = SQRT(Y_NUM)
K :   #Y_DEN = SQRT(Y_DEN)
K :   #IF ( Y_DEN GT 1E-10 ) THEN #Y_RAT = Y_NUM/Y_DEN
C :
K :   #DELETE  Y_NUM Y_DEN
C :
K : #SKIPIT1:
C :
C :----------------------------------------------------------------------
C :
K :   #IF ( Y_CONFLG NE 0 ) THEN GOTO NOMESS3
K :   /LOG MESSAGE  "Displacement Calculation Complete"
C :
K : #NOMESS3:
C :                                                 
C :*********************************************************************
C :---------------------------------------------------------------------
C :
K :   #IF ( Y_CONFLG EQ 0 ) THEN GOTO GAP6
C :
C :     Compute contact element displacements, ECDISP.
C :
K :   #IF ( Y_GNL EQ 0 ) THEN GOTO GAP6A
C :
C :     Geometric nonlinear contact solution:
C :     Boolean contact element displacement increments, ECDISP, from
C :     nodal displacement increment since last contact element update,
C :     ECIDIS = (DISP - DISCRC) and add to current ECDISP (only
C :     if purely geometric nonlinear). 
C :
K :   #IF ( Y_PLAS EQ 1 OR Y_MODTYP EQ 0 ) THEN /HMPACK MATRIX KILL  ECDISP
K :   /HMPACK MATRIX KILL  DISDUM ECIDIS
K :   /HMPACK MATRIX COPY  DISP DISDUM
K :   /HMPACK MATRIX INQUIRE  DISCRC "Z_EXIST"
K :   #IF (Z_EXIST EQ 1)THEN /HMPACK MATHEMATICS SUBTRACT  DISCRC DISDUM
K :   /HMPACK BOOLEAN_RHS EXPANSION  EBCDIS DISDUM ECDISP
K :   /HMPACK BOOLEAN_RHS EXPANSION  EBCDIS DISDUM ECIDIS
K :   /HMPACK MATRIX KILL  DISDUM
C :
K :   #GOTO GAP6
C :
K : #GAP6A:
C :
C :     Geometric linear contact solution:
C :     Boolean contact element displacements, ECDISP, from nodal
C :     displacements, DISP.
C :
K :   /HMPACK MATRIX KILL  ECDISP ECIDIS
K :   /HMPACK BOOLEAN_RHS EXPANSION  EBCDIS DISP ECDISP
K :   /HMPACK MATRIX COPY  ECDISP ECIDIS
C : 
K : #GAP6:
C :
C :---------------------------------------------------------------------
C :*********************************************************************
C :
C :---------------------------------------------------------------------
C : 
C :     Compute plastic strain increment if required. For stability
C :     do not calculate plastic strains on the first equilibrium
C :     iteration of a geometric nonlinear or contact solution.
C :
K :   #IF ( Y_MTNL EQ 0 ) THEN GOTO NWLDST
K :   #IF ( Y_GNL NE 0 AND Y_CIEQ LT 2 ) THEN GOTO NWLDST 
K :   #IF ( Y_CONFLG NE 0 AND Y_CIEQ LT 2 ) THEN GOTO NWLDST
C : 
K :   #IF ( Y_CONFLG NE 0 ) THEN GOTO NOMESS4
K :   /LOG MESSAGE  "Begin Plastic Strain Increment Calculation"
C :
K : #NOMESS4:
C :
K :   /FO EL PI  EOFF ELTPOF TNVOFF ELDISP ELNTMP EPSPLA EPSCRE ALPHA,
K :              DEPSPL DALPHA SCVOFF KAPPA DKAPPA ELCFOF ELDSPC FEQTOT,
K :              ELDIPR ELLDPR ELEDIS
C :
K :   #IF ( Y_CONFLG NE 0 ) THEN GOTO NWLDST
K :   /LOG MESSAGE  "End Plastic Strain Increment Calculation"
C :
K : #NWLDST:
C :
C :*********************************************************************
C :---------------------------------------------------------------------
C :
C :     End of contact iteration loop, set contact convergence flag,
C :     Y_CNVI and decide if another contact iteration is required.
C :                                           
K :   #IF ( Y_CONFLG EQ 0 ) THEN GOTO GAP8
C :
K :   #IF ( Y_NCS2 EQ 0 AND Y_NCS3 EQ 0 ) THEN GOTO GAP8
C :
K :   /HMPACK MATRIX KILL  TLAMDA
K :   /HMPACK MATRIX COPY  ELAMDA TLAMDA
K :   /HMPACK MATHEMATICS ADD  DLAMDA TLAMDA
C :
K :   #IF ( Y_CITI LE 1 ) THEN GOTO GAP7
C :
K :   #Y_INNR = 1
K :   #Y_CMFX = 0.0D0
K :   /INSPECT CONTACT_CONVERGENCE  TLAMDA TLAMDO "Y_INNR" "Y_CMFX",
K :                                 "Y_FRAT" "Y_CNVI" "Y_MODTYP"
C :
K : #GAP7:
C :
K :   /HMPACK MATRIX COPY  TLAMDA TLAMDO
K :   /HMPACK MATRIX KILL  TLAMDA
C :
C :     If geoemtric changes are not too large and the user specified 
C :     number of maximum inner loop iterations is not exceeded, go to 
C :     the next inner contact loop.
C :
K :   #IF ( Y_NCS2+Y_NCS3 EQ 0 ) THEN GOTO GAP8
K :   #IF ( Y_GNL EQ 1 AND Y_RAT GT Y_CONTOL ) THEN GOTO GAP8
K :   #IF ( Y_CITI LT Y_MAXI ) THEN GOTO DOGAPI
K :   #IF ( Y_CITI LT 2 ) THEN GOTO GAP8
K :   #IF ( Y_CNVI EQ 1 ) THEN GOTO INRCNV
C :
C :   inner loop failed to converge in user-specifed no. of iterations
C :
K :   #Y_ERNO = 21646
K :   #Y_NMAR = 2
K :   #Y_SVER = 2
K :   /LOG ERROR "Y_ERNO" "Y_SVER" "Y_NMAR" "Y_CIEQ" "Y_CITI"
C :
K :   #INRCNV:
K :   #X_FRAT = "Contact force convergence ratio:  "+Y_FRAT
K :   /LOG MESSAGE  X_FRAT
C :
K : #GAP8:
C :
C :---------------------------------------------------------------------
C :*********************************************************************
C :
C :     Check strain convergence if required.
C :
K :   #IF ( Y_NSTN EQ 0 OR Y_CIEQ EQ 1 ) THEN GOTO NSTRN
C :
C :     Compute strain increment, OLDSTN is the negative of the strain
C :     from the last equilibrium iteration.
C : 
K :   /HMPACK MATHEMATICS ADD  ELSTRN OLDSTN
C :
K :   /INSPECT NL_VECT_CONVERGENCE  "Y_STN" "Y_CIEQ" OLDSTN ELSTRN,
K :                                 "Y_CVG" "Y_RATS" "Y_DMIN" "Y_IAUSTP"
K :   #IF ( Y_CVG LT 0 ) THEN GOTO DIVRGE
K :   #IF ( Y_CVG EQ 1 ) THEN #Y_GLBC = Y_GLBC+1
C :
K : #NSTRN:
C :
C :----------------------------------------------------------------------
C :
C :     Check stress convergence if required.
C :
K :   #IF ( Y_NSTS EQ 0 OR Y_CIEQ EQ 1 ) THEN GOTO NSTRS
C :
C :     Compute stress increment.
C :
K :   /HMPACK MATHEMATICS ADD  ELSTRS OLDSTS
C :
K :   /INSPECT NL_VECT_CONVERGENCE  "Y_STS" "Y_CIEQ" OLDSTS ELSTRS,
K :                                 "Y_CVG" "Y_RATT" "Y_DMIN" "Y_IAUSTP"
K :   #IF ( Y_CVG LT 0 ) THEN GOTO DIVRGE
K :   #IF ( Y_CVG EQ 1 ) THEN #Y_GLBC = Y_GLBC+1
C :
K : #NSTRS:
C :
C :----------------------------------------------------------------------
C :
C :     Check force convergence if required.
C :
K :   #IF ( Y_RSCH EQ 0 OR Y_CIEQ EQ 1 ) THEN GOTO NRESID 
K :   /INSPECT NL_VECT_CONVERGENCE  "Y_FOR" "Y_CIEQ" FRESCK FEXT,
K :                                 "Y_CVG" "Y_RATF" "Y_DMIN" "Y_IAUSTP"
K :   #IF ( Y_CVG LT 0 ) THEN GOTO DIVRGE
K :   #IF ( Y_CVG EQ 1 ) THEN #Y_GLBC = Y_GLBC+1
C :
K : #NRESID:
C :
C :----------------------------------------------------------------------
C : 
C :     Check displacement convergence if required, skip on first 
C :     equilibrium iteration.
C :
K :   #IF ( Y_DSPL EQ 0 ) THEN GOTO NDISP
K :   #IF ( Y_CIEQ EQ 1 ) THEN GOTO NDISP2 
K :   /HMPACK MATRIX COPY  DISP DISDUM
K :   /HMPACK MATHEMATICS SUBTRACT  DISPLS DISDUM
K :   /INSPECT NL_VECT_CONVERGENCE  "Y_DIS" "Y_CIEQ" DISDUM DISP,
K :                                 "Y_CVG" "Y_RATD" "Y_DMIN" "Y_IAUSTP"
K :   #IF ( Y_CVG LT 0 ) THEN GOTO DIVRGE
K :   #IF ( Y_CVG EQ 1 ) THEN #Y_GLBC = Y_GLBC+1
C :
K : #NDISP2:
K :   /HMPACK MATRIX COPY  DISP DISPLS
C :
K : #NDISP:
C :                                        
C :----------------------------------------------------------------------
C : 
C :     Check energy convergence if required.
C : 
K :   #IF ( Y_ENRG EQ 0 ) THEN GOTO NENRG 
K :   #IF ( Y_CIEQ NE 1 ) THEN GOTO ENRG1
C :
C :     On first equilibrium iteration compute strain energy, TOTENO.
C :
K :   /HMPACK MATRIX KILL  TOTENO
K :   /HMPACK MATHEMATICS TRANSPOSE_MULTIPLY  FEXTSV DISP TOTENO
K :   #GOTO NENRG 
C :
K : #ENRG1:
C :
C :     Compute strain energy increment, DELEN.
C :
K :   /HMPACK MATRIX KILL  TOTEN DELEN
K :   /HMPACK MATHEMATICS TRANSPOSE_MULTIPLY  FEXTSV DISP TOTEN
K :   /HMPACK MATRIX COPY  TOTEN DELEN
K :   /HMPACK MATHEMATICS SUBTRACT  TOTENO DELEN
K :   /INSPECT NL_VECT_CONVERGENCE  "Y_ENG" "Y_CIEQ" DELEN TOTEN,
K :                                 "Y_CVG" "Y_RATE" "Y_DMIN" "Y_IAUSTP"
K :   /HMPACK MATRIX COPY  TOTEN TOTENO
K :   /HMPACK MATRIX KILL  DELEN TOTEN
K :   #IF ( Y_CVG LT 0 ) THEN GOTO DIVRGE
K :   #IF ( Y_CVG EQ 1 ) THEN #Y_GLBC = Y_GLBC+1
C :
K : #NENRG:
C :
C :----------------------------------------------------------------------
C :
C :     End of equilibrium iteration loop, check for errors, determine 
C :     stiffness parameter, check for convergence.
C :
K :   /LOG MESSAGE  "----- End of Iteration -----"
K :   /LOG MESSAGE " "
C :
K :   /GLOBAL_PARAMETER GET  "Y_PMOD"
K :   #IF ( Y_PMOD EQ 0 ) THEN GOTO QUIT
C :
C :     Determine if stiffness must be reformed based on stiffness 
C :     change, set the flag, Y_LUPK. Compute the current stiffness
C :     parameter, Y_CURSTF. Set the limit point flag, Y_LATLIM, if
C :     a limit point is detected.
C :
K :   /HMPACK PROGRAM_FILE RUN
K :                                "stfpar.prg"
C :
C :     Update the previous iteration stiffness parameter, Y_PRESTF, if
C :     the stiffness matrix is to be reformed, Y_LUPK = 1.
C :
K :   #IF ( Y_LUPK EQ 1 ) THEN #Y_PRESTF = Y_CURSTF
C :
C :     Check global convergence, if not converged do another iteration.
C :
K :   /HMPACK PROGRAM_FILE RUN
K :                                "cnvgrat.prg"
C :
K :   #IF ( Y_CIEQ LT 3 ) THEN GOTO ITERLP
K :   #IF ( Y_GLBC GE Y_GLBM AND Y_NCSC EQ 0 ) THEN GOTO CONVRG
C :
K :   #GOTO ITERLP
C :
C :----------------------------------------------------------------------
C :
C :     Programability error
C :
K : #PERROR:
C :
K :   #OUTPUT "A program file error has occurred.. Error #", Z_ERROR
K :   #Y_ERNO = 17907
K :   #Y_NMAR = 1
K :   #Y_SVER = 4
K :   /LOG ERROR  "Y_ERNO" "Y_SVER" "Y_NMAR" "Z_ERROR"
C :
C :     Set to verification mode Y_PMOD=0.
C :
K :   #Y_PMOD = 0
K :   /GLOBAL_PARAMETER PUT  "Y_PMOD"
K :   #GOTO QUIT
C :
C :----------------------------------------------------------------------
C :
K : #DIVRGE:
C :
K :   /HMPACK PROGRAM_FILE RUN
K :                             "cnvgrat.prg"
C :
K :   #Y_IDIVRG = 1
K :   /LOG MESSAGE  "Solution point terminated due to divergence"
K :   /LOG MESSAGE  "    Use smaller load increments"
K :   #GOTO QUIT
C :
C :----------------------------------------------------------------------
C :
K : #CONVRG:
C :
C :     The solution has converged, so set the current stiffness 
C :     parameter for the step, Y_CURK, equal to the current iteration
C :     stiffness parameter, Y_CURSTF. Also set the current step length, 
C :     Y_CURL, equal to the time increment, Y_DELTT.
C :
K :   #Y_CURK = Y_CURSTF
K :   #Y_CURL = Y_DELTT
C :
C :*********************************************************************
C :---------------------------------------------------------------------
C :
K :   #IF ( Y_CONFLG EQ 0 ) THEN GOTO GAP9
C : 
C :     For contact analysis, update contact traction forces, ELAMDA, 
C :     create element contact forces, ELCFOR, contact tractions, ELTRCT, 
C :     nodal contact forces, CTFORC, and update old contact displacements, 
C :     ECDISO. Note: ELTRCT = ELAMDA / AREA.
C :
K :   #Y_CONV = 1
K :   /FORM ELEMENT STATIC CONTACT_FORCE  ECTOFF ECSTAT ECDISP ECDISO,
K :                            ELAMDA DLAMDA ELCFOR "Y_CONV" "Y_CIEQ",
K :                            "Y_MODTYP" ECIDIS
K :   /HMPACK MATRIX KILL  ELCFOR ELTRCT CTFORC
K :   /FORM ELEMENT STATIC TRACTION_FORCE  ECTOFF ECSTAT ELAMDA,
K :                                        ELCFOR ELCFRT ELTRCT "Y_MODTYP"
K :   /HMPACK BOOLEAN_RHS PARTITION  EBOCFU ELCFOR CTFORC
C :
K :   /HMPACK MATRIX COPY  DISP DISPRV
C :
K : #GAP9:
C :
C :---------------------------------------------------------------------
C :*********************************************************************
C :
C :     Backup applied loads, FEXT, and displacements, DISP, for 
C :     stiffness parameter calculation.
C :
K :   /HMPACK MATRIX KILL  FEXTB2 DISPB2
K :   /HMPACK MATRIX INQUIRE  FEXT "Z_EXIST"
K :   #IF ( Z_EXIST EQ 1 ) THEN /HMPACK MATRIX COPY  FEXT FEXTB2
K :   /HMPACK MATRIX INQUIRE  DISP "Z_EXIST"
K :   #IF ( Z_EXIST EQ 1 ) THEN /HMPACK MATRIX COPY  DISP DISPB2
C :
C :     Update plastic strains, EPSPLA, back stress, ALPHA, and isotropic
C :     hardening variable, KAPPA. Also save the plastic strain and back
C :     stress prior to update for use in the calculation of the
C :     stiffness matrix at the beginning of the next load step, 
C :
K :   /HMPACK MATRIX COPY  EPSPLA OEPSPL
K :   /HMPACK MATRIX COPY  ALPHA  OALPHA
K :   /HMPACK MATHEMATICS ADD  DEPSPL EPSPLA
K :   /HMPACK MATHEMATICS ADD  DALPHA ALPHA
K :   /HMPACK MATHEMATICS ADD  DKAPPA KAPPA
K :   /HMPACK MATRIX KILL  DKAPPA
C :
C :     Compute the displacement increment from the last user specified 
C :     solution point if arc length is inactive.
C :
K :   #IF ( Y_ARCSRC EQ 1 ) THEN GOTO QUIT
K :   /HMPACK MATRIX KILL  DELDIS
K :   /HMPACK MATRIX COPY  DISP DELDIS
K :   /HMPACK MATHEMATICS SUBTRACT DISPBK DELDIS
C :
C :----------------------------------------------------------------------
C :
K : #QUIT:                   
C :
K :   #Y_NITRN = Y_CIEQ
E :