Ronald D. Kriz
Associate Professor
Engineering Science and Mechanics
Virginia Polytechnic Institute and State University
Blacksburg, Virginia 24061
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.
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.
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.
MAIN Program: Linear Finite Element Model with In-core equation solver (LFEMI)
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.
eqnum
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.
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 |
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