Free-Edge FEM Program Development


Development of a Laminate Free-Edge
Finite Element Computer Program

by

Ronald D. Kriz
Associate Professor
Engineering Science and Mechanics
Virginia Polytechnic Institute and State University
Blacksburg, Virginia 24061


Prerequisites for reader: Introductory course in Finite Element Modeling.


1.0 Motivation

Development of this program was motivated largely by working with a finite element computer program, developed at Virginia Tech in the late 1970s, that predicted stresses near the "Free-Edge" of a laminated composite. This program was mostly developed by graduate students who worked on a variety of projects funded by the Co-operative NASA Langley Virginia Tech Composite Materials Program, under the direction of Professor Carl Herakovich. Several VPI&SU reports are referenced below that document this development, see Refs.[1-3]. Over the course of five years many additions where made to the code, where not any one of the contributors had a comprehensive knowledge of the entire FEM Free-Edge program. The program also grew to a massive size, largely due to the fact that Virginia Tech at that time (1973-1979) had one of the first IBM360 virtual memory computers that allowed programmers to take advantage of this particular large memory architecture. As a result the code could not be ported to other large computers that had not embraced the virtual memory architecture in the late 1970s. Consequently the program developed here was the result of an attempt to port the Virginia Tech program to a CDC 7600 computer at the National Bureau of Standards at Boulder Colorado in 1980. At that time CDC computers, the precursor computers to the Cray, had fixed "In-Core"* memory, as we called it at that time, of about 120,000 double precision floating point numbers where as the Virginia Tech code required over 6,000,000 double precision floating point numbers.


* NOTE: "In-Core" memory was actually very small, 10-3m, ferro-magnetic disk with a hole in the center through which wires, that carried currents, could change the on-off bit, zero or one.

2.0 Introduction: Background

After close examination of the existing Virginia Tech FEM "Free-Edge" computer program, it was decided that necessary modifications were not possible. Therefore the author chose to not only rewrite the program to run on nonvirtual memory computer architectures, but also decided to carefully document the creation of this particular Finite Element Model "Free-Edge" program so that future students and interested users could more easily use and modify the program for future development. This effort began in January 1980 as part of the first of a two year National Research Post Doctoral Research Program at the National Bureau of Standards (now NIST) at Boulder Colorado. Because of the emphasis on developing a well documented code the author chose to use a popular code called SAPIV which was developed at the Earthquake Engineering Research Center at the University of California at Berkeley. The code was freely available to the public. This particular code formulation was very well documented by the authors, K.L. Bathe and E.L. Wilson, who also wrote a popular textbook, "Numerical Methods in Finite Element Analysis" and reports from the Earthquake Engineering Research Center, see Ref.[4,5]. The other popular finite element text book in the late 1970s was written by Zienkiewicz, Ref.[6], but only code fragments existed in the Appendices and no code related to Zienkiewicz's text book was commercially available. Commercially available codes such as NASTRAN and ASKA were expensive and were difficult to add new functionality to the existing code. SAPV followed by SAPIV, had an additional attractive feature at that time which allowed nonvirtual memory computers to solve for large FEM models by using a modified Gauss elimination "out-of-core" storage scheme for solving large simultaneous equations that exceeded "In-Core" memory. SAPV was rewritten to move blocks of data in and out of the Gauss Elimination equation solver to "scratch files" (external data storage devices such as tape drives at that time) which allowed solving large FEM models on small "In-Core" computers. Of course this added feature resulted in a trade-off towards longer overall computation time, but this did allow researchers to solve large FEM models on small "In-Core" memory computers. This proved eventful with the emerging limited memory desktop computers of the early 1980s and the advent of expensive but affordable hard disk drives as external "out-of-core" data storage devices.

3.0 Introduction: Program Description

Based on the rationale described in the previous section the author proceeded to study the text book by Bathe and Wilson and the source code of SAPIV and SAPV from which a new FEM "Free-Edge" program was constructed. Of course SAPIV and SAPV did not have the necessary "Generalized Out-of-Plane Strain" finite element needed for the laminate "Free-Edge" problem, so this element was incorporated into a "SAPIV-like" program. The finite element formulation for the "Free-Edge" problem is outlined prior to this section. Because of the limited computer memory available on NBS computers at that time, the existing SAPIV program was not modified but instead a much smaller memory program was written from "scratch" that would not only run on NBS computers but also the new desktop PC computers of the early 1980s. Even today this strategy has proven useful for students who wish to use this same computer program on their personal computers. Because of the small size of the computer program and liberal use of comment statements through out the program, this program is particularly useful to introduce students who have only had an introductory finite element analysis course. Most of the documentation on the construction of this finite element computer program is given in the reference text book, Ref.[4]. This program uses the same notation as defined in Ref.[4] and SAPIV and SAPV. The documentation provided below together with Ref.[4,5] should provide the reader the necessary description that would allow program modification for future development.

PROGRAM OBJECTIVES:
  1. Calculate the stress gradients near the free edge of a multi-layered composite laminate.

  2. Calculate this stress gradient using the minimum amount of "In-Core" storage.

  3. Structure and organize the coded from an instructive view-point for future development.


4.0 Main Program and Subroutine Description


MAIN Program: Linear Finite Element Model with In-core equation solver (LFEMI)

LFEMI uses the idea of one shared memory variable called SHARE that is stored in a FORTRAN COMMON memory block. Before each subroutine is called, SHARE is partitioned into segments and assigned to different variables (see comments in main program prior to each subroutine call). When possible data is stored on scratch files and only arrayed variables used in the next subroutine are passed through the SHARE variable common block. Comments embedded throughout this computer program should demonstrate how this SHARE variable is implemented as a shared memory variable.

SUBROUTINES: The following subroutines are called in sequence. A brief description is given for each subroutine. It is important to note how the shared memory variable SHARE is repartitioned prior to each subroutine call.

NEDATA: Called by MAIN. Reads Node-E lement DATA from tape unit LR and calculates element areas and geometric parameters which are stored for reference on tapes (1,2,3).

eqnum

EQNUM Called by MAIN. Assigns EQuation NUMbers to nodal degrees-of-freedon (DOF).

CIJ Called by MAIN. Calculates elastic constants, Cij for each element which is stored on tape (4).

ADDRS Called by MAIN. Calculates the diagonal ADDReSses for the global stiffness array.

ELMNT Called by MAIN. Calculates: 1) stiffness array for each ELeMeNT, 2) Right-hand-side load vector for each element, 3) Thermal-moisture strains for each element.

ASSEMB Called by ELMNT. ASSEMBles: 1) Element stiffness arrays into global stiffness array, 2) Element load vector into global load vector.

NDLOAD Called by MAIN. Adds prescribed NoDal LOADs to global load vector.

NDDISP Called by MAIN. Adds large stiffnesses to diagonals of global stiffness array so as to simulate a prescribed NoDal DISPlacements.

COLSOL Called by MAIN. COLumn linear equation SOLver calculates nodal displacements

NDFORC Called by MAIN. Recovers NoDal FORCes corresponding to prescribed displacements.

DISP Called by MAIN. Recovers: 1) NoDal DISPlacements, 2) Thermal-moisture strain normal to the plane of elements.

STR Called by MAIN. From the displacements calculates STResses and STRains.


5.0 Assignment of Data-Tape Unit Numbers (local scratch files)

If your computer has a large memory architecture it is easy to store these arrays as variables in a labeled common block. If your computer has limited memory you may want to continue to store these arrays as scratch files. The later is slower but you can solve larger FEM grids. An estimate of the size of these arrays is given in the note below the table.

Tape Unit Number Stored Variables Range No. of Words Stored
(Single Precision)
1 ND1,ND2, ND3
(Node no.s assoc. each element)
1->NE 3*NE
2 THRAD
(Fiber angle (radians) each element
1->NE NE
3 AA,BB,CC,DD,EE,GG,AREA
(Geometric parameters each element)
1->NE 7*NE
4 D11,D12,D13,D16,D22, ...,D66
(Stiffnesses each element)
1->NE 13*NE
5
6
Data read from TAPE5="filename"
Data read from TAPE6="filename"
. .
7 ENO,EHO,ENHO,EZO
(Thermal-moisture strains each element)
1->NE 4*NE
8 A(NCOFS)
(Global Stiffnesses)
1->NCOFS NCOFS
9 NOD,IDIRN,DISP
(Prescribed displs.;Node no., Magnitude)
1->NDISP 3*NDISP

Total no. of single precision words stored on local scratch files

               Total=28*NE+NCOFS+3*NDISP

For the CDC 7600 computer a total of 53,248 PRUs are allowed. There are 64 words for one PRU. Therefore a total of 3,599,872 words can be stored on local scratch files. On present day desktop workstations the size of scratch files is only limited by the available hard disk storage.


6.0 Output File Data Output: KEY Parameter Defined

The Print KEY parameter was designed to selectively print out results from different portions of the code to help diagnose errors. If KEY is 7 all results are printed out.

Data sent to output file -- IF --> Value of KEY Parameter is
Heading
(First one to five lines of comments from input file)
Always Printed
Grid Parameters
(No. of Nodes, Elements, etc.)
Always Printed
Material Properties
(Elastic Prop.s, Expansion Coeff.s)
Always Printed
Node-Element Geometry
(Assignment of node no.s to elements, coordinates of nodes, element areas and geometric parameters)
KEY = 2 or 7
Load Types and Magnituds
(Mechanical, Thermal and Moisture)
Always Printed
Equation Numbers Assigned to Nodal DOFs KEY = 2 or 7
Material Elastic Properties KEY = 3 or 7
Element Material Stiffnesses Each Element KEY = 3 or 7
Element Stiffnesses Each Element KEY = 3 or 7
Column Heights of Global Stiffness Array KEY = 4 or 7
Band Width of Global Stiffness Array Always Printed
Number of Terms in Global Stiffness Array Always Printed
Diagonal Addresses of Global Stiffness Array KEY = 2 or 7
Total Size of Arrays Stored in Common Blocks Always Printed
Connectivity Array Each Element KEY = 3 or 7
Stiffness Array and Force Vector Each Element KEY = 3 or 7
Global Stiffness Array and Load Vector KEY = 4 or 7
Nodal Load Information Always Printed
Nodal Displacement Information Always Printed
Nodal Forces Resulting from Prescribed Nodal Displacements Always Printed
Stresses Each Element Always Printed
Strains for each element and nodal displacements KEY = 5 or 7


7.0 References:

  1. G.D. Renieri and C.T. Herakovich, "Nonlinear Analysis of Laminated Fibrous Composites", Virginia Polytechnic Institute and State University, VPI-E-76-10, 1976.

  2. P.W. Hsu and C.T. Herakovich, "A limiting Analysis for Edge Effects in Angle-Ply Laminates", Virginia Virginia Polytechnic Institute and State University, Report VPI-E-76-18, September 1976.

  3. A.P. Nagarkar and C.T. Herakovich, "Nonlinear Temperature Dependent Failure Analysis of Finite Width Composite Laminates,", Virginia Polytechnic Institute and State University, Report VPI-E-79-36, September, 1979.

  4. K.J. Bathe and E.L. Wilson, Numerical Methods in Finite Element Analysis, Prentice-Hall Inc., New Jersey, 1976.

  5. K.J. Bathe, E.L. Wilson, and F.E. Peterson, "SAPIV-Structural Analysis Program for Static and Dynamic Response of Linear Systems," Earthquake Engineering Research Center, Report No. UCB/EERC-73/11, University of California, Berkeley, June 1973, revised April 1974.

  6. O.C. Zienkiewicz, "The Finite Element Method in Engineering Science, " Second Edition, McGraw-Hill, London, 1971.


Appendix A: SOURCE CODE with Embedded Notes



NOTE: This is identically the same code that is compiled and executed on the NPIB form.

PROGRAM edge 
C 
C23456789012345678901234567890123456789012345678901234567890123456789012
C        1         2         3         4         5         6         7
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C  *                     M A I N   P R O G R A M                       *
C  *                          NAME: ( LFEMI )                          *
C  *       L  INEAR   F  INITE   E  LEMENT   M  ODEL   I  N-CORE       *
C  *                                                                   *
C  *          USED TO SOLVE THE LAMINATE "FRRE-EDGE" PROBLEM           *
C  *                                                                   *
C  *  CORE STORAGE MINIMIZED                                           *
C  *    1. STORES ELEMENT-NODE DATA ON SCRATCH FILES (TAPE 1-4,7-9).   *
C  *    2. DIMENSIONED VARIABLES SHARE HIGH SPEED CORE STORAGE.        *
C  *    3. USES THE SKYLINE STORAGE SCHEME, DEVELOPED BY BATHE AND     *
C  *       WILSON ( NUMERICAL METHODS IN FINITE ELEMENT ANALYSIS), TO  * 
C  *       STORE THE STIFFNESS ARRAY AS A SINGLE DIMENSIONED ARRAY     *
C  *       SUCH THAT A COLUMN ORIENTED GAUSS ELIMINATION TECHNIQUE IS  *
C  *       USED IN SUBROUTINE COLSOL TO SOLVE THE SET OF LINEAR EQNS.  *
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C 
      REAL NU12,NU13,NU23 
C 
      COMMON/RDK1/LW,LR,KEY 
      COMMON/RDK2/E11,E22,E33,G12,G13,G23,NU12,NU13,NU23
      COMMON/RDK3/AL1,AL2,AL3,BT1,BT2,BT3 
C 
      DIMENSION TITLE(20,5) 
C 
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C  *  THE FOLLOWING TWO CARDS DETERMINE THE MAXIMUM HIGH SPEED CORE    *
C  *  STORAGE SHARED BY DIMENSIONED VARIABLES IN VARIOUS SUBROUTINES.  *
C  *     TO CHANGE STORAGE, CHANGE BOTH MTOT AND SHARE(MTOT).          *
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C 
      COMMON SHARE(75021)
      MTOT=75021 
C 
C   ASSIGN TAPE NO.S 5 AND 6 TO INPUT AND OUTPUT FILES RESPECTIVELY 
C 
      OPEN(5,FILE='edge.dat',STATUS='OLD',ERR=8888)
      OPEN(6,FILE='edge.out',STATUS='UNKNOWN',ERR=8888)
      OPEN(7,FILE='edge_stress.out',STATUS='UNKNOWN',ERR=8888)
C 
C           H  E  A  D  I  N  G   
C 
      WRITE(6,99)
   99 FORMAT(1H1) 
      WRITE(6,999) 
  999 FORMAT(/,80(1H*),/) 
      READ(5,10)((TITLE(I,J),I=1,20),J=1,5)
   10 FORMAT(20A4)
      WRITE(6,11)((TITLE(I,J),I=1,20),J=1,5) 
   11 FORMAT(1X,20A4) 
      WRITE(6,999) 
C 
C         G  R  I  D     P  A  R  A  M  E  T  E  R  S 
C 
      READ(5,20)NE,NDS,NLOAD,NDISP,KEY 
   20 FORMAT(5I6) 
      READ(5,21)SMH,SMZ
   21 FORMAT(2F12.0)
C 
      WRITE(6,22)NE,NDS,NLOAD,NDISP,SMH,SMZ,KEY
   22 FORMAT(//,5X,'G R I D   P A R A M E T E R S',/, 
     1' NUMBER OF ELEMENTS .........................',I4,/, 
     2' NUMBER OF NODES ............................',I4,/, 
     3' NUMBER OF NODE LOADS .......................',I4,/, 
     4' NO. OF PRESCRIBED NODE DISPLACEMENTS .......',I4,/, 
     5' SCALE FOR HORIZONTAL COORD.S  H(I) .........',F12.5,/, 
     6' SCALE FOR VERTICAL   COORD.S  Z(I) .........',F12.5,/, 
     7' OUTPUT FILE PRINT KEY ......................',I4) 
C 
C         M  A  T  E  R  I  A  L    P  R  O  P  E  R  I  T  I  E  S 
C 
      WRITE(6,30)
   30 FORMAT(//,10X,'M A T E R I A L   P R O P E R T I E S',/)       
C 
C     READ(5,40)E11,E22,E33
      READ(5,*)E11,E22,E33
C     READ(5,40)G12,G13,G23
      READ(5,*)G12,G13,G23
C     READ(5,40)NU12,NU13,NU23 
      READ(5,*)NU12,NU13,NU23 
C  40 FORMAT(3(6X,E10.3)) 
      WRITE(6,41)E11,E22,E33 
   41 FORMAT(/,15X,'YOUNGS MODULI',/,' E11=',E10.3,' E22=',E10.3,     
     *' E33=',E10.3,/)
      WRITE(6,42)G12,G13,G23 
   42 FORMAT(/,15X,'SHEAR MODULI',/,' G12=',E10.3,' G13=',E10.3,
     *' G23=',E10.3,/)
      WRITE(6,43)NU12,NU13,NU23
   43 FORMAT(/,15X,'POISSONS RATIO',/,' NU12=',E10.3,' NU13=',E10.3,
     *' NU23=',E10.3,//)
C 
C     READ(5,50)AL1,AL2,AL3
      READ(5,*)AL1,AL2,AL3
C     READ(5,50)BT1,BT2,BT3
      READ(5,*)BT1,BT2,BT3
C  50 FORMAT(3(8X,E10.3)) 
C 
      WRITE(6,51)AL1,AL2,AL3
   51 FORMAT(/,15X,'THERMAL EXPANSION COEFFICIENTS',/,
     1' ALPHA-1=',E10.3,' ALPHA-2=',E10.3,' ALPHA-3=',E10.3,/)
      WRITE(6,52)BT1,BT2,BT3 
   52 FORMAT(/,15X,'MOISTURE EXPANSION COEFFECIENTS',/, 
     2' BETA-1=',E10.3,' BETA-2=',E10.3,' BETA-3=',E10.3,//)
C 
C              L O A D    I N F O R M A T I O N 
C  O T H E R   T H E N   L O A D S   A T   N O D E   P O I N T S
C 
      WRITE(6,60)
   60 FORMAT(1H1,10X,'L O A D   I N F O R M A T I O N',/,2X,
     *'O T H E R   T H E N   L O A D S   A T   N O D E   P O I N T S',
     *//) 
C 
C     READ(5,61)LTYPE
      READ(5,*)LTYPE
C  61 FORMAT(I6)
C     READ(5,62)SLOAD
      READ(5,*)SLOAD
C  62 FORMAT(F12.0) 
C     READ(5,63)TLOAD,WLOAD
      READ(5,*)TLOAD,WLOAD
C  63 FORMAT(2F12.0)
C 
      WRITE(6,64)LTYPE 
   64 FORMAT(/,10X,'LOAD TYPE = ',I1,' ( = 0, NODE LOADS ONLY)',
     */,23X,' ( = 1, STRAIN APPLIED NORMAL TO ELEMENTS ONLY)',
     */,23X,' ( = 2, CHANGE IN THERMAL-MOISTURE CONDITIONS ONLY)',//)
      WRITE(6,65)LTYPE 
   65 FORMAT(5X,'LOAD TYPE =',I2,///) 
      IF(LTYPE.LE.1)WRITE(6,66)SLOAD 
   66 FORMAT(' STRAIN APPLIED NORMAL TO H-Z PLANE OF ELEMENTS',/,
     1'         SLOAD=',E11.4,//) 
      IF(LTYPE.EQ.2)WRITE(6,67)TLOAD,WLOAD 
   67 FORMAT(' STRAINS NORMAL TO PLANE OF ELEMENTS',/, 
     4'          TO BE CALCULATED',/, 
     5'    FOR NONZERO  TEMPERATURE CHANGE, TLOAD=',E11.4,/,
     6'                 MOISTURE  CHANGE,   WLOAD=',E11.4,/)
C 
C      N O D E - E L E M E N T    I N F O R M A T I O N 
C 
C     ADDRESSES 
      N1=1
      N2=N1+3*NDS 
      N3=N2+NDS 
C 
C          NODE BOUNDARY CONDITION CODES  ID(NDS,3) 
C          NODE HORIZONTAL CO-ORDINATES    H(NDS) 
C          NODE VERTICAL CO-ORDINATES      Z(NDS) 
C 
C                 ID(NDS,3)  H(NDS)    Z(NDS)
      CALL NEDATA(SHARE(N1),SHARE(N2),SHARE(N3),NDS,NE,SMH,SMZ) 
C 
C  C A L C U L A T E   M A T E R I A L   S T I F F N E S S E S
C             F O R   E A C H   E L E M E N T 
C 
      CALL CIJ(NE)
C 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *      F.  E.  M.     S  O  L  U  T  I  O  N        *
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
C  A S S I G N   E Q U A T I O N   N U M B E R S   TO 
C      E A C H   D E G R E E   O F   F R E E D O M
C            STORED IN ARRAY ID(NDS,3) 
C 
C                ID(NDS,3)
      CALL EQNUM(SHARE(N1),NDS,NEQ,NEQF,LTYPE)
C 
C C A L C U L A T E   A D D R E S S E S    U S E D   T O    S T O R E 
C      T H E   G L O B A L   S T I F F N E S S   A R R A Y
C 
C  GLOBAL ARRAY SIZE, NCOFS (NO. OF COEFFICIENTS)
C  GLOBAL ADDRESSES OF DIAGONAL TERMS, MAXA(NEQ1) 
C  GLOBAL HALF BANDWIDTH INCLUDING DIAGONAL, MBAND
C  ELEMENT ARRAY SIZE, NDOFE (DEGREES OF FREEDOM PER ELEMENT) 
C  THERMAL-MOISTURE LOAD EQN. NO., NEQF 
C 
C          CALCULATE NO. EQN.S PLUS ONE, NEQ1 USED TO 
C        STORE DIAGONAL TERMS OF GLOBAL STIFFNESS ARRAY 
      NEQ1=NEQ+1
C 
C                ID(NDS,3) MAXA(NEQ1) 
      CALL ADDRS(SHARE(N1),SHARE(N2),NE,NDS,NEQ1,LTYPE,NEQF,MBAND,NCOFS,
     *                               NDOFE) 
C 
C   C A L C U L A T E   T O T A L   S I Z E   O F   A R R A Y S 
C         T O   B E   S T O R E D   I N   C O M M O N
C 
C   ARRAYS USED IN SUBROUTINE ELMNT AND SHARED BY MAIN PROG. COMMON 
C       
C    ARRAY NAMES     COMMON   ADDRESSES 
C     ID( NDS,3 ) = SHARE(N1)  N1=1 
C   MAXA( NEQ1  ) = SHARE(N2)  N2=N1+3*NDS
C      A( NCOFS ) = SHARE(N3)  N3=N2+NEQ1 
C      V(  NEQ  ) = SHARE(N4)  N4=N3+NCOFS
C 
      NTOT1=3*NDS+NEQ1+NCOFS+NEQ 
C 
C   ARRAYS USED IN SUBROUTINE STR AND SHARED BY MAIN PROG. COMMON 
C 
C    ARRAY NAMES     COMMON   ADDRESSES 
C     ID( NDS,3 ) = SHARE(N1)  N1=1   SPACE ALLOCATED FOR ID ARRAY  
C                                            PRIOR TO SUBPROG. STR 
C      U( NDS,3 ) = SHARE(N2)  N2=N1+NDS*3 
C     EN(   NE  ) = SHARE(N3)  N3=N2+NDS*3
C     EH(   NE  ) = SHARE(N4)  N4=N3+NE 
C     EZ(   NE  ) = SHARE(N5)  N5=N4+NE 
C    EHZ(   NE  ) = SHARE(N6)  N6=N5+NE 
C    ENZ(   NE  ) = SHARE(N7)  N7=N6+NE 
C    ENH(   NE  ) = SHARE(N8)  N8=N7+NE 
C 
      NTOT2=6*NDS+6*NE
C 
C    CHOOSE LARGEST SHARED ARRAY SIZE 
C 
      NTOT=NTOT1
      IF(NTOT2.GT.NTOT1)NTOT=NTOT2
C 
      WRITE(6,99)
      WRITE(6,999) 
      WRITE(6,70)MTOT,NTOT1,NTOT2 
   70 FORMAT(' TOTAL SIZE OF ARRAYS STORED IN COMMON',/,
     1' AVAILBLE(',I6,')    REQUIRED BY SUBROUTINE ELMNT  (',I6,')',/,
     233X,'SUBROUTINE   STR  (',I6,')') 
      WRITE(6,999) 
C 
      IF(NTOT.GT.MTOT)STOP
C 
C   C A L C U L A T E   V A R I O U S   E L E M E N T   A R R A Y S 
C        CONNECTIVITY       STIFFNESS           FORCE 
C        LM( NDOFE )    AK( NDOFE,NDOFE )    F( NDOFE ) 
C 
C                 ----------  AND  ---------- 
C       A S S E M B L E   E L E M E N T   A R R A Y S   I N T O 
C                    G L O B A L    A R R A Y S 
C              STIFFNESS                    LOAD
C              A( NCOFS )                  V( NEQ ) 
C 
C     ADDRESSES 
      N3=N2+NEQ1
      N4=N3+NCOFS 
C 
C                ID(NDS,3) MAXA(NEQ1) A(NCOFS)  V(NEQ)
      CALL ELMNT(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4),NE,NDS,
     *           LTYPE,SLOAD,TLOAD,WLOAD,NCOFS,NEQ,NEQ1,NEQF,NDOFE) 
C 
C    I F   A V A I L A B L E
C        S U P E R P O S E   N O D E   F O R C E S
C 
C                               ID(NDS,3)  V(NEQ) 
      IF(NLOAD.GT.0)CALL NDLOAD(SHARE(N1),SHARE(N4),NDS,NEQ,NLOAD)
C 
C    I F   A V A I L A B L E    I N C O R P O R A T E 
C        P R E S C R I B E D   N O D E   D I S P L A C E M E N T S
C 
C                               ID(NDS,3) MAXA(NEQ1) A(NCOFS)  V(NEQ) 
      IF(NDISP.GT.0)CALL NDDISP(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4),
     *                          NDS,NEQ1,NCOFS,NEQ,NDISP) 
C 
C**********  S O L V E   S Y S T E M   O F   E Q N. S  *************
C 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   *                                                                 * 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
C   CALL COLSOL TWICE; 1ST TO TRIANGULARIZE  A(NCOFS) 
C                      2ND TO SOLVE FOR DISPLACEMENTS
C                      LOADS STORED IN V(NEQ) REPLACED BY DISPLACEMENTS 
C 
C                 MAXA(NEQ1) A(NCOFS)  V(NEQ) 
      CALL COLSOL(SHARE(N2),SHARE(N3),SHARE(N4),NEQ,NCOFS,NEQ1,1) 
      CALL COLSOL(SHARE(N2),SHARE(N3),SHARE(N4),NEQ,NCOFS,NEQ1,2) 
C 
C*************  S Y S T E M   O F   E Q N. S   S O L V E D  ************ 
C 
C   I F   D I S P L A C E M E N T S   W H E R E   P R E S C R I B E D 
C   R E C O V E R   C O R R E S P O N D I N G   N O D E   F O R C E S
C 
C                               ID(NDS,3) MAXA(NEQ1) A(NCOFS)  V(NEQ) 
      IF(NDISP.GT.0)CALL NDFORC(SHARE(N1),SHARE(N2),SHARE(N3),SHARE(N4),
     *                          NDS,NEQ1,NCOFS,NEQ,NDISP,MBAND) 
C 
C            R E C O V E R    D I S P L A C E M E N T S  U(NDS,3) 
C  F R O M  COLSOL  S O L T I O N S   S T O R E D   I N   A R R A Y  V(NEQ) 
C 
C               ID(NDS,3) U(NDS,3)   V(NEQ) 
      CALL DISP(SHARE(N1),SHARE(N2),SHARE(N4),NDS,NEQ,NEQF,LTYPE,UKEEP) 
C 
C      C A L C U L A T E   S T R E S S E S   A N D   S T R A I N S
C 
C     ADDRESSES 
      N2=N1+3*NDS 
      N3=N2+3*NDS
      N4=N3+NE
      N5=N4+NE
      N6=N5+NE
      N7=N6+NE
      N8=N7+NE
C 
C     SYMMETRIC STRAIN TENSOR COMPONENTS
C        W.R.T. LAMINATE AXIS (N,H,Z) 
C     ( EN , EH , EZ , EHZ , ENZ , ENH )
C 
C              U(NDS,3)   EN(NE)    EH(NE)    EZ(NE)    EHZ(NE) 
      CALL STR(SHARE(N2),SHARE(N3),SHARE(N4),SHARE(N5),SHARE(N6), 
     *         SHARE(N7),SHARE(N8),NE,NDS,LTYPE,SLOAD,UKEEP)
C               ENZ(NE)    ENH(NE)
C 
      REWIND 1
      REWIND 2
      REWIND 3
      REWIND 4

      CLOSE(1)
      CLOSE(2)
      CLOSE(3)
      CLOSE(4)
 8888 STOP
      END 


C***********************************************************************
      SUBROUTINE NEDATA(ID,H,Z,NDS,NE,SMH,SMZ) 
C         CALLING PROGRAM 
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C  *   P R O G R A M                                         *
C  *     READS NODE - ELEMENT DATA FROM TAPE UNIT LR.        *
C  *     CALCULATES ELEMENT AREAS AND GEOMETRIC PARAMETERS   *
C  *     STORES PARAMENTERS FOR REFERENCE ON TAPES (1,2,3)   *
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C 
      COMMON/RDK1/LW,LR,KEY 
C 
C  DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON
C 
      DIMENSION ID(NDS,3),H(NDS),Z(NDS) 
C 
C                  R E A D   N O D E    D A T A
C 
C   WRITE HEADER ON OUTPUT FILE 
C 
      IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(6,10)
   10 FORMAT(1H1,5X,'N O D E    I N F O R M A T I O N',//,
     16X,'(FIXED=1',/,8X,'FREE=0)',/,' NODE  MOTION     CO-ORDINATES',/,
     2'  NO. N  H  Z  ',6X,'H',8X,'Z')
C 
C        L O O P    T H R U    A L L   N O D E S
C 
      DO 40 I=1,NDS
C 
C   READ NODE DATA FROM INPUT FILE
C 
      READ(5,20)ID(I,1),ID(I,2),ID(I,3),H(I),Z(I)
   20 FORMAT(6X,3I3,2F12.0) 
C 
C    WRITE NODE DATA ON OUTPUT FILE 
C 
      IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(6,30)I,ID(I,1),ID(I,2),ID(I,3), 
     *                                        H(I),Z(I)
   30 FORMAT(1X,I4,3(1X,I1,1X),1X,2F9.5)
C 
   40 CONTINUE
C 
C       E N D   N O D E   L O O P 
C 
C 
C        R E A D  -  S T O R E   E L E M E N T   D A T A
C 
C         C A L C U L A T E :  RESCALED NODE CO-ORDINATES
C                              ASSOCIATED WITH EACH ELEMENT 
C 
C    C A L C U L A T E  -  S T O R E  1. ELEMENT FIBER ANGLE ORIENTATION
C                                     2. ELEMENT GEOMETRIC PARAMETERS 
C                                                   AND AREAS 
C 
C  INITIALIZE TOTAL GRID AREA 
C 
      TAREA=0 
C 
C    WRITE HEADER ON OUTPUT FILE
C 
      IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(6,50)
   50 FORMAT(1H1,5X,'E L E M E N T   I N F O R M A T I O N',//, 
     1' ELMT------NODES----- PLY',/,
     2'  NO.  1ST  2ND  3RD  ANGLE   AREA')
C 
C        CALCULATE PI (RADIANS TO DEGREES CONVERSION) 
C 
      PI=ACOS(-1.)
C 
C    L O O P   T H R U   A L L   E L E M E N T S
C 
      DO 90 K=1,NE
C 
C    READ ELEMENT DATA FROM INPUT FILE
C 
      READ(5,60)ND1,ND2,ND3,THETA
   60 FORMAT(6X,3I6,F8.2)
C 
C    STORE NODES ASSOCIATED WITH EACH ELEMENT ON TAPE 1 
C 
      WRITE(1)ND1,ND2,ND3 
C 
C   STORE FIBER ANGLE (IN RADIANS) FOR EACH ELEMENT ON TAPE 2 
C      ANGLE ORIENTATION (CLOCKWISE ROTATION, POSITIVE) 
C         WITH RESPECT TO LAMINATE NORMAL N-AXIS
C 
      THRAD=-THETA*PI/180. 
      WRITE(2)THRAD 
C 
C    NODE CO-ORDINATES ASSOCIATED WITH EACH ELEMENT RESCALED
C 
      H1=H(ND1)*SMH 
      H2=H(ND2)*SMH 
      H3=H(ND3)*SMH 
      Z1=Z(ND1)*SMZ 
      Z2=Z(ND2)*SMZ 
      Z3=Z(ND3)*SMZ 
C 
C    ELEMENT AREA CALCULATED
C 
      TERM1=-(-H1+H2)*(Z1+Z2)/2.0 
      TERM2=-(-H2+H3)*(Z2+Z3)/2.0 
      TERM3=-(-H3+H1)*(Z3+Z1)/2.0 
      AREA=TERM1+TERM2+TERM3
C 
C    ELEMENT GEOMETRIC PARAMETERS CALCULATED
C 
      AA=(Z2-Z3)/2.0
      BB=(H3-H2)/2.0
      CC=(Z3-Z1)/2.0
      DD=(H1-H3)/2.0
      EE=(Z1-Z2)/2.0
      GG=(H2-H1)/2.0
C 
C    STORE AREA AND GEOMETRIC PARAMETERS ON TAPE 3
C 
      WRITE(3)AA,BB,CC,DD,EE,GG,AREA
C 
C    WRITE ELEMENT DATA ON OUTPUT FILE
C 
      IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(6,70)K,ND1,ND2,ND3,THETA,AREA 
   70 FORMAT(4(1X,I4),1X,F6.2,E11.4)
C 
C    CALCULATE TOTAL GRID AREA
C 
      TAREA=TAREA+AREA
C 
   90 CONTINUE
C 
C      E N D   E L E M E N T   L O O P
C 
C    WRITE GRID AREA ON OUTPUT FILE 
C 
      WRITE(6,100)TAREA 
  100 FORMAT(//,5X,'TOTAL GRID AREA=',E10.3) 
C 
      RETURN
      END 


C***********************************************************************
      SUBROUTINE EQNUM(ID,NDS,NEQ,NEQF,LTYPE)
C         CALLING PROGRAM - MAIN
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  
C  *   P R O G R A M                                             *
C  *     GENERATES EQUATION NUMBERS BY REPLACING THE FREE DISPL  *
C  *     BOUNDARY CONDITION CODE (ID=0) WITH AN EQUATION NUMBER. * 
C  *     EQUATION NUMBERS ASSOCIATED WITH EACH FREE NODE DEGREE  *
C  *     OF FREEDOM (DOF) ARE STORED IN ARRAY ID(NDS,3) WHICH IS * 
C  *     SHARED WITH MAIN PROGRAM CONNOM SHARE(MTOT)             *
C  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C 
      COMMON/RDK1/LW,LR,KEY 
C 
C   DIMENSIONED VARIABLE SHARED WITH MAIN PROGRAM COMMON
C 
      DIMENSION ID(NDS,3) 
C 
C   INITIALIZE VARIABLES DEFINING EQUATION NUMBERS
C 
      IEQN=0 
      NEQF=0
C 
C   WRITE HEADER ON OUTPUT FILE 
C 
      IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(6,10)
   10 FORMAT(1H1,13X,'EQUATION NUMBERS',//,5(12X,'NODE MOTION  '),/,
     *           5(' # NODE  NORM  HORIZ VERT'))
C 
C   SEQUENTIALLY TEST THE DISPLACEMENT BOUNDARY CONDITION 
C   CODES( FIXED=1 OR FREE=0 ) FOR EACH NODE.  THERE ARE  
C   THREE POSSIBLE DEGREES OF FREEDOM ASSOCIATED WITH EACH NODE.
C      ID(I,1)=0, I-TH NODE IS FREE TO MOVE NORMAL TO THE H-Z PLANE 
C      ID(I,2)=0, I-TH NODE IS FREE TO MOVE IN THE HORIZ. H-DIRECTION 
C      ID(I,3)=0, I-TH NODE IS FREE TO MOVE IN THE VERT.  Z-DIRECTION 
C   IF A NODE IS FREE TO MOVE THEN AN EQUATION NUMBER (IEQN) IS 
C   ASSOCIATED WITH THAT NODE DEGREE OF FREEDOM.
C 
      DO 40 I=1,NDS 
      DO 20 J=1,3 
      IF(ID(I,J))6,5,6
    5 IEQN=IEQN+1 
      ID(I,J)=IEQN 
      GO TO 20
    6 ID(I,J)=0 
   20 CONTINUE
   40 CONTINUE
C 
C    WRITE EQN. NO.S ON OUTPUT FILE 
C 
      IF((KEY.EQ.2).OR.(KEY.EQ.7))WRITE(6,3)(I,(ID(I,J),J=1,3),I=1,NDS) 
    3 FORMAT(5(' #',I5,3I6)) 
C 
C      LAST VALUE OF IEQN IS EQUAL TO THE TOTAL NO. OF EQUATIONS
C 
      NEQ=IEQN
C 
C    WRITE TOTAL NUMBER OF EQUATIONS DUE TO NODAL D.O.F. S
C                ON OUTPUT FILE 
C 
      WRITE(6,45)NEQ 
   45 FORMAT(//,' TOTAL NO. OF EQN.S DUE TO NODAL DEGREES OF FREEDOM', 
     *I5,//) 
C 
C        IF LOAD TYPE IS MECHANICAL LOADS ONLY
C                    THEN RETURN
C 
      IF(LTYPE.LE.1)RETURN
C 
C   WHEN LOAD TYPE (LTYPE=2) IS EXCLUSIVELY THERMAL AND/OR MOISTURE 
C     ADD ONE ADDITIONAL EQN TO THE SET OF SIMULTANEOUS EQUATIONS 
C 
      NEQ=NEQ+1 
C 
C            NEQ = TOTAL NO. OF EQN.S TO BE SOLVED
C           NEQF = EQN. NO. WHICH RELATES UNKNOWN THERMAL OR MOISTURE 
C                 STRAIN (UKEEP) TO THE THERMAL OR MOISTURE LOAD B(NEQF)
C 
C 
C          ADDITIONAL EQN. NO. NEQF IS LOCATED AS THE LAST EQN.
C        OF THE ASSEMBLED (GLOBAL) STIFFNESS ARRAY SUCH THAT THE 
C    COLUMN ORIENTED STORAGE OF THE NO. OF COEFF.S (NCOFS) IS MINIMIZED 
C 
      NEQF=NEQ
C 
C    WRITE ADDITIONAL EQUATION NUMBER ON OUTPUT FILE
C 
      WRITE(6,50)NEQF 
   50 FORMAT(' ADDITIONAL EQN. FOR THERMAL-MOISTURE LOAD, NEQF=',I6) 
C 
      RETURN
      END 



C***********************************************************************
      SUBROUTINE CIJ(NE)
C         CALLING PROGRAM - MAIN
C   *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
C   *  P R O G R A M                                         *
C   *      CALCULATE AND STORE MATERIAL STIFFNESSES FOR      *
C   *            EACH ELEMENT OF THE GRID                    *
C   *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
C 
      REAL NU12,NU13,NU23,NU21,NU31,NU32
C 
      COMMON/RDK1/LW,LR,KEY 
      COMMON/RDK2/E11,E22,E33,G12,G13,G23,NU12,NU13,NU23
C 
C             CALCULATE MATERIAL STIFFNESSES
C 
      NU21=NU12*E22/E11 
      NU31=NU13*E33/E11 
      NU32=NU23*E33/E22 
C 
      DEL=1-NU12*NU21-NU23*NU32-NU31*NU13-2.0*NU21*NU32*NU13 
C 
      C11=(1.0-NU23*NU32)*E11/DEL 
      C12=(NU21+NU31*NU23)*E11/DEL 
      C13=(NU31+NU21*NU32)*E11/DEL 
      C22=(1.0-NU13*NU31)*E22/DEL 
      C23=(NU32+NU12*NU31)*E22/DEL 
      C33=(1.0-NU12*NU21)*E33/DEL 
      C44=G23 
      C55=G13 
      C66=G12 
C 
C            WRITE MATERIAL STIFFNESSES ON OUTPUT FILE 
C 
      IF((KEY.EQ.3).OR.(KEY.EQ.7))WRITE(6,20)C11,C12,C13,C22,C23,C33, 
     *                                        C44,C55,C66 
   20 FORMAT(1H1,/,15X,'MATERIAL STIFFNESSES',//,
     1' C11=',E11.4,' C12=',E11.4,' C13=',E11.4,/,
     216X,' C22=',E11.4,' C23=',E11.4,/,32X,' C33=',E11.4,/,
     348X,' C44=',E11.4,/,64X,' C55=',E11.4,/,80X,' C66=',E11.4)
C 
C      STIFFNESS ARRAY TRANSFORMATIONS USING THRAD=THETA(RADIANS)
C                     FOR THE K-TH ELEMENT
C                                                                     
      REWIND 2
C 
C      WRITE HEADER ON OUTPUT FILE
C 
      IF((KEY.EQ.3).OR.(KEY.EQ.7))WRITE(6,25)
   25 FORMAT(//,15X,'ELEMENT MATERIAL STIFFNESSES')
C 
      DO 10 K=1,NE
      READ(2)THRAD
      C=COS(THRAD)
      SN=SIN(THRAD) 
      C2=C*C 
      C4=C2*C2 
      S2=SN*SN 
      S4=S2*S2 
      D11=C11*C4+2.*C2*S2*(C12+2.*C66)+C22*S4 
      D12=C2*S2*(C11+C22-4.*C66)+C12*(S4+C4) 
      D13=C13*C2+C23*S2 
      D16=-C*SN*(C11*C2-C22*S2-(C12+2.*C66)*(C2-S2)) 
      D22=C11*S4+2.*C2*S2*(C12+2.*C66)+C22*C4 
      D23=C13*S2+C23*C2 
      D26=-C*SN*(C11*S2-C22*C2+(C12+2.*C66)*(C2-S2)) 
      D33=C33 
      D36=C*SN*(C23-C13)
      D44=C44*C2+C55*S2
      D45=C*SN*(C44-C55)
      D55=C55*C2+C44*S2
      D66=(C11+C22-2.*C12)*C2*S2+C66*(C2-S2)**2 
C 
C        STORE ELEMENT STIFFNESSES ON TAPE 4
C 
      WRITE(4)D11,D12,D13,D16,D22,D23,D26,D33,D36,D44,D45,D55,D66 
C 
C               WRITE ELEMENT STIFFNESSES ON OUTPUT FILE 
C 
      IF((KEY.EQ.3).OR.(KEY.EQ.7))WRITE(6,35)K,D11,D12,D13,D16,D22,D23,
     *D26,D33,D36,D44,D45,D55,D66 
   35 FORMAT(' ELEMENT NO.=',I4,/,' D11=',E11.4,' D12=',E11.4,
     1' D13=',E11.4,33X,'D16=',E11.4,/,16X,' D22=',E11.4,' D23=',E11.4,
     2,32X,' D26=',E11.4,/,32X,' D33=',E11.4,32X,' D36=',E11.4,/,
     348X,' D44=',E11.4,' D45=',E11.4,/,64X,' D55=',E11.4,/,80X,
     4' D66=',E11.4)
C 
   10 CONTINUE 
C 
      RETURN
      END 



C***********************************************************************
      SUBROUTINE ADDRS(ID,MAXA,NE,NDS,NEQ1,LTYPE,NEQF,MBAND,NCOFS,NDOFE)
C        CALLING PROGRAM - MAIN 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *   P R O G R A M                                                 * 
C   *    1. DETERMINES COLUMN HEIGHTS, STORED TEMPORARILY             * 
C   *                  IN ARRAY MAXA(NEQ).                            * 
C   *    2. DETERMINES ADDRESSES = ( I-TH SUBSCRIPT IN A(I) ) FOR     * 
C   *          THE DIAGONAL TERMS IN THE GLOBAL ARRAY A(I)            * 
C   *               ADDRESSES STORED IN MAXA(NEQ+1).                  * 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
      COMMON/RDK1/LW,LR,KEY 
C 
C   DIMENSIONED VARIABLES ID AND MAXA SHARED WITH MAIN PROGRAM COMMON 
C 
      DIMENSION ID(NDS,3),MAXA(NEQ1),LM(10)
C 
C        CLEAR LM AND MAXA ARRAYS 
C 
      DO 5 IEQN=1,NEQ1
    5 MAXA(IEQN)=0
      DO 10 I=1,10 
   10 LM(I)=0 
C 
C   USING EQN. NO.S STORED IN ID(NDS,3) ARRAY TOGETHER WITH NODE
C   NO.S (ND1,ND2,ND3) STORED ON TAPE 1  CALCULATE CONNECTIVITY 
C                  LM ARRAYS FOR EACH ELEMENT.
C 
      REWIND 1
C 
C      L O O P   T H R U   A L L   E L E M E N T S
C 
      DO 40 K=1,NE
C 
C  READ NODE NO.S ASSOCIATED WITH K-TH ELEMENT
C 
      READ(1)ND1,ND2,ND3
C 
C       ASSIGN EQUATION NUMBERS TO ELEMENT CONNECTIVITY ARRAY 
C 
C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 1ST NODE)=ID(ND1,1)=LM(1)
C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 2ND NODE)=ID(ND2,1)=LM(2)
C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 3RD NODE)=ID(ND3,1)=LM(3)
C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 1ST NODE)=ID(ND1,2)=LM(4)
C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 2ND NODE)=ID(ND2,2)=LM(5)
C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 3RD NODE)=ID(ND3,2)=LM(6)
C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 1ST NODE)=ID(ND1,3)=LM(7)
C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 2ND NODE)=ID(ND2,3)=LM(8)
C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 3RD NODE)=ID(ND3,3)=LM(9)
C 
      LM(1)=ID(ND1,1) 
      LM(2)=ID(ND2,1) 
      LM(3)=ID(ND3,1) 
      LM(4)=ID(ND1,2) 
      LM(5)=ID(ND2,2) 
      LM(6)=ID(ND3,2) 
      LM(7)=ID(ND1,3) 
      LM(8)=ID(ND2,3) 
      LM(9)=ID(ND3,3) 
C 
C             ******* EXCLUSIVELY *******
C  FOR THERMAL - MOISTURE LOADS ONLY INCLUDE EXTRA DOF
C  BY ASSINGING EQN NO. NEQF TO LM(10) FOR EACH ELEMENT
C 
      IF(LTYPE.EQ.2)LM(10)=NEQF
C 
C  CHOOSE NO. OF DEGREES OF FREEDOM PER ELEMENT, NDOFE
C 
C       DEFINE NDOFE=9  1. FOR A STRAIN APPLIED NORMAL TO THE 
C                               ELEMENT H-Z PLANE, SLOAD
C                       2. FOR NODE LOADS OR PRESCRIBED DISPL.S 
      NDOFE=9 
C 
C    EXCLUSIVELY DEFINE 
C         NDOFE=10      1. FOR A CHANGE IN TEMPERATURE, TLOAD 
C                       2. FOR A CHANGE IN ABSORBED MOISTURE, WLOAD 
C 
      IF(LTYPE.EQ.2)NDOFE=10 
C 
C         CALCULATE THE COLUMN HEIGHTS ABOVE THE DIAGONAL 
C            OF THE ASSEMBLED GLOBAL STIFFNESS ARRAY
C       (COLUMN HEIGHTS TEMPORARILY STORED IN MAXA ARRAY) 
C 
C   THIS IS DONE BY: 
C         COMPARING VALUES IN LM ARRAYS FOR ELEMENTS
C             WHICH SHARE A DEGREE OF FREEDOM(DOF)
C     COLUMN HEIGHT = MAXA( SHARED DOF ) - LOWEST DOF IN SHARED LM ARRAY
C 
      DO 30 I=1,NDOFE 
      L=LM(I) 
      IF(L.EQ.0)GO TO 30
      DO 20 J=1,NDOFE 
      M=LM(J) 
      IF(M.EQ.0)GO TO 20
      IF(M.GT.L)GO TO 20
      N=L-M+1 
      IF(N.GT.MAXA(L))MAXA(L)=N 
   20 CONTINUE
   30 CONTINUE
   40 CONTINUE
C 
C       E N D   E L E M E N T   L O O P 
C 
C   CALCUALTE WIDTH OF HALFBAND PLUS DIAGONAL, MBAND
C    WRITE COLUMN HEIGHTS AND MBAND ON OUTPUT FILE
C 
      MBAND=0 
      NEQ=NEQ1-1
      DO 55 I=1,NEQ 
      IF(MAXA(I).GT.MBAND)MBAND=MAXA(I) 
   55 CONTINUE
      IF((KEY.EQ.4).OR.(KEY.EQ.7))WRITE(6,50)(I,MAXA(I),I=1,NEQ1)
   50 FORMAT(4(' # COLUMN=',I6,' HEIGHT=',I6,' #')) 
      WRITE(6,56)MBAND 
   56 FORMAT('  WIDTH OF HALFBAND PLUS DIAGONAL, MBAND=',I4)
C 
C   USING COLUMN HEIGHTS STORED IN MAXA( NEQ ) 
C      DETERMINE ADDRESSES = ( I-TH SUBSCRIPT IN A(I) )
C      FOR DIAGONAL TERMS IN THE GLOBAL ARRAY A(I)
C      AND STORE VALUES IN MAXA( NEQ + 1 ). 
C 
C  FIRST COLUMN DIAGONAL ADDRESS BY DEFINITION IS ONE 
C 
      NN=1
      MAXA(1)=1 
      MAXM1=1 
C 
C   ADD COLUMN HEIGHTS TO PREVIOUS DIAGONAL ADDRESS 
C 
      DO 60 I=2,NEQ1
      NN=NN+MAXM1 
      MAXM1=MAXA(I) 
      MAXA(I)=NN
   60 CONTINUE
C 
C  WRITE GLOBAL ARRAY SIZE AND DIAGONAL ADDRESSES ON THE OUTPUT FILE 
C 
      NCOFS=MAXA(NEQ1)
      WRITE(6,70)NCOFS 
   70 FORMAT(/,' NO. OF COEFF.S IN GLOBAL STIFFNESS ARRAY, NCOFS=',I6,/) 
C 
      IF((KEY.EQ.6).OR.(KEY.EQ.7))WRITE(6,75)(I,MAXA(I),I=1,NEQ) 
   75 FORMAT(3(' # COLUMN=',I6,' DIAGONAL ADDRESS=',I6,' #')) 
C 
      RETURN
      END 



C***********************************************************************
      SUBROUTINE ELMNT(ID,MAXA,A,V,NE,NDS,LTYPE,SLOAD,TLOAD,WLOAD,
     *                 NCOFS,NEQ,NEQ1,NEQF,NDOFE) 
C        CALLING PROGRAM - MAIN 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *   P R O G R A M                                                 * 
C   *     CALCULATES THERMAL-MOISTURE STRAINS.                        * 
C   *     CALCULATES ELEMENT ARRAYS LM(I),AK(I,J),AND F(I).           * 
C   *     ASSEMBLES GLOBAL ARRAYS A(I) AND V(I) FROM ELEMENT ARRAYS.  * 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
      COMMON/RDK1/LW,LR,KEY 
      COMMON/RDK3/AL1,AL2,AL3,BT1,BT2,BT3 
C 
C  DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON
C 
      DIMENSION ID(NDS,3),MAXA(NEQ1),A(NCOFS),V(NEQ)
C 
C  DIMENSION ELEMENT ARRAYS 
C 
      DIMENSION AK(10,10),F(10),LM(10)
C 
C        CLEAR A(NCOFS) AND V(NEQ) ARRAYS 
C 
      DO 1 I=1,NCOFS
    1 A(I)=0.0
      DO 2 I=1,NEQ
    2 V(I)=0.0
C 
C  C A L C U L A T E   T H E R M A L - M O I S T U R E   S T R A I N S
C 
C      STRN1,  STRAIN PARALLEL TO PLY FIBER DIRECTION
C      STRN2,  STRAIN NORMAL TO PLY FIBER DIRECTION
C      STRN12, SHEAR STRAIN IN THE THIN PLY 1-2 PLANE 
C      STRN3,  STRAIN NORMAL TO THE PLY 1-2 PLANE 
C 
      IF(LTYPE.LE.1)GO TO 5 
C 
      STRN1=AL1*TLOAD+BT1*WLOAD 
      STRN2=AL2*TLOAD+BT2*WLOAD 
      STRN12=TLOAD*(AL2-AL1)+WLOAD*(BT2-BT1)
      STRN3=AL3*TLOAD+BT3*WLOAD 
C 
C    R E W I N D   S T O R E D   E L E M E N T   D A T A
C 
C       TAPE 1  ND1,ND2,ND3 ....... NODE NO.S ASSIGNED TO EACH ELEMENT 
    5 REWIND 1  
C 
C       TAPE 2  THRAD=THETA(RADIANS) .. ELEMENT FIBER AXIS ORIENTATION
      REWIND 2  
C 
C       TAPE 3  AA,BB,...,GG,ARR ...... ELEMENT GEOMETRIC PARAMETERS
      REWIND 3  
C 
C       TAPE 4  D11,D12,...,D66 ....... ELEMENT STIFFNESSES 
      REWIND 4  
C 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *         L O O P   T H R U   A L L   E L E M E N T S             * 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
      DO 70 K=1,NE
C 
C   F O R M A T I O N   O F   C O N N E C T I V I T Y   A R R A Y 
C                          LM ( NDOFE ) 
C 
C     READ NODE NO.S ASSOCIATED WITH THE K-TH ELEMENT 
C 
      READ(1)ND1,ND2,ND3
C 
C   ASSIGN EQUATION NUMBERS TO ELEMENT CONNECTIVITY ARRAY 
C 
C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 1ST NODE)=ID(ND1,1)=LM(1)
C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 2ND NODE)=ID(ND2,1)=LM(2)
C(EQN NO ASSOCIATED WITH DOF IN N-DIRECTION OF 3RD NODE)=ID(ND3,1)=LM(3)
C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 1ST NODE)=ID(ND1,2)=LM(4)
C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 2ND NODE)=ID(ND2,2)=LM(5)
C(EQN NO ASSOCIATED WITH DOF IN H-DIRECTION OF 3RD NODE)=ID(ND3,2)=LM(6)
C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 1ST NODE)=ID(ND1,3)=LM(7)
C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 2ND NODE)=ID(ND2,3)=LM(8)
C(EQN NO ASSOCIATED WITH DOF IN Z-DIRECTION OF 3RD NODE)=ID(ND3,3)=LM(9)
C 
      LM(1)=ID(ND1,1) 
      LM(2)=ID(ND2,1) 
      LM(3)=ID(ND3,1) 
      LM(4)=ID(ND1,2) 
      LM(5)=ID(ND2,2) 
      LM(6)=ID(ND3,2) 
      LM(7)=ID(ND1,3) 
      LM(8)=ID(ND2,3) 
      LM(9)=ID(ND3,3) 
C 
C             ******* EXCLUSIVELY *******
C  FOR A THERMAL - MOISTURE LOADS ONLY INCLUDE EXTRA DOF  
C  BY ASSIGNING EQN. NO. NEQF TO LM(10) FOR EACH ELEMENT 
C 
      IF(LTYPE.EQ.2)LM(10)=NEQF
C 
C      WRITE LM ARRAY FOR EACH ELEMENT ON OUTPUT FILE 
C 
      IF((KEY.EQ.3).OR.(KEY.EQ.7))GO TO 9 
      GO TO 15
C 
    9 IF(NDOFE.EQ.9)WRITE(6,10)K,NDOFE,(LM(I),I=1,NDOFE) 
      IF(NDOFE.EQ.10)WRITE(6,11)K,NDOFE,(LM(I),I=1,NDOFE) 
   10 FORMAT(1H1,' ELEMENT NO.',I4,'   CONNECTIVITY ARRAY LM(',I2,')',
     *9(I4,1X))
   11 FORMAT(1H1,' ELEMENT NO.',I4,'   CONNECTIVITY ARRAY LM(',I2,')',
     *10(I4,1X)) 
C 
C                 F O R M A T I O N   O F 
C         E L E M E N T A L    S T I F F N E S S
C              A R R A Y    AK ( NDOFE,NDOFE )
C                 -------  AND  ------- 
C     E L E M E N T A L   R I G H T - H A N D - S I D E 
C             F O R C E   V E C T O R   F ( NDOFE ) 
C 
C       READ ELEMENT DATA 
C 
   15 READ(3)AA,BB,CC,DD,EE,GG,ARR
      READ(4)D11,D12,D13,D16,D22,D23,D26,D33,D36,D44,D45,D55,D66
C 
C   CALCULATE STIFFNESS ARRAY FOR EACH ELEMENT
C 
      AK(1,1)=(D55*BB*BB+D66*AA*AA)/ARR
      AK(1,2)=(D55*BB*DD+D66*AA*CC)/ARR
      AK(1,3)=(D55*BB*GG+D66*AA*EE)/ARR
      AK(1,4)=(D26*AA*AA+D45*BB*BB)/ARR
      AK(1,5)=(D26*CC*AA+D45*BB*DD)/ARR
      AK(1,6)=(D26*EE*AA+D45*BB*GG)/ARR
      AK(1,7)=(D36*BB*AA+D45*BB*AA)/ARR
      AK(1,8)=(D36*DD*AA+D45*BB*CC)/ARR
      AK(1,9)=(D36*GG*AA+D45*BB*EE)/ARR
      AK(2,2)=(D55*DD*DD+D66*CC*CC)/ARR
      AK(2,3)=(D55*DD*GG+D66*CC*EE)/ARR
      AK(2,4)=(D26*AA*CC+D45*DD*BB)/ARR
      AK(2,5)=(D26*CC*CC+D45*DD*DD)/ARR
      AK(2,6)=(D26*EE*CC+D45*DD*GG)/ARR
      AK(2,7)=(D36*BB*CC+D45*DD*AA)/ARR
      AK(2,8)=(D36*DD*CC+D45*DD*CC)/ARR
      AK(2,9)=(D36*GG*CC+D45*DD*EE)/ARR
      AK(3,3)=(D55*GG*GG+D66*EE*EE)/ARR
      AK(3,4)=(D26*AA*EE+D45*GG*BB)/ARR
      AK(3,5)=(D26*CC*EE+D45*GG*DD)/ARR
      AK(3,6)=(D26*EE*EE+D45*GG*GG)/ARR
      AK(3,7)=(D36*BB*EE+D45*GG*AA)/ARR
      AK(3,8)=(D36*DD*EE+D45*GG*CC)/ARR
      AK(3,9)=(D36*GG*EE+D45*GG*EE)/ARR
      AK(4,4)=(D22*AA*AA+D44*BB*BB)/ARR
      AK(4,5)=(D22*AA*CC+D44*BB*DD)/ARR
      AK(4,6)=(D22*AA*EE+D44*BB*GG)/ARR
      AK(4,7)=(D44*BB*AA+D23*AA*BB)/ARR
      AK(4,8)=(D44*BB*CC+D23*AA*DD)/ARR
      AK(4,9)=(D44*BB*EE+D23*AA*GG)/ARR
      AK(5,5)=(D22*CC*CC+D44*DD*DD)/ARR
      AK(5,6)=(D22*CC*EE+D44*DD*GG)/ARR
      AK(5,7)=(D44*DD*AA+D23*CC*BB)/ARR
      AK(5,8)=(D44*DD*CC+D23*CC*DD)/ARR
      AK(5,9)=(D44*DD*EE+D23*CC*GG)/ARR
      AK(6,6)=(D22*EE*EE+D44*GG*GG)/ARR
      AK(6,7)=(D44*GG*AA+D23*EE*BB)/ARR
      AK(6,8)=(D44*GG*CC+D23*EE*DD)/ARR
      AK(6,9)=(D44*GG*EE+D23*EE*GG)/ARR
      AK(7,7)=(D33*BB*BB+D44*AA*AA)/ARR
      AK(7,8)=(D33*BB*DD+D44*AA*CC)/ARR
      AK(7,9)=(D33*BB*GG+D44*AA*EE)/ARR
      AK(8,8)=(D33*DD*DD+D44*CC*CC)/ARR
      AK(8,9)=(D33*DD*GG+D44*CC*EE)/ARR
      AK(9,9)=(D33*GG*GG+D44*EE*EE)/ARR
C 
C   CALCULATE THE RIGHT-HAND-SIDE FORCE VECTOR DUE TO A 
C   UNIFORMLY APPLIED STRAIN (SLOAD) NORMAL TO THE H-Z PLANE
C 
      F(1)=-D16*AA*SLOAD
      F(2)=-D16*CC*SLOAD
      F(3)=-D16*EE*SLOAD
      F(4)=-D12*AA*SLOAD
      F(5)=-D12*CC*SLOAD
      F(6)=-D12*EE*SLOAD
      F(7)=-D13*BB*SLOAD
      F(8)=-D13*DD*SLOAD
      F(9)=-D13*GG*SLOAD
C 
C   IF LOAD TYPE IS A MECHANICALLY APPLIED STRAIN NORMAL TO 
C   THE ELEMENT PLANE THEN LOOP AROUND ADDITIONAL EQUATIONS 
C 
      IF(LTYPE.LE.1)GO TO 20
C 
C   FOR THERMAL-MOISTURE LOADING (EXCLUSIVELY) 
C     1. INCLUDE TERMS IN THE TENTH ROW/COLUMN OF STIFFNESS ARRAY 
C 
      AK(1,10)=D16*AA 
      AK(2,10)=D16*CC 
      AK(3,10)=D16*EE 
      AK(4,10)=D12*AA
      AK(5,10)=D12*CC
      AK(6,10)=D12*EE
      AK(7,10)=D13*BB
      AK(8,10)=D13*DD
      AK(9,10)=D13*GG
      AK(10,10)=D11*ARR 
C 
C     2. CALCULATE RIGHT-HAND-SIDE FORCE VECTORS
C 
C                          2.1 TRANSFORM 
C                     THERMAL-MOISTURE STRAINS 
C         FROM MATERIAL AXIS(1,2,3) TO LAMINATE AXIS(N,H,Z) 
C 
      READ(2)THRAD
C 
      C=COS(THRAD)
      S=SIN(THRAD)
      ENO=C*C*STRN1+S*S*STRN2 
      EHO=S*S*STRN1+C*C*STRN2 
      ENHO=2.0*S*C*STRN12 
      EZO=STRN3 
C 
C           2.2 STORE THERMAL-MOISTURE STRAINS ON TAPE 7
C 
      WRITE(7)ENO,EHO,ENHO,EZO
C 
C       2.3 CALCULATE LOADS DUE TO THERMAL-MOISTURE STRAINS 
C 
      F(1)=(D16*ENO+D26*EHO+D36*EZO+D66*ENHO)*AA 
      F(2)=(D16*ENO+D26*EHO+D36*EZO+D66*ENHO)*CC 
      F(3)=(D16*ENO+D26*EHO+D36*EZO+D66*ENHO)*EE 
      F(4)=(D12*ENO+D22*EHO+D23*EZO+D26*ENHO)*AA 
      F(5)=(D12*ENO+D22*EHO+D23*EZO+D26*ENHO)*CC 
      F(6)=(D12*ENO+D22*EHO+D23*EZO+D26*ENHO)*EE 
      F(7)=(D13*ENO+D23*EHO+D33*EZO+D36*ENHO)*BB 
      F(8)=(D13*ENO+D23*EHO+D33*EZO+D36*ENHO)*DD 
      F(9)=(D13*ENO+D23*EHO+D33*EZO+D36*ENHO)*GG 
      F(10)=(D11*ENO+D12*EHO+D13*EZO+D16*ENHO)*ARR 
C 
C   ASSIGN VALUES TO SYMMETRIC PORTION OF ELEMENTAL 
C         STIFFNESS ARRAY  AK(I,J)=AK(J,I)
C 
   20 JJ=1
      DO 40 II=1,NDOFE
      DO 30 J=JJ,NDOFE
   30 AK(J,II)=AK(II,J) 
   40 JJ=JJ+1 
C 
C    WRITE STIFFNESS ARRAY AND FORCE VECTOR ON OUTPUT FILE
C 
      IF((KEY.EQ.3).OR.(KEY.EQ.7))GO TO 45 
      GO TO 60
C 
   45 IF(NDOFE.EQ.10)GO TO 55
C 
      WRITE(6,50)K 
   50 FORMAT(//,15X,'ELEMENT NO.',I4,/,12X,'STIFFNESS ARRAY AK(I,J)', 
     */,' ROW/COL    1',9X,'2',9X,'3',9X,'4',9X,'5',9X,'6',9X,'7',9X,'8' 
     *,9X,'9',5X,'FORCE') 
      DO 51 I=1,NDOFE 
   51 WRITE(6,52)I,(AK(I,J),J=1,NDOFE),F(I)
   52 FORMAT(1X,I4,3X,10(E10.3)) 
      GO TO 60
C 
   55 WRITE(6,56)K 
   56 FORMAT(//,15X,'ELEMENT NO.',I4,/,12X,'STIFFNESS ARRAY AK(I,J)', 
     */,' ROW/COL    1',9X,'2',9X,'3',9X,'4',9X,'5',9X,'6',9X,'7',9X,'8' 
     ,9X,'9',9X,'10',5X,'FORCE')
      DO 57 I=1,NDOFE 
   57 WRITE(6,58)I,(AK(I,J),J=1,NDOFE),F(I)
   58 FORMAT(1X,I4,3X,11(E10.3)) 
C 
C                     A  S  S  E  M  B  L  E
C    ELEMENTAL ARRAYS ( AK , F ) INTO GLOBAL ARRAYS ( A , V ) 
C USING THE CONNECTIVITY AND DIAGONAL ADDRESS ARRAYS ( LM , MAXA )
C 
   60 CALL ASSEMB(NDOFE,AK,F,LM,MAXA,A,V,NEQ1,NCOFS,NEQ)
C 
   70 CONTINUE
C 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *           E N D   E L E M E N T   L O O P                       * 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
C   WRITE ASSEMBLED  1. GLOBAL STIFFNESS ARRAY ON OUTPUT FILE 
C                    2. RIGHT-HAND-SIDE LOAD VECTOR ON OUTPUT FILE
C 
      IF((KEY.EQ.4).OR.(KEY.EQ.7))GO TO 74 
      RETURN
C 
C    WRITE HEADER AND FIRST DIAGONAL TERM ON OUTPUT FILE
C 
   74 WRITE(6,75)A(1)
   75 FORMAT(1H1,//,' GLOBAL STIFFNESS ARRAY',//,'   COLUMN NO.    1',
     */,7X,' A(   1)=',E11.4,'  1 RST DIAGONAL POSITION') 
C 
C   START LOOP AT 2 ND DIAGONAL TERM
C 
      IDIAG=2 
C 
C    LOOP THRU A(I) TERMS 
C 
      DO 80 I=2,NCOFS 
      IF(I.EQ.MAXA(IDIAG))GO TO 77
      WRITE(6,76)I,A(I)
   76 FORMAT(8X,'A(',I4,')=',E14.7) 
      GO TO 80
   77 IF(I.NE.NCOFS)WRITE(6,78)IDIAG,I,A(I),IDIAG
   78 FORMAT(/,'   COLUMN NO.',I4,/,8X,'A(',I4,')=',E14.7,
     *' DIAGONAL NO.',I4) 
C 
      IF(I.EQ.NCOFS)WRITE(6,79)I,A(I)
   79 FORMAT(//,8X,'A(',I4,')=',E14.7,' LAST STORED VALUE IN A(I)') 
C 
      IDIAG=IDIAG+1 
C 
   80 CONTINUE
C 
C   WRITE RIGHT-HAND-SIDE LOAD VECTOR ON OUTPUT FILE
C 
      WRITE(6,81)
   81 FORMAT(1H1,//,'  RIGHT-HAND-SIDE LOAD VECTOR',//) 
C 
      DO 82 I=1,NEQ
   82 WRITE(6,83)I,V(I)
   83 FORMAT(7X,'V(',I4,')=',E14.7) 
C 
      RETURN
      END 



C***********************************************************************
      SUBROUTINE ASSEMB(NDOFE,AK,F,LM,MAXA,A,V,NEQ1,NCOFS,NEQ)
C        CALLING PROGRAM - ELMNT
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *   P R O G R A M                                                 * 
C   *     ASSEMBLES : 1. EACH ELEMENT STIFFNESS ARRAY AK(I,J)         *
C   *                    ONTO THE GLOBAL ARRAY A(I).                  * 
C   *                 2. EACH ELEMENT FORCE VECTOR F(I) ONTO THE      * 
C   *                    RIGHT-HAND-SIDE LOAD VECTOR V(I).            * 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
      DIMENSION AK(10,10),F(10),LM(10)
      DIMENSION MAXA(NEQ1),A(NCOFS),V(NEQ)
C 
      DO 20 I=1,NDOFE 
      L=LM(I) 
      IF(L.EQ.0)GO TO 20
      NN=MAXA(L)
      DO 10 J=I,NDOFE 
      M=LM(J) 
      IF(M.EQ.0)GO TO 10
      MM=MAXA(M)
      KK=M-L
      MM=MM+KK
      IF(KK.LT.0)MM=NN-KK 
C 
      A(MM)=A(MM)+AK(I,J) 
   10 CONTINUE
C 
      V(L)=V(L)+F(I)
   20 CONTINUE
C 
      RETURN
      END 



C***********************************************************************
      SUBROUTINE NDLOAD(ID,V,NDS,NEQ,NLOAD) 
C        CALLING PROGRAM - MAIN 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *   P R O G R A M                                                 * 
C   *     ADDS NODAL LOADS TO THE EXISTING ASSEMBLED GLOBAL           * 
C   *            RIGHT-HAND-SIDE LOAD VECTOR V(NEQ)                   * 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
      COMMON/RDK1/LW,LR,KEY 
C 
C   DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON 
C 
      DIMENSION ID(NDS,3),V(NEQ)
C 
C   USING EQN NO.S STORED IN ID(NDS,2) ARRAY
C      RECOVER EQN NO.(IEQN) WHICH CORRESPONDS TO NODAL LOAD
C      DIRECTION (IDIRN) AND ADD NODE LOAD (FLOAD) TO THE 
C      RIGHT-HAND-SIDE LOAD VECTOR (V(IEQN)). 
C 
C   WRITE HEADER ON OUTPUT FILE 
C 
      WRITE(6,10)NLOAD
   10 FORMAT(1H1,//,'  NUMBER OF APPLIED LOADS=',I4,/)
C 
C        L O O P   T H R U   L O A D S
C 
      DO 60 ILOAD=1,NLOAD 
C 
C   READ 1. NOD  , NODE NO. WHERE I-TH LOAD IS APPLIED
C        2. IDIRN, DIRECTION OF I-TH LOAD 
C        3. FLOAD, MAGNITUDE OF I-TH LOAD 
C 
      READ(5,20)NOD,IDIRN,FLOAD
   20 FORMAT(2I6,F12.0) 
C 
C     RECOVER EQN NO. 
C 
      IEQN=ID(NOD,IDIRN)
C 
      WRITE(6,30)ILOAD,NOD,IDIRN,IEQN,FLOAD,IEQN,V(IEQN) 
   30 FORMAT(/,5X,'ILOAD=',I4,'  NODE NO.=',I4,'  LOAD DIRECTION=', 
     *I3,'  EQN NO.=',I5,/,10X,'LOAD MAGNITUDE=',E10.3,/, 
     *15X,'OLD RIGHT-HAND-SIDE LOAD VECTOR, V(',I4,')=',E11.4)
C 
C   ADD NODAL LOAD (FLOAD) TO EXISTING LOAD VECTOR
C 
      V(IEQN)=V(IEQN)+FLOAD 
C 
      WRITE(6,50)IEQN,V(IEQN)
   50 FORMAT(15X,'NEW RIGHT-HAND-SIDE LOAD VECTOR,V(',I4,')=',E11.4,/)
C 
   60 CONTINUE
C 
      RETURN
      END 



C***********************************************************************
      SUBROUTINE NDDISP(ID,MAXA,A,V,NDS,NEQ1,NCOFS,NEQ,NDISP) 
C        CALLING PROGRAM - MAIN 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *   P R O G R A M                                                 
C   *     ADDS A HUGE SPRING STIFFNESS (HUGEK) TO THE ASSEMBLED       * 
C   *     STIFFNESS (A(MAXA(IEQN))) AND LOAD (V(IEQN)) COMPONENTS     * 
C   *     SO AS TO SIMULATE A PRESCRIBED NODAL DISPLACEMENT.          * 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
      COMMON/RDK1/LW,LR,KEY 
C 
C   DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON 
C 
      DIMENSION ID(NDS,3),MAXA(NEQ1),A(NCOFS),V(NEQ)
C 
C   USING EQN NO.S STORED IN ID(NDS,3) ARRAY
C      RECOVER EQN NO. (IEQN) CORRESPONDING TO NODE DISPLACEMENT (DISP).
C      THEN   1. ADD LARGE SPRING STIFFNESS (HUGEK) TO DIAGONAL 
C                TERN ( A(MAXA(IEQN)) ) OF THE GLOBAL STIFFNESS ARRAY.
C             2. ADD LARGE SPRING FORCE ( HUGEK*DISP ) TO THE 
C                RIGHT-HAND-SIDE LOAD VECTOR ( V(IEQN) ). 
C 
      HUGEK=1.0E+20 
C 
C     PRIOR TO ADDITION OF HUGE SPRING STIFFNESS
C     STORE GLOBAL STIFFNESS ARRAY ON TAPE 8
C 
      WRITE(8)(A(I),I=1,NCOFS)
C 
C   WRITE HEADER ON OUTPUT FILE 
C 
      WRITE(LW,10)NDISP 
   10 FORMAT(1H1,//,'  NUMBER OF PRESCRIBED DISPL.S=',I4,/) 
C 
C     L O O P   T H R U   P R E S C R I B E D   D I S P L. S
C 
      DO 60 IDISP=1,NDISP 
C 
C   READ 1. NOD  , NODE NO. WHERE I-TH DISPL. IS PRESCRIBED
C        2. IDIRN, DIRECTION OF I-TH DISPL. 
C        3. DISP , MAGNITUDE OF I-TH DISPL. 
C 
      READ(LR,20)NOD,IDIRN,DISP 
   20 FORMAT(2I6,F12.0) 
C 
C   STORE PRIESCRIBED DISPL. DATA ON TAPE 9 
C 
      WRITE(9)NOD,IDIRN,DISP
C 
C    RECOVER EQN NO. (IEQN) AND DIAGONAL ADDRESSES MAXA(IEQN) 
C 
      IEQN=ID(NOD,IDIRN)
      II=MAXA(IEQN) 
C 
      WRITE(LW,30)IDISP,NOD,IDIRN,IEQN,DISP,II,A(II),IEQN,V(IEQN) 
   30 FORMAT(/,5X,' IDISP=',I4,'  NODE NO.=',I4,'  DISPL. DIRECTION=',
     *I3,/,6X,'EQN NO.=',I5,14X,'DISPL. MAGNITUDE=',E10.3,//, 
     *5X,'OLD DIAGONAL A(',I4,')=',E11.4, 
     *'  OLD RIGHT-HAND-SIDE VECTOR V(',I4,')=',E11.4)
C 
C   ADD SPRING STIFFNESS AND SPRING FORCE 
C 
      A(II)=A(II)+HUGEK 
      V(IEQN)=V(IEQN)+HUGEK*DISP
C 
      WRITE(LW,50)II,A(II),IEQN,V(IEQN) 
   50 FORMAT(5X,'NEW DIAGONAL A(',I4,')=',E11.4,
     *'  NEW RIGHT-HAND-SIDE VECTOR V(',I4,')=',E11.4,/)
C 
   60 CONTINUE
C 
      RETURN
      END 
 



C***********************************************************************
      SUBROUTINE COLSOL(MAXA,A,V,NN,NWK,NNM,KKK)
C         CALLING PROGRAM - MAIN
C *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
C *    THIS SUBROUTINE IS CALLED TWICE                                 *
C *                                                                    * 
C * THE FIRST TIME                                                     *
C *     * TO PERFORM L D L TRANSPOSE OPERATION                         * 
C *          ON THE GLOBAL STIFFNESS MATRIX                            * 
C *                                                                    * 
C * THE SECOND TIME                                                    * 
C *       TO SOLVE THE SYSTEM OF EQUATIONS                             *
C *                                                                    *
C * PROGRAM                                                            *
C *    TO SOLVE FINITE ELEMENT STATIC EQUILIBRIUM EQUATIONS IN         *
C *    CORE, USING COMPACTED STORAGE AND COLUMN REDUCTION SCHEME       *
C * - - INPUT VARIABLES - -                                            *
C *    A(NWK)      = STIFFNESS MATRIX STORED IN COMPACTED FORM         *
C *    V(NN)       = RIGHT-HAND-SIDE LOAD VECTOR                       *
C *    MAXA(NNM)   = VECTOR CONTAINING ADDRESSES OF DIAGONAL           *
C *                  ELEMENTS OF STIFFNESS MATRIX IN A                 *
C *    NN = NEQ    = NUMBER OF EQUATIONS                               *
C *  NWK (NOT USED)= NUMBER OF ELEMENTS BELOW SKYLINE OF MATRIX        *
C *  NNM (NOT USED)=  NN  +  1                                         *
C *    KKK         = INPUT FLAG                                        *
C *        EQ. 1 TRIANLGULARIZATION OF STIFFNESS MATRIX                *
C *        EQ. 2 REDUCTION AND BACK-SUBSTITUTION OF LOAD VECTOR        *
C *    IOUT        = NUMBER OF OUTPUT DEVICE                           *
C * - - OUTPUT - -                                                     *
C *    A(NWK)      = D AND L - FACTORS OF STIFFNESS MATRIX             *
C *    V(NN)       = DISPLACEMENT VECTOR                               *
C *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
C 
C     IMPLICIT REAL*8(A-H,0-Z)                                          
C 
C *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
C *    THIS PROGRAM IS USED IN SINGLE PRECISION ARITHMETIC ON          *
C *    CDC EQUIPMENT AND DOUBLE PRECISION ARITHMETIC ON IBM            *
C *    OR UNIVAC MACHINES.  ACTIVATE, DEACTIVATE OR ADJUST ABOVE       *
C *    CARD FOR SINGLE OR DOUBLE PRECISION ARITHMETIC.                 *
C *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 
C 
      COMMON/RDK1/IOUT,LR,KEY 
C 
      DIMENSION MAXA(NNM),A(NWK),V(NN)
C 
C *** PERFORM L*D*L(T) FACTORIZATION OF STIFFNESS MATRIX                
C                                                                       
      IF (KKK-2) 40,150,150                                             
   40 DO 140 N=1,NN                                                     
      KN=MAXA(N)                                                        
      KL=KN+1                                                           
      KU=MAXA(N+1)-1                                                    
      KH=KU-KL                                                          
      IF (KH) 110,90,50                                                 
   50 K=N-KH                                                            
      IC=0                                                              
      KLT=KU                                                            
      DO 80 J=1,KH                                                      
      IC=IC+1                                                           
      KLT=KLT-1                                                         
      KI=MAXA(K)                                                        
      ND=MAXA(K+1)-KI-1                                                 
      IF (ND) 80,80,60                                                  
   60 KK=IC                                                             
      IF (ND.LT.IC)KK=ND                                                
      C=0.0                                                             
      DO 70 L=1,KK                                                      
   70 C=C+A(KI+L)*A(KLT+L)                                              
      A(KLT)=A(KLT)-C                                                   
   80 K=K+1                                                             
   90 K=N                                                               
      B=0.0                                                             
      DO 100 KK=KL,KU                                                   
      K=K-1                                                             
      KI=MAXA(K)                                                        
      C=A(KK)/A(KI)                                                     
      B=B+C*A(KK)                                                       
  100 A(KK)=C                                                           
      A(KN)=A(KN)-B                                                     
  110 IF (A(KN)) 120,120,140                                            
  120 WRITE (IOUT,2000) N,A(KN)                                         
      STOP                                                              
  140 CONTINUE                                                          
      RETURN                                                            
C                                                                       
C *** REDUCE RIGHT-HAND-SIDE LOAD VECTOR                                
C                                                                       
  150 DO 180 N=1,NN                                                     
      KL=MAXA(N)+1                                                      
      KU=MAXA(N+1)-1                                                    
      IF (KU-KL) 180,160,160                                            
  160 K=N                                                               
      C=0.0                                                             
      DO 170 KK=KL,KU                                                   
      K=K-1                                                             
  170 C=C+A(KK)*V(K)                                                    
      V(N)=V(N)-C                                                       
  180 CONTINUE                                                          
C                                                                       
C *** BACK SUBSTITUTE                                                   
C                                                                       
      DO 200 N=1,NN                                                     
      K=MAXA(N)                                                         
  200 V(N)=V(N)/A(K)                                                    
      IF (NN.EQ.1) RETURN                                               
      N=NN                                                              
      DO 230 L=2,NN                                                     
      KL=MAXA(N)+1                                                      
      KU=MAXA(N+1)-1                                                    
      IF (KU-KL) 230,210,210                                            
  210 K=N                                                               
      DO 220 KK=KL,KU                                                   
      K=K-1                                                             
  220 V(K)=V(K)-A(KK)*V(N)                                              
  230 N=N-1                                                             
      RETURN                                                            
C 
 2000 FORMAT (//48H STOP - STIFFNESS MATRIX NOT POSITIVE DEFINITE     ,/
     */  32H NONPOSITIVE PIVOT FOR EQUATION   ,I4//                     
     *   10H PIVOT =   ,E20.12)                                         
      END 



C***********************************************************************
      SUBROUTINE NDFORC(ID,MAXA,A,V,NDS,NEQ1,NCOFS,NEQ,NDISP,MBAND) 
C        CALLING PROGRAM - MAIN 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *   P R O G R A M                                                 * 
C   *     RECOVER FORCES AT NODES CORRESPONDING TO PRESCRIBED         * 
C   *     NODAL DISPLACEMENTS.  FORCES CALCULATED FROM COMPONENTS     * 
C   *     OF GLOBAL ARRAY A(I) STORED ON TAPE 8 PRIOR TO SUPERPOSING  * 
C   *     A LARGE SPRING STIFFNESS (HUGEK) ON DIAGONAL TERMS.         * 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
      COMMON/RDK1/LW,LR,KEY 
C 
C   DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON 
C 
      DIMENSION ID(NDS,3),MAXA(NEQ1),A(NCOFS),V(NEQ)
C 
C   CALCULATE BANDWITH (MB) NOT INCLUDING DIAGONAL TERM
C   CALCULATE (LIMIT) EQN NO. WHERE TRUNCATION OF TERMS TO
C     RIGHT OF THE DIAGONAL OCCURS FOR (IEQN.GT.LIMIT)
C 
      MB=MBAND-1
      LIMIT=NEQ1-MBAND
C 
C   REPLACE GLOBAL STIFFNESS ARRAY TERMS A(NCOFS) USED IN COLSOL WITH 
C   VALUES STORED ON TAPE 8 PRIOR TO THE ADDITION OF HUGEK STIFFNESS. 
C 
      REWIND 8
      READ(8)(A(I),I=1,NCOFS) 
C 
C         USING EQN NO.S STORED IN ID(NDS,3) AND MAXA(NEQ+1) ARRAYS
C      RECOVER EQN NO. (IEQN) CORRESPONDING TO PRESCRIBED NODAL DISPL.
C 
C   TO CALCULATE THE FORCE IN THE DIRECTION OF THE PRESCRIBED 
C   DISPLACEMENT, SUM UP THE PRODUCT OF KNOWN DISPLACEMENTS STORED
C   IN ARRAY V(I) WITH THE ROW OF STIFFNESSES A(I) TERMS STORED 
C   ON TAPE 8, WHICH COMPRISES EQN NO. (IEQN).
C 
C     WRITE HEADER ON OUTPUT FILE 
C 
      WRITE(6,10)NDISP 
   10 FORMAT(1H1,//,'  NUMBER OF FORCES CORRESPONDING TO PRESCRIBED NODA
     *L DISPLACEMENTS=',I4)
      WRITE(6,999) 
  999 FORMAT(/,80(1H*),/) 
C 
      REWIND 9
C 
C     L O O P    T H R U    P R E S C R I B E D   D I S P L. S
C 
      DO 70 IDISP=1,NDISP 
C 
C   RECOVER PRESCRIBED DISPL. DATA STORED ON TAPE 9 
C      NOD  , NODE NO. WHERE I-TH DISPL. IS PRESCRIBED
C      IDIRN, DIRECTION OF I-TH DISPL.
C      DISP , MAGNITUDE OF I-TH DISPL.
C 
      READ(9)NOD,IDIRN,DISP 
C 
      IEQN=ID(NOD,IDIRN)
      II=MAXA(IEQN) 
C 
C   INITIALIZE VARIABLES USED TO POSITION SUMMATION OF TERMS
C               AT LEFT SIDE OF EQN NO. (IEQN) 
C 
C      N , NO. OF COMPONENTS TO THE LEFT AND INCLUDING DIAGONAL 
C      IN, SUBSCRIPT OF FIRST DISPL. IN EQN NO. (IEQN)
C      IM, SUBSCRIPT OF FIRST STIFFNESS COMPONENT IN EQN NO. (IEQN) 
C 
      N=MAXA(IEQN+1)-II 
      IN=IEQN-N+1 
      IM=MAXA(IEQN+1)-1 
C 
C   SUM UP PRODUCTS OF DISPLACEMENTS V(I) AND STIFFNESSES A(I)
C          FOR TERMS TO THE LEFT AND INCLUDING DIAGONAL 
C 
      FORCE=0.0 
C 
      DO 30 K=1,N 
      FORCE=FORCE+V(IN)*A(IM) 
C 
      IF(KEY.EQ.7)WRITE(6,20)IN,V(IN),IM,A(IM),FORCE 
   20 FORMAT(/,'  DISPLACEMENT, V(',I4,')=',E11.4,
     *'   STIFFNESS, A(',I4,')=',E11.4,'   FORCE=',E11.4) 
C 
      IM=IM-1 
      IN=IN+1 
   30 CONTINUE
C 
C     SUM UP REMAINING PRODUCT OF TERMS TO RIGHT OF DIAGONAL.
C   BUT FIRST RECALCULATE TRUNCATED BANDWIDTH (MB) TO THE RIGHT  
C  OF THE DIAGONAL IF THE EQN NO. (IEQN) IS GREATER THAN (LIMIT)
C 
      IF(IEQN.GT.LIMIT)MB=MBAND-1-IEQN+LIMIT
      IF(MB.EQ.0)GO TO 55 
C 
      DO 50 IB=1,MB 
      IM=MAXA(IN)+IB
      IF(IM.GE.MAXA(IN+1))GO TO 40
      FORCE=FORCE+V(IN)*A(IM) 
      IF(KEY.EQ.7)WRITE(6,20)IN,V(IN),IM,A(IM),FORCE 
   40 IN=IN+1 
   50 CONTINUE
C 
   55 WRITE(6,60)NOD,IDIRN,IEQN,DISP,FORCE
   60 FORMAT(///,'  NODE NO.=',I5,'  DIRECTION=',I2,'   EQN. NO.=',I5,/,
     *10X,'PRESCRIBED DISPLACEMENT=',E11.4,'  RESULTING FORCE=',E11.4)
      WRITE(6,999) 
   70 CONTINUE
C 
C     E N D   P R E S C R I B E D   D I S P L. S   L O O P
C 
      RETURN
      END 



C***********************************************************************
      SUBROUTINE DISP(ID,U,V,NDS,NEQ,NEQF,LTYPE,UKEEP)
C        CALLING PROGRAM - MAIN 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *   P R O G R A M                                                 * 
C   *     1. RECOVERS NODAL DISPL.S STORED IN V(NEQ) BY COLSOL.       * 
C   *         AND ASSIGNS NODE DISPL.S TO ARRAY U(NDS,3).             * 
C   *     2. RECOVERS THERMAL-MOISTURE STRAIN (UKEEP).                * 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
      COMMON/RDK1/LW,LR,KEY 
C 
C   DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON 
C 
      DIMENSION ID(NDS,3),U(NDS,3),V(NEQ) 
C 
C    RECOVER TOTAL STRAIN NORMAL TO ELEMENT H-Z PLANE
C     DUE TO A CHANGE IN THERMAL-MOISTURE CONDITION
C 
      IF(LTYPE.EQ.2)UKEEP=V(NEQF) 
      IF(LTYPE.EQ.2)WRITE(6,1)UKEEP
    1 FORMAT(1H1,' FOR CHANGES IN THERMAL-MOISTURE CONDITIONS',/, 
     *' TOTAL STRAIN NORMAL TO H-Z ELEMENT PLANE=',E14.7) 
C 
C   RECOVER NODAL DISPLACEMENTS 
C 
      IEQN=1
      DO 10 I=1,NDS 
      DO 10 J=1,3 
      U(I,J)=0.0
      IF(IEQN.EQ.NEQF)IEQN=IEQN+1 
      IF(ID(I,J).EQ.0)GO TO 10
      U(I,J)=V(IEQN)
      IEQN=IEQN+1 
   10 CONTINUE
C 
      RETURN
      END 



C***********************************************************************
      SUBROUTINE STR(U,EN,EH,EZ,EHZ,ENZ,ENH,NE,NDS,LTYPE,SLOAD,UKEEP)
C        CALLING PROGRAM - MAIN 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C   *   P R O G R A M                                                 * 
C   *     CALCULATES 1. STR AINS FROM DISPLACEMENTS                   * 
C   *                2. STR ESSESES FROM STRAINS                      * 
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
C 
      COMMON/RDK1/LW,LR,KEY 
C 
C   DIMENSIONED VARIABLES SHARED WITH MAIN PROGRAM COMMON 
C 
      DIMENSION U(NDS,3),EN(NE),EH(NE),EZ(NE),EHZ(NE),ENZ(NE),ENH(NE) 
C 
C   R E W I N D   S T O R E D   E L E M E N T   D A T A 
C 
C       TAPE 1  ND1,ND2,ND3 ....... NODE NO.S ASSIGNED TO EACH ELEMENT
      REWIND 1
C   
C       TAPE 3  AA,BB,...GG,ARR ... ELEMENT GEOMETRIC PARAMETERS
      REWIND 3
C 
C       TAPE 4  D11,D12,...,D66 ... ELEMENT STIFFNESSES 
      REWIND 4
C 
C       TAPE 7  ENO,EHO,ENHO,EZO .. ELEMENT THERMAL-MOISTURE STRAINS
      REWIND 7
C 
C   WRITE HEADER ON OUTPUT FILE 
C 
      WRITE(6,10)
   10 FORMAT(1H1,//,T10,14H** STRESSES **,//,T2,8HELE. NO.,T15,4HSIGN,
     *T30,4HSIGH,T45,4HSIGZ,T60,5HTAUHZ,T75,5HTAUNZ,T90,5HTAUNH,//) 
C 
C       L O O P   T H R U   E L E M E N T S 
C 
      DO 40 K=1,NE
C 
C   READ ELEMENT DATA 
      READ(1)ND1,ND2,ND3
      READ(3)AA,BB,CC,DD,EE,GG,ARR
      READ(4)D11,D12,D13,D16,D22,D23,D26,D33,D36,D44,D45,D55,D66
C 
C   CALCULATE TOTAL STRAINS FROM DISPLACEMENTS
C 
      EN(K)=SLOAD 
      IF(LTYPE.EQ.2)EN(K)=UKEEP 
      EH(K)=(AA*U(ND1,2)+CC*U(ND2,2)+EE*U(ND3,2))/ARR
      EZ(K)=(BB*U(ND1,3)+DD*U(ND2,3)+GG*U(ND3,3))/ARR
      EHZ(K)=(BB*U(ND1,2)+DD*U(ND2,2)+GG*U(ND3,2)+AA*U(ND1,3)+
     1CC*U(ND2,3)+EE*U(ND3,3))/ARR
      ENZ(K)=(BB*U(ND1,1)+DD*U(ND2,1)+GG*U(ND3,1))/ARR
      ENH(K)=(AA*U(ND1,1)+CC*U(ND2,1)+EE*U(ND3,1))/ARR
C 
C   IF LOAD CASE IS A CHANGE IN THERMAL-MOISTURE CONDITION
C     THEN RECOVER ELASTIC STRAINS FROM TOTAL STRAINS 
C             OTHERWISE LOOP TO STATEMENT 20
C 
      IF(LTYPE.LE.1)GO TO 20
C 
C   READ THERMAL-MOISTURE STRAINS 
C 
      READ(7)ENO,EHO,ENHO,EZO 
C 
C        RECOVER MECHANICAL (ELASTIC) STRAINS FROM TOTAL STRAINS
C 
      EN(K)=EN(K)-ENO 
      EH(K)=EH(K)-EHO 
      ENH(K)=ENH(K)-ENHO
      EZ(K)=EZ(K)-EZO 
C 
C   CALCULATE STRESSES FROM ELASTIC STRAINS 
C 
   20 SIGN=D11*EN(K)+D12*EH(K)+D13*EZ(K)+D16*ENH(K) 
      SIGH=D12*EN(K)+D22*EH(K)+D23*EZ(K)+D26*ENH(K) 
      SIGZ=D13*EN(K)+D23*EH(K)+D33*EZ(K)+D36*ENH(K) 
      SIGHZ=D44*EHZ(K)+D45*ENZ(K) 
      SIGNZ=D45*EHZ(K)+D55*ENZ(K) 
      SIGNH=D16*EN(K)+D26*EH(K)+D36*EZ(K)+D66*ENH(K)
C 
C   WRITE STRESSES ON OUTPUT FILE 
C 
      WRITE(6,30)K,SIGN,SIGH,SIGZ,SIGHZ,SIGNZ,SIGNH
      WRITE(7,31)K,SIGN,SIGH,SIGZ,SIGHZ,SIGNZ,SIGNH
   30 FORMAT(T2,I5,6(1X,E14.7)) 
   31 FORMAT(I5,6(1X,E12.5))
C 
   40 CONTINUE
C 
C           E N D   E L E M E N T   L O O P 
C 
C   WRITE STRAINS AND DISPLACEMENTS ON OUTPUT FILE
C 
      IF((KEY.EQ.5).OR.(KEY.EQ.7))GO TO 45 
      RETURN
C 
C          STRAINS
C 
   45 WRITE(6,50)
   50 FORMAT(1H1,//,T10,13H** STRAINS **,//,T2,8HELE. NO.,T15,2HEN,T30, 
     *2HEH,T45,2HEZ,T60,3HEHZ,T75,3HENZ,T90,3HENH,//) 
C 
      DO 60 K=1,NE
   60 WRITE(6,70)K,EN(K),EH(K),EZ(K),EHZ(K),ENZ(K),ENH(K)
   70 FORMAT(T2,I5,6(1X,E14.7)) 
C 
C          DISPLACEMENTS
C 
      WRITE(6,80)
   80 FORMAT(1H1,//,16X,29H          NODAL DISPLACEMENTS,/,10X,65HNODE
     * NO.    U-DISPLACEMENT    V-DISPLACEMENT     W-DISPLACEMENT)
C 
      DO 90 I=1,NDS 
   90 WRITE(6,100)I,U(I,1),U(I,2),U(I,3) 
  100 FORMAT(13X,I5,3(7X,E12.5))
      RETURN
      END 


Comments /Concerns contact information:
Course content in this section: email Dr Ron Kriz
Interactive Courseware Development in this section: email Dr Ron Kriz
Information created 1/30/00 last updated 2/2/00
http://www.jwave.vt.edu/crcd/kriz/modules/module06/fem_prog/fem_prog_dev.html